MueLu_EminPFactory_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_EMINPFACTORY_DEF_HPP
#define MUELU_EMINPFACTORY_DEF_HPP
 
#include <Xpetra_Matrix.hpp>
#include <Xpetra_StridedMapFactory.hpp>
 
#include "MueLu_EminPFactory_decl.hpp"
#include "MueLu_ProductOperator.hpp"
#include "MueLu_FlatOperator.hpp"
 
#include "MueLu_Constraint.hpp"
#include "MueLu_MasterList.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_PerfUtils.hpp"
#include "MueLu_Behavior.hpp"
 
#ifdef HAVE_MUELU_BELOS
#include <BelosLinearProblem.hpp>
#include <BelosSolverFactory.hpp>
#include <BelosXpetraAdapter.hpp>
#endif
 
namespace MueLu {
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
typename Teuchos::ScalarTraits<Scalar>::magnitudeType
EminPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
    ComputeProlongatorEnergyNorm(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
                                 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& P,
                                 Teuchos::FancyOStream& out) {
  using MagnitudeType = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
  RCP<Matrix> AP;
  AP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *P, false, AP, out, true, true);
  RCP<Matrix> PtAP;
  PtAP             = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *AP, false, PtAP, out, true, true);
  RCP<Vector> diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(PtAP->getRowMap());
  PtAP->getLocalDiagCopy(*diag);
  MagnitudeType norm = diag->norm1();
  norm               = Teuchos::ScalarTraits<MagnitudeType>::squareroot(norm);
  return norm;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const ParameterList> EminPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
  RCP<ParameterList> validParamList = rcp(new ParameterList());
 
#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
  SET_VALID_ENTRY("emin: num iterations");
  SET_VALID_ENTRY("emin: num reuse iterations");
  SET_VALID_ENTRY("emin: iterative method");
  {
    validParamList->getEntry("emin: iterative method").setValidator(rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("cg", "sd", "gmres"))));
  }
#undef SET_VALID_ENTRY
 
  validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory for the matrix A used during internal iterations");
  validParamList->set<RCP<const FactoryBase>>("P", Teuchos::null, "Generating factory for the initial guess");
  validParamList->set<RCP<const FactoryBase>>("Constraint", Teuchos::null, "Generating factory for constraints");
 
  validParamList->set<RCP<Matrix>>("P0", Teuchos::null, "Initial guess at P");
  validParamList->set<bool>("Keep P0", false, "Keep an initial P0 (for reuse)");
 
  validParamList->set<RCP<Constraint>>("Constraint0", Teuchos::null, "Initial Constraint");
  validParamList->set<bool>("Keep Constraint0", false, "Keep an initial Constraint (for reuse)");
 
  return validParamList;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void EminPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& fineLevel, Level& coarseLevel) const {
  Input(fineLevel, "A");
 
  static bool isAvailableP0          = false;
  static bool isAvailableConstraint0 = false;
 
  // Here is a tricky little piece of code
  // We don't want to request (aka call Input) when we reuse and P0 is available
  // However, we cannot run something simple like this:
  //   if (!coarseLevel.IsAvailable("P0", this))
  //     Input(coarseLevel, "P");
  // The reason is that it works fine during the request stage, but fails
  // in the release stage as we _construct_ P0 during Build process. Therefore,
  // we need to understand whether we are in Request or Release mode
  // NOTE: This is a very unique situation, please try not to propagate the
  // mode check any further
 
  // #define externalSuppliedP0
  if (coarseLevel.GetRequestMode() == Level::REQUEST) {
#ifdef externalSuppliedP0
    isAvailableP0 = coarseLevel.IsAvailable("P0");
#else
    isAvailableP0 = coarseLevel.IsAvailable("P0", this);
#endif
    isAvailableConstraint0 = coarseLevel.IsAvailable("Constraint0", this);
  }
 
  if (isAvailableP0 == false)
    Input(coarseLevel, "P");
 
  if (isAvailableConstraint0 == false)
    Input(coarseLevel, "Constraint");
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void EminPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& fineLevel, Level& coarseLevel) const {
  BuildP(fineLevel, coarseLevel);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void EminPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level& fineLevel, Level& coarseLevel) const {
  FactoryMonitor m(*this, "Prolongator minimization", coarseLevel);
 
  const ParameterList& pL = GetParameterList();
 
  // Get the matrix
  RCP<Matrix> A = Get<RCP<Matrix>>(fineLevel, "A");
 
  if (restrictionMode_) {
    SubFactoryMonitor m2(*this, "Transpose A", coarseLevel);
 
    A = Utilities::Transpose(*A, true);
  }
 
  // Get/make constraint operator
  RCP<Constraint> X;
  if (coarseLevel.IsAvailable("Constraint0", this)) {
    // Reuse data
    X = coarseLevel.Get<RCP<Constraint>>("Constraint0", this);
    GetOStream(Runtime0) << "Reusing Constraint0" << std::endl;
 
  } else {
    // Construct data
    X = Get<RCP<Constraint>>(coarseLevel, "Constraint");
  }
 
  // Get/make initial guess
  // NOTE: the main assumption here that P0 satisfies both constraints:
  //   - nonzero pattern
  //   - nullspace preservation
  RCP<Matrix> P0;
  int numIts;
#ifdef externalSuppliedP0
  if (coarseLevel.IsAvailable("P0")) {
#else
  if (coarseLevel.IsAvailable("P0", this)) {
#endif
    // Reuse data
#ifdef externalSuppliedP0
    P0 = coarseLevel.Get<RCP<Matrix>>("P0");
#else
    P0 = coarseLevel.Get<RCP<Matrix>>("P0", this);
#endif
    numIts = pL.get<int>("emin: num reuse iterations");
    GetOStream(Runtime0) << "Reusing P0" << std::endl;
 
  } else {
    GetOStream(Runtime0) << "Getting P from coarseLevel" << std::endl;
    // Construct data
    P0     = Get<RCP<Matrix>>(coarseLevel, "P");
    numIts = pL.get<int>("emin: num iterations");
  }
 
  if (Behavior::debug() && IsPrint(Statistics1)) {
    SubFactoryMonitor m2(*this, "Statistics", coarseLevel);
    GetOStream(Statistics1) << "Energy norm of P0: " << ComputeProlongatorEnergyNorm(A, P0, GetOStream(Statistics0)) << std::endl;
    GetOStream(Statistics1) << "Constraint residual norm of P0: " << X->ResidualNorm(P0) << std::endl;
  }
 
  std::string solverType = pL.get<std::string>("emin: iterative method");
  if (solverType == "cg")
    solverType = "Pseudo Block CG";
  else if (solverType == "gmres")
    solverType = "Pseudo Block GMRES";
 
  RCP<Matrix> P;
#ifdef HAVE_MUELU_BELOS
  if (numIts > 0) {
    // Construct diagonal preconditioner
    RCP<const Vector> invDiagonal = Utilities::GetMatrixDiagonalInverse(*A);
    auto invD                     = MatrixFactory::Build(invDiagonal);
    auto flatInvD                 = rcp(new FlatOperator(invD, X));
 
    // Construct projection of A
    // In principle we should be using the operator proj * A * proj.
    // But since the initial guess already satisfies the constraints and all subsequent iterates do too,
    // the application of proj on the right side is a no-op.
    auto flatA                     = rcp(new FlatOperator(A, X));
    std::vector<RCP<Operator>> ops = {X, flatA};  // Instead of ops = {X, flatA, X};
    RCP<Operator> projectedFlatA   = rcp(new ProductOperator(ops));
 
    // Vectors for rhs and solution
    auto B    = MultiVectorFactory::Build(projectedFlatA->getRangeMap(), 1);
    auto vecP = MultiVectorFactory::Build(projectedFlatA->getDomainMap(), 1);
    // initialize solution vector to P0
    X->AssignMatrixEntriesToVector(*P0, *vecP);
 
    auto belosProjectedFlatA = rcp(new Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(projectedFlatA));
    auto belosFlatInvDiag    = rcp(new Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(flatInvD));
 
    {
      using MV = Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
      using OP = Belos::OperatorT<MV>;
 
      auto problem = rcp(new Belos::LinearProblem<Scalar, MV, OP>());
      problem->setOperator(belosProjectedFlatA);
      problem->setRightPrec(belosFlatInvDiag);
      problem->setLHS(vecP);
      problem->setRHS(B);
      problem->setLabel("Emin");
      TEUCHOS_ASSERT(problem->setProblem());
 
      auto belosList = rcp(new Teuchos::ParameterList());
      belosList->set("Timer Label", "Emin");
      belosList->set("Maximum Iterations", numIts);
      belosList->set("Implicit Residual Scaling", "None");
      belosList->set("Verbosity", Belos::Errors + Belos::Warnings + Belos::StatusTestDetails);
      belosList->set("Output Frequency", 1);
      belosList->set("Output Style", Belos::Brief);
      auto out = GetMueLuOStream();
      belosList->set("Output Stream", out->getOStream());
 
      Belos::SolverFactory<Scalar, MV, OP> solverFactory;
      auto solver = solverFactory.create(solverType, belosList);
      solver->setProblem(problem);
 
      if (numIts > 0) solver->solve();
    }
    // Convert from vector to matrix
    P = X->GetMatrixWithEntriesFromVector(*vecP);
  } else {
    P = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(P0);
  }
  if (P0->IsView("stridedMaps"))
    P->CreateView("stridedMaps", P0);
#else
  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Energy minimization multigrid requires Belos to be enabled.");
#endif
 
  if (Behavior::debug() && IsPrint(Statistics1)) {
    SubFactoryMonitor m2(*this, "Statistics", coarseLevel);
    // Compute constraint residual norm if we are not in a debug build
    GetOStream(Statistics1) << "Energy norm of P: " << ComputeProlongatorEnergyNorm(A, P, GetOStream(Statistics0)) << std::endl;
    GetOStream(Statistics1) << "Norm of constraint residual: " << X->ResidualNorm(P) << std::endl;
  }
 
  // NOTE: EXPERIMENTAL and FRAGILE
  if (!P->IsView("stridedMaps")) {
    if (A->IsView("stridedMaps") == true) {
      GetOStream(Runtime1) << "Using A to fillComplete P" << std::endl;
 
      // FIXME: X->GetPattern() actually returns a CrsGraph.
      // CrsGraph has no knowledge of Xpetra's sup/Matrix views. As such,
      // it has no idea about strided maps. We create one, which is
      // most likely incorrect for many use cases.
      std::vector<size_t> stridingInfo(1, 1);
      RCP<const StridedMap> dMap = StridedMapFactory::Build(X->GetPattern()->getDomainMap(), stridingInfo);
 
      P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), dMap);
    }
  }
 
  // Level Set
  if (!restrictionMode_) {
    // The factory is in prolongation mode
    Set(coarseLevel, "P", P);
 
    if (pL.get<bool>("Keep P0")) {
      // NOTE: we must do Keep _before_ set as the Needs class only sets if
      //  a) data has been requested (which is not the case here), or
      //  b) data has some keep flag
      coarseLevel.Keep("P0", this);
      Set(coarseLevel, "P0", P);
    }
    if (pL.get<bool>("Keep Constraint0")) {
      // NOTE: we must do Keep _before_ set as the Needs class only sets if
      //  a) data has been requested (which is not the case here), or
      //  b) data has some keep flag
      coarseLevel.Keep("Constraint0", this);
      Set(coarseLevel, "Constraint0", X);
    }
 
    if (IsPrint(Statistics2)) {
      RCP<ParameterList> params = rcp(new ParameterList());
      params->set("printLoadBalancingInfo", true);
      params->set("printCommInfo", true);
      GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*P, "P", params);
    }
 
  } else {
    // The factory is in restriction mode
    RCP<Matrix> R;
    {
      SubFactoryMonitor m2(*this, "Transpose P", coarseLevel);
 
      R = Utilities::Transpose(*P, true);
    }
 
    Set(coarseLevel, "R", R);
 
    if (IsPrint(Statistics2)) {
      RCP<ParameterList> params = rcp(new ParameterList());
      params->set("printLoadBalancingInfo", true);
      params->set("printCommInfo", true);
      GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*R, "R", params);
    }
  }
}
 
}  // namespace MueLu
 
#endif  // MUELU_EMINPFACTORY_DEF_HPP