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 "HYPRE_IJ_mv.h"
#include "HYPRE_parcsr_ls.h"
#include "krylov.h"
#include "_hypre_parcsr_mv.h"
#include "_hypre_IJ_mv.h"
#include "HYPRE_parcsr_mv.h"
#include "HYPRE.h"
#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 "Epetra_MpiComm.h"

#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{
  public:
    FunctionParameter(Hypre_Chooser chooser, int (*funct_name)(HYPRE_Solver, int), int param1) :
      chooser_(chooser),
      option_(0),
      int_func_(funct_name),
      int_param1_(param1) {}

    FunctionParameter(Hypre_Chooser chooser, int (*funct_name)(HYPRE_Solver, double), double param1):
      chooser_(chooser),
      option_(1),
      double_func_(funct_name),
      double_param1_(param1) {}

    FunctionParameter(Hypre_Chooser chooser, int (*funct_name)(HYPRE_Solver, double, int), double param1, int param2):
      chooser_(chooser),
      option_(2),
      double_int_func_(funct_name),
      int_param1_(param2),
      double_param1_(param1) {}

    FunctionParameter(Hypre_Chooser chooser, int (*funct_name)(HYPRE_Solver, int, int), int param1, int param2):
      chooser_(chooser),
      option_(3),
      int_int_func_(funct_name),
      int_param1_(param1),
      int_param2_(param2) {}

    FunctionParameter(Hypre_Chooser chooser, int (*funct_name)(HYPRE_Solver, int*), int *param1):
      chooser_(chooser),
      option_(4),
      int_star_func_(funct_name),
      int_star_param_(param1) {}

    FunctionParameter(Hypre_Chooser chooser, int (*funct_name)(HYPRE_Solver, double*), double* param1):
      chooser_(chooser),
      option_(5),
      double_star_func_(funct_name),
      double_star_param_(param1) {}

    int CallFunction(HYPRE_Solver solver, HYPRE_Solver precond){
      if(chooser_ == Solver){
        if(option_ == 0){
          return int_func_(solver, int_param1_);
        } else if(option_ == 1){
          return double_func_(solver, double_param1_);
        } else if(option_ == 2){
          return double_int_func_(solver, double_param1_, int_param1_);
        } else if (option_ == 3){
          return int_int_func_(solver, int_param1_, int_param2_);
        } else if (option_ == 4){
          return int_star_func_(solver, int_star_param_);
        } else {
          return double_star_func_(solver, double_star_param_);
        }
      } else {
        if(option_ == 0){
          return int_func_(precond, int_param1_);
        } else if(option_ == 1){
          return double_func_(precond, double_param1_);
        } else if(option_ == 2){
          return double_int_func_(precond, double_param1_, int_param1_);
        } else if(option_ == 3) {
          return int_int_func_(precond, int_param1_, int_param2_);
        } else if(option_ == 4) {
          return int_star_func_(precond, int_star_param_);
        } else {
          return double_star_func_(precond, double_star_param_);
        }
      }
    }

  private:
    Hypre_Chooser chooser_;
    int option_;
    int (*int_func_)(HYPRE_Solver, int);
    int (*double_func_)(HYPRE_Solver, double);
    int (*double_int_func_)(HYPRE_Solver, double, int);
    int (*int_int_func_)(HYPRE_Solver, int, int);
    int (*int_star_func_)(HYPRE_Solver, int*);
    int (*double_star_func_)(HYPRE_Solver, double*);
    int int_param1_;
    int int_param2_;
    double double_param1_;
    int *int_star_param_;
    double *double_star_param_;
};

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, 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, Hypre_Solver Solver);


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


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

    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 *MySimpleMap_;}

  const Epetra_Map& OperatorRangeMap() const{ return *MySimpleMap_;}
  
  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_);}

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

  virtual ostream& Print(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 AddFunToList(Teuchos::RCP<FunctionParameter> NewFun);

  int Hypre_BoomerAMGCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_BoomerAMGCreate(solver);}

  int Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_ParaSailsCreate(comm, solver);}

  int Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_EuclidCreate(comm, solver);}

  int Hypre_AMSCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_AMSCreate(solver);}

  int Hypre_ParCSRHybridCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_ParCSRHybridCreate(solver);}

  int Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_ParCSRPCGCreate(comm, solver);}

  int Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_ParCSRGMRESCreate(comm, solver);}

  int Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_ParCSRFlexGMRESCreate(comm, solver);}

  int Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_ParCSRLGMRESCreate(comm, solver);}

  int Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
    { return HYPRE_ParCSRBiCGSTABCreate(comm, solver);}

  // @}
  // @{ Internal data
  
  Teuchos::RefCountPtr<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_;
  mutable HYPRE_IJVector XHypre_;
  mutable HYPRE_IJVector YHypre_;
  mutable HYPRE_ParVector ParX_;
  mutable HYPRE_ParVector ParY_;
  mutable hypre_ParVector *XVec_;
  mutable hypre_ParVector *YVec_;
  mutable hypre_Vector *XLocal_;
  mutable hypre_Vector *YLocal_;
  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 *IsSolverSetup_;
  bool *IsPrecondSetup_;
  Hypre_Chooser SolveOrPrec_;
  Teuchos::RefCountPtr<Epetra_Map> MySimpleMap_;
  int NumFunsToCall_;
  Hypre_Solver SolverType_;
  Hypre_Solver PrecondType_;
  bool UsePreconditioner_;
  std::vector<Teuchos::RCP<FunctionParameter> > FunsToCall_;
  bool NiceRowMap_;
};

#endif // HAVE_HYPRE
#endif /* IFPACK_HYPRE_H */