solver.hpp

Go to the documentation of this file.

#ifndef PARGEMSLR_SOLVER_H
#define PARGEMSLR_SOLVER_H
 
#include "../utils/utils.hpp"
#include "../vectors/vector.hpp"
#include "../vectors/sequential_vector.hpp"
#include "../vectors/parallel_vector.hpp"
#include "../matrices/matrix.hpp"
#include "../matrices/csr_matrix.hpp"
#include "../matrices/dense_matrix.hpp"
#include "../matrices/parallel_csr_matrix.hpp"
 
namespace pargemslr
{
   
   enum SolverTypeEnum
   {
      kSolverUndefined = -1,
      kSolverFgmres,
      kSolverIlu,
      kSolverBJ,
      kSolverFloat,
      kSolverGemslr,
      kSolverParGemslr
   };
   
    template <class MatrixType, class VectorType, typename DataType>
   class SolverClass
    {
   protected:
      
      PrecisionEnum                                      _solver_precision;
      
      SolverTypeEnum                                     _solver_type;
      
      bool                                               _is_mixed;
      
      MatrixType                                         *_matrix;
      
      bool                                               _own_matrix;
      
      SolverClass<MatrixType, VectorType, DataType>      *_preconditioner;
      
      bool                                               _own_preconditioner;
      
      VectorType                                         *_solution;
      
      VectorType                                         *_right_hand_side;
      
      bool                                               _ready;
      
      int                                                _print_option;
      
      int CheckReadySetups(const char *str) const
      {
         if(this->_ready)
         {
            char str1[1024];
            snprintf( str1, 1024, "Change setup %s after preconditioner is built is invalid.", str );
            PARGEMSLR_ERROR(str1);
            return PARGEMSLR_ERROR_FUNCTION_CALL_ERR;
         }
         return PARGEMSLR_SUCCESS;
      }
      
    public:
      
      SolverClass()
      {
         this->_solver_precision = kUnknownPrecision;
         this->_solver_type = kSolverUndefined;
         /* if vector type not equal to data type, this is a mixed precision preconditioner */
         this->_is_mixed = PargemslrIsDoublePrecision<VectorType>::value != PargemslrIsDoublePrecision<DataType>::value;
         this->_matrix = NULL;
         this->_own_matrix = false;
         this->_preconditioner = NULL;
         this->_own_preconditioner = false;
         this->_solution = NULL;
         this->_right_hand_side = NULL;
         this->_ready = false;
         this->_print_option = 0;
      }
      
      SolverClass(const SolverClass<MatrixType, VectorType, DataType> &solver)
      {
         this->_solver_precision = solver._solver_precision;
         this->_solver_type = solver._solver_type;
         this->_is_mixed = solver._is_mixed;
         this->_matrix = solver._matrix;
         this->_own_matrix = false;
         this->_preconditioner = solver._preconditioner;
         this->_own_preconditioner = false;
         this->_right_hand_side = solver._right_hand_side;
         this->_ready = solver._ready;
         this->_print_option = solver._print_option;
      }
      
      SolverClass(SolverClass<MatrixType, VectorType, DataType> &&solver)
      {
         this->_solver_precision = solver._solver_precision;
         solver._solver_precision = kUnknownPrecision;
         this->_solver_type = solver._solver_type;
         solver._solver_type = kSolverUndefined;
         this->_is_mixed = solver._is_mixed;
         this->_matrix = solver._matrix;
         solver._matrix = NULL;
         this->_own_matrix = solver._own_matrix;
         solver._own_matrix = false;
         this->_preconditioner = solver._preconditioner;
         solver._preconditioner = NULL;
         this->_own_preconditioner = solver._own_preconditioner;
         solver._own_preconditioner = false;
         this->_solution = solver._solution;
         solver._solution = NULL;
         this->_right_hand_side = solver._right_hand_side;
         solver._right_hand_side = NULL;
         this->_ready = solver._ready;
         solver._ready = false;
         this->_print_option = solver._print_option;
         solver._print_option = 0;
      }
      
      SolverClass<MatrixType, VectorType, DataType>& operator= (const SolverClass<MatrixType, VectorType, DataType> &solver)
      {
         this->Clear();
         this->_solver_precision = solver._solver_precision;
         this->_solver_type = solver._solver_type;
         this->_is_mixed = solver._is_mixed;
         this->_own_matrix = false;
         this->_preconditioner = solver._preconditioner;
         this->_own_preconditioner = false;
         this->_solution = solver._solution;
         this->_right_hand_side = solver._right_hand_side;
         this->_ready = solver._ready;
         this->_print_option = solver._print_option;
         return *this;
      }
      
      SolverClass<MatrixType, VectorType, DataType>& operator= (SolverClass<MatrixType, VectorType, DataType> &&solver)
      {
         this->Clear();
         this->_solver_precision = solver._solver_precision;
         solver._solver_precision = kUnknownPrecision;
         this->_solver_type = solver._solver_type;
         solver._solver_type = kSolverUndefined;
         this->_is_mixed = solver._is_mixed;
         this->_matrix = solver._matrix;
         solver._matrix = NULL;
         this->_own_matrix = solver._own_matrix;
         solver._own_matrix = false;
         this->_preconditioner = solver._preconditioner;
         solver._preconditioner = NULL;
         this->_own_preconditioner = solver._own_preconditioner;
         solver._own_preconditioner = false;
         this->_solution = solver._solution;
         solver._solution = NULL;
         this->_right_hand_side = solver._right_hand_side;
         solver._right_hand_side = NULL;
         this->_ready = solver._ready;
         solver._ready = false;
         this->_print_option = solver._print_option;
         solver._print_option = 0;
         return *this;
      }
      
      virtual int Clear()
      {
         
         this->_solver_precision = kUnknownPrecision;
         this->_solver_type = kSolverUndefined;
         
         if(this->_matrix && this->_own_matrix)
         {
            this->_matrix->Clear();
            delete this->_matrix;
         }
         this->_matrix = NULL;
         
         if(this->_preconditioner && this->_own_preconditioner)
         {
            this->_preconditioner->Clear();
            delete this->_preconditioner;
         }
         this->_preconditioner = NULL;
         this->_solution = NULL;
         this->_right_hand_side = NULL;
         
         this->_own_matrix = false;
         this->_own_preconditioner = false;
         this->_ready = false;
         this->_print_option = 0;
         
         return PARGEMSLR_SUCCESS;
      }
     
      virtual ~SolverClass()
      {
         /* call base clear function */
         SolverClass<MatrixType, VectorType, DataType>::Clear();
      }
      
      virtual int Setup( VectorType &x, VectorType &rhs) = 0;
      
      virtual int Solve( VectorType &x, VectorType &rhs) = 0;
      
      virtual int Solve( DenseMatrixClass<DataType> &x, DenseMatrixClass<DataType> &rhs)
      {
         PARGEMSLR_ERROR("Error: Solve with multiple right-hand-side unimplemented for this class");
         return PARGEMSLR_ERROR_FUNCTION_CALL_ERR;
      }
      
      virtual long int  GetNumNonzeros()
      {
         PARGEMSLR_ERROR("Error: GetNumNonzeros unimplemented for this class");
         return PARGEMSLR_ERROR_FUNCTION_CALL_ERR;
      }
      
      MatrixType* GetMatrix()
      {
         return _matrix;
      }
      
      SolverClass<MatrixType, VectorType, DataType>* GetPreconditioner()
      {
         return _preconditioner;
      }
      
      virtual int SetWithParameterArray(double *params)
      {
         PARGEMSLR_ERROR("Error: SetWithParameterArray unimplemented for this class");
         return PARGEMSLR_ERROR_FUNCTION_CALL_ERR;
      }
      
      virtual int SetSolveLocation( const int &location) = 0;
      
      int SetMatrix( MatrixType &matrix)
      {
         return this->SetMatrixP(&matrix);
      }
      
      int SetMatrixP( MatrixType *matrix)
      {
         if(matrix == this->_matrix)
         {
            /* set to the same matrix, do nothing */
            return PARGEMSLR_SUCCESS;
         }
         
         if(this->_matrix && this->_own_matrix)
         {
            this->_matrix->Clear();
            delete this->_matrix;
         }
         
         /* set a new matrix */
         this->_matrix = matrix;
         
         if(this->_preconditioner && this->_own_preconditioner)
         {
            this->_preconditioner->Clear();
            delete this->_preconditioner;
         }
         
         /* also destroy the preconditioner */
         this->_preconditioner = NULL;
         this->_solution = NULL;
         this->_right_hand_side = NULL;
         
         this->_ready = false;
         
         this->_own_matrix = false;
         this->_own_preconditioner = false;
         
         return PARGEMSLR_SUCCESS;
      }
      
      int SetPreconditioner( SolverClass<MatrixType, VectorType, DataType> &precond)
      {
         return this->SetPreconditionerP(&precond);
      }
      
      int SetPreconditionerP( SolverClass<MatrixType, VectorType, DataType>* precond)
      {
         if(this->_preconditioner && this->_own_preconditioner)
         {
            this->_preconditioner->Clear();
            delete this->_preconditioner;
         }
         this->_preconditioner = precond;
         this->_own_preconditioner = false;
         return PARGEMSLR_SUCCESS;
      }
      
      bool IsReady() const
      {
         return this->_ready;
      }
      
      int SetOwnMatrix(bool own_matrix)
      {
         this->_own_matrix = own_matrix;
         return PARGEMSLR_SUCCESS;
      }
      
      int SetOwnPreconditioner(bool own_preconditioner)
      {
         this->_own_preconditioner = own_preconditioner;
         return PARGEMSLR_SUCCESS;
      }
      
      int SetPrintOption(int print_option)
      {
         this->_print_option = print_option;
         return PARGEMSLR_SUCCESS;
      }
      
      PrecisionEnum GerPrecision() const
      {
         return this->_solver_precision;
      }
      
      SolverTypeEnum GetSolverType() const
      {
         return this->_solver_type;
      }
      
    };
   
   typedef SolverClass<CsrMatrixClass<float>, SequentialVectorClass<float>, float>           solver_csr_seq_float;
   typedef SolverClass<CsrMatrixClass<double>, SequentialVectorClass<double>, float>         solver_csr_seq_mix;
   typedef SolverClass<CsrMatrixClass<double>, SequentialVectorClass<double>, double>        solver_csr_seq_double;
   typedef SolverClass<CsrMatrixClass<complexs>, SequentialVectorClass<complexs>, complexs>  solver_csr_seq_complexs;
   typedef SolverClass<CsrMatrixClass<complexd>, SequentialVectorClass<complexd>, complexs>  solver_csr_seq_complexmix;
   typedef SolverClass<CsrMatrixClass<complexd>, SequentialVectorClass<complexd>, complexd>  solver_csr_seq_complexd;
   
   typedef SolverClass<ParallelCsrMatrixClass<float>, ParallelVectorClass<float>, float>           solver_csr_par_float;
   typedef SolverClass<ParallelCsrMatrixClass<double>, ParallelVectorClass<double>, double>        solver_csr_par_double;
   typedef SolverClass<ParallelCsrMatrixClass<double>, ParallelVectorClass<double>, double>        solver_csr_par_double;
   typedef SolverClass<ParallelCsrMatrixClass<complexs>, ParallelVectorClass<complexs>, complexs>  solver_csr_par_complexs;
   typedef SolverClass<ParallelCsrMatrixClass<complexd>, ParallelVectorClass<complexd>, complexs>  solver_csr_par_complexmix;
   typedef SolverClass<ParallelCsrMatrixClass<complexd>, ParallelVectorClass<complexd>, complexd>  solver_csr_par_complexd;
}
 
#endif