MueLu_CreateXpetraPreconditioner.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 PACKAGES_MUELU_ADAPTERS_XPETRA_MUELU_CREATEXPETRAPRECONDITIONER_HPP_
#define PACKAGES_MUELU_ADAPTERS_XPETRA_MUELU_CREATEXPETRAPRECONDITIONER_HPP_
 
 
#include <Teuchos_XMLParameterListHelpers.hpp>
#include <Xpetra_CrsMatrix.hpp>
#include <Xpetra_MultiVector.hpp>
 
#include <MueLu.hpp>
 
#include <MueLu_Exceptions.hpp>
#include <MueLu_Hierarchy.hpp>
#include <MueLu_MasterList.hpp>
#include <MueLu_ParameterListInterpreter.hpp>
#include <MueLu_Utilities.hpp>
#include <MueLu_HierarchyUtils.hpp>
#include <MueLu_ML2MueLuParameterTranslator.hpp>
#include <stdlib.h>
 
#include <MueLu_KokkosTuningInterface.hpp>
 
namespace MueLu {
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
CreateXpetraPreconditioner(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> op,
                           const Teuchos::ParameterList& inParamList) {
  using SC = Scalar;
  using LO = LocalOrdinal;
  using GO = GlobalOrdinal;
  using NO = Node;
 
  using HierarchyManager         = MueLu::HierarchyManager<SC, LO, GO, NO>;
  using HierarchyUtils           = MueLu::HierarchyUtils<SC, LO, GO, NO>;
  using Hierarchy                = MueLu::Hierarchy<SC, LO, GO, NO>;
  using ParameterListInterpreter = ParameterListInterpreter<SC, LO, GO, NO>;
 
  bool hasParamList = inParamList.numParams();
 
  RCP<HierarchyManager> mueLuFactory;
 
  // Rip off non-serializable data before validation
  Teuchos::ParameterList nonSerialList, paramList;
  MueLu::ExtractNonSerializableData(inParamList, paramList, nonSerialList);
 
  std::string label;
  if (hasParamList && paramList.isParameter("hierarchy label")) {
    label = paramList.get<std::string>("hierarchy label");
  } else
    label = op->getObjectLabel();
 
  RCP<Teuchos::Time> tm;
  std::string timerName;
 
  // Do Kokkos tuning, if requested
  if (paramList.isSublist("kokkos tuning: muelu parameter mapping") &&
      paramList.sublist("kokkos tuning: muelu parameter mapping").isParameter("kokkos context id")) {
    // Time tuning separately
    if (label != "")
      timerName = "MueLu tuning time (" + label + ")";
    else
      timerName = "MueLu tuning time";
    tm = Teuchos::TimeMonitor::getNewTimer(timerName);
    tm->start();
 
    MueLu::KokkosTuningInterface KokkosTuner(op->getMap()->getComm());
    KokkosTuner.SetParameterList(paramList);
    KokkosTuner.SetMueLuParameters(paramList);
    tm->stop();
    tm->incrementNumCalls();
    tm = Teuchos::null;
  }
 
  // Setup Timer
  if (label != "")
    timerName = "MueLu setup time (" + label + ")";
  else
    timerName = "MueLu setup time";
 
  tm = Teuchos::TimeMonitor::getNewTimer(timerName);
  tm->start();
 
  std::string syntaxStr = "parameterlist: syntax";
  if (hasParamList && paramList.isParameter(syntaxStr) && paramList.get<std::string>(syntaxStr) == "ml") {
    paramList.remove(syntaxStr);
    paramList = *MueLu::ML2MueLuParameterTranslator::translate(paramList, "");
  }
  //  Need to check if Muelu option is inconsistent with user data provided
  bool Minv_Supplied = false, M_Supplied = false, MinvA_Supplied = false;
  if (inParamList.isSublist("user data")) {
    const Teuchos::ParameterList& userList = inParamList.sublist("user data");
    if (userList.isParameter("M")) M_Supplied = true;
    if (userList.isParameter("Minv")) Minv_Supplied = true;
    if (userList.isParameter("MinvA")) MinvA_Supplied = true;
  }
  if (inParamList.isParameter("aggregation: strength-of-connection: matrix") && (inParamList.get<std::string>("aggregation: strength-of-connection: matrix") == "MinvA")) {
    if (inParamList.isSublist("project auxiliary matrices")) {
      auto projectList = inParamList.sublist("project auxiliary matrices");
      TEUCHOS_TEST_FOR_EXCEPTION(projectList.isParameter("M") && !M_Supplied, Exceptions::Incompatible, "MueLu_CreateXpetraPreconditioner: Must supply M as it is listed in the project auxiliary matrices sublist");
      TEUCHOS_TEST_FOR_EXCEPTION(projectList.isParameter("Minv") && (projectList.get("Minv", "") == "NoFactory") && !Minv_Supplied, Exceptions::Incompatible,
                                 "MueLu_CreateXpetraPreconditioner: Must supply Minv  as NoFactory is listed as supplier of Minv  in the project auxiliary matrices sublist");
      TEUCHOS_TEST_FOR_EXCEPTION(projectList.isParameter("MinvA") && (projectList.get("MinvA", "") == "NoFactory") && !MinvA_Supplied, Exceptions::Incompatible,
                                 "MueLu_CreateXpetraPreconditioner: Must supply MinvA as NoFactory is listed as supplier of MinvA in the project auxiliary matrices sublist");
    } else {  // default behavior if sublist("project auxiliary matrices") not user-supplied requires "M" to be user-supplied.
      TEUCHOS_TEST_FOR_EXCEPTION(!M_Supplied && !Minv_Supplied, Exceptions::Incompatible, "MueLu_CreateXpetraPreconditioner: Must supply M or Minv when 'aggregation: strength-of-connection: matrix'= MinvA and sublist('project auxiliary matrices') not supplied.");
    }
  }
 
  mueLuFactory = rcp(new ParameterListInterpreter(paramList, op->getDomainMap()->getComm()));
 
  // Create Hierarchy
  RCP<Hierarchy> H = mueLuFactory->CreateHierarchy(label);
  H->setlib(op->getDomainMap()->lib());
 
  // Set fine level operator
  H->GetLevel(0)->Set("A", op);
  H->SetProcRankVerbose(op->getDomainMap()->getComm()->getRank());
 
  // Stick the non-serializible data on the hierarchy.
  HierarchyUtils::AddNonSerializableDataToHierarchy(*mueLuFactory, *H, nonSerialList);
 
  mueLuFactory->SetupHierarchy(*H);
 
  tm->stop();
  tm->incrementNumCalls();
 
  if (H->GetVerbLevel() & Statistics0) {
    const bool alwaysWriteLocal = true;
    const bool writeGlobalStats = true;
    const bool writeZeroTimers  = false;
    const bool ignoreZeroTimers = true;
    const std::string filter    = timerName;
    Teuchos::TimeMonitor::summarize(op->getRowMap()->getComm().ptr(), H->GetOStream(Statistics0), alwaysWriteLocal, writeGlobalStats,
                                    writeZeroTimers, Teuchos::Union, filter, ignoreZeroTimers);
  }
 
  tm->reset();
 
  return H;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
CreateXpetraPreconditioner(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> op,
                           const std::string& xmlFileName) {
  Teuchos::ParameterList paramList;
  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *op->getDomainMap()->getComm());
  return CreateXpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(op, paramList);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
CreateXpetraPreconditioner(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> op) {
  Teuchos::ParameterList paramList;
  return CreateXpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(op, paramList);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void ReuseXpetraPreconditioner(const Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
                               Teuchos::RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& H) {
  std::string label = H->GetLevel(0)->getObjectLabel();
 
  std::string timerName;
  if (label != "")
    timerName = "MueLu setup time (" + label + ")";
  else
    timerName = "MueLu setup time";
  RCP<Teuchos::Time> tm = Teuchos::TimeMonitor::getNewTimer(timerName);
  tm->start();
 
  typedef Scalar SC;
  typedef LocalOrdinal LO;
  typedef GlobalOrdinal GO;
  typedef Node NO;
 
  typedef Xpetra::Matrix<SC, LO, GO, NO> Matrix;
  typedef Xpetra::Operator<SC, LO, GO, NO> Operator;
 
  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), Exceptions::RuntimeError,
                             "MueLu::ReuseXpetraPreconditioner: Hierarchy has no levels in it");
  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), Exceptions::RuntimeError,
                             "MueLu::ReuseXpetraPreconditioner: Hierarchy has no fine level operator");
  RCP<Level> level0 = H->GetLevel(0);
 
  RCP<Operator> O0 = level0->Get<RCP<Operator>>("A");
  RCP<Matrix> A0   = Teuchos::rcp_dynamic_cast<Matrix>(O0);
 
  if (!A0.is_null()) {
    // If a user provided a "number of equations" argument in a parameter list
    // during the initial setup, we must honor that settings and reuse it for
    // all consequent setups.
    A->SetFixedBlockSize(A0->GetFixedBlockSize());
  }
  level0->Set("A", A);
 
  H->SetupRe();
 
  tm->stop();
  tm->incrementNumCalls();
 
  if (H->GetVerbLevel() & Statistics0) {
    const bool alwaysWriteLocal = true;
    const bool writeGlobalStats = true;
    const bool writeZeroTimers  = false;
    const bool ignoreZeroTimers = true;
    const std::string filter    = timerName;
    Teuchos::TimeMonitor::summarize(A->getRowMap()->getComm().ptr(), H->GetOStream(Statistics0), alwaysWriteLocal, writeGlobalStats,
                                    writeZeroTimers, Teuchos::Union, filter, ignoreZeroTimers);
  }
 
  tm->reset();
}
 
}  // namespace MueLu
 
#define XPETRA_CREATEXPETRAPRECONDITIONER_SHORT
 
#endif /* PACKAGES_MUELU_ADAPTERS_XPETRA_MUELU_CREATEXPETRAPRECONDITIONER_HPP_ */