HypreSolver.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
 
/// This HYPRE library interface has been taken originally from MFEM and modified
/// to suit the needs for the MOAB library.
/// Modified by: Vijay Mahadevan
 
#ifndef MOAB_HYPRESOLVER
#define MOAB_HYPRESOLVER
 
#include "moab/MOABConfig.h"
#include "moab/Core.hpp"
 
#ifdef MOAB_HAVE_MPI
#include <mpi.h>
 
// Enable internal hypre timing routines
#ifndef HYPRE_TIMING
#define HYPRE_TIMING
#endif
 
#include "HypreParVector.hpp"
#include "HypreParMatrix.hpp"
 
namespace moab
{
 
/// Base class for solvers
class AbstractSolver
{
 protected:
 const Eigen::SparseMatrix< double >* m_operator;
 
 public:
 /// If true, use the second argument of Mult as an initial guess
 bool iterative_mode;
 
 /// Initialize a Solver
 explicit AbstractSolver( bool iter_mode )
 {
 iterative_mode = iter_mode;
 }
 
 /// Set/update the solver for the given operator
 virtual void SetOperator( const HypreParMatrix& op ) = 0;
};
 
/// Abstract class for hypre's solvers and preconditioners
class HypreSolver : public AbstractSolver
{
 public:
 typedef HypreParVector Vector;
 
 protected:
 /// The linear system matrix
 HypreParMatrix* A;
 
 /// Right-hand side and solution vector
 mutable HypreParVector *B, *X;
 
 /// Was hypre's Setup function called already?
 mutable int setup_called;
 
 public:
 HypreSolver( bool iterative = true );
 
 HypreSolver( HypreParMatrix* _A, bool iterative = true );
 
 /// Typecast to HYPRE_Solver -- return the solver
 virtual operator HYPRE_Solver() const = 0;
 
 virtual void Verbosity( int /*print_level*/ )
 { /* nothing to do */
 }
 
 /// Set the hypre solver to be used as a preconditioner
 virtual void SetPreconditioner( HypreSolver& /*precond*/ )
 { /* Nothing to do */
 }
 
 /// hypre's internal Setup function
 virtual HYPRE_PtrToParSolverFcn SetupFcn() const = 0;
 /// hypre's internal Solve function
 virtual HYPRE_PtrToParSolverFcn SolveFcn() const = 0;
 
 virtual void SetOperator( const HypreParMatrix& /*op*/ )
 {
 MB_SET_ERR_RET( "HypreSolvers do not support SetOperator!" );
 }
 
 virtual void GetFinalResidualNorm( double& normr ) = 0;
 
 virtual void GetNumIterations( int& num_iterations ) = 0;
 
 /// Solve the linear system Ax=b
 virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
 
 virtual ~HypreSolver();
};
 
/// PCG solver in hypre
class HyprePCG : public HypreSolver
{
 private:
 HYPRE_Solver pcg_solver;
 
 public:
 HyprePCG( HypreParMatrix& _A );
 
 void SetTol( double tol );
 void SetMaxIter( int max_iter );
 void SetLogging( int logging );
 virtual void Verbosity( int print_lvl );
 
 /// Set the hypre solver to be used as a preconditioner
 virtual void SetPreconditioner( HypreSolver& precond );
 
 /** Use the L2 norm of the residual for measuring PCG convergence, plus
 (optionally) 1) periodically recompute true residuals from scratch; and
 2) enable residual-based stopping criteria. */
 void SetResidualConvergenceOptions( int res_frequency = -1, double rtol = 0.0 );
 
 /// non-hypre setting
 void SetZeroInintialIterate()
 {
 iterative_mode = false;
 }
 
 virtual void GetFinalResidualNorm( double& normr )
 {
 HYPRE_ParCSRPCGGetFinalRelativeResidualNorm( pcg_solver, &normr );
 }
 
 virtual void GetNumIterations( int& num_iterations )
 {
 HYPRE_Int num_it;
 HYPRE_ParCSRPCGGetNumIterations( pcg_solver, &num_it );
 num_iterations = internal::to_int( num_it );
 }
 
 /// The typecast to HYPRE_Solver returns the internal pcg_solver
 virtual operator HYPRE_Solver() const
 {
 return pcg_solver;
 }
 
 /// PCG Setup function
 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRPCGSetup;
 }
 /// PCG Solve function
 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRPCGSolve;
 }
 
 /// Solve Ax=b with hypre's PCG
 virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
 using HypreSolver::Solve;
 
 virtual ~HyprePCG();
};
 
/// GMRES solver in hypre
class HypreGMRES : public HypreSolver
{
 private:
 HYPRE_Solver gmres_solver;
 
 public:
 HypreGMRES( HypreParMatrix& _A );
 
 void SetTol( double atol, double rtol = 1e-20 );
 void SetMaxIter( int max_iter );
 void SetKDim( int dim );
 void SetLogging( int logging );
 virtual void Verbosity( int print_lvl );
 
 /// Set the hypre solver to be used as a preconditioner
 virtual void SetPreconditioner( HypreSolver& precond );
 
 /** Use the L2 norm of the residual for measuring PCG convergence, plus
 (optionally) 1) periodically recompute true residuals from scratch; and
 2) enable residual-based stopping criteria. */
 void SkipRealResidualCheck()
 {
 HYPRE_GMRESSetSkipRealResidualCheck( gmres_solver, 1 );
 }
 
 /// non-hypre setting
 void SetZeroInintialIterate()
 {
 iterative_mode = false;
 }
 
 virtual void GetFinalResidualNorm( double& normr )
 {
 HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm( gmres_solver, &normr );
 }
 
 virtual void GetNumIterations( int& num_iterations )
 {
 HYPRE_Int num_it;
 HYPRE_GMRESGetNumIterations( gmres_solver, &num_it );
 num_iterations = internal::to_int( num_it );
 }
 
 /// The typecast to HYPRE_Solver returns the internal gmres_solver
 virtual operator HYPRE_Solver() const
 {
 return gmres_solver;
 }
 
 /// GMRES Setup function
 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRGMRESSetup;
 }
 /// GMRES Solve function
 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRGMRESSolve;
 }
 
 /// Solve Ax=b with hypre's GMRES
 virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
 using HypreSolver::Solve;
 
 virtual ~HypreGMRES();
};
 
/// The identity operator as a hypre solver
class HypreIdentity : public HypreSolver
{
 public:
 virtual operator HYPRE_Solver() const
 {
 return NULL;
 }
 
 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)hypre_ParKrylovIdentitySetup;
 }
 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)hypre_ParKrylovIdentity;
 }
 
 virtual void GetNumIterations( int& num_iterations )
 {
 num_iterations = 1;
 }
 
 virtual void GetFinalResidualNorm( double& normr )
 {
 normr = 0.0;
 }
 
 virtual ~HypreIdentity() {}
};
 
/// Jacobi preconditioner in hypre
class HypreDiagScale : public HypreSolver
{
 public:
 HypreDiagScale() : HypreSolver() {}
 explicit HypreDiagScale( HypreParMatrix& A ) : HypreSolver( &A ) {}
 virtual operator HYPRE_Solver() const
 {
 return NULL;
 }
 
 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRDiagScaleSetup;
 }
 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRDiagScale;
 }
 
 virtual void GetNumIterations( int& num_iterations )
 {
 num_iterations = 1;
 }
 
 virtual void GetFinalResidualNorm( double& normr )
 {
 normr = 0.0;
 }
 
 HypreParMatrix* GetData()
 {
 return A;
 }
 virtual ~HypreDiagScale() {}
};
 
/// The ParaSails preconditioner in hypre
class HypreParaSails : public HypreSolver
{
 private:
 HYPRE_Solver sai_precond;
 
 public:
 HypreParaSails( HypreParMatrix& A );
 
 void SetSymmetry( int sym );
 
 virtual void GetNumIterations( int& num_iterations )
 {
 num_iterations = 1;
 }
 
 virtual void GetFinalResidualNorm( double& normr )
 {
 normr = 0.0;
 }
 
 /// The typecast to HYPRE_Solver returns the internal sai_precond
 virtual operator HYPRE_Solver() const
 {
 return sai_precond;
 }
 
 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_ParaSailsSetup;
 }
 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_ParaSailsSolve;
 }
 
 virtual ~HypreParaSails();
};
 
/// The BoomerAMG solver in hypre
class HypreBoomerAMG : public HypreSolver
{
 private:
 HYPRE_Solver amg_precond;
 
 /// Default, generally robust, BoomerAMG options
 void SetDefaultOptions();
 
 // If amg_precond is NULL, allocates it and sets default options.
 // Otherwise saves the options from amg_precond, destroys it, allocates a new
 // one, and sets its options to the saved values.
 void ResetAMGPrecond();
 
 public:
 HypreBoomerAMG();
 
 HypreBoomerAMG( HypreParMatrix& A );
 
 virtual void SetOperator( const HypreParMatrix& op );
 
 virtual void GetNumIterations( int& num_iterations )
 {
 num_iterations = 1;
 }
 
 virtual void GetFinalResidualNorm( double& normr )
 {
 normr = 0.0;
 }
 
 /** More robust options for systems, such as elasticity. Note that BoomerAMG
 assumes Ordering::byVDIM in the finite element space used to generate the
 matrix A. */
 void SetSystemsOptions( int dim );
 
 virtual void Verbosity( int print_level )
 {
 HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
 }
 
 /// The typecast to HYPRE_Solver returns the internal amg_precond
 virtual operator HYPRE_Solver() const
 {
 return amg_precond;
 }
 
 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_BoomerAMGSetup;
 }
 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
 {
 return (HYPRE_PtrToParSolverFcn)HYPRE_BoomerAMGSolve;
 }
 
 virtual ~HypreBoomerAMG();
};
 
} // namespace moab
 
#endif // MOAB_HAVE_MPI
 
#endif