Commit b2fde510 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Hypre solvers: refactor the reset and setMatrix functions to allow reusing the... 

Hypre solvers: refactor the reset and setMatrix functions to allow reusing the preconditioner in time-evolving problems
parent 22fe7cd2

src/Examples/Hypre/tnl-hypre-ex5.cpp

+3 −11
Original line number Diff line number Diff line
@@ -218,14 +218,12 @@ main( int argc, char* argv[] ) 
      TNL::Solvers::Linear::HypreBoomerAMG solver;



      // Set the matrix of the linear system

      // WARNING: setMatrix resets the preconditioner, including setting

      //          default options.

      // NOTE: The wrapper class sets its own default options that are

      //       different from Hypre. The overriding settings below result in

      //       the same state as the hypre-ex5.c example.

      solver.setMatrix( parcsr_A );



      // Set some parameters (See Reference Manual for more parameters)

      // NOTE: The wrapper class sets its own default options that are

      //       different from Hypre. The overriding settings below result in

      //       the same state as the hypre-ex5.c example.

      HYPRE_BoomerAMGSetPrintLevel( solver, 3 );    // Print solve info + parameters

      HYPRE_BoomerAMGSetOldDefault( solver );       // Falgout coarsening with modified classical interpolation

      HYPRE_BoomerAMGSetRelaxType( solver, 6 );     // Sym. G-S/Jacobi hybrid relaxation
@@ -271,8 +269,6 @@ main( int argc, char* argv[] ) 
      solver.setPreconditioner( precond );



      // Set the matrix of the linear system

      // WARNING: setMatrix resets the preconditioner, including setting

      //          default options.

      solver.setMatrix( parcsr_A );



      // Set some parameters (See Reference Manual for more parameters)
@@ -303,8 +299,6 @@ main( int argc, char* argv[] ) 
      solver.setPreconditioner( precond );



      // Set the matrix of the linear system

      // WARNING: setMatrix resets the preconditioner, including setting

      //          default options.

      solver.setMatrix( parcsr_A );



      // Set some parameters (See Reference Manual for more parameters)
@@ -333,8 +327,6 @@ main( int argc, char* argv[] ) 
      solver.setPreconditioner( precond );



      // Set the matrix of the linear system

      // WARNING: setMatrix resets the preconditioner, including setting

      //          default options.

      solver.setMatrix( parcsr_A );



      // Set some parameters (See Reference Manual for more parameters)

src/TNL/Solvers/Linear/Hypre.h

+31 −39
Original line number Diff line number Diff line
@@ -72,16 +72,38 @@ public: 
   /**

    * \brief Set the matrix of the linear system to be solved.

    *

    * This function also resets the current state of the solver causing the

    * Hypre's internal setup function to be called before the solve function.

    * This function also resets the internal flag indicating whether the Hypre

    * setup function was called for the current matrix.

    *

    * \param reuse_setup When true, the result of the previous setup phase will

    *                    be preserved, i.e., the solver (and preconditioner)

    *                    will not be updated for the new matrix when calling

    *                    the \ref solve method.

    */

   virtual void setMatrix( const Matrices::HypreParCSRMatrix& op )

   virtual void setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup = false )

   {

      A = &op;

      if( setup_called && reuse_setup )

         setup_called = true;

      else

         setup_called = false;

   }



   //! Solve the linear system Ax=b

   /**

    * \brief Setup the solver for solving the linear system Ax=b.

    *

    * Calling this function repeatedly has no effect until the internal flag is

    * reset via the \ref setMatrix function.

    */

   virtual void

   setup( const Containers::HypreParVector& b, Containers::HypreParVector& x ) const;



   /**

    * \brief Solve the linear system Ax=b.

    *

    * This function checks if the \ref setup function was already called and

    * calls it if necessary.

    */

   virtual void

   solve( const Containers::HypreParVector& b, Containers::HypreParVector& x ) const;
@@ -125,7 +147,7 @@ public: 
   ~HyprePCG() override;



   void

   setMatrix( const Matrices::HypreParCSRMatrix& op ) override;

   setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup = false ) override;



   //! Set the Hypre solver to be used as a preconditioner

   void
@@ -219,7 +241,7 @@ public: 
   ~HypreBiCGSTAB() override;



   void

   setMatrix( const Matrices::HypreParCSRMatrix& op ) override;

   setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup = false ) override;



   //! Set the Hypre solver to be used as a preconditioner

   void
@@ -302,7 +324,7 @@ public: 
   ~HypreGMRES() override;



   void

   setMatrix( const Matrices::HypreParCSRMatrix& op ) override;

   setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup = false ) override;



   //! Set the Hypre solver to be used as a preconditioner

   void
@@ -390,7 +412,7 @@ public: 
   ~HypreFlexGMRES() override;



   void

   setMatrix( const Matrices::HypreParCSRMatrix& op ) override;

   setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup = false ) override;



   //! Set the Hypre solver to be used as a preconditioner

   void
@@ -518,11 +540,6 @@ public: 
 */

class HypreParaSails : public HypreSolver

{

private:

   //! \brief Reset the preconditioner to the default state

   void

   reset( MPI_Comm comm );



public:

   explicit HypreParaSails( MPI_Comm comm );
@@ -530,9 +547,6 @@ public: 


   ~HypreParaSails() override;



   void

   setMatrix( const Matrices::HypreParCSRMatrix& op ) override;



   HYPRE_PtrToParSolverFcn

   setupFcn() const override

   {
@@ -567,11 +581,6 @@ public: 
 */

class HypreEuclid : public HypreSolver

{

private:

   //! \brief Reset the preconditioner to the default state

   void

   reset( MPI_Comm comm );



public:

   explicit HypreEuclid( MPI_Comm comm );
@@ -579,9 +588,6 @@ public: 


   ~HypreEuclid() override;



   void

   setMatrix( const Matrices::HypreParCSRMatrix& op ) override;



   HYPRE_PtrToParSolverFcn

   setupFcn() const override

   {
@@ -616,10 +622,6 @@ private: 
   void

   setDefaultOptions();



   //! \brief Reset the preconditioner to the default state (including default options).

   void

   reset();



public:

   HypreILU();
@@ -634,9 +636,6 @@ public: 
      HYPRE_ILUSetPrintLevel( solver, print_level );

   }



   void

   setMatrix( const Matrices::HypreParCSRMatrix& op ) override;



   HYPRE_PtrToParSolverFcn

   setupFcn() const override

   {
@@ -676,10 +675,6 @@ private: 
   void

   setDefaultOptions();



   //! \brief Reset the preconditioner to the default state (including default options).

   void

   reset();



public:

   HypreBoomerAMG();
@@ -687,9 +682,6 @@ public: 


   ~HypreBoomerAMG() override;



   void

   setMatrix( const Matrices::HypreParCSRMatrix& op ) override;



   //! \brief More robust options for systems of PDEs

   void

   setSystemsOptions( int dim, bool order_bynodes = false );

src/TNL/Solvers/Linear/Hypre.hpp

+47 −97
Original line number Diff line number Diff line
@@ -22,15 +22,15 @@ HypreSolver::HypreSolver( const Matrices::HypreParCSRMatrix& A ) 
}



void

HypreSolver::solve( const Containers::HypreParVector& b, Containers::HypreParVector& x ) const

HypreSolver::setup( const Containers::HypreParVector& b, Containers::HypreParVector& x ) const

{

   HYPRE_Int err_flag;

   if( A == nullptr ) {

      throw std::runtime_error( "HypreSolver::solve(...) : HypreParCSRMatrix A is missing" );

   }

   if( A == nullptr )

      throw std::runtime_error( "HypreSolver::setup(...) : HypreParCSRMatrix A is missing" );



   if( setup_called )

      return;



   if( ! setup_called ) {

      err_flag = setupFcn()( *this, *A, b, x );

   const HYPRE_Int err_flag = setupFcn()( *this, *A, b, x );

   if( err_flag != 0 ) {

      if( error_mode == WARN_HYPRE_ERRORS )

         std::cout << "Error during setup! Error code: " << err_flag << std::endl;
@@ -38,10 +38,20 @@ HypreSolver::solve( const Containers::HypreParVector& b, Containers::HypreParVec 
         throw std::runtime_error( "Error during setup! Error code: " + std::to_string( err_flag ) );

   }

   hypre_error_flag = 0;



   setup_called = true;

}



   err_flag = solveFcn()( *this, *A, b, x );

void

HypreSolver::solve( const Containers::HypreParVector& b, Containers::HypreParVector& x ) const

{

   if( A == nullptr )

      throw std::runtime_error( "HypreSolver::solve(...) : HypreParCSRMatrix A is missing" );



   if( ! setup_called )

      setup( b, x );



   const HYPRE_Int err_flag = solveFcn()( *this, *A, b, x );

   if( err_flag != 0 ) {

      if( error_mode == WARN_HYPRE_ERRORS )

         std::cout << "Error during solve! Error code: " << err_flag << std::endl;
@@ -70,11 +80,11 @@ HyprePCG::~HyprePCG() 
}



void

HyprePCG::setMatrix( const Matrices::HypreParCSRMatrix& op )

HyprePCG::setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup )

{

   HypreSolver::setMatrix( op );

   HypreSolver::setMatrix( op, reuse_setup );

   if( precond != nullptr )

      precond->setMatrix( *A );

      precond->setMatrix( *A, reuse_setup );

}



void
@@ -132,11 +142,11 @@ HypreBiCGSTAB::~HypreBiCGSTAB() 
}



void

HypreBiCGSTAB::setMatrix( const Matrices::HypreParCSRMatrix& op )

HypreBiCGSTAB::setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup )

{

   HypreSolver::setMatrix( op );

   HypreSolver::setMatrix( op, reuse_setup );

   if( precond != nullptr )

      precond->setMatrix( *A );

      precond->setMatrix( *A, reuse_setup );

}



void
@@ -182,11 +192,11 @@ HypreGMRES::~HypreGMRES() 
}



void

HypreGMRES::setMatrix( const Matrices::HypreParCSRMatrix& op )

HypreGMRES::setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup )

{

   HypreSolver::setMatrix( op );

   HypreSolver::setMatrix( op, reuse_setup );

   if( precond != nullptr )

      precond->setMatrix( *A );

      precond->setMatrix( *A, reuse_setup );

}



void
@@ -234,11 +244,11 @@ HypreFlexGMRES::~HypreFlexGMRES() 
}



void

HypreFlexGMRES::setMatrix( const Matrices::HypreParCSRMatrix& op )

HypreFlexGMRES::setMatrix( const Matrices::HypreParCSRMatrix& op, bool reuse_setup )

{

   HypreSolver::setMatrix( op );

   HypreSolver::setMatrix( op, reuse_setup );

   if( precond != nullptr )

      precond->setMatrix( *A );

     precond->setMatrix( *A, reuse_setup );

}



void
@@ -271,12 +281,12 @@ HypreFlexGMRES::postSolveHook() const 


HypreParaSails::HypreParaSails( MPI_Comm comm )

{

   reset( comm );

   HYPRE_ParaSailsCreate( comm, &solver );

}



HypreParaSails::HypreParaSails( const Matrices::HypreParCSRMatrix& A ) : HypreSolver( A )

{

   reset( A.getCommunicator() );

   HYPRE_ParaSailsCreate( A.getCommunicator() , &solver );

}



HypreParaSails::~HypreParaSails()
@@ -284,22 +294,6 @@ HypreParaSails::~HypreParaSails() 
   HYPRE_ParaSailsDestroy( solver );

}



void

HypreParaSails::reset( MPI_Comm comm )

{

   if( solver != nullptr )

      HYPRE_ParaSailsDestroy( solver );

   HYPRE_ParaSailsCreate( comm, &solver );

}



void

HypreParaSails::setMatrix( const Matrices::HypreParCSRMatrix& op )

{

   if( A != nullptr )

      reset( A->getCommunicator() );

   HypreSolver::setMatrix( op );

}



HypreEuclid::HypreEuclid( MPI_Comm comm )

{

   HYPRE_EuclidCreate( comm, &solver );
@@ -315,30 +309,16 @@ HypreEuclid::~HypreEuclid() 
   HYPRE_EuclidDestroy( solver );

}



void

HypreEuclid::reset( MPI_Comm comm )

{

   if( solver != nullptr )

      HYPRE_EuclidDestroy( solver );

   HYPRE_EuclidCreate( comm, &solver );

}



void

HypreEuclid::setMatrix( const Matrices::HypreParCSRMatrix& op )

{

   if( A != nullptr )

      reset( A->getCommunicator() );

   HypreSolver::setMatrix( op );

}



HypreILU::HypreILU()

{

   reset();

   HYPRE_ILUCreate( &solver );

   setDefaultOptions();

}



HypreILU::HypreILU( const Matrices::HypreParCSRMatrix& A ) : HypreSolver( A )

{

   reset();

   HYPRE_ILUCreate( &solver );

   setDefaultOptions();

}



HypreILU::~HypreILU()
@@ -365,30 +345,16 @@ HypreILU::setDefaultOptions() 
   HYPRE_ILUSetLocalReordering( solver, 1 );

}



void

HypreILU::reset()

{

   if( solver != nullptr )

      HYPRE_ILUDestroy( solver );

   HYPRE_ILUCreate( &solver );

   setDefaultOptions();

}



void

HypreILU::setMatrix( const Matrices::HypreParCSRMatrix& op )

{

   reset();

   HypreSolver::setMatrix( op );

}



HypreBoomerAMG::HypreBoomerAMG()

{

   reset();

   HYPRE_BoomerAMGCreate( &solver );

   setDefaultOptions();

}



HypreBoomerAMG::HypreBoomerAMG( const Matrices::HypreParCSRMatrix& A ) : HypreSolver( A )

{

   reset();

   HYPRE_BoomerAMGCreate( &solver );

   setDefaultOptions();

}



HypreBoomerAMG::~HypreBoomerAMG()
@@ -450,22 +416,6 @@ HypreBoomerAMG::setDefaultOptions() 
   HYPRE_BoomerAMGSetTol( solver, 0.0 );

}



void

HypreBoomerAMG::reset()

{

   if( solver != nullptr )

      HYPRE_BoomerAMGDestroy( solver );

   HYPRE_BoomerAMGCreate( &solver );

   setDefaultOptions();

}



void

HypreBoomerAMG::setMatrix( const Matrices::HypreParCSRMatrix& op )

{

   reset();

   HypreSolver::setMatrix( op );

}



void

HypreBoomerAMG::setSystemsOptions( int dim, bool order_bynodes )

{