Ifpack_Hypre.h

/*@HEADER
// ***********************************************************************
//
//       Ifpack: Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2002) 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 IFPACK_HYPRE_H
#define IFPACK_HYPRE_H
 
#include "Ifpack_ConfigDefs.h"
#ifdef HAVE_HYPRE
 
#include "Ifpack_Preconditioner.h"
#include "Ifpack_Condest.h"
#include "Ifpack_ScalingType.h"
#include "Epetra_CompObject.h"
#include "Epetra_MultiVector.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsGraph.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_BlockMap.h"
#include "Epetra_Map.h"
#include "Epetra_Object.h"
#include "Epetra_Comm.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Time.h"
#include "Teuchos_RefCountPtr.hpp"
#include "Teuchos_ArrayRCP.hpp"
#include "Epetra_MpiComm.h"
 
 
#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;
 
// Will only work if Hypre is built with HYPRE_BIGINT=OFF
typedef int HYPRE_Int;
 
typedef HYPRE_Int (*HYPRE_PtrToParSolverFcn)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
 
#ifndef HYPRE_ENUMS
#define HYPRE_ENUMS
enum Hypre_Solver{
    BoomerAMG,
    ParaSails,
    Euclid,
    AMS,
    Hybrid,
    PCG,
    GMRES,
    FlexGMRES,
    LGMRES,
    BiCGSTAB
};

enum Hypre_Chooser{
    Solver,
    Preconditioner
};
#endif //HYPRE_ENUMS
 
class FunctionParameter;
 
namespace Teuchos {
  class ParameterList;
}


 
class Ifpack_Hypre: public Ifpack_Preconditioner {
 
public:
  // @{ Constructors and destructors.
  Ifpack_Hypre(Epetra_RowMatrix* A);

  ~Ifpack_Hypre(){ Destroy();}
 
  // @}
  // @{ Construction methods

  int Initialize();

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


  int 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.
  */
  int SetParameters(Teuchos::ParameterList& parameterlist);


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


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


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


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


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


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


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


    int SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int**), int** 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<Epetra_MultiVector> coords);

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

    int CallFunctions() const;


  int SetUseTranspose(bool UseTranspose_in) {UseTranspose_ = UseTranspose_in; return(0);};
 
  // @}
 
  // @{ Mathematical functions.
  // Applies the matrix to X, returns the result in Y.
  int Apply(const Epetra_MultiVector& X,
               Epetra_MultiVector& Y) const{ return(Multiply(false,X,Y));}


  int Multiply(bool Trans, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;


  int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;

  double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
                 const int MaxIters = 1550,
                 const double Tol = 1e-9,
                 Epetra_RowMatrix* Matrix_in = 0);

  double Condest() const{ return(Condest_);}
 
  // @}
  // @{ Query methods

  const char* Label() const {return(Label_);}

  int SetLabel(const char* Label_in)
  {
    strcpy(Label_,Label_in);
    return(0);
  }

  const Epetra_Map& OperatorDomainMap() const{ return *GloballyContiguousRowMap_;}

  const Epetra_Map& OperatorRangeMap() const{ return *GloballyContiguousRowMap_;}

  double NormInf() const {return(0.0);};

  bool HasNormInf() const {return(false);};

  bool UseTranspose() const {return(UseTranspose_);};

  const Epetra_Comm & Comm() const{return(A_->Comm());};

  const Epetra_RowMatrix& Matrix() const{ return(*A_);}

#if 0
  const HYPRE_IJMatrix& HypreMatrix()
  {
    if(IsInitialized() == false)
      Initialize();
    return(HypreA_);
  }
#endif

  virtual std::ostream& Print(std::ostream& os) const;

  virtual int NumInitialize() const{ return(NumInitialize_);}

  virtual int NumCompute() const{ return(NumCompute_);}

  virtual int NumApplyInverse() const{ return(NumApplyInverse_);}

  virtual double InitializeTime() const{ return(InitializeTime_);}

  virtual double ComputeTime() const{ return(ComputeTime_);}

  virtual double ApplyInverseTime() const{ return(ApplyInverseTime_);}

  virtual double InitializeFlops() const{ return(0.0);}

  virtual double ComputeFlops() const{ return(ComputeFlops_);}

  virtual double ApplyInverseFlops() const{ return(ApplyInverseFlops_);}
 
private:
 
  // @}
  // @{ Private methods

  Ifpack_Hypre(const Ifpack_Hypre& RHS) : Time_(RHS.Comm()){}

  Ifpack_Hypre& operator=(const Ifpack_Hypre& /*RHS*/){ return(*this);}

  void Destroy();

  MPI_Comm GetMpiComm() const
    { return (dynamic_cast<const Epetra_MpiComm*>(&A_->Comm()))->GetMpiComm();}


  int Solve(bool Trans, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
 

  int NumGlobalRows() const {return(A_->NumGlobalRows());};

  int NumGlobalCols() const {return(A_->NumGlobalCols());};

  int NumMyRows() const {return(A_->NumMyRows());};

  int NumMyCols() const {return(A_->NumMyCols());};

  int SetSolverType(Hypre_Solver solver);

  int SetPrecondType(Hypre_Solver precond);

  int CreateSolver();

  int CreatePrecond();

  int CopyEpetraToHypre();

  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 Epetra_Map> MakeContiguousColumnMap(Teuchos::RCP<const Epetra_RowMatrix> &Matrix) const;
  // @}
  // @{ Internal data

  Teuchos::RCP<Epetra_RowMatrix> A_;
  Teuchos::ParameterList List_;
  bool UseTranspose_;
  double Condest_;
  bool IsInitialized_;
  bool IsComputed_;
  char Label_[160];
  int NumInitialize_;
  int NumCompute_;
  mutable int NumApplyInverse_;
  double InitializeTime_;
  double ComputeTime_;
  mutable double ApplyInverseTime_;
  double ComputeFlops_;
  mutable double ApplyInverseFlops_;
  mutable Epetra_Time Time_;

  mutable HYPRE_IJMatrix HypreA_;
  mutable HYPRE_ParCSRMatrix ParMatrix_;

  Teuchos::RCP<const Epetra_CrsMatrix> 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<Epetra_MultiVector> 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 (Ifpack_Hypre::*SolverCreatePtr_)(MPI_Comm, HYPRE_Solver*);
  int (*SolverDestroyPtr_)(HYPRE_Solver);
  int (*SolverSetupPtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
  int (*SolverSolvePtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
  int (*SolverPrecondPtr_)(HYPRE_Solver, HYPRE_PtrToParSolverFcn, HYPRE_PtrToParSolverFcn, HYPRE_Solver);
  int (Ifpack_Hypre::*PrecondCreatePtr_)(MPI_Comm, HYPRE_Solver*);
  int (*PrecondDestroyPtr_)(HYPRE_Solver);
  int (*PrecondSetupPtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
  int (*PrecondSolvePtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
 
  bool IsSolverCreated_;
  bool IsPrecondCreated_;
  Hypre_Chooser SolveOrPrec_;
  Teuchos::RCP<const Epetra_Map> GloballyContiguousRowMap_;
  Teuchos::RCP<const Epetra_Map> GloballyContiguousColMap_;
  Teuchos::RCP<const Epetra_Map> GloballyContiguousNodeRowMap_;
  Teuchos::RCP<const Epetra_Map> GloballyContiguousNodeColMap_;
  int NumFunsToCall_;
  Hypre_Solver SolverType_;
  Hypre_Solver PrecondType_;
  bool UsePreconditioner_;
  std::vector<Teuchos::RCP<FunctionParameter> > FunsToCall_;
  bool Dump_;
  mutable Teuchos::ArrayRCP<double> VectorCache_;
 
};
 
#endif // HAVE_HYPRE
#endif /* IFPACK_HYPRE_H */