context: // <copyright file="Ilutp.cs" company="Math.NET">
  // Math.NET Numerics, part of the Math.NET Project
  // http://numerics.mathdotnet.com
  // http://github.com/mathnet/mathnet-numerics
  //
  // Copyright (c) 2009-2010 Math.NET
  //
  // Permission is hereby granted, free of charge, to any person
  // obtaining a copy of this software and associated documentation
  // files (the "Software"), to deal in the Software without
  // restriction, including without limitation the rights to use,
  // copy, modify, merge, publish, distribute, sublicense, and/or sell
  // copies of the Software, and to permit persons to whom the
  // Software is furnished to do so, subject to the following
  // conditions:
  //
  // The above copyright notice and this permission notice shall be
  // included in all copies or substantial portions of the Software.
  //
  // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  // WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  // OTHER DEALINGS IN THE SOFTWARE.
  // </copyright>
  
  using System;
  using System.Collections.Generic;
  using Nequeo.Science.Math.LinearAlgebra.Solvers;
  using Nequeo.Science.Math.Properties;
  
  namespace Nequeo.Science.Math.LinearAlgebra.Complex32.Solvers
  {
      using Nequeo.Science.Math;
  
      /// <summary>
      /// This class performs an Incomplete LU factorization with drop tolerance
      /// and partial pivoting. The drop tolerance indicates which additional entries
      /// will be dropped from the factorized LU matrices.
      /// </summary>
      /// <remarks>
      /// The ILUTP-Mem algorithm was taken from: <br/>
      /// ILUTP_Mem: a Space-Efficient Incomplete LU Preconditioner
      /// <br/>
      /// Tzu-Yi Chen, Department of Mathematics and Computer Science, <br/>
      /// Pomona College, Claremont CA 91711, USA <br/>
      /// Published in: <br/>
      /// Lecture Notes in Computer Science <br/>
      /// Volume 3046 / 2004 <br/>
      /// pp. 20 - 28 <br/>
      /// Algorithm is described in Section 2, page 22
      /// </remarks>
      public sealed class ILUTPPreconditioner : IPreconditioner<Complex32>
      {
          /// <summary>
          /// The default fill level.
          /// </summary>
          public const double DefaultFillLevel = 200.0;
  
          /// <summary>
          /// The default drop tolerance.
          /// </summary>
          public const double DefaultDropTolerance = 0.0001;
  
          /// <summary>
          /// The decomposed upper triangular matrix.
          /// </summary>
          SparseMatrix _upper;
  
          /// <summary>
          /// The decomposed lower triangular matrix.
          /// </summary>
          SparseMatrix _lower;
  
          /// <summary>
          /// The array containing the pivot values.
          /// </summary>
          int[] _pivots;
  
          /// <summary>
          /// The fill level.
          /// </summary>
          double _fillLevel = DefaultFillLevel;
  
          /// <summary>
          /// The drop tolerance.
          /// </summary>
          double _dropTolerance = DefaultDropTolerance;
  
          /// <summary>
          /// The pivot tolerance.
          /// </summary>
          double _pivotTolerance;
  
          /// <summary>
          /// Initializes a new instance of the <see cref="ILUTPPreconditioner"/> class with the default settings.
          /// </summary>
          public ILUTPPreconditioner()
          {
          }
  
          /// <summary>
          /// Initializes a new instance of the <see cref="ILUTPPreconditioner"/> class with the specified settings.
          /// </summary>
          /// <param name="fillLevel">
          /// The amount of fill that is allowed in the matrix. The value is a fraction of
          /// the number of non-zero entries in the original matrix. Values should be positive.
          /// </param>
          /// <param name="dropTolerance">
          /// The absolute drop tolerance which indicates below what absolute value an entry
          /// will be dropped from the matrix. A drop tolerance of 0.0 means that no values
          /// will be dropped. Values should always be positive.
          /// </param>
          /// <param name="pivotTolerance">
          /// The pivot tolerance which indicates at what level pivoting will take place. A
          /// value of 0.0 means that no pivoting will take place.
          /// </param>
          public ILUTPPreconditioner(double fillLevel, double dropTolerance, double pivotTolerance)
          {
              if (fillLevel < 0)
              {
                  throw new ArgumentOutOfRangeException("fillLevel");
              }
  
              if (dropTolerance < 0)
              {
                  throw new ArgumentOutOfRangeException("dropTolerance");
              }
  
              if (pivotTolerance < 0)
              {
                  throw new ArgumentOutOfRangeException("pivotTolerance");
              }
  
              _fillLevel = fillLevel;
              _dropTolerance = dropTolerance;
              _pivotTolerance = pivotTolerance;
          }
  
          /// <summary>
          /// Gets or sets the amount of fill that is allowed in the matrix. The
          /// value is a fraction of the number of non-zero entries in the original
          /// matrix. The standard value is 200.
          /// </summary>
          /// <remarks>
          /// <para>
          /// Values should always be positive and can be higher than 1.0. A value lower
          /// than 1.0 means that the eventual preconditioner matrix will have fewer
          /// non-zero entries as the original matrix. A value higher than 1.0 means that
          /// the eventual preconditioner can have more non-zero values than the original
          /// matrix.
          /// </para>
          /// <para>
          /// Note that any changes to the <b>FillLevel</b> after creating the preconditioner
          /// will invalidate the created preconditioner and will require a re-initialization of
          /// the preconditioner.
          /// </para>
          /// </remarks>
          /// <exception cref="ArgumentOutOfRangeException">Thrown if a negative value is provided.</exception>
          public double FillLevel
          {
              get { return _fillLevel; }
              set
              {
                  if (value < 0)
                  {
                      throw new ArgumentOutOfRangeException("value");
                  }
  
                  _fillLevel = value;
              }
          }
  
          /// <summary>
          /// Gets or sets the absolute drop tolerance which indicates below what absolute value
          /// an entry will be dropped from the matrix. The standard value is 0.0001.
          /// </summary>
          /// <remarks>
          /// <para>
          /// The values should always be positive and can be larger than 1.0. A low value will
          /// keep more small numbers in the preconditioner matrix. A high value will remove
          /// more small numbers from the preconditioner matrix.
          /// </para>
          /// <para>
          /// Note that any changes to the <b>DropTolerance</b> after creating the preconditioner
          /// will invalidate the created preconditioner and will require a re-initialization of
          /// the preconditioner.
          /// </para>
          /// </remarks>
          /// <exception cref="ArgumentOutOfRangeException">Thrown if a negative value is provided.</exception>
          public double DropTolerance
          {
              get { return _dropTolerance; }
              set
              {
                  if (value < 0)
                  {
                      throw new ArgumentOutOfRangeException("value");
                  }
  
                  _dropTolerance = value;
              }
          }
  
          /// <summary>
          /// Gets or sets the pivot tolerance which indicates at what level pivoting will
          /// take place. The standard value is 0.0 which means pivoting will never take place.
          /// </summary>
          /// <remarks>
          /// <para>
          /// The pivot tolerance is used to calculate if pivoting is necessary. Pivoting
          /// will take place if any of the values in a row is bigger than the
          /// diagonal value of that row divided by the pivot tolerance, i.e. pivoting
          /// will take place if <b>row(i,j) > row(i,i) / PivotTolerance</b> for
          /// any <b>j</b> that is not equal to <b>i</b>.
          /// </para>
          /// <para>
          /// Note that any changes to the <b>PivotTolerance</b> after creating the preconditioner
          /// will invalidate the created preconditioner and will require a re-initialization of
          /// the preconditioner.
          /// </para>
          /// </remarks>
          /// <exception cref="ArgumentOutOfRangeException">Thrown if a negative value is provided.</exception>
          public double PivotTolerance
          {
              get { return _pivotTolerance; }
              set
              {
                  if (value < 0)
                  {
                      throw new ArgumentOutOfRangeException("value");
                  }
  
                  _pivotTolerance = value;
              }
          }
  
          /// <summary>
          /// Returns the upper triagonal matrix that was created during the LU decomposition.
          /// </summary>
          /// <remarks>
          /// This method is used for debugging purposes only and should normally not be used.
          /// </remarks>
          /// <returns>A new matrix containing the upper triagonal elements.</returns>
          internal Matrix<Complex32> UpperTriangle()
          {
              return _upper.Clone();
          }
  
          /// <summary>
          /// Returns the lower triagonal matrix that was created during the LU decomposition.
          /// </summary>
          /// <remarks>
          /// This method is used for debugging purposes only and should normally not be used.
          /// </remarks>
          /// <returns>A new matrix containing the lower triagonal elements.</returns>
          internal Matrix<Complex32> LowerTriangle()
          {
              return _lower.Clone();
          }
  
          /// <summary>
          /// Returns the pivot array. This array is not needed for normal use because
          /// the preconditioner will return the solution vector values in the proper order.
          /// </summary>
          /// <remarks>
          /// This method is used for debugging purposes only and should normally not be used.
          /// </remarks>
          /// <returns>The pivot array.</returns>
          internal int[] Pivots()
          {
              var result = new int[_pivots.Length];
              for (var i = 0; i < _pivots.Length; i++)
              {
                  result[i] = _pivots[i];
              }
  
              return result;
          }
  
          /// <summary>
          /// Initializes the preconditioner and loads the internal data structures.
          /// </summary>
          /// <param name="matrix">
          /// The <see cref="Matrix"/> upon which this preconditioner is based. Note that the
          /// method takes a general matrix type. However internally the data is stored
          /// as a sparse matrix. Therefore it is not recommended to pass a dense matrix.
          /// </param>
          /// <exception cref="ArgumentNullException"> If <paramref name="matrix"/> is <see langword="null" />.</exception>
          /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
          public void Initialize(Matrix<Complex32> matrix)
          {
              if (matrix == null)
              {
                  throw new ArgumentNullException("matrix");
              }
  
              if (matrix.RowCount != matrix.ColumnCount)
              {
                  throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix");
              }
  
              var sparseMatrix = matrix as SparseMatrix ?? SparseMatrix.OfMatrix(matrix);
  
              // The creation of the preconditioner follows the following algorithm.
              // spaceLeft = lfilNnz * nnz(A)
              // for i = 1, .. , n
              // {
              //    w = a(i,*)
              //    for j = 1, .. , i - 1
              //    {
              //        if (w(j) != 0)
              //        {
              //            w(j) = w(j) / a(j,j)
              //            if (w(j) < dropTol)
              //            {
              //                w(j) = 0;
              //            }
              //            if (w(j) != 0)
              //            {
              //                w = w - w(j) * U(j,*)
              //            }
              //        }
              //    }
              //
              //    for j = i, .. ,n
              //    {
              //        if w(j) <= dropTol * ||A(i,*)||
              //        {
              //            w(j) = 0
              //        }
              //    }
              //
              //    spaceRow = spaceLeft / (n - i + 1) // Determine the space for this row
              //    lfil = spaceRow / 2  // space for this row of L
              //    l(i,j) = w(j) for j = 1, .. , i -1 // only the largest lfil elements
              //
              //    lfil = spaceRow - nnz(L(i,:))  // space for this row of U
              //    u(i,j) = w(j) for j = i, .. , n // only the largest lfil - 1 elements
              //    w = 0
              //
              //    if max(U(i,i + 1: n)) > U(i,i) / pivTol then // pivot if necessary
              //    {
              //        pivot by swapping the max and the diagonal entries
              //        Update L, U
              //        Update P
              //    }
              //    spaceLeft = spaceLeft - nnz(L(i,:)) - nnz(U(i,:))
              // }
              // Create the lower triangular matrix
              _lower = new SparseMatrix(sparseMatrix.RowCount);
  
              // Create the upper triangular matrix and copy the values
              _upper = new SparseMatrix(sparseMatrix.RowCount);
  
              // Create the pivot array
              _pivots = new int[sparseMatrix.RowCount];
              for (var i = 0; i < _pivots.Length; i++)
              {
                  _pivots[i] = i;
              }
  
              var workVector = new DenseVector(sparseMatrix.RowCount);
              var rowVector = new DenseVector(sparseMatrix.ColumnCount);
              var indexSorting = new int[sparseMatrix.RowCount];
  
              // spaceLeft = lfilNnz * nnz(A)
              var spaceLeft = (int) _fillLevel*sparseMatrix.NonZerosCount;
  
              // for i = 1, .. , n
              for (var i = 0; i < sparseMatrix.RowCount; i++)
              {
                  // w = a(i,*)
                  sparseMatrix.Row(i, workVector);
  
                  // pivot the row
                  PivotRow(workVector);
                  var vectorNorm = workVector.InfinityNorm();
  
                  // for j = 1, .. , i - 1)
                  for (var j = 0; j < i; j++)
                  {
                      // if (w(j) != 0)
                      // {
                      //     w(j) = w(j) / a(j,j)
                      //     if (w(j) < dropTol)
                      //     {
                      //         w(j) = 0;
                      //     }
                      //     if (w(j) != 0)
                      //     {
                      //         w = w - w(j) * U(j,*)
                      //     }
                      if (workVector[j] != 0.0f)
                      {
                          // Calculate the multiplication factors that go into the L matrix
                          workVector[j] = workVector[j]/_upper[j, j];
                          if (workVector[j].Magnitude < _dropTolerance)
                          {
                              workVector[j] = 0.0f;
                          }
  
                          // Calculate the addition factor
                          if (workVector[j] != 0.0f)
                          {
                              // vector update all in one go
                              _upper.Row(j, rowVector);
  
                              // zero out columnVector[k] because we don't need that
                              // one anymore for k = 0 to k = j
                              for (var k = 0; k <= j; k++)
                              {
                                  rowVector[k] = 0.0f;
                              }
  
                              rowVector.Multiply(workVector[j], rowVector);
                              workVector.Subtract(rowVector, workVector);
                          }
                      }
                  }
  
                  // for j = i, .. ,n
                  for (var j = i; j < sparseMatrix.RowCount; j++)
                  {
                      // if w(j) <= dropTol * ||A(i,*)||
                      // {
                      //     w(j) = 0
                      // }
                      if (workVector[j].Magnitude <= _dropTolerance*vectorNorm)
                      {
                          workVector[j] = 0.0f;
                      }
                  }
  
                  // spaceRow = spaceLeft / (n - i + 1) // Determine the space for this row
                  var spaceRow = spaceLeft/(sparseMatrix.RowCount - i + 1);
  
                  // lfil = spaceRow / 2  // space for this row of L
  
gt:                 var fillLevel = spaceRow/2;