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 

// 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