new linear solver k3biilu2 is added

examples/MatSolve/main.cpp

@@ -15,7 +15,7 @@ int main(int argc, char ** argv)
 	int rank,procs;
 	if( argc < 3 )
 	{
-		std::cout << "Usage: " << argv[0] << " method_number<0:INNER_ILU2,1:INNER_MLILUC,2:PETSc,3:Trilinos_Aztec,4:Trilinos_Belos,5:Trilinos_Ifpack,6:Trilinos_ML,7:ANI> matrix.mtx [right_hand_side.rhs] [solver_options.txt]" << std::endl;
+		std::cout << "Usage: " << argv[0] << " method_number<0:INNER_ILU2,1:INNER_DDPQILUC,2:INNER_MPTILUC,3:INNER_MPTILU2,4:Trilinos_Aztec,5:Trilinos_Belos,6:Trilinos_ML,7:Trilinos_Ifpack,8:PETSc,9:ANI,10:FCBIILU2,11:K3BIILU2> matrix.mtx [right_hand_side.rhs] [solver_options.txt]" << std::endl;
 		return -1;
 	}
 	Solver::Type type;
@@ -23,14 +23,16 @@ int main(int argc, char ** argv)
 	{
 		case  0: type = Solver::INNER_ILU2;      break;
 		case  1: type = Solver::INNER_DDPQILUC;  break;
-		case 2: type = Solver::PETSc; break;
-		case 3: type = Solver::Trilinos_Aztec; break;
-		case 4: type = Solver::Trilinos_Belos; break;
-		case 5: type = Solver::Trilinos_Ifpack; break;
+		case  2: type = Solver::INNER_MPTILUC;   break;
+		case  3: type = Solver::INNER_MPTILU2;   break;
+		case  4: type = Solver::Trilinos_Aztec;  break;
+		case  5: type = Solver::Trilinos_Belos;  break;
 		case  6: type = Solver::Trilinos_ML;     break;
-    case 7: type = Solver::ANI; break;
-		case 8: type = Solver::INNER_MPTILUC; break;
-		case 9: type = Solver::INNER_MPTILU2; break;
+		case  7: type = Solver::Trilinos_Ifpack; break;
+		case  8: type = Solver::PETSc;           break;
+		case  9: type = Solver::ANI;             break;
+		case 10: type = Solver::FCBIILU2;        break;
+		case 11: type = Solver::K3BIILU2;        break;
 	}
 	Solver::Initialize(&argc,&argv,argc > 4 ? argv[4] : NULL); // Initialize the linear solver in accordance with args
 	{

inmost_solver.h

@@ -51,8 +51,8 @@ namespace INMOST
 			Trilinos_Ifpack,///< external Solver AztecOO with Ifpack preconditioner.
 			PETSc,          ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/
 			ANI,            ///< external Solver from ANI3D based on ILU2 (sequential Fortran version), @see http://ani3d.sourceforge.net/
-			FCBIILU2,       //   external FCBIILU2 Solver (BIILU2 parallel F2C version).
-			K3BIILU2        //   inner    K3BIILU2 Solver (BIILU2 parallel version).
+			FCBIILU2,       ///< external FCBIILU2 Solver (BIILU2 parallel F2C version).
+			K3BIILU2        ///< inner    K3BIILU2 Solver (BIILU2 parallel version).
 		};
 
 		static INMOST_MPI_Type & GetRowEntryType() {return RowEntryType;}

k3d.h

+//------------------------------------------------------------------------------------------------
+// File: k3d.h
+//------------------------------------------------------------------------------------------------
+
+#ifndef __K3D_H
+#define __K3D_H
+
+#include <iostream>
+#include <fstream>
+//#include <cstring>
+//#include <string>
+#include <cmath>
+#include <vector>
+
+using namespace std;
+
+namespace k3d {
+
+struct SSolverParams // Class that contains solver parameters
+{
+   int prec_float; ///< The type of preconditioner floating data to be used: prec_float == 1 (default) - use float type, prec_float == 2 - use double type
+   int ncycle;     ///< The number of sparsity extension cycles when determining backward overlap for blocks, 0 < ncycle, default value ncycle = 1
+   int ordtype;    ///< The ordering type when performing extended factorization: ordtype == 0 - no ordering, ordtype = 1 (default) - RCM ordering for initial, main and separator parts
+   int collap;     ///< The parameter that controls print sparsity condensing: collap == -1 (defalut) - no print, 0 < collap - print sparsity as .ps files for each block with condensing into collap times
+   int sctype;     ///< The type of scaling: sctype == -1 - no scaling, sctype == 0 (default) - usual scaling via diagonal value, sctype == 1 - balancing scaling is performed via rows/colums norms computations
+   int nitersc;    ///< The number of balancing scaling iterations, 0 < nitersc, default value nitersc = 3
+   int fcttype;    ///< The structural factorization type parameter, fcttype == -1 - no structural restrictions onto sparsity (default), 0 <= fcttype - elements of order fcttype at most are remaind in fct
+   double pivmin;  ///< The minimal pivot value, -1 < pivmin < 1.0, for negative value pivoting is ignored, default value pivmin = -1.0e0
+   double tau1;    ///< First  threshold value, 0 <= \tau1 < 1.0, default value \tau1 = 0.01e0
+   double tau2;    ///< Second threshold value, 0 <= \tau2 <= \tau1, default value \tau2 = 0.0001e0
+   double theta;   ///< Diagonal correction parameter, 0 <= \theta <= 1.0, default value \theta = 0.1e0
+// Constructors and destructor
+      /// Epmty constructor
+   SSolverParams () {
+      prec_float = 1; ncycle = 1; ordtype = 1; collap = -1; sctype = 0; nitersc = 3; 
+      fcttype = -1; pivmin = -1.0e0; tau1 = 1.0e-2; tau2 = 1.0e-4; theta = 0.1e0;
+   }; 
+      /// Set defaluts
+   void SetDefaults () {
+      prec_float = 1; ncycle = 1; ordtype = 1; collap = -1; sctype = 0; nitersc = 3; 
+      fcttype = -1; pivmin = -1.0e0; tau1 = 1.0e-2; tau2 = 1.0e-4; theta = 0.1e0;
+   };
+};
+
+   /// Perform binary search
+int BinarySearch (long long _isup, int _nblks, long long *_blks, int _iblkprev);
+
+   /// Print array
+void PrintArray (ostream &_fout, const char *_name, int _isize, int *_iarr);
+void PrintArray (ostream &_fout, const char *_name, int _isize, long long *_iarr);
+void PrintArray (ostream &_fout, const char *_name, int _isize, float *_farr);
+void PrintArray (ostream &_fout, const char *_name, int _isize, double *_farr);
+
+ostream &SetPw (ostream &stream); ///< Set ouput manip
+void Round (double &dx, double &dy); ///< Round values
+
+class CExchange
+{
+// Data
+public:
+// Control functions
+   static void Synchronize (void *_comm); /// Barrier synchronization of the CPUs
+   static int GetNproc (void *_comm);     /// Get the number of CPUs
+   static int GetMyid (void *_comm);      /// Get the ID number of the CPU
+// Exchange functions
+      /// Collective array operation
+   static void ExchangeArray (void *_comm, char _Datatype, char _Operation,
+                              int _Length, void *_Array);
+      /// Collective array operation (implementation)
+   static void ExchangeArray_impl (void *_comm, char _Datatype, char _Operation,
+                              int _Length, void *_Array);
+      /// Major data exchange function
+   static void DataExchange (void *_comm, vector<int> &_CpuIDSend, vector<vector<char> > &_ObjSend,
+                              vector<int> &_CpuIDRecv, vector<vector<char> > &_ObjRecv);
+      /// Major data exchange function (implementation)
+   static void DataExchange_impl (void *_comm, vector<int> &_CpuIDSend, vector<vector<char> > &_ObjSend,
+                                    vector<int> &_CpuIDRecv, vector<vector<char> > &_ObjRecv);
+      /// Data exchange function (implementation)
+   static void DataExchange_impl2 (void *_comm,
+                                    int _NObjSend, int *_CpuIDSend, int *_ObjSizeSend, char **_ObjSend,
+                                    int &_NObjRecv, int *&_CpuIDRecv, int *&_ObjSizeRecv, char **&_ObjRecv);
+      /// Exchange data requests
+   static void DataExchangeRequest (void *_comm,
+                                    int _NObjSend, int *_CpuIDSend, int *_ObjSizeSend,
+                                    int &_NObjRecv, int *&_CpuIDRecv, int *&_ObjSizeRecv);
+// Send/recv/test/wait/probe functions
+      /// Asynchronous send data to the other CPU
+   static void ISend (void *_comm, int _rank, int _msgtag, int _length, char *_buffer, int _indrecv, void *_recv_arr);
+      /// Asynchronous receive data from the other CPU
+   static void IRecv (void *_comm, int _rank, int _msgtag, int _length, char *_buffer, int _indrecv, void *_recv_arr);
+      /// Wait for completion of exchanges
+   static void WaitAll  (int _count, void *_recvarr, void *_statarr);
+      /// Allocate send/receive requests
+   static void AllocateRecvs (int _nrecv, void *&_recvarr);
+      /// Allocate statuses
+   static void AllocateStats (int _nstat, void *&_statarr);
+      /// Delete send/receive requests
+   static void DeleteRecvs (void *_recvarr);
+      /// Delete statuses
+   static void DeleteStats (void *_statarr);
+};
+
+class CSortInt /// Support sorting of int data
+{
+public:
+// Data
+   int ival;  // ivalue contains integer value to be sorted
+   int i2val; // i2val contains secondary integer value to be sorted
+// Functions
+// Constructors and destructor
+   CSortInt () {}; // Default constructor
+   CSortInt (int _ival, int _i2val) { // Memory allocation constructor with zero data
+      ival = _ival;
+      i2val = _i2val;
+   };
+   ~CSortInt () {}; // Destructor
+// Operator functions
+   CSortInt &operator= (const CSortInt &_ii2) { // Equality operator
+      ival = _ii2.ival;
+      i2val = _ii2.i2val;
+      return *this;
+   };
+// Comparison function
+   inline bool operator< (const CSortInt& _obj) const {
+      return ival < _obj.ival;
+   };
+};
+
+class CSortInt2 /// Support sorting of double int data
+{
+public:
+// Data
+   int ixval; // ixval contains x integer value to be sorted
+   int iyval; // iyval contains y integer value to be sorted
+   int itail; // itail contains secondary integer value to be sorted
+// Functions
+// Constructors and destructor
+   CSortInt2 () {}; // Default constructor
+   CSortInt2 (int _ixval, int _iyval, int _itail) { // Memory allocation constructor with zero data
+      ixval = _ixval;
+      iyval = _iyval;
+      itail = _itail;
+   };
+   ~CSortInt2 () {}; // Destructor
+// Operator functions
+   CSortInt2 &operator= (const CSortInt2 &_ii2) { // Equality operator
+      ixval = _ii2.ixval;
+      iyval = _ii2.iyval;
+      itail = _ii2.itail;
+      return *this;
+   };
+// Comparison function
+   inline bool operator< (const CSortInt2& _obj) const {
+      if(ixval < _obj.ixval)
+         return true;
+      else if(ixval > _obj.ixval)
+         return false;
+      return iyval < _obj.iyval;
+   };
+};
+
+template <typename _FltVect>
+class CVect /// Class that supports vector operations
+{
+public:
+   static void SetByZeroes (int _n, _FltVect *_x); ///< Set vector data by zeroes
+   static _FltVect ScProd (int _n, _FltVect *_x, _FltVect *_y); ///< Compute scalar product
+   static void CopyVector (int _n, _FltVect *_x, _FltVect *_y); ///< Copy array
+   static void AddReplaceVector (int _n, _FltVect *_x1, _FltVect *_x1plus2); ///< Add vector data
+   static void SubtractReplaceVector (int _n, _FltVect *_x1, _FltVect *_x2minus1); ///< Subtract vector data
+   static void UpdateVector (int _n, _FltVect *_value, _FltVect *_arr_x, _FltVect *_arr_y); ///< Update array (axpy): y=a*x+y
+   static void UpdateVectorMinus (int _n, _FltVect *_value, _FltVect *_arr_x, _FltVect *_arr_y); ///< Update array for minus alpha (axpy)
+   static void UpdateVectorReversed (int _n, _FltVect *_value, _FltVect *_arr_x, _FltVect *_arr_y); ///< Update array reversed (aypx)
+   static void OrderVector (int _n, const vector<int> &_order, 
+                              const _FltVect *_x, _FltVect *_x_ord); ///< Order vector data
+   static void InvOrderVector (int _n, const vector<int> &_order, 
+                              const _FltVect *_x, _FltVect *_x_ord); ///< Inverse order vector data
+};
+
+template <typename _Int, typename _Flt>
+class CIlu2_impl // Class that supports factorization computations
+{
+public:
+      /// Compute optimal ordering for the matrix
+   static void ComputeOptimalOrder (int _ordtype,
+                                    const int _n, const vector<_Int> &_ia_alu, const vector<_Int> &_ja_alu,
+                                    vector<int> &_order);
+      /// Compute optimal ordering for the matrix via splitting into 3 blocks: 1-st order Schur complement for third block, which is the separator between first block and remaining data
+   static void ComputeOptimalOrderSchur (int _ordtype, 
+                                          const int _n, int _n1, const vector<_Int> &_ia_alu, const vector<_Int> &_ja_alu,
+                                          int &_n2, vector<int> &_order);
+      /// Compute ordered matrix (sparsity only)
+   static void ReorderMatrixSp (int _n, const vector<int> &_order,
+                              vector<_Int> &_ia_alu, vector<_Int> &_ja_alu,
+                              vector<_Int> &_ia_alu_ord, vector<_Int> &_ja_alu_ord);
+      /// Compute ordered matrix
+   static void ReorderMatrix (int _n, const vector<int> &_order,
+                              vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu,
+                              vector<_Int> &_ia_alu_ord, vector<_Int> &_ja_alu_ord, vector<_Flt> &_a_alu_ord);
+      /// Perform ILU2 point factorization of the block with future diagonal modification
+   static void Ilu2BlockTransform (int _sctype, int _nitersc, int _fcttype, double _pivmin, double _tau1, double _tau2, double _theta,
+                                    int _n, 
+                                    vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu,
+                                    vector<_Int> &_ia_lu, vector<_Int> &_ida_lu, vector<_Int> &_ja_lu, vector<_Flt> &_a_lu,
+                                    double &_eigmin_att, double &_eigmax_att);
+      /// Perform ILU2 point factorization of the block with future diagonal modification
+   static void Ilu2Block (int _sctype, int _nitersc, int _fcttype, double _pivmin, double _tau1, double _tau2, double _theta,
+                           int _n, 
+                           vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu,
+                           vector<_Int> &_ia_l, vector<_Int> &_ja_l, vector<_Flt> &_a_l,
+                           vector<_Int> &_ia_u, vector<_Int> &_ja_u, vector<_Flt> &_a_u,
+                           double &_eigmin_att, double &_eigmax_att);
+      /// Compute scaling
+   static void ComputeScaling (int _sctype, int _nitersc,
+                                 int _n, 
+                                 vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu,
+                                 vector<_Flt> &_sclL, vector<_Flt> &_sclU);
+      /// Perform explicit scaling
+   static void MatrixScale (int _n, vector<_Flt> &_sclL, vector<_Flt> &_sclU,
+                              vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu);
+      /// Symmetrize sparsity
+   static void SymmetrizeSparsity (int _n, const vector<_Int> &_ia_alu, const vector<_Int> &_ja_alu,
+                                    vector<_Int> &_ia_symm, vector<_Int> &_ja_symm);
+      /// Split matrix data into L and U parts for sparsity only
+   static void SplitLUSp (int _n, 
+                           const vector<_Int> &_ia_alu, const vector<_Int> &_ja_alu,
+                           vector<_Int> &_ia_l, vector<_Int> &_ja_l,
+                           vector<_Int> &_ia_u, vector<_Int> &_ja_u);
+      /// Split matrix data into L and U parts
+   static void SplitLU (int _n, 
+                        vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu,
+                        vector<_Int> &_ia_l, vector<_Int> &_ja_l, vector<_Flt> &_a_l,
+                        vector<_Int> &_ia_u, vector<_Int> &_ja_u, vector<_Flt> &_a_u);
+      /// Transpose square matrix for sparsity only
+   static void TransposeSp (int _n, 
+                           vector<_Int> &_ia_a, vector<_Int> &_ja_a,
+                           vector<_Int> &_ia_at, vector<_Int> &_ja_at);
+      /// Transpose square matrix
+   static void Transpose (int _n, 
+                           vector<_Int> &_ia_a, vector<_Int> &_ja_a, vector<_Flt> &_a_a,
+                           vector<_Int> &_ia_at, vector<_Int> &_ja_at, vector<_Flt> &_a_at);
+      /// Combine sparsity of L and U
+   static void CombineLUSp (int _n, 
+                              vector<_Int> &_ia_al, vector<_Int> &_ja_al,
+                              vector<_Int> &_ia_au, vector<_Int> &_ja_au,
+                              vector<_Int> &_ia_alu, vector<_Int> &_ja_alu);
+      /// Combine L and U data into extended pairs
+   static void CombinePairs (int _n, 
+                              vector<_Int> &_ia_al, vector<_Int> &_ja_al, vector<_Flt> &_a_al,
+                              vector<_Int> &_ia_au, vector<_Int> &_ja_au, vector<_Flt> &_a_au,
+                              vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu);
+      /// Split pairs fct data into L and U parts with post filtering
+   static void SplitPairsFilter (double _tau1, int _n, 
+                                 vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu,
+                                 vector<_Int> &_ia_l, vector<_Int> &_ja_l, vector<_Flt> &_a_l,
+                                 vector<_Int> &_ia_u, vector<_Int> &_ja_u, vector<_Flt> &_a_u);
+      /// Combine L and U data into extended pairs
+   static void CombineRowsLU (int _n, 
+                              vector<_Int> &_ia_al, vector<_Int> &_ja_al, vector<_Flt> &_a_al,
+                              vector<_Int> &_ia_au, vector<_Int> &_ja_au, vector<_Flt> &_a_au,
+                              vector<_Int> &_ia_alu, vector<_Int> &_ida_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu);
+      /// Perform ILU2 point factorization of the block with future diagonal modification (no structural control)
+   static void Ilu2BlockIlu2 (double _pivmin, double _tau1, double _tau2, double _theta,
+                              int _n, int _n_ini, 
+                              vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<_Flt> &_a_alu,
+                              vector<_Int> &_ia_lu, vector<_Int> &_ja_lu, vector<_Flt> &_a_lu,
+                              double &_eigmin_att, double &_eigmax_att);
+      /// Perform ILU2 point factorization of the block with future diagonal modification (with structural control)
+   static void Ilu2BlockIlu2 (int _fcttype, double _pivmin, double _tau1, double _tau2, double _theta,
+                              int _n, int _n_ini, 
+                              vector<_Int> &_ia_alu, vector<_Int> &_ja_alu, vector<char> &_ja_char_alu, vector<_Flt> &_a_alu,
+                              vector<_Int> &_ia_lu, vector<_Int> &_ja_lu, vector<char> &_ja_char_lu, vector<_Flt> &_a_lu,
+                              double &_eigmin_att, double &_eigmax_att);
+      /// Compute inverse scaling
+   static void InverseDiag (int _n, vector<_Flt> &_sclU, vector<_Flt> &_invsclU);
+      /// Rescale factor back
+   static void RescaleU (int _n, vector<_Flt> &_sclU, vector<_Flt> &_invsclU,
+                          vector<_Int> &_ia_u, vector<_Int> &_ja_u, vector<_Flt> &_a_u);
+      /// Balance diagonal
+   static void BalanceDiag (int _n, vector<_Int> &_ia_a, vector<_Int> &_ja_a, vector<_Flt> &_a_a,
+                              vector<_Int> &_ia_l, vector<_Int> &_ja_l, vector<_Flt> &_a_l,
+                              vector<_Int> &_ia_u, vector<_Int> &_ja_u, vector<_Flt> &_a_u,
+                              double &_diacorr_min, double &_diacorr_max);
+};
+
+template <typename _Int, typename _Flt>
+class CMatrix // Class that supports matrix data
+{
+   int n_list;               ///< The number of elements in the List and Ia arrays
+   int n_list2;              ///< The number of elements in the List2 array
+   int nz_ja;                ///< The number of elements in Ja array
+   int nz_ja2;               ///< The number of elements in Ja2 array
+   int nz_a;                 ///< The number of elements in A array
+   vector<_Int> list_matr;   ///< List part of a matrix
+   vector<_Int> list2_matr;  ///< List2 part of a matrix
+   vector<_Int> ia_matr;     ///< Ia part of a matrix
+   vector<_Int> ja_matr;     ///< Ja part of a matrix
+   vector<_Int> ja2_matr;    ///< Ja2 part of a matrix
+   vector<_Flt> a_matr;      ///< A  part of a matrix
+// Get template names
+   typedef _Int   TInt_type; ///< Reference to int type
+   typedef _Flt TFloat_type; ///< Reference to float type
+// Constructors and destructor
+public:
+      /// Epmty constructor
+   CMatrix () { n_list = 0; list_matr.resize(1); list_matr[0] = 0; list2_matr.resize(1); list2_matr[0] = 0; 
+                  ia_matr.resize(1); ia_matr[0] = 0; ja_matr.resize(1); ja_matr[0] = 0; ja2_matr.resize(1); ja2_matr[0] = 0; 
+                  a_matr.resize(1); a_matr[0] = (_Flt)0.0e0; n_list2 = 0; nz_ja = 0; nz_ja2 = 0; nz_a = 0; }; 
+      /// Copy constructor
+   CMatrix (const CMatrix<_Int,_Flt> &_aa); 
+      /// Equality operator
+   CMatrix<_Int,_Flt> &operator= (const CMatrix<_Int,_Flt> &_aa); 
+      /// Init by data
+   CMatrix (int _n, _Int *_ia, _Int *_ja) {
+      n_list = _n; n_list2 = 0; int nzja = (int)_ia[_n]; nz_ja = nzja; nz_ja2 = 0; nz_a = 0; 
+      list_matr.resize(1); list_matr[0] = 0; list2_matr.resize(1); list2_matr[0] = 0; 
+      ia_matr.resize(_n+1); ja_matr.resize(nzja+1); _Int *pia = &ia_matr[0]; _Int *pja = &ja_matr[0];
+      int i; for (i=0;i<=_n;i++) pia[i] = _ia[i]; for (i=0;i<nzja;i++) pja[i] = _ja[i];
+      ja2_matr.resize(1); ja2_matr[0] = 0; this->SetIdentityList();
+   };
+      /// Init by data
+   CMatrix (int _nlist, _Int *_list, _Int *_ia, _Int *_ja) {
+      n_list = _nlist; n_list2 = 0; int nzja = (int)_ia[_nlist]; nz_ja = nzja; nz_ja2 = 0; nz_a = 0; 
+      list_matr.resize(_nlist+1); list2_matr.resize(1); list2_matr[0] = 0; 
+      ia_matr.resize(_nlist+1); ja_matr.resize(nzja+1);
+      _Int *plist = &list_matr[0]; _Int *pia = &ia_matr[0]; _Int *pja = &ja_matr[0];
+      int i; for (i=0;i<_nlist;i++) plist[i] = _list[i]; for (i=0;i<=_nlist;i++) pia[i] = _ia[i]; for (i=0;i<nzja;i++) pja[i] = _ja[i];
+      ja2_matr.resize(1); ja2_matr[0] = 0; a_matr.resize(1); a_matr[0] = (_Flt)0.0e0; 
+   };
+      /// Init by data
+   CMatrix (int _nlist, _Int *_list, _Int *_ia, _Int *_ja, _Flt *_a) {
+      n_list = _nlist; n_list2 = 0; int nzja = (int)_ia[_nlist]; nz_ja = nzja; nz_ja2 = 0; nz_a = 0; 
+      list_matr.resize(_nlist+1); list2_matr.resize(1); list2_matr[0] = 0;
+      ia_matr.resize(_nlist+1); ja_matr.resize(nzja+1); ja2_matr.resize(1); ja2_matr[0] = 0; a_matr.resize(nzja+1);
+      _Int *plist = &list_matr[0]; _Int *pia = &ia_matr[0]; _Int *pja = &ja_matr[0]; _Flt *pa = &a_matr[0];
+      int i; for (i=0;i<_nlist;i++) plist[i] = _list[i]; for (i=0;i<=_nlist;i++) pia[i] = _ia[i]; 
+      for (i=0;i<nzja;i++) pja[i] = _ja[i]; for (i=0;i<nzja;i++) pa[i] = _a[i];
+   };
+      /// Init by data
+   CMatrix (int _n, _Int *_ia, _Int *_ja, _Flt *_a) {
+      n_list = _n; n_list2 = 0; int nzja = (int)_ia[_n]; nz_ja = nzja; nz_ja2 = 0; nz_a = nzja; 
+      list_matr.resize(1); list_matr[0] = 0; list2_matr.resize(1); list2_matr[0] = 0; 
+      ia_matr.resize(_n+1); ja_matr.resize(nzja+1); a_matr.resize(nzja+1);
+      _Int *pia = &ia_matr[0]; _Int *pja = &ja_matr[0]; _Flt *pa = &a_matr[0];
+      int i; for (i=0;i<=_n;i++) pia[i] = _ia[i]; for (i=0;i<nzja;i++) pja[i] = _ja[i]; for (i=0;i<nzja;i++) pa[i] = _a[i];
+      this->SetIdentityList(); ja2_matr.resize(1); ja2_matr[0] = 0; 
+   };
+      /// Destructor
+   ~CMatrix () {};
+public:
+// External functions
+// Get/set functions
+      /// Get N
+   int GetN () {return n_list;}; 
+   const int GetN () const {return n_list;}; 
+      /// Get Nlist
+   int GetNlist () {return n_list;}; 
+   const int GetNlist () const {return n_list;}; 
+      /// Get Nlist2
+   int GetNlist2 () {return n_list2;}; 
+   const int GetNlist2 () const {return n_list2;}; 
+      /// Get Nzja
+   int GetNzja () {return nz_ja;}; 
+   const int GetNzja () const {return nz_ja;}; 
+      /// Get Nzja2
+   int GetNzja2 () {return nz_ja2;}; 
+   const int GetNzja2 () const {return nz_ja2;}; 
+      /// Get Nza
+   int GetNza () {return nz_a;}; 
+   const int GetNza () const {return nz_a;}; 
+      /// Get List
+   vector<_Int> * GetList () {return &list_matr;}; 
+   const vector<_Int> * GetList () const {return &list_matr;}; 
+      /// Get List array
+   _Int * GetListArr () {return &list_matr[0];}; 
+   const _Int * GetListArr () const {return &list_matr[0];}; 
+      /// Get List2
+   vector<_Int> * GetList2 () {return &list2_matr;}; 
+   const vector<_Int> * GetList2 () const {return &list2_matr;}; 
+      /// Get List2 array
+   _Int * GetList2Arr () {return &list2_matr[0];}; 
+   const _Int * GetList2Arr () const {return &list2_matr[0];}; 
+      /// Get Ia
+   vector<_Int> * GetIa () {return &ia_matr;}; 
+   const vector<_Int> * GetIa () const {return &ia_matr;}; 
+      /// Get Ia array
+   _Int * GetIaArr () {return &ia_matr[0];}; 
+   const _Int * GetIaArr () const {return &ia_matr[0];}; 
+      /// Get Ja
+   vector<_Int> * GetJa () {return &ja_matr;}; 
+   const vector<_Int> * GetJa () const {return &ja_matr;}; 
+      /// Get Ja array
+   _Int * GetJaArr () {return &ja_matr[0];}; 
+   const _Int * GetJaArr () const {return &ja_matr[0];}; 
+      /// Get Ja2
+   vector<_Int> * GetJa2 () {return &ja2_matr;}; 
+   const vector<_Int> * GetJa2 () const {return &ja2_matr;}; 
+      /// Get Ja2 array
+   _Int * GetJa2Arr () {return &ja2_matr[0];}; 
+   const _Int * GetJa2Arr () const {return &ja2_matr[0];}; 
+      /// Get A
+   vector<_Flt> * GetA () {return &a_matr;}; 
+   const vector<_Flt> * GetA () const {return &a_matr;}; 
+      /// Get A array
+   _Flt * GetAArr () {return &a_matr[0];}; 
+   const _Flt * GetAArr () const {return &a_matr[0];}; 
+      /// Set N
+   void SetN (int _n) {n_list = _n;}; 
+      /// Set Nlist
+   void SetNlist (int _n) {n_list = _n;}; 
+      /// Set Nzja
+   void SetNzja (int _nzja) {nz_ja = _nzja;}; 
+      /// Set Nza
+   void SetNza (int _nza) {nz_a = _nza;}; 
+// Other
+      /// Resize sp function
+   void ResizeAndSetAllSp (int _nlist, int _nlist2, int _nzja, int _nzja2) {
+      this->n_list = _nlist; this->n_list2 = _nlist2; this->nz_ja = _nzja; this->nz_ja2 = _nzja2; this->nz_a = 0; 
+      this->list_matr.resize (_nlist+1); this->list2_matr.resize (_nlist2+1); this->ia_matr.resize (_nlist+1); this->ia_matr[0] = 0;
+      this->ja_matr.resize (_nzja+1); this->ja2_matr.resize (_nzja2+1); this->a_matr.resize (1);
+   };
+      /// Resize function
+   void ResizeAndSetAll (int _nlist, int _nlist2, int _nzja, int _nzja2, int _nza) {
+      this->n_list = _nlist; this->n_list2 = _nlist2; this->nz_ja = _nzja; this->nz_ja2 = _nzja2; this->nz_a = _nza; 
+      this->list_matr.resize (_nlist+1); this->list2_matr.resize (_nlist2+1); this->ia_matr.resize (_nlist+1); this->ia_matr[0] = 0;
+      this->ja_matr.resize (_nzja+1); this->ja2_matr.resize (_nzja2+1); this->a_matr.resize (_nza+1);
+   };
+      /// Replace function
+   void ReplaceFree (CMatrix<_Int,_Flt> &_aa) {
+      this->n_list = _aa.n_list; this->n_list2 = _aa.n_list2; this->nz_ja = _aa.nz_ja; this->nz_ja2 = _aa.nz_ja2; this->nz_a = _aa.nz_a; 
+      this->list_matr.swap (_aa.list_matr); this->list2_matr.swap (_aa.list2_matr); this->ia_matr.swap (_aa.ia_matr); 
+      this->ja_matr.swap (_aa.ja_matr); this->ja2_matr.swap (_aa.ja2_matr); this->a_matr.swap (_aa.a_matr); _aa.Clean();
+   };
+      /// Clean matrix data
+   void Clean () {
+      vector<_Int> list_dummy(1); vector<_Int> list2_dummy(1); vector<_Int> ia_dummy(1); 
+      vector<_Int> ja_dummy(1); vector<_Int> ja2_dummy(1); vector<_Flt> a_dummy(1);
+      list_matr.swap (list_dummy); list2_matr.swap (list2_dummy); ia_matr.swap (ia_dummy); 
+      ja_matr.swap (ja_dummy); ja2_matr.swap (ja2_dummy); a_matr.swap (a_dummy); 
+      n_list = 0; n_list2 = 0; nz_ja = 0; nz_ja2 = 0; nz_a = 0; ia_matr[0] = 0;
+   };
+      /// Set identity list
+   void SetIdentityList () {list_matr.resize(n_list+1); int i;  for (i=0;i<n_list;i++) list_matr[i] = (_Int)i;}; 
+      /// Get sparsity
+   void GetSparsity (const CMatrix<_Int,_Flt> &_a_sp);
+      /// Compute transposed sparsity for incomplete list
+   void TransposedSparsityListSp (int &_icycle, int *_imask, int *_indarr, int *_iptr, int *_listloc, int *_ialoc,
+                                    CMatrix<_Int,_Flt> &_at) const;
+      /// Compute transposed matrix for incomplete list
+   void TransposedSparsityList (int &_icycle, int *_imask, int *_indarr, int *_iptr, int *_listloc, int *_ialoc,
+                                 CMatrix<_Int,_Flt> &_at) const;
+      /// Extend sparsity by zeroes
+   void ExtendSparsity (const CMatrix<_Int,_Flt> &_a_sp,
+                         CMatrix<_Int,_Flt> &_a_ext) const;
+      /// Perform multiplication of matrices stored by rows, resulting matrix is also stored by rows
+   static void MatrixByMatrixMultiply (int &_icycle, int *_imask, int *_imask1, int *_indarr, int *_listloc, _Flt *_fmask,
+                                          CMatrix<_Int,_Flt> &_a, CMatrix<_Int,_Flt> &_b,
+                                          CMatrix<_Int,_Flt> &_a_times_b);
+      /// Compute the packed size
+   int GetPackedSize () const;
+      /// Fill char array by the packed data
+   void FillPacked (int _length, char *_obj) const;
+      /// Unpack data
+   void UnPack (int _length, char *_obj);
+};
+
+template <typename _Int, typename _Flt, typename _FltVect>
+class CMvmSlv_impl // Class that supports multiplication computations
+{
+public:
+      /// Multiply by super sparse matrix by rows and add result into prescribed positions
+   static void MvmA (int _nlist, const vector<_Int> &_list_alu, 
+                        const vector<_Int> &_ia_alu, const vector<_Int> &_ja_alu, const vector<_Flt> &_a_alu,
+                        const _FltVect *_x, _FltVect *_ax);
+      /// Multiply by super sparse matrix by columns add result into prescribed positions
+   static void MvmAT (int _nlist, const vector<_Int> &_list_alu, 
+                        const vector<_Int> &_ia_alu, const vector<_Int> &_ja_alu, const vector<_Flt> &_a_alu,
+                        const _FltVect *_x, _FltVect *_ax);
+      /// Solve with L, L is stored by columns (diag is inverted)
+   static void SolveL (int _n, const vector<_Int> &_ia_l, const vector<_Int> &_ja_l, const vector<_Flt> &_a_l,
+                        const _FltVect *_x, _FltVect *_lx);
+      /// Solve with U, U is stored by rows (diag is inverted)
+   static void SolveU (int _n, const vector<_Int> &_ia_u, const vector<_Int> &_ja_u, const vector<_Flt> &_a_u,
+                        const _FltVect *_x, _FltVect *_ux);
+};
+
+template <typename _Int, typename _Flt, typename _FltVect>
+class CMvmSlv // Class that supports multiplication computations
+{
+public:
+      /// Multiply by super sparse matrix by rows and add result into prescribed positions
+   static void MvmA (const CMatrix<_Int,_Flt> &_a_matr,
+                     const _FltVect *_x, _FltVect *_ax) {
+      CMvmSlv_impl<_Int, _Flt, _FltVect>::MvmA (_a_matr.GetNlist(), *(_a_matr.GetList()), *(_a_matr.GetIa()), *(_a_matr.GetJa()), *(_a_matr.GetA()), 
+                                                            _x, _ax);
+   };
+      /// Multiply by super sparse matrix by rows and add result into prescribed positions
+   static void MvmAT (const CMatrix<_Int,_Flt> &_a_matr,
+                     const _FltVect *_x, _FltVect *_ax) {
+      CMvmSlv_impl<_Int, _Flt, _FltVect>::MvmAT (_a_matr.GetNlist(), *(_a_matr.GetList()), *(_a_matr.GetIa()), *(_a_matr.GetJa()), *(_a_matr.GetA()), 
+                                                            _x, _ax);
+   };
+      /// Solve with L, L is stored by columns (diag is inverted)
+   static void SolveL (const CMatrix<_Int,_Flt> &_l_matr,
+                        const _FltVect *_x, _FltVect *_lx) {
+      CMvmSlv_impl<_Int, _Flt, _FltVect>::SolveL (_l_matr.GetN(), *(_l_matr.GetIa()), *(_l_matr.GetJa()), *(_l_matr.GetA()), 
+                                                            _x, _lx);
+   };
+      /// Solve with U, U is stored by rows (diag is inverted)
+   static void SolveU (const CMatrix<_Int,_Flt> &_u_matr,
+                        const _FltVect *_x, _FltVect *_ux) {
+      CMvmSlv_impl<_Int, _Flt, _FltVect>::SolveU (_u_matr.GetN(), *(_u_matr.GetIa()), *(_u_matr.GetJa()), *(_u_matr.GetA()), 
+                                                            _x, _ux);
+   };
+};
+
+template <typename _Int, typename _Flt>
+class CIlu2 // Class that supports factorization computations
+{
+public:
+// Functions that support matrix structure as a whole
+      /// Compute optimal ordering for the matrix
+   static void ComputeOptimalOrder (const CMatrix<_Int,_Flt> &_a_matr, int _ordtype, vector<int> &_order) {
+      CIlu2_impl<_Int,_Flt>::ComputeOptimalOrder (_ordtype, _a_matr.GetNlist(), *(_a_matr.GetIa()), *(_a_matr.GetJa()), _order);
+   };
+      /// Compute ordered matrix