belosSolver_decl.hh

Go to the documentation of this file.

 
#ifndef BELOSSOLVER_DECL_HH_
#define BELOSSOLVER_DECL_HH_
 
#include "belosSolver.hh"
#include <BelosLinearProblem.hpp>
#include <BelosBlockGmresSolMgr.hpp>
#include <BelosPseudoBlockCGSolMgr.hpp>
#include <BelosBlockCGSolMgr.hpp>
#include <BelosMinresSolMgr.hpp>
#include <BelosPseudoBlockGmresSolMgr.hpp>
#include <Teuchos_ArrayView.hpp>
#include <Ifpack2_Factory.hpp>
#include "ctype.h"
 
namespace concepts {
 
  template<class T>
  Teuchos::RCP<
      Belos::SolverManager<T, Tpetra::MultiVector<T, int>, Tpetra::Operator<T> > > createBelosSolverMgr(
      std::string solverType,
      const Teuchos::RCP<Teuchos::ParameterList>& solverParam,
      Teuchos::RCP<
          Belos::LinearProblem<T, Tpetra::MultiVector<T, int>,
              Tpetra::Operator<T> > > linearProblem)
 
      {
    typedef Tpetra::MultiVector<T, int> MV;
    typedef Tpetra::Operator<T> OP;
 
    std::transform(solverType.begin(), solverType.end(), solverType.begin(),
        ::tolower);
 
    if (solverType == "blockgmres") {
      return Teuchos::RCP<Belos::BlockGmresSolMgr<T, MV, OP> >(
          new Belos::BlockGmresSolMgr<T, MV, OP>(linearProblem, solverParam));
    } else if (solverType == "gmres") {
      return Teuchos::RCP<Belos::PseudoBlockGmresSolMgr<T, MV, OP> >(
          new Belos::PseudoBlockGmresSolMgr<T, MV, OP>(linearProblem,
              solverParam));
 
//    } else if (solverType == "minres") {
//      return Teuchos::RCP<Belos::MinresSolMgr<T, MV, OP> >(
//          new Belos::MinresSolMgr<T, MV, OP>(linearProblem, solverParam));
//
    } else if (solverType == "blockcg") {
      return Teuchos::rcp(
          new Belos::BlockCGSolMgr<T, MV, OP>(linearProblem, solverParam));
    }
 
    else if (solverType == "pseudocg") {
      return Teuchos::rcp(
          new Belos::PseudoBlockCGSolMgr<T, MV, OP>(linearProblem, solverParam));
    } else {
      std::cerr << "ERROR: solver type '" << solverType << "' unknown."
          << std::endl;
      return Teuchos::null;
    }
 
  }
 
  template<class T>
  Teuchos::RCP<Ifpack2::Preconditioner<T> > createIfpackPrec(
      std::string precType, Teuchos::RCP<Teuchos::ParameterList> precParam,
      const Teuchos::RCP<const Tpetra::CrsMatrix<T, int> > A) {
 
    Ifpack2::Factory factory;
    try {
      // Set up the preconditioner of that type.
      typename Teuchos::RCP<Ifpack2::Preconditioner<T, int> > prec =
          factory.create(precType, A);
      prec->setParameters(*precParam);
      return prec;
    } catch (std::exception &e) {
      std::cerr << "ERROR: preconditioner type '" << precType << "' unknown."
          << std::endl;
    }
    return Teuchos::null;
 
  }
 
  template<class T>
  bool BelosSolver<T>::createPrec_() {
    prec_ = createIfpackPrec<T>(precType_, precParam_, lp_->getCrsMat());
    return !Teuchos::is_null(prec_);
  }
 
  template<class T>
  bool BelosSolver<T>::createSolver_() {
    solverManager_ = createBelosSolverMgr(solverType_, solverParam_,
        Teuchos::RCP<Belos::LinearProblem<T, MV, OP> >(lp_));
    solverManager_->setProblem(lp_);
    return !Teuchos::is_null(solverManager_);
  }
 
  template<class T>
  bool BelosSolver<T>::phasePrecInit_(bool verbose) {
    Teuchos::Time timer("precond");
//    prec_->setParam(precParam);
    prec_->initialize();
    timer.start();
    prec_->compute();
    timer.stop();
    lp_->setLeftPrec(prec_);
//    lp_->setRightPrec(prec_);
    if (verbose) {
      std::cout << "... Finished Computing Ifpack2 preconditioner (time: "
          << timer.totalElapsedTime() << "s)" << std::endl;
    }
 
    typename Teuchos::ScalarTraits<T>::magnitudeType condest =
        prec_->computeCondEst(Ifpack2::Cheap);
    if (verbose) {
      std::cout << "Condition estimate(cheap) for preconditioner: " << condest;
      std::cout << std::endl;
    }
    return true;
  }
  ;
 
  template<class T>
  bool BelosSolver<T>::phaseSolve_(bool verbose) {
    verbose = (Comm_->getRank() == 0);
 
    // Signal that we are done setting up the linear problem
    bool ierr = lp_->setProblem();
 
    // Check the return from setProblem(). If this is true, there was NO
    // error. This probably means we did not specify enough information for
    // the eigenproblem.
    conceptsAssert(ierr == true, concepts::Assertion());
 
    Belos::ReturnType solverRet = solverManager_->solve();
 
    // Check return code of the solver: Unconverged, Failed, or OK
    switch (solverRet) {
 
      // UNCONVERGED
      case Belos::Unconverged:
        std::cout << "##### Iterative solver did not converge!" << std::endl;
        throw std::logic_error("Belos  did not converge!");
 
        // CONVERGED
      case Belos::Converged:
        if (verbose)
          std::cout << "Iterative solver converged!" << std::endl;
        return true;
    }
    conceptsAssert(false, concepts::Assertion());
    return false;
  }
  ;
 
  template<class T>
  void BelosSolver<T>::apply_(const Vector<T>& fncY, Vector<T>& fncX) {
    lp_->setConceptsRHS(fncY, Comm_);
    //solve
    try {
      phaseSolve_();
    } catch (std::logic_error& e) {
      std::cerr << "***** ERROR during solving process: " << e.what()
          << std::endl;
    }
    lp_->writeSolution(fncX, Comm_);
 
  }
 
  template<class T>
  void BelosSolver<T>::apply_() {
    lp_->setConceptsRHS(Comm_);
    try {
      phaseSolve_();
    } catch (std::logic_error& e) {
      std::cerr << "***** ERROR during solving process: " << e.what()
          << std::endl;
    }
    lp_->writeSolution(Comm_);
 
  }
 
} // namespace
#endif /* BELOSSOLVER_DECL_HH_ */