MueLu_ShiftedLaplacian_def.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//        MueLu: A package for multigrid based preconditioning
//
// Copyright 2012 NTESS and the MueLu contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef MUELU_SHIFTEDLAPLACIAN_DEF_HPP
#define MUELU_SHIFTEDLAPLACIAN_DEF_HPP
 
#include "MueLu_ShiftedLaplacian_decl.hpp"
 
#if defined(HAVE_MUELU_IFPACK2) and defined(HAVE_MUELU_TPETRA)
 
#include <MueLu_AmalgamationFactory.hpp>
#include <MueLu_CoalesceDropFactory.hpp>
#include <MueLu_CoarseMapFactory.hpp>
#include <MueLu_CoupledRBMFactory.hpp>
#include <MueLu_DirectSolver.hpp>
#include <MueLu_GenericRFactory.hpp>
#include <MueLu_Hierarchy.hpp>
#include <MueLu_Ifpack2Smoother.hpp>
#include <MueLu_PFactory.hpp>
#include <MueLu_PgPFactory.hpp>
#include <MueLu_RAPFactory.hpp>
#include <MueLu_RAPShiftFactory.hpp>
#include <MueLu_SaPFactory.hpp>
#include <MueLu_ShiftedLaplacian.hpp>
#include <MueLu_ShiftedLaplacianOperator.hpp>
#include <MueLu_SmootherFactory.hpp>
#include <MueLu_SmootherPrototype.hpp>
#include <MueLu_TentativePFactory.hpp>
#include <MueLu_TransPFactory.hpp>
#include <MueLu_UncoupledAggregationFactory.hpp>
#include <MueLu_Utilities.hpp>
 
namespace MueLu {
 
// Destructor
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~ShiftedLaplacian() {}
 
// Input
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameters(Teuchos::RCP<Teuchos::ParameterList> paramList) {
  // Parameters
  coarseGridSize_ = paramList->get("MueLu: coarse size", 1000);
  numLevels_      = paramList->get("MueLu: levels", 3);
  int stype       = paramList->get("MueLu: smoother", 8);
  if (stype == 1) {
    Smoother_ = "jacobi";
  } else if (stype == 2) {
    Smoother_ = "gauss-seidel";
  } else if (stype == 3) {
    Smoother_ = "symmetric gauss-seidel";
  } else if (stype == 4) {
    Smoother_ = "chebyshev";
  } else if (stype == 5) {
    Smoother_ = "krylov";
  } else if (stype == 6) {
    Smoother_ = "ilut";
  } else if (stype == 7) {
    Smoother_ = "riluk";
  } else if (stype == 8) {
    Smoother_ = "schwarz";
  } else if (stype == 9) {
    Smoother_ = "superilu";
  } else if (stype == 10) {
    Smoother_ = "superlu";
  } else {
    Smoother_ = "schwarz";
  }
  smoother_sweeps_     = paramList->get("MueLu: sweeps", 5);
  smoother_damping_    = paramList->get("MueLu: relax val", 1.0);
  ncycles_             = paramList->get("MueLu: cycles", 1);
  iters_               = paramList->get("MueLu: iterations", 500);
  solverType_          = paramList->get("MueLu: solver type", 1);
  restart_size_        = paramList->get("MueLu: restart size", 100);
  recycle_size_        = paramList->get("MueLu: recycle size", 25);
  isSymmetric_         = paramList->get("MueLu: symmetric", true);
  ilu_leveloffill_     = paramList->get("MueLu: level-of-fill", 5);
  ilu_abs_thresh_      = paramList->get("MueLu: abs thresh", 0.0);
  ilu_rel_thresh_      = paramList->get("MueLu: rel thresh", 1.0);
  ilu_diagpivotthresh_ = paramList->get("MueLu: piv thresh", 0.1);
  ilu_drop_tol_        = paramList->get("MueLu: drop tol", 0.01);
  ilu_fill_tol_        = paramList->get("MueLu: fill tol", 0.01);
  schwarz_overlap_     = paramList->get("MueLu: overlap", 0);
  schwarz_usereorder_  = paramList->get("MueLu: use reorder", true);
  int combinemode      = paramList->get("MueLu: combine mode", 1);
  if (combinemode == 0) {
    schwarz_combinemode_ = Tpetra::ZERO;
  } else {
    schwarz_combinemode_ = Tpetra::ADD;
  }
  tol_ = paramList->get("MueLu: tolerance", 0.001);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setProblemMatrix(RCP<Matrix>& A) {
  A_ = A;
  if (A_ != Teuchos::null)
    TpetraA_ = toTpetra(A_);
  if (LinearProblem_ != Teuchos::null)
    LinearProblem_->setOperator(TpetraA_);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setProblemMatrix(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraA) {
  TpetraA_ = TpetraA;
  if (LinearProblem_ != Teuchos::null)
    LinearProblem_->setOperator(TpetraA_);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setPreconditioningMatrix(RCP<Matrix>& P) {
  P_                  = P;
  GridTransfersExist_ = false;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setPreconditioningMatrix(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraP) {
  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraP));
  P_                                                                      = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
  GridTransfersExist_                                                     = false;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setstiff(RCP<Matrix>& K) {
  K_ = K;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setstiff(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraK) {
  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraK));
  K_                                                                      = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setmass(RCP<Matrix>& M) {
  M_ = M;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setmass(RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> >& TpetraM) {
  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Atmp = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraM));
  M_                                                                      = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Atmp));
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setcoords(RCP<MultiVector>& Coords) {
  Coords_ = Coords;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setNullSpace(RCP<MultiVector> NullSpace) {
  NullSpace_ = NullSpace;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setLevelShifts(std::vector<Scalar> levelshifts) {
  levelshifts_ = levelshifts;
  numLevels_   = levelshifts_.size();
}
 
// initialize
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::initialize() {
  TentPfact_     = rcp(new TentativePFactory);
  Pfact_         = rcp(new SaPFactory);
  PgPfact_       = rcp(new PgPFactory);
  TransPfact_    = rcp(new TransPFactory);
  Rfact_         = rcp(new GenericRFactory);
  Acfact_        = rcp(new RAPFactory);
  Acshift_       = rcp(new RAPShiftFactory);
  Amalgfact_     = rcp(new AmalgamationFactory);
  Dropfact_      = rcp(new CoalesceDropFactory);
  UCaggfact_     = rcp(new UncoupledAggregationFactory);
  CoarseMapfact_ = rcp(new CoarseMapFactory);
  Manager_       = rcp(new FactoryManager);
  Manager_->SetFactory("UnAmalgamationInfo", Amalgfact_);
  Teuchos::ParameterList params;
  params.set("lightweight wrap", true);
  params.set("aggregation: drop scheme", "classical");
  Dropfact_->SetParameterList(params);
  Manager_->SetFactory("Graph", Dropfact_);
  Manager_->SetFactory("Aggregates", UCaggfact_);
  Manager_->SetFactory("CoarseMap", CoarseMapfact_);
  Manager_->SetFactory("Ptent", TentPfact_);
  if (isSymmetric_ == true) {
    Manager_->SetFactory("P", Pfact_);
    Manager_->SetFactory("R", TransPfact_);
  } else {
    Manager_->SetFactory("P", PgPfact_);
    Manager_->SetFactory("R", Rfact_);
    solverType_ = 10;
  }
 
  // choose smoother
  if (Smoother_ == "jacobi") {
    precType_ = "RELAXATION";
    precList_.set("relaxation: type", "Jacobi");
    precList_.set("relaxation: sweeps", smoother_sweeps_);
    precList_.set("relaxation: damping factor", smoother_damping_);
  } else if (Smoother_ == "gauss-seidel") {
    precType_ = "RELAXATION";
    precList_.set("relaxation: type", "Gauss-Seidel");
    precList_.set("relaxation: sweeps", smoother_sweeps_);
    precList_.set("relaxation: damping factor", smoother_damping_);
  } else if (Smoother_ == "symmetric gauss-seidel") {
    precType_ = "RELAXATION";
    precList_.set("relaxation: type", "Symmetric Gauss-Seidel");
    precList_.set("relaxation: sweeps", smoother_sweeps_);
    precList_.set("relaxation: damping factor", smoother_damping_);
  } else if (Smoother_ == "chebyshev") {
    precType_ = "CHEBYSHEV";
  } else if (Smoother_ == "krylov") {
    precType_ = "KRYLOV";
    precList_.set("krylov: iteration type", krylov_type_);
    precList_.set("krylov: number of iterations", krylov_iterations_);
    precList_.set("krylov: residual tolerance", 1.0e-8);
    precList_.set("krylov: block size", 1);
    precList_.set("krylov: preconditioner type", krylov_preconditioner_);
    precList_.set("relaxation: sweeps", 1);
    solverType_ = 10;
  } else if (Smoother_ == "ilut") {
    precType_ = "ILUT";
    precList_.set("fact: ilut level-of-fill", ilu_leveloffill_);
    precList_.set("fact: absolute threshold", ilu_abs_thresh_);
    precList_.set("fact: relative threshold", ilu_rel_thresh_);
    precList_.set("fact: drop tolerance", ilu_drop_tol_);
    precList_.set("fact: relax value", ilu_relax_val_);
  } else if (Smoother_ == "riluk") {
    precType_ = "RILUK";
    precList_.set("fact: iluk level-of-fill", ilu_leveloffill_);
    precList_.set("fact: absolute threshold", ilu_abs_thresh_);
    precList_.set("fact: relative threshold", ilu_rel_thresh_);
    precList_.set("fact: drop tolerance", ilu_drop_tol_);
    precList_.set("fact: relax value", ilu_relax_val_);
  } else if (Smoother_ == "schwarz") {
    precType_ = "SCHWARZ";
    precList_.set("schwarz: overlap level", schwarz_overlap_);
    precList_.set("schwarz: combine mode", schwarz_combinemode_);
    precList_.set("schwarz: use reordering", schwarz_usereorder_);
    //    precList_.set("schwarz: filter singletons", true); // Disabled due to issues w/ Ifpack2/Zoltan2 w.r.t. Issue #560 - CMS 8/26/16
    precList_.set("order_method", schwarz_ordermethod_);
    precList_.sublist("schwarz: reordering list").set("order_method", schwarz_ordermethod_);
    precList_.sublist("schwarz: subdomain solver parameters").set("fact: ilut level-of-fill", ilu_leveloffill_);
    precList_.sublist("schwarz: subdomain solver parameters").set("fact: absolute threshold", ilu_abs_thresh_);
    precList_.sublist("schwarz: subdomain solver parameters").set("fact: relative threshold", ilu_rel_thresh_);
    precList_.sublist("schwarz: subdomain solver parameters").set("fact: drop tolerance", ilu_drop_tol_);
    precList_.sublist("schwarz: subdomain solver parameters").set("fact: relax value", ilu_relax_val_);
  } else if (Smoother_ == "superilu") {
    precType_ = "superlu";
    precList_.set("RowPerm", ilu_rowperm_);
    precList_.set("ColPerm", ilu_colperm_);
    precList_.set("DiagPivotThresh", ilu_diagpivotthresh_);
    precList_.set("ILU_DropRule", ilu_drop_rule_);
    precList_.set("ILU_DropTol", ilu_drop_tol_);
    precList_.set("ILU_FillFactor", ilu_leveloffill_);
    precList_.set("ILU_Norm", ilu_normtype_);
    precList_.set("ILU_MILU", ilu_milutype_);
    precList_.set("ILU_FillTol", ilu_fill_tol_);
    precList_.set("ILU_Flag", true);
  } else if (Smoother_ == "superlu") {
    precType_ = "superlu";
    precList_.set("ColPerm", ilu_colperm_);
    precList_.set("DiagPivotThresh", ilu_diagpivotthresh_);
  }
  // construct smoother
  smooProto_ = rcp(new Ifpack2Smoother(precType_, precList_));
  smooFact_  = rcp(new SmootherFactory(smooProto_));
#if defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLU)
  coarsestSmooProto_ = rcp(new DirectSolver("Superlu", coarsestSmooList_));
#elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_KLU2)
  coarsestSmooProto_ = rcp(new DirectSolver("Klu", coarsestSmooList_));
#elif defined(HAVE_MUELU_AMESOS2) and defined(HAVE_AMESOS2_SUPERLUDIST)
  coarsestSmooProto_ = rcp(new DirectSolver("Superludist", coarsestSmooList_));
#else
  coarsestSmooProto_ = rcp(new Ifpack2Smoother(precType_, precList_));
#endif
  coarsestSmooFact_ = rcp(new SmootherFactory(coarsestSmooProto_, Teuchos::null));
 
  // For setupSlowRAP and setupFastRAP, the prolongation/restriction matrices
  // are constructed with the stiffness matrix. These matrices are kept for future
  // setup calls; this is achieved by calling Hierarchy->Keep(). It is particularly
  // useful for multiple frequency problems - when the frequency/preconditioner
  // changes, you only compute coarse grids (RAPs) and setup level smoothers when
  // you call Hierarchy->Setup().
  if (K_ != Teuchos::null) {
    Manager_->SetFactory("Smoother", Teuchos::null);
    Manager_->SetFactory("CoarseSolver", Teuchos::null);
    Hierarchy_ = rcp(new Hierarchy(K_));
    if (NullSpace_ != Teuchos::null)
      Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
    if (isSymmetric_ == true) {
      Hierarchy_->Keep("P", Pfact_.get());
      Hierarchy_->Keep("R", TransPfact_.get());
      Hierarchy_->SetImplicitTranspose(true);
    } else {
      Hierarchy_->Keep("P", PgPfact_.get());
      Hierarchy_->Keep("R", Rfact_.get());
    }
    Hierarchy_->Keep("Ptent", TentPfact_.get());
    Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
    Hierarchy_->Setup(*Manager_, 0, numLevels_);
    GridTransfersExist_ = true;
  }
  // Use preconditioning matrix to setup prolongation/restriction operators
  else {
    Manager_->SetFactory("Smoother", smooFact_);
    Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
    Hierarchy_ = rcp(new Hierarchy(P_));
    if (NullSpace_ != Teuchos::null)
      Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
    if (isSymmetric_ == true)
      Hierarchy_->SetImplicitTranspose(true);
    Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
    Hierarchy_->Setup(*Manager_, 0, numLevels_);
    GridTransfersExist_ = true;
  }
 
  // Belos Linear Problem and Solver Manager
  BelosList_ = rcp(new Teuchos::ParameterList("GMRES"));
  BelosList_->set("Maximum Iterations", iters_);
  BelosList_->set("Convergence Tolerance", tol_);
  BelosList_->set("Verbosity", Belos::Errors + Belos::Warnings + Belos::StatusTestDetails);
  BelosList_->set("Output Frequency", 1);
  BelosList_->set("Output Style", Belos::Brief);
  BelosList_->set("Num Blocks", restart_size_);
  BelosList_->set("Num Recycled Blocks", recycle_size_);
}
 
// setup coarse grids for new frequency
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupFastRAP() {
  int numLevels = Hierarchy_->GetNumLevels();
  Manager_->SetFactory("Smoother", smooFact_);
  Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
  Hierarchy_->GetLevel(0)->Set("A", P_);
  Hierarchy_->Setup(*Manager_, 0, numLevels);
  setupSolver();
}
 
// setup coarse grids for new frequency
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSlowRAP() {
  int numLevels = Hierarchy_->GetNumLevels();
  Acshift_->SetShifts(levelshifts_);
  Manager_->SetFactory("Smoother", smooFact_);
  Manager_->SetFactory("CoarseSolver", coarsestSmooFact_);
  Manager_->SetFactory("A", Acshift_);
  Manager_->SetFactory("K", Acshift_);
  Manager_->SetFactory("M", Acshift_);
  Hierarchy_->GetLevel(0)->Set("A", P_);
  Hierarchy_->GetLevel(0)->Set("K", K_);
  Hierarchy_->GetLevel(0)->Set("M", M_);
  Hierarchy_->Setup(*Manager_, 0, numLevels);
  setupSolver();
}
 
// setup coarse grids for new frequency
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupNormalRAP() {
  // Only setup hierarchy again if preconditioning matrix has changed
  if (GridTransfersExist_ == false) {
    Hierarchy_ = rcp(new Hierarchy(P_));
    if (NullSpace_ != Teuchos::null)
      Hierarchy_->GetLevel(0)->Set("Nullspace", NullSpace_);
    if (isSymmetric_ == true)
      Hierarchy_->SetImplicitTranspose(true);
    Hierarchy_->SetMaxCoarseSize(coarseGridSize_);
    Hierarchy_->Setup(*Manager_, 0, numLevels_);
    GridTransfersExist_ = true;
  }
  setupSolver();
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setupSolver() {
  // Define Preconditioner and Operator
  MueLuOp_ = rcp(new MueLu::ShiftedLaplacianOperator<SC, LO, GO, NO>(Hierarchy_, A_, ncycles_, subiters_, option_, tol_));
  // Belos Linear Problem
  if (LinearProblem_ == Teuchos::null)
    LinearProblem_ = rcp(new LinearProblem);
  LinearProblem_->setOperator(TpetraA_);
  LinearProblem_->setRightPrec(MueLuOp_);
  if (SolverManager_ == Teuchos::null) {
    std::string solverName;
    SolverFactory_ = rcp(new SolverFactory());
    if (solverType_ == 1) {
      solverName = "Block GMRES";
    } else if (solverType_ == 2) {
      solverName = "Recycling GMRES";
    } else {
      solverName = "Flexible GMRES";
    }
    SolverManager_ = SolverFactory_->create(solverName, BelosList_);
    SolverManager_->setProblem(LinearProblem_);
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::resetLinearProblem() {
  LinearProblem_->setOperator(TpetraA_);
}
 
// Solve phase
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Belos::ReturnType ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::solve(const RCP<TMV> B, RCP<TMV>& X) {
  // Set left and right hand sides for Belos
  LinearProblem_->setProblem(X, B);
  // iterative solve
  return SolverManager_->solve();
}
 
// Solve phase
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multigrid_apply(const RCP<MultiVector> B,
                                                                                  RCP<MultiVector>& X) {
  // Set left and right hand sides for Belos
  Hierarchy_->Iterate(*B, *X, 1, true, 0);
}
 
// Solve phase
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multigrid_apply(const RCP<Tpetra::MultiVector<SC, LO, GO, NO> > B,
                                                                                  RCP<Tpetra::MultiVector<SC, LO, GO, NO> >& X) {
  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > XpetraX = Teuchos::rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(X));
  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > XpetraB = Teuchos::rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(B));
  // Set left and right hand sides for Belos
  Hierarchy_->Iterate(*XpetraB, *XpetraX, 1, true, 0);
}
 
// Get most recent iteration count
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
int ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetIterations() {
  int numiters = SolverManager_->getNumIters();
  return numiters;
}
 
// Get most recent solver tolerance achieved
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
typename Teuchos::ScalarTraits<Scalar>::magnitudeType
ShiftedLaplacian<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetResidual() {
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
  MT residual = SolverManager_->achievedTol();
  return residual;
}
 
}  // namespace MueLu
 
#define MUELU_SHIFTEDLAPLACIAN_SHORT
 
#endif  // if defined(HAVE_MUELU_IFPACK2)
#endif  // MUELU_SHIFTEDLAPLACIAN_DEF_HPP