Ifpack2_Hypre_decl.hpp

/*@HEADER
// ***********************************************************************
//
//       Ifpack2:  Templated Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2009) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//@HEADER
*/
 
#ifndef IFPACK2_HYPRE_DECL_HPP
#define IFPACK2_HYPRE_DECL_HPP
 
#include "Ifpack2_ConfigDefs.hpp"
#if defined(HAVE_IFPACK2_HYPRE) && defined(HAVE_IFPACK2_MPI)
 
#include "Ifpack2_Preconditioner.hpp"
#include "Ifpack2_Details_CanChangeMatrix.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_Vector.hpp"
#include "Tpetra_CrsGraph.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_CrsMatrix.hpp"
 
#include "Teuchos_RefCountPtr.hpp"
#include "Teuchos_ArrayRCP.hpp"
 
#include "Ifpack2_Hypre_FunctionParameters.hpp"
 
#include <map>
 
// Hypre forward declarations (to avoid downstream header pollution)
struct hypre_IJMatrix_struct;
typedef struct hypre_IJMatrix_struct *HYPRE_IJMatrix;
struct hypre_IJVector_struct;
typedef struct hypre_IJVector_struct *HYPRE_IJVector;
struct hypre_ParCSRMatrix_struct;
typedef struct hypre_ParCSRMatrix_struct* HYPRE_ParCSRMatrix;
struct hypre_ParVector_struct;
typedef struct hypre_ParVector_struct * HYPRE_ParVector;
struct hypre_Solver_struct;
typedef struct hypre_Solver_struct *HYPRE_Solver;
struct hypre_ParVector_struct;
typedef struct hypre_ParVector_struct hypre_ParVector;
//struct hypre_Vector;
 
namespace Ifpack2 {
 
 
#ifndef HYPRE_ENUMS
#define HYPRE_ENUMS
  enum Hypre_Solver{
    BoomerAMG,
    ParaSails,
    Euclid,
    AMS,
    Hybrid,
    PCG,
    GMRES,
    FlexGMRES,
    LGMRES,
    BiCGSTAB
  };
  
  enum Hypre_Chooser{
    Hypre_Is_Solver,
    Hypre_Is_Preconditioner
  };
#endif //HYPRE_ENUMS
 


 
template<class MatrixType>
class Hypre: 
    virtual public Ifpack2::Preconditioner<typename MatrixType::scalar_type,
                                           typename MatrixType::local_ordinal_type,
                                           typename MatrixType::global_ordinal_type,
                                           typename MatrixType::node_type>,
    virtual public Ifpack2::Details::CanChangeMatrix<Tpetra::RowMatrix<typename MatrixType::scalar_type,
                                                                       typename MatrixType::local_ordinal_type,
                                                                       typename MatrixType::global_ordinal_type,
                                                                       typename MatrixType::node_type> >
{
public:
public:


  typedef MatrixType matrix_type;

  typedef typename MatrixType::scalar_type scalar_type;

  typedef typename MatrixType::local_ordinal_type local_ordinal_type;

  typedef typename MatrixType::global_ordinal_type global_ordinal_type;

  typedef typename MatrixType::node_type::device_type device_type;

  typedef typename MatrixType::node_type node_type;

  typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;

  typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type,
                            global_ordinal_type, node_type> row_matrix_type;
 
  static_assert (std::is_same<MatrixType, row_matrix_type>::value,
                 "Ifpack2::Hypre: MatrixType must be a Tpetra::RowMatrix "
                 "specialization.  Don't use Tpetra::CrsMatrix here.");
 

  typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
                            global_ordinal_type, node_type> crs_matrix_type;

  typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;

  typedef Tpetra::Vector<scalar_type, local_ordinal_type,
                         global_ordinal_type, node_type> vector_type;

  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
                         global_ordinal_type, node_type> multivector_type;
  
  // Hypre Specs
  // This will need to be either int or long long depending on how Hypre was built
  //    typedef global_ordinal_type global_ordinal_type;
 
  typedef global_ordinal_type (*HYPRE_PtrToParSolverFcn)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);

  // \name Constructors and destructors

  explicit Hypre(const Teuchos::RCP<const row_matrix_type>& A);

  ~Hypre();
 
  // @}
  // @{ Construction methods
  void initialize();

  bool isInitialized() const{ return(IsInitialized_);}


  void compute();

  bool isComputed() const{ return(IsComputed_);}
 
  /* This method is only available if the Teuchos package is enabled.
     This method recognizes six parameter names: Solver,
     Preconditioner, SolveOrPrecondition, SetPreconditioner, NumFunctions and Functions. These names are
     case sensitive. Solver requires an enumerated parameter of type Hypre_Solver. Preconditioner is similar
     except requires the type be a preconditioner. The options are listed below:
                       Solvers                            Preconditioners
                       BoomerAMG                          BoomerAMG
                       AMS                                ParaSails
                       Hybrid                             AMS
                       PCG (Default)                      Euclid (Default)
                       GMRES
                       FlexGMRES
                       LGMRES
                       BiCGSTAB
     SolveOrPrecondition takes enumerated type Hypre_Chooser, Solver will solve the system, Preconditioner will apply the preconditioner.
     SetPreconditioner takes a boolean, true means the solver will use the preconditioner.
     NumFunctions takes an int that describes how many parameters will be passed into Functions. (This needs to be correct.)
     Functions takes an array of Ref Counted Pointers to an object called FunctionParameter. This class is implemented in Ifpack_Hypre.h.
     The object takes whether it is Solver or Preconditioner that we are setting a parameter for.
     The function in Hypre that sets the parameter, and the parameters for that function. An example is below:
 
     RCP<FunctionParameter> functs[2];
     functs[0] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetMaxIter, 1000)); // max iterations
     functs[1] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTol, 1e-7)); // conv. tolerance
     list.set("NumFunctions", 2);
     list.set<RCP<FunctionParameter>*>("Functions", functs);
     NOTE: SetParameters() must be called to use ApplyInverse(), the solvers will not be created otherwise. An empty list is acceptable to use defaults.
  */
  void setParameters(const Teuchos::ParameterList& parameterlist);


    int SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type), global_ordinal_type parameter);


    int SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, double), double parameter);


    int SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, double, global_ordinal_type), double parameter1, global_ordinal_type parameter2);


    int SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type, double), global_ordinal_type parameter1, double parameter2);


    int SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type, global_ordinal_type), global_ordinal_type parameter1, global_ordinal_type parameter2);


    int SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, double*), double* parameter);


    int SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type*), global_ordinal_type* parameter);


    int SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type**), global_ordinal_type** parameter);


 
    int SetParameter(Hypre_Chooser chooser, Hypre_Solver Solver);


 
    int SetParameter(bool UsePreconditioner){ UsePreconditioner_ = UsePreconditioner; return 0;}


    int SetParameter(Hypre_Chooser chooser) { SolveOrPrec_ = chooser; return 0;}

  int SetCoordinates(Teuchos::RCP<multivector_type> coords);

  int SetDiscreteGradient(Teuchos::RCP<const crs_matrix_type> G);

    int CallFunctions() const;



  virtual void
  setMatrix (const Teuchos::RCP<const row_matrix_type>& A);

  void
  apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
         Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
         Teuchos::ETransp mode = Teuchos::NO_TRANS,
         scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
         scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero()) const;

  Teuchos::RCP<const map_type> getDomainMap() const;

  Teuchos::RCP<const map_type> getRangeMap() const;

  bool hasTransposeApply() const;

  void
  applyMat (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
            Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
            Teuchos::ETransp mode = Teuchos::NO_TRANS) const;



  Teuchos::RCP<const Teuchos::Comm<int> > getComm() const;

  Teuchos::RCP<const row_matrix_type> getMatrix() const;

  Teuchos::RCP<const Tpetra::CrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> >
  getCrsMatrix() const;

  int getNumInitialize() const;

  int getNumCompute() const;

  int getNumApply() const;

  double getInitializeTime() const;

  double getComputeTime() const;

  double getApplyTime() const;



  std::string description() const;

  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const;

 
private:
  // @{ Private methods

  typedef Teuchos::ScalarTraits<typename MatrixType::scalar_type> STS;

  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> MV;

  Hypre(const Hypre<MatrixType>&);

  Hypre<MatrixType>& operator= (const Hypre<MatrixType>&);

  int SetSolverType(Hypre_Solver solver);

  int SetPrecondType(Hypre_Solver precond);

  int CreateSolver();

  int CreatePrecond();

  int CopyTpetraToHypre();

  int AddFunToList(Teuchos::RCP<FunctionParameter> NewFun);

  int Hypre_BoomerAMGCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_AMSCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_ParCSRHybridCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver);

  int Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver);
  
  Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > 
  MakeContiguousColumnMap(Teuchos::RCP<const crs_matrix_type> &Matrix) const;

  void Destroy();
 
  // @}
  // @{ Internal data
  Teuchos::RCP<const MatrixType> A_;
  Teuchos::ParameterList List_;

  bool IsInitialized_;
  bool IsComputed_;
  int NumInitialize_;
  int NumCompute_;
  mutable int NumApply_;
  double InitializeTime_;
  double ComputeTime_;
  mutable double ApplyTime_;
  double ComputeFlops_;
  mutable double ApplyFlops_;  
 

  mutable HYPRE_IJMatrix HypreA_;
  mutable HYPRE_ParCSRMatrix ParMatrix_;

  Teuchos::RCP<const crs_matrix_type> G_;
  mutable HYPRE_IJMatrix HypreG_;
  mutable HYPRE_ParCSRMatrix ParMatrixG_;

  mutable HYPRE_IJVector XHypre_;
  mutable HYPRE_IJVector YHypre_;
  mutable HYPRE_ParVector ParX_;
  mutable HYPRE_ParVector ParY_;
  mutable Teuchos::RCP<hypre_ParVector> XVec_;
  mutable Teuchos::RCP<hypre_ParVector> YVec_;
 
  Teuchos::RCP<multivector_type> Coords_;
  mutable HYPRE_IJVector xHypre_;
  mutable HYPRE_IJVector yHypre_;
  mutable HYPRE_IJVector zHypre_;
  mutable HYPRE_ParVector xPar_;
  mutable HYPRE_ParVector yPar_;
  mutable HYPRE_ParVector zPar_;

  mutable HYPRE_Solver Solver_;
  mutable HYPRE_Solver Preconditioner_;
  //  The following are pointers to functions to use the solver and preconditioner.
  int (Hypre::*SolverCreatePtr_)(MPI_Comm, HYPRE_Solver*);
  global_ordinal_type (*SolverDestroyPtr_)(HYPRE_Solver);
  global_ordinal_type (*SolverSetupPtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
  global_ordinal_type (*SolverSolvePtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
  global_ordinal_type (*SolverPrecondPtr_)(HYPRE_Solver, HYPRE_PtrToParSolverFcn, HYPRE_PtrToParSolverFcn, HYPRE_Solver);
  int (Hypre::*PrecondCreatePtr_)(MPI_Comm, HYPRE_Solver*);
  global_ordinal_type (*PrecondDestroyPtr_)(HYPRE_Solver);
  global_ordinal_type (*PrecondSetupPtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
  global_ordinal_type (*PrecondSolvePtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
 
  bool IsSolverCreated_;
  bool IsPrecondCreated_;
  Hypre_Chooser SolveOrPrec_;
  Teuchos::RCP<const map_type> GloballyContiguousRowMap_;
  Teuchos::RCP<const map_type> GloballyContiguousColMap_;
  Teuchos::RCP<const map_type> GloballyContiguousNodeRowMap_;
  Teuchos::RCP<const map_type> GloballyContiguousNodeColMap_;
  int NumFunsToCall_;
  Hypre_Solver SolverType_;
  Hypre_Solver PrecondType_;
  bool UsePreconditioner_;
  std::vector<Teuchos::RCP<FunctionParameter> > FunsToCall_;
  bool Dump_;
  mutable Teuchos::ArrayRCP<double> VectorCache_;
 
};
 
 
}//end Ifpack2 namespace
 
#endif // HAVE_IFPACK2_HYPRE && HAVE_IFPACK2_MPI
#endif /* IFPACK2_HYPRE_DECL_HPP */