Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
HypreSolver.hpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.org.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 /// This HYPRE library interface has been taken originally from MFEM and modified
00013 /// to suit the needs for the MOAB library.
00014 /// Modified by: Vijay Mahadevan
00015 
00016 #ifndef MOAB_HYPRESOLVER
00017 #define MOAB_HYPRESOLVER
00018 
00019 #include "moab/MOABConfig.h"
00020 #include "moab/Core.hpp"
00021 
00022 #ifdef MOAB_HAVE_MPI
00023 #include <mpi.h>
00024 
00025 // Enable internal hypre timing routines
00026 #ifndef HYPRE_TIMING
00027 #define HYPRE_TIMING
00028 #endif
00029 
00030 #include "HypreParVector.hpp"
00031 #include "HypreParMatrix.hpp"
00032 
00033 namespace moab
00034 {
00035 
00036 /// Base class for solvers
00037 class AbstractSolver
00038 {
00039   protected:
00040     const Eigen::SparseMatrix< double >* m_operator;
00041 
00042   public:
00043     /// If true, use the second argument of Mult as an initial guess
00044     bool iterative_mode;
00045 
00046     /// Initialize a Solver
00047     explicit AbstractSolver( bool iter_mode )
00048     {
00049         iterative_mode = iter_mode;
00050     }
00051 
00052     /// Set/update the solver for the given operator
00053     virtual void SetOperator( const HypreParMatrix& op ) = 0;
00054 };
00055 
00056 /// Abstract class for hypre's solvers and preconditioners
00057 class HypreSolver : public AbstractSolver
00058 {
00059   public:
00060     typedef HypreParVector Vector;
00061 
00062   protected:
00063     /// The linear system matrix
00064     HypreParMatrix* A;
00065 
00066     /// Right-hand side and solution vector
00067     mutable HypreParVector *B, *X;
00068 
00069     /// Was hypre's Setup function called already?
00070     mutable int setup_called;
00071 
00072   public:
00073     HypreSolver( bool iterative = true );
00074 
00075     HypreSolver( HypreParMatrix* _A, bool iterative = true );
00076 
00077     /// Typecast to HYPRE_Solver -- return the solver
00078     virtual operator HYPRE_Solver() const = 0;
00079 
00080     virtual void Verbosity( int /*print_level*/ )
00081     { /* nothing to do */
00082     }
00083 
00084     /// Set the hypre solver to be used as a preconditioner
00085     virtual void SetPreconditioner( HypreSolver& /*precond*/ )
00086     { /* Nothing to do */
00087     }
00088 
00089     /// hypre's internal Setup function
00090     virtual HYPRE_PtrToParSolverFcn SetupFcn() const = 0;
00091     /// hypre's internal Solve function
00092     virtual HYPRE_PtrToParSolverFcn SolveFcn() const = 0;
00093 
00094     virtual void SetOperator( const HypreParMatrix& /*op*/ )
00095     {
00096         MB_SET_ERR_RET( "HypreSolvers do not support SetOperator!" );
00097     }
00098 
00099     virtual void GetFinalResidualNorm( double& normr ) = 0;
00100 
00101     virtual void GetNumIterations( int& num_iterations ) = 0;
00102 
00103     /// Solve the linear system Ax=b
00104     virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
00105 
00106     virtual ~HypreSolver();
00107 };
00108 
00109 /// PCG solver in hypre
00110 class HyprePCG : public HypreSolver
00111 {
00112   private:
00113     HYPRE_Solver pcg_solver;
00114 
00115   public:
00116     HyprePCG( HypreParMatrix& _A );
00117 
00118     void SetTol( double tol );
00119     void SetMaxIter( int max_iter );
00120     void SetLogging( int logging );
00121     virtual void Verbosity( int print_lvl );
00122 
00123     /// Set the hypre solver to be used as a preconditioner
00124     virtual void SetPreconditioner( HypreSolver& precond );
00125 
00126     /** Use the L2 norm of the residual for measuring PCG convergence, plus
00127         (optionally) 1) periodically recompute true residuals from scratch; and
00128         2) enable residual-based stopping criteria. */
00129     void SetResidualConvergenceOptions( int res_frequency = -1, double rtol = 0.0 );
00130 
00131     /// non-hypre setting
00132     void SetZeroInintialIterate()
00133     {
00134         iterative_mode = false;
00135     }
00136 
00137     virtual void GetFinalResidualNorm( double& normr )
00138     {
00139         HYPRE_ParCSRPCGGetFinalRelativeResidualNorm( pcg_solver, &normr );
00140     }
00141 
00142     virtual void GetNumIterations( int& num_iterations )
00143     {
00144         HYPRE_Int num_it;
00145         HYPRE_ParCSRPCGGetNumIterations( pcg_solver, &num_it );
00146         num_iterations = internal::to_int( num_it );
00147     }
00148 
00149     /// The typecast to HYPRE_Solver returns the internal pcg_solver
00150     virtual operator HYPRE_Solver() const
00151     {
00152         return pcg_solver;
00153     }
00154 
00155     /// PCG Setup function
00156     virtual HYPRE_PtrToParSolverFcn SetupFcn() const
00157     {
00158         return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRPCGSetup;
00159     }
00160     /// PCG Solve function
00161     virtual HYPRE_PtrToParSolverFcn SolveFcn() const
00162     {
00163         return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRPCGSolve;
00164     }
00165 
00166     /// Solve Ax=b with hypre's PCG
00167     virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
00168     using HypreSolver::Solve;
00169 
00170     virtual ~HyprePCG();
00171 };
00172 
00173 /// GMRES solver in hypre
00174 class HypreGMRES : public HypreSolver
00175 {
00176   private:
00177     HYPRE_Solver gmres_solver;
00178 
00179   public:
00180     HypreGMRES( HypreParMatrix& _A );
00181 
00182     void SetTol( double atol, double rtol = 1e-20 );
00183     void SetMaxIter( int max_iter );
00184     void SetKDim( int dim );
00185     void SetLogging( int logging );
00186     virtual void Verbosity( int print_lvl );
00187 
00188     /// Set the hypre solver to be used as a preconditioner
00189     virtual void SetPreconditioner( HypreSolver& precond );
00190 
00191     /** Use the L2 norm of the residual for measuring PCG convergence, plus
00192         (optionally) 1) periodically recompute true residuals from scratch; and
00193         2) enable residual-based stopping criteria. */
00194     void SkipRealResidualCheck()
00195     {
00196         HYPRE_GMRESSetSkipRealResidualCheck( gmres_solver, 1 );
00197     }
00198 
00199     /// non-hypre setting
00200     void SetZeroInintialIterate()
00201     {
00202         iterative_mode = false;
00203     }
00204 
00205     virtual void GetFinalResidualNorm( double& normr )
00206     {
00207         HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm( gmres_solver, &normr );
00208     }
00209 
00210     virtual void GetNumIterations( int& num_iterations )
00211     {
00212         HYPRE_Int num_it;
00213         HYPRE_GMRESGetNumIterations( gmres_solver, &num_it );
00214         num_iterations = internal::to_int( num_it );
00215     }
00216 
00217     /// The typecast to HYPRE_Solver returns the internal gmres_solver
00218     virtual operator HYPRE_Solver() const
00219     {
00220         return gmres_solver;
00221     }
00222 
00223     /// GMRES Setup function
00224     virtual HYPRE_PtrToParSolverFcn SetupFcn() const
00225     {
00226         return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRGMRESSetup;
00227     }
00228     /// GMRES Solve function
00229     virtual HYPRE_PtrToParSolverFcn SolveFcn() const
00230     {
00231         return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRGMRESSolve;
00232     }
00233 
00234     /// Solve Ax=b with hypre's GMRES
00235     virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
00236     using HypreSolver::Solve;
00237 
00238     virtual ~HypreGMRES();
00239 };
00240 
00241 /// The identity operator as a hypre solver
00242 class HypreIdentity : public HypreSolver
00243 {
00244   public:
00245     virtual operator HYPRE_Solver() const
00246     {
00247         return NULL;
00248     }
00249 
00250     virtual HYPRE_PtrToParSolverFcn SetupFcn() const
00251     {
00252         return (HYPRE_PtrToParSolverFcn)hypre_ParKrylovIdentitySetup;
00253     }
00254     virtual HYPRE_PtrToParSolverFcn SolveFcn() const
00255     {
00256         return (HYPRE_PtrToParSolverFcn)hypre_ParKrylovIdentity;
00257     }
00258 
00259     virtual void GetNumIterations( int& num_iterations )
00260     {
00261         num_iterations = 1;
00262     }
00263 
00264     virtual void GetFinalResidualNorm( double& normr )
00265     {
00266         normr = 0.0;
00267     }
00268 
00269     virtual ~HypreIdentity() {}
00270 };
00271 
00272 /// Jacobi preconditioner in hypre
00273 class HypreDiagScale : public HypreSolver
00274 {
00275   public:
00276     HypreDiagScale() : HypreSolver() {}
00277     explicit HypreDiagScale( HypreParMatrix& A ) : HypreSolver( &A ) {}
00278     virtual operator HYPRE_Solver() const
00279     {
00280         return NULL;
00281     }
00282 
00283     virtual HYPRE_PtrToParSolverFcn SetupFcn() const
00284     {
00285         return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRDiagScaleSetup;
00286     }
00287     virtual HYPRE_PtrToParSolverFcn SolveFcn() const
00288     {
00289         return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRDiagScale;
00290     }
00291 
00292     virtual void GetNumIterations( int& num_iterations )
00293     {
00294         num_iterations = 1;
00295     }
00296 
00297     virtual void GetFinalResidualNorm( double& normr )
00298     {
00299         normr = 0.0;
00300     }
00301 
00302     HypreParMatrix* GetData()
00303     {
00304         return A;
00305     }
00306     virtual ~HypreDiagScale() {}
00307 };
00308 
00309 /// The ParaSails preconditioner in hypre
00310 class HypreParaSails : public HypreSolver
00311 {
00312   private:
00313     HYPRE_Solver sai_precond;
00314 
00315   public:
00316     HypreParaSails( HypreParMatrix& A );
00317 
00318     void SetSymmetry( int sym );
00319 
00320     virtual void GetNumIterations( int& num_iterations )
00321     {
00322         num_iterations = 1;
00323     }
00324 
00325     virtual void GetFinalResidualNorm( double& normr )
00326     {
00327         normr = 0.0;
00328     }
00329 
00330     /// The typecast to HYPRE_Solver returns the internal sai_precond
00331     virtual operator HYPRE_Solver() const
00332     {
00333         return sai_precond;
00334     }
00335 
00336     virtual HYPRE_PtrToParSolverFcn SetupFcn() const
00337     {
00338         return (HYPRE_PtrToParSolverFcn)HYPRE_ParaSailsSetup;
00339     }
00340     virtual HYPRE_PtrToParSolverFcn SolveFcn() const
00341     {
00342         return (HYPRE_PtrToParSolverFcn)HYPRE_ParaSailsSolve;
00343     }
00344 
00345     virtual ~HypreParaSails();
00346 };
00347 
00348 /// The BoomerAMG solver in hypre
00349 class HypreBoomerAMG : public HypreSolver
00350 {
00351   private:
00352     HYPRE_Solver amg_precond;
00353 
00354     /// Default, generally robust, BoomerAMG options
00355     void SetDefaultOptions();
00356 
00357     // If amg_precond is NULL, allocates it and sets default options.
00358     // Otherwise saves the options from amg_precond, destroys it, allocates a new
00359     // one, and sets its options to the saved values.
00360     void ResetAMGPrecond();
00361 
00362   public:
00363     HypreBoomerAMG();
00364 
00365     HypreBoomerAMG( HypreParMatrix& A );
00366 
00367     virtual void SetOperator( const HypreParMatrix& op );
00368 
00369     virtual void GetNumIterations( int& num_iterations )
00370     {
00371         num_iterations = 1;
00372     }
00373 
00374     virtual void GetFinalResidualNorm( double& normr )
00375     {
00376         normr = 0.0;
00377     }
00378 
00379     /** More robust options for systems, such as elasticity. Note that BoomerAMG
00380         assumes Ordering::byVDIM in the finite element space used to generate the
00381         matrix A. */
00382     void SetSystemsOptions( int dim );
00383 
00384     virtual void Verbosity( int print_level )
00385     {
00386         HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
00387     }
00388 
00389     /// The typecast to HYPRE_Solver returns the internal amg_precond
00390     virtual operator HYPRE_Solver() const
00391     {
00392         return amg_precond;
00393     }
00394 
00395     virtual HYPRE_PtrToParSolverFcn SetupFcn() const
00396     {
00397         return (HYPRE_PtrToParSolverFcn)HYPRE_BoomerAMGSetup;
00398     }
00399     virtual HYPRE_PtrToParSolverFcn SolveFcn() const
00400     {
00401         return (HYPRE_PtrToParSolverFcn)HYPRE_BoomerAMGSolve;
00402     }
00403 
00404     virtual ~HypreBoomerAMG();
00405 };
00406 
00407 }  // namespace moab
00408 
00409 #endif  // MOAB_HAVE_MPI
00410 
00411 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines