ilu.hpp

Go to the documentation of this file.

#ifndef PARGEMSLR_ILU_H
#define PARGEMSLR_ILU_H
 
#include <iostream>
 
#include "../utils/utils.hpp"
#include "../vectors/int_vector.hpp"
#include "../vectors/sequential_vector.hpp"
#include "../matrices/csr_matrix.hpp"
#include "../solvers/solver.hpp"
 
/* IMPORTANT NOTE:
 * this file contains some functions from the package HYPRE 
 * parcsr_ls/par_ilu_setup.c
 * parcsr_ls/par_ilu_solve.c
 * 
 * Modifications are made to fit the ParGEMSLR data structure
 * the following is the license statement of HYPRE
 * 
 * Include the following functions:
 * 
 * ----------------------------------------------------------------------------------------------------------------
 * 
 * Intellectual Property Notice
 * ----------------------------
 * HYPRE is licensed under the Apache License, Version 2.0 (see LICENSE-APACHE or
 * http://www.apache.org/licenses/LICENSE-2.0) or the MIT license (see LICENSE-MIT
 * or http://opensource.org/licenses/MIT), at your option.
 * 
 * Copyrights and patents in the HYPRE project are retained by contributors.
 * No copyright assignment is required to contribute to HYPRE.
 * 
 * HYPRE's original open source license was LGPL-2.1.  Consent from contributors
 * was received and documented before relicensing to Apache-2.0/MIT.
 * 
 * SPDX usage
 * ----------
 * 
 * Individual files contain SPDX tags instead of the full license text.
 * This enables machine processing of license information based on the SPDX
 * License Identifiers that are available here: https://spdx.org/licenses/
 * 
 * Files that are dual-licensed as Apache-2.0 OR MIT contain the following
 * text in the license header:
 * 
 *     SPDX-License-Identifier: (Apache-2.0 OR MIT)
 * 
 * 
 * External Software
 * -----------------
 * 
 * External software in hypre is covered by various permissive licenses.  A summary
 * listing follows.  See the license included with the software for full details.
 * 
 * Directory: src/blas
 * License: University of Tennessee
 * 
 * Directory: src/lapack
 * License: University of Tennessee
 * 
 * ----------------------------------------------------------------------------------------------------------------
 * 
 */
 
#ifdef PARGEMSLR_CUDA
#include "cusparse.h"
#endif
 
namespace pargemslr
{
   
   enum IluOptionEnum
   {
      kIluOptionILUT,
      kIluOptionPartialILUT,
      kIluOptionILUK
   };
   
   enum IluReorderingOptionEnum
   {
      kIluReorderingNo,
      kIluReorderingRcm,
      kIluReorderingAmd,
      kIluReorderingNd,
      kIluReorderingDdpq,
      kIluReorderingInput
   };
   
   enum IluOpenMPOptionEnum
   {
      kIluOpenMPNo,
      kIluOpenMPLevelScheduling,
      kIluOpenMPPoly
   };
   
   template <class MatrixType, class VectorType, typename DataType>
   class IluClass: public SolverClass<MatrixType, VectorType, DataType>
   {
   private:
      
      DataType                            _diag_shift_milu;
      
      VectorType                          _x_temp;
      
#ifdef PARGEMSLR_CUDA
 
      CsrMatrixClass<DataType>            _LDU;
      
      csrsv2Info_t                        _matL_info;
      
      csrsv2Info_t                        _matU_info;
      
#endif
      
      int                                 _n;
      
      int                                 _nnz;
      
      bool                                _complex_shift;
      
      CsrMatrixClass<DataType>            _L;
      
      SequentialVectorClass<DataType>     _D;
      
      CsrMatrixClass<DataType>            _U;
      
      int                                 _nB;
      
      CsrMatrixClass<DataType>            _E;
      
      CsrMatrixClass<DataType>            _F;
      
      CsrMatrixClass<DataType>            _S;
      
#ifdef PARGEMSLR_OPENMP
      
      VectorType                          _y_temp;
      
      VectorType                          _z_temp;
      
      CsrMatrixClass<DataType>            _L_poly;
      
      CsrMatrixClass<DataType>            _D_poly;
      
      CsrMatrixClass<DataType>            _U_poly;
      
      CsrMatrixClass<DataType>            _L_level;
      
      SequentialVectorClass<DataType>     _D_level;
      
      CsrMatrixClass<DataType>            _U_level;
      
      vector_int                          _level_ptr_l;
      
      vector_int                          _levels_l;
      
      vector_int                          _level_ptr_u;
      
      vector_int                          _levels_u;
      
      int                                 _levels_l_start;
      
      int                                 _levels_l_end;
      
      int                                 _levels_u_start;
      
      int                                 _levels_u_end;
      
#endif
      
      int                                 _fill_level;
      
      typename std::conditional<PargemslrIsDoublePrecision<DataType>::value, 
                                 double, 
                                 float>::type _droptol;
      
      typename std::conditional<PargemslrIsDoublePrecision<DataType>::value, 
                                 double, 
                                 float>::type _droptol_ef;
      
      typename std::conditional<PargemslrIsDoublePrecision<DataType>::value, 
                                 double, 
                                 float>::type _droptol_s;
      
      int                                 _max_row_nnz;
      
      int                                 _max_row_nnz_s;
      
      int                                 _option;
      
      int                                 _omp_option; 
      
      int                                 _poly_order;
      
      int                                 _perm_option;
      
      bool                                _modified;
      
      IntVectorClass<int>                 _row_perm_vec;
      
      IntVectorClass<int>                 _col_perm_vec;
      
      int                                 _location;
      
      int                                 _matrix_location;
      
      int                                 _cusparse_ready;
      
      template <typename T>
      int ApplyShift( T &w, int nz, T norm);
      
      template <typename T>
      int ApplyShift( ComplexValueClass<T> &w, int nz, T norm);
      
      template <typename T>
      int         Qsplit(T *a, int *ind, int n, int Ncut); 
 
      int SetupILUKSymbolic( int n, int *A_i, int *A_j, int lfil, int *perm, int *qperm, int *rqperm, int *iw, int *L_i, int *U_i, int **L_jp, int **U_jp);
      
      virtual int SolveHost( VectorType &x, VectorType &rhs);
 
      virtual int SolveLHost( VectorType &x, VectorType &rhs);
      
      virtual int SolveUHost( VectorType &x, VectorType &rhs);
      
#ifdef PARGEMSLR_CUDA
      
      virtual int SolveDevice( VectorType &x, VectorType &rhs);
      
      int SolveLDevice( VectorType &x, VectorType &rhs);
      
      int SolveUDevice( VectorType &x, VectorType &rhs);
      
#endif
 
#ifdef PARGEMSLR_OPENMP
 
      int         BuildLevels();
      
      int         BuildPoly();
      
      int         SolveHostOmp( VectorType &x, VectorType &rhs);
      
      int         SolveHostOmpPoly( VectorType &x, VectorType &rhs, int option);
      
#endif
 
    public:
 
      IluClass();
      
      IluClass(const IluClass<MatrixType, VectorType, DataType> &precond);
      
      IluClass(IluClass<MatrixType, VectorType, DataType> &&precond);
      
      IluClass<MatrixType, VectorType, DataType>& operator=(const IluClass<MatrixType, VectorType, DataType> &precond);
      
      IluClass<MatrixType, VectorType, DataType>& operator=(IluClass<MatrixType, VectorType, DataType> &&precond);
      
      virtual int Clear();
     
      virtual ~IluClass();
      
      virtual int Setup( VectorType &x, VectorType &rhs);
      
      int Setup(MatrixType &L_mat, VectorType &D_vec, MatrixType &U_mat, int *pperm, int *qperm);
      
      int SetupPermutation( VectorType &x, VectorType &rhs);
      
      int SetupILUT( VectorType &x, VectorType &rhs);
      
      int SetupILUK( VectorType &x, VectorType &rhs);
      
      template <typename T>
      static void Swap( T *v, int i, int j)
      {
         T temp;
         temp = v[i];
         v[i] = v[j];
         v[j] = temp;
      }
      
      template <typename T1, typename T2>
      static void Swap(T1 *v, T2 *w, int i, int j)
      {
         T1 temp;
         T2 temp2;
 
         temp = v[i];
         v[i] = v[j];
         v[j] = temp;
         temp2 = w[i];
         w[i] = w[j];
         w[j] = temp2;
      }
 
      template <typename T1, typename T2>
      static void AbsMinHeapAdd_RN(T1 *heap, T2 *v, int len)
      {
         /* parent, left, right */
         int p;
         len--;/* now len is the current index */
         while(len > 0)
         {
            /* get the parent index */
            p = (len-1)/2;
            if(PargemslrAbs(heap[-p]) < PargemslrAbs(heap[-len]))
            {
               /* this is smaller */
               IluClass<MatrixType, VectorType, DataType>::Swap<T2, T1>( v, heap, -p, -len);
               len = p;
            }
            else
            {
               break;
            }
         }
      }
      
      template <typename T1, typename T2>
      static void AbsMinHeapRemove_RN(T1 *heap, T2 *v, int len)
      {
         /* parent, left, right */
         int p,l,r;
         len--;/* now len is the max index */
         /* swap the first element to last */
         IluClass<MatrixType, VectorType, DataType>::Swap<T2, T1>( v, heap, 0, -len);
         p = 0;
         l = 1;
         /* while I'm still in the heap */
         while(l < len)
         {
            r = 2*p+2;
            /* two childs, pick the smaller one */
            l = r >= len || PargemslrAbs(heap[-l])>PargemslrAbs(heap[-r]) ? l : r;
            if(PargemslrAbs(heap[-l])>PargemslrAbs(heap[-p]))
            {
               IluClass<MatrixType, VectorType, DataType>::Swap<T2, T1>( v, heap, -l, -p);
               p = l;
               l = 2*p+1;
            }
            else
            {
               break;
            }
         }
      }
 
      template <typename T1, typename T2>
      static void MinHeapAdd_NNR(int *heap, T1 *v1, T2 *v2, int len)
      {
         /* parent, left, right */
         int p;
         len--;/* now len is the current index */
         while(len > 0)
         {
            /* get the parent index */
            p = (len-1)/2;
            if(heap[p] > heap[len])
            {
               /* this is smaller */
               IluClass<MatrixType, VectorType, DataType>::Swap<T2>( v2, heap[p], heap[len]);
               IluClass<MatrixType, VectorType, DataType>::Swap<int, T1>( heap, v1, p, len);
               len = p;
            }
            else
            {
               break;
            }
         }
      }
      
      template <typename T1, typename T2>
      static void MinHeapRemovd_NNR(int *heap, T1 *v1, T2 *v2, int len)
      {
         /* parent, left, right */
         int p,l,r;
         len--;/* now len is the max index */
         /* swap the first element to last */
         IluClass<MatrixType, VectorType, DataType>::Swap<T2>(v2, heap[0], heap[len]);
         IluClass<MatrixType, VectorType, DataType>::Swap<int, T1>(heap, v1, 0, len);
         p = 0;
         l = 1;
         /* while I'm still in the heap */
         while(l < len)
         {
            r = 2*p+2;
            /* two childs, pick the smaller one */
            l = r >= len || heap[l]<heap[r] ? l : r;
            if(heap[l]<heap[p])
            {
               IluClass<MatrixType, VectorType, DataType>::Swap<T2>(v2,heap[p],heap[l]);
               IluClass<MatrixType, VectorType, DataType>::Swap<int, T1>( heap, v1, l, p);
               p = l;
               l = 2*p+1;
            }
            else
            {
               break;
            }
         }
      }
      
      template <typename T>
      static void Qsplit(T *a, int *ind, int left, int bound, int right)
      {
         int i, last;
         if (left >= right)
         {
            return;
         }
         IluClass<MatrixType, VectorType, DataType>::Swap<int, T>( ind, a, left, (left+right)/2);
         last = left;
         for(i = left + 1 ; i <= right ; i ++)
         {
            if(PargemslrAbs(a[i]) > PargemslrAbs(a[left]))
            {
               IluClass<MatrixType, VectorType, DataType>::Swap<int, T>( ind, a, ++last, i);
            }
         }
         IluClass<MatrixType, VectorType, DataType>::Swap<int, T>( ind, a, left, last);
         IluClass<MatrixType, VectorType, DataType>::Qsplit<T>( a, ind, left, bound, last-1);
         if(bound > last)
         {
             IluClass<MatrixType, VectorType, DataType>::Qsplit<T>( a, ind, last+1, bound, right);
         }
      }
      
      int SetupPartialILUT(VectorType &x, VectorType &rhs);
      
      virtual int Solve( VectorType &x, VectorType &rhs);
      
      int SolveL( VectorType &x, VectorType &rhs);
      
      int SolveU( VectorType &x, VectorType &rhs);
      
      virtual int SetSolveLocation( const int &location);
      
      virtual int MoveData( const int &location);
      
      /* --------- SETS and GETS ----------- */
      
      virtual int SetWithParameterArray(double *params)
      {
         if(this->_ready)
         {
            PARGEMSLR_ERROR("Change setup after preconditioner is built.");
            return PARGEMSLR_ERROR_FUNCTION_CALL_ERR;
         }
         
         this->_droptol = params[PARGEMSLR_IO_ILU_DROPTOL_B_LOCAL];
         this->_droptol_ef = params[PARGEMSLR_IO_ILU_DROPTOL_EF_LOCAL];
         this->_droptol_s = params[PARGEMSLR_IO_ILU_DROPTOL_S_LOCAL];
         this->_max_row_nnz = params[PARGEMSLR_IO_ILU_ROWNNZ_B_LOCAL];
         this->_max_row_nnz_s = params[PARGEMSLR_IO_ILU_ROWNNZ_S_LOCAL];
         this->_perm_option = params[PARGEMSLR_IO_ILU_PERM_OPTION_LOCAL];
         this->_omp_option = params[PARGEMSLR_IO_ILU_OMP_OPTION_LOCAL];
         this->_poly_order = params[PARGEMSLR_IO_POLY_ORDER];
         this->_print_option = params[PARGEMSLR_IO_GENERAL_PRINT_LEVEL];
         this->_complex_shift = params[PARGEMSLR_IO_ADVANCED_USE_COMPLEX_SHIFT] != 0.0;
         this->_diag_shift_milu = params[PARGEMSLR_IO_ADVANCED_DIAG_SHIFT_MODIFIED];
         
         return PARGEMSLR_SUCCESS;
      }
      
      int         SetLevelOfFill( int lfil)
      {
         this->_fill_level = lfil;
         return PARGEMSLR_SUCCESS;
      }
      
      int         SetIluComplexShift( bool complex_shift)
      {
         this->_complex_shift = complex_shift;
         return PARGEMSLR_SUCCESS;
      }
      
      template <typename T>
      int         SetDropTolerance( T droptol)
      {
         this->_droptol = droptol;
         return PARGEMSLR_SUCCESS;
      }
      
      template <typename T>
      int         SetDropToleranceEF( T droptol_ef)
      {
         if(this->_ready)
         {
            PARGEMSLR_WARNING("Change setting after Setup is not going to change the preconditioner.");
         }
         this->_droptol_ef = droptol_ef;
         return PARGEMSLR_SUCCESS;
      }
      
      template <typename T>
      int         SetDropToleranceS( T droptol_s)
      {
         if(this->_ready)
         {
            PARGEMSLR_WARNING("Change setting after Setup is not going to change the preconditioner.");
         }
         this->_droptol_s = droptol_s;
         return PARGEMSLR_SUCCESS;
      }
      
      int         SetMaxNnzPerRow( int max_row_nnz);
      
      int         SetMaxNnzPerRowSPart( int max_row_nnz_s);
      
      int         SetNB(int nB);
      
      int         SetOption( int option);
      
      int         SetOpenMPOption( int omp_option);
      
      int         SetPolyOrder( int order);
      
      int         SetPermutationOption( int perm_option);
      
      int         SetModified( bool modified);
      
      template <typename T>
      int         SetModifiedShift( T milu_shift)
      {
         if(this->_ready)
         {
            PARGEMSLR_WARNING("Change setting after Setup is not going to change the preconditioner.");
         }
         this->_diag_shift_milu = milu_shift;
         return PARGEMSLR_SUCCESS;
      }
      
      int         GetSize();
      
      virtual long int  GetNumNonzeros();
      
      CsrMatrixClass<DataType>& GetL();
      
      SequentialVectorClass<DataType>& GetD();
      
      CsrMatrixClass<DataType>& GetU();
      
      CsrMatrixClass<DataType>& GetE();
      
      CsrMatrixClass<DataType>& GetF();
      
      CsrMatrixClass<DataType>& GetS();
      
      int GetNB();
      
      IntVectorClass<int>& GetRowPermutationVector();
      
      IntVectorClass<int>& GetColPermutationVector();
      
    };
   
   typedef IluClass<CsrMatrixClass<float>, SequentialVectorClass<float>, float>        precond_ilu_csr_seq_float;
   typedef IluClass<CsrMatrixClass<double>, SequentialVectorClass<double>, double>     precond_ilu_csr_seq_double;
   typedef IluClass<CsrMatrixClass<complexs>, SequentialVectorClass<complexs>, complexs>  precond_ilu_csr_seq_complexs;
   typedef IluClass<CsrMatrixClass<complexd>, SequentialVectorClass<complexd>, complexd> precond_ilu_csr_seq_complexd;
   
}
 
#endif