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");
    paramList.remove("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);
    std::string paramXML = MueLu::ML2MueLuParameterTranslator::translate(paramList, "");
    paramList            = *Teuchos::getParametersFromXmlString(paramXML);
  }
  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_ */