PardisoSolver.h

// =======================================================================
// @(#)ParDirectSolver.h
//
// Class:      TPardisoSolver
// Purpose:    Solve equation system by OpenMP based direct
//         solver PARDISO          
//
// Author:     
//
// History:    
//
// =======================================================================

#ifndef __PARDIRECTSOLVER__
#define __PARDIRECTSOLVER__

#include <Database.h>
#include <MooNMD_Io.h>
#include <DirectSolver.h>

#ifdef __2D__
  #include <SquareMatrix2D.h>
  #include <Matrix2D.h>
#else
  #include <SquareMatrix3D.h>
  #include <Matrix3D.h>
#endif

#ifdef _OPENMP
  #include <omp.h>
#endif

class TPardisoSolver : public TDirectSolver
{
  protected:
    void **pt;

    int *iparam;
    
    double *dparam;

    int *KCol;
    int KCol_size;

    int *RowPtr;
    int RowPtr_size;

    double *Entries;
    int Entries_size;

    int N_Eq;
  
    int N_Entries;
    

    int num_prc;

    int mtype;

    bool pp;
    int rhs_index;
    
    std::ofstream *dat;
    double time_a;
    double time_f;
    double time_s;
    int runs_a;
    int runs_f;
    int runs_s;

  protected:
    void ShiftIndicies();
    void PrintIparam();
//     void PrintDparam();
    void ErrorMsg(int);
    void Sort();
    void FreeMemory();
    void AllocMemory();

    int SetMsgLvl()
    { return TDatabase::ParamDB->SC_VERBOSE > 0 ? 1 : 0;};

#ifdef __2D__
    // nstype 1
    void FillRowPtr(int N_U, int N_P, int N_Active,
            int *RowPtrA, int *RowPtrB, int *ColPtrB);

    void GetTransposedArrays(TMatrix *B, int *&KRowB,
                 int *&ColPtrB, int *&MapB);
#endif

  public:
    TPardisoSolver();

    ~TPardisoSolver();

#ifdef __2D__

    void SetMatrix(TSquareMatrix2D *Matrix);
    
    void SetMatrix(TSquareMatrix2D *A, TMatrix *B1, TMatrix *B2);

    void SetMatrixPar(TSquareMatrix2D *A, TMatrix *B1, TMatrix *B2);

   void SetMatrix(TSquareMatrix2D *sqmatrixA,
          TMatrix2D *matrixB1T, TMatrix2D *matrixB2T,
          TMatrix2D *matrixB1,  TMatrix2D *matrixB2);

   void SetMatrix(TSquareMatrix2D *sqmatrixA11, TSquareMatrix2D *sqmatrixA12,
          TSquareMatrix2D *sqmatrixA21, TSquareMatrix2D *sqmatrixA22,
          TMatrix2D *matrixB1T, TMatrix2D *matrixB2T,
          TMatrix2D *matrixB1,  TMatrix2D *matrixB2);
#endif

#ifdef __3D__

    void SetMatrix(TSquareMatrix3D *sqmatrixA,
           TMatrix3D *matrixB1T, TMatrix3D *matrixB2T, TMatrix3D *matrixB3T,
           TMatrix3D *matrixB1, TMatrix3D *matrixB2, TMatrix3D *matrixB3);
        
    void SetMatrix(TSquareMatrix3D *sqmatrixA11, TSquareMatrix3D *sqmatrixA12,
              TSquareMatrix3D *sqmatrixA13, TSquareMatrix3D *sqmatrixA21,
              TSquareMatrix3D *sqmatrixA22, TSquareMatrix3D *sqmatrixA23,
              TSquareMatrix3D *sqmatrixA31, TSquareMatrix3D *sqmatrixA32,
              TSquareMatrix3D *sqmatrixA33,
              TMatrix3D *matrixB1T, TMatrix3D *matrixB2T, TMatrix3D *matrixB3T,
              TMatrix3D *matrixB1, TMatrix3D *matrixB2, TMatrix3D *matrixB3);
#endif

    void Analyse();
    void Factorize();
    void Solve(double *sol, double *rhs);
    void FactorizeSolve(double *sol, double *rhs);

    void BenchReset();

};

#endif