belosSolver.hh

Go to the documentation of this file.

 
#ifndef BELOSSOLVER_HH_
#define BELOSSOLVER_HH_
 
#include "belosLinProb.hh"
#include "sparseMatrix.hh"
#include <BelosSolverManager.hpp>
#include <BelosTpetraAdapter.hpp>
#include <Ifpack2_Preconditioner.hpp>
#include <Tpetra_DefaultPlatform.hpp>
#include <Tpetra_MultiVector.hpp>
#include <Teuchos_RCP.hpp>
#include <Tpetra_CrsMatrix.hpp>
 
namespace concepts {
 
  // *********************************************************** BelosSolver **
 
  template<class T>
  class BelosSolver: public VecOperator<T> {
  public:
 
    //typenames for better reading
    typedef Tpetra::MultiVector<T, int> MV;
    typedef Tpetra::Operator<T> OP;
    typedef Tpetra::CrsMatrix<T, int, int> CRSmat;
    typedef BelosLinProb<T, MV, OP> LP;
    typedef Belos::SolverManager<T, MV, OP> solverManager;
 
    BelosSolver(concepts::SparseMatrix<T>& sparse,
        Teuchos::RCP<const Teuchos::Comm<int> > comm);
 
    BelosSolver(Teuchos::RCP<const Teuchos::Comm<int> > comm);
    inline virtual ~BelosSolver() {
    }
    ;
 
    inline virtual void operator()() {
      apply_();
    }
    ;
 
    inline virtual void operator()(const Vector<T>& fncY, Vector<T>& fncX) {
      apply_(fncY, fncX);
    }
    ;
 
    inline void setSolverParam(Teuchos::RCP<Teuchos::ParameterList> param) {
      solverParam_ = param;
    }
    ;
 
    inline void setSolverParam(int maxIter, int maxRestarts, double tol = 1e-10,
        int blockSize = 1, int outFreq = 1) {
      int verbosity = Belos::Warnings + Belos::Errors + Belos::FinalSummary;
      solverParam_ = Teuchos::RCP<Teuchos::ParameterList>(
          new Teuchos::ParameterList());
      solverParam_->set("Verbosity", verbosity);
      solverParam_->set("Block Size", blockSize);
      solverParam_->set("Maximum Iterations", maxIter);
//      if (solverType_ != "minres")
//        solverParam_->set("Maximum Restarts", maxRestarts);
      solverParam_->set("Convergence Tolerance", tol);
      solverParam_->set("Output Frequency", outFreq);
      solverParam_->set("Output Stream",
          Teuchos::rcpFromRef<std::ostream>(std::cout));
    }
    ;
 
    inline Teuchos::RCP<Teuchos::ParameterList> getSolverParam() {
      return solverParam_;
    }
    ;
 
    inline void setPrecParam(Teuchos::RCP<Teuchos::ParameterList> param) {
      precParam_ = param;
    }
    ;
 
    //TODO right expl. for fill??
    inline void setPrecParam(double fill, double tol = 1e-5, int absTresh = 0,
        int relTresh = 1, double relax = 0) {
      precParam_ = Teuchos::RCP<Teuchos::ParameterList>(
          new Teuchos::ParameterList("IFPACK2"));
      precParam_->set("fact: level-of-fill", fill);
      precParam_->set("fact: ilut level-of-fill", fill);
      precParam_->set("fact: riluk level-of-fill", fill);
      precParam_->set("fact: absolute threshold", absTresh);
      precParam_->set("fact: relative threshold", relTresh);
      precParam_->set("fact: relax value", relax);
      precParam_->set("fact: drop tolerance", tol);
    }
    ;
 
    inline Teuchos::RCP<Teuchos::ParameterList> getPrecParam() {
      return precParam_;
    }
    ;
 
    inline void setSolverManager(Teuchos::RCP<solverManager> manager) {
      solverType_ = "";
      solverManager_ = manager;
    }
    ;
 
    inline Teuchos::RCP<solverManager> getSolverManager() {
      return solverManager_;
    }
    ;
 
    inline Teuchos::RCP<LP> getLinearProblem() {
      return lp_;
    }
 
    inline void setSolverType(std::string type) {
      solverType_ = type;
    }
    ;
 
    inline std::string getSolverType() {
      return solverType_;
    }
    ;
 
    inline void setPrecType(std::string type) {
      precType_ = type;
    }
    ;
 
    //Getter for the preconditioner type
    inline std::string getPrecType() {
      return precType_;
    }
    ;
 
    inline void prepare() {
      //type and param should be setted
      createPrec_();
 
      //type and param should be setted
      createSolver_();
 
 
      //precond
      phasePrecInit_((Comm_->getRank() == 0));
 
      prepared_ = true;
    }
 
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
 
    virtual void apply_(const Vector<T>& fncY, Vector<T>& fncX);
 
    virtual void apply_();
 
    bool createSolver_();
 
    bool createPrec_();
 
    bool phasePrecInit_(bool verbose = true);
 
    bool phaseSolve_(bool verbose = true);
 
    Teuchos::RCP<const Teuchos::Comm<int> > Comm_;
 
    Teuchos::RCP<LP> lp_;
 
    Teuchos::RCP<Teuchos::ParameterList> solverParam_;
    Teuchos::RCP<Teuchos::ParameterList> precParam_;
 
    Teuchos::RCP<solverManager> solverManager_;
 
    Teuchos::RCP<Ifpack2::Preconditioner<T, int> > prec_;
 
    //type of solver/preconditioner as string
    std::string solverType_;
    std::string precType_;
 
    bool prepared_;
  }
  ;
 
  //Constructors
  template<class T>
  BelosSolver<T>::BelosSolver(concepts::SparseMatrix<T>& sparse,
      Teuchos::RCP<const Teuchos::Comm<int> > Comm) :
      VecOperator<T>(sparse.dimY(), sparse.dimX()) {
    Comm_ = Comm;
    prepared_ = false;
    lp_ = Teuchos::RCP<LP>(new LP(sparse, Comm));
  }
  ;
 
  template<class T>
  BelosSolver<T>::BelosSolver(Teuchos::RCP<const Teuchos::Comm<int> > Comm) :
      VecOperator<T>(1, 1) {
    Comm_ = Comm;
    prepared_ = false;
    lp_ = Teuchos::RCP<LP>(new LP(Comm));
  }
  ;
 
} /* namespace concepts */
#endif /* BELOSSOLVER_HH_ */