MueLu_RebalanceAcFactory_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_REBALANCEACFACTORY_DEF_HPP
#define MUELU_REBALANCEACFACTORY_DEF_HPP
 
#include <Xpetra_Matrix.hpp>
#include <Xpetra_CrsMatrix.hpp>
#include <Xpetra_CrsMatrixWrap.hpp>
#include <Xpetra_MatrixFactory.hpp>
 
#include "MueLu_RebalanceAcFactory_decl.hpp"
 
#include "MueLu_MasterList.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_PerfUtils.hpp"
#include "MueLu_RAPFactory.hpp"
 
namespace MueLu {
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RebalanceAcFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::RebalanceAcFactory() = default;
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RebalanceAcFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~RebalanceAcFactory() = default;
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const ParameterList> RebalanceAcFactory<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("repartition: use subcommunicators");
  SET_VALID_ENTRY("repartition: use subcommunicators in place");
#undef SET_VALID_ENTRY
 
  validParamList->set<std::string>("Matrix name", "A", "Name of the matrix that will be transferred on the coarse grid (level key)");
 
  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A for rebalancing");
  validParamList->set<RCP<const FactoryBase> >("Importer", Teuchos::null, "Generating factory of the importer");
  validParamList->set<RCP<const FactoryBase> >("InPlaceMap", Teuchos::null, "Generating factory of the InPlaceMap");
 
  return validParamList;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void RebalanceAcFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level & /* fineLevel */, Level &coarseLevel) const {
  const Teuchos::ParameterList &pL = GetParameterList();
  std::string matrixName           = pL.get<std::string>("Matrix name");
  coarseLevel.DeclareInput(matrixName, GetFactory("A").get(), this);
  if (pL.isParameter("repartition: use subcommunicators in place") && pL.get<bool>("repartition: use subcommunicators in place") == true) {
    Input(coarseLevel, "InPlaceMap");
  } else
    Input(coarseLevel, "Importer");
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void RebalanceAcFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level & /* fineLevel */, Level &coarseLevel) const {
  const Teuchos::ParameterList &pL = GetParameterList();
  auto matrixName                  = pL.get<std::string>("Matrix name");
  std::string coarseMatrixName;
  if (matrixName.size() == 1)
    coarseMatrixName = matrixName + "c";
  else
    coarseMatrixName = matrixName + "_coarse";
 
  FactoryMonitor m(*this, "Computing " + coarseMatrixName, coarseLevel);
 
  RCP<Matrix> originalAc = coarseLevel.Get<RCP<Matrix> >(matrixName, GetFactory("A").get());
 
  // This is a short-circuit for if we want to leave A where it is, but restrict its communicator
  // to something corresponding to a given map.  Maxwell1 is the prime customer for this.
  bool inPlace = pL.get<bool>("repartition: use subcommunicators in place");
  if (inPlace) {
    SubFactoryMonitor subM(*this, "Rebalancing existing " + coarseMatrixName + " in-place", coarseLevel);
    RCP<const Map> newMap = Get<RCP<const Map> >(coarseLevel, "InPlaceMap");
 
    originalAc->removeEmptyProcessesInPlace(newMap);
 
    // The "in place" still leaves a dummy matrix here.  That needs to go
    if (newMap.is_null()) originalAc = Teuchos::null;
 
    Set(coarseLevel, matrixName, originalAc);
    return;
  }
 
  RCP<const Import> rebalanceImporter = Get<RCP<const Import> >(coarseLevel, "Importer");
 
  if (rebalanceImporter != Teuchos::null) {
    RCP<Matrix> rebalancedAc;
    {
      SubFactoryMonitor subM(*this, "Rebalancing existing " + coarseMatrixName, coarseLevel);
      RCP<const Map> targetMap = rebalanceImporter->getTargetMap();
 
      ParameterList XpetraList;
      if (pL.get<bool>("repartition: use subcommunicators") == true) {
        GetOStream(Runtime0) << "Replacing maps with a subcommunicator" << std::endl;
        XpetraList.set("Restrict Communicator", true);
      }
      // NOTE: If the communicator is restricted away, Build returns Teuchos::null.
      XpetraList.set("Timer Label", "MueLu::RebalanceAc-" + Teuchos::toString(coarseLevel.GetLevelID()));
      {
        SubFactoryMonitor subM2(*this, "Rebalancing existing " + coarseMatrixName + ": MatrixFactory::Build", coarseLevel);
        rebalancedAc = MatrixFactory::Build(originalAc, *rebalanceImporter, *rebalanceImporter, targetMap, targetMap, rcp(&XpetraList, false));
      }
 
      if (!rebalancedAc.is_null()) {
        if (originalAc->IsFixedBlockSizeSet())
          rebalancedAc->SetFixedBlockSize(originalAc->GetFixedBlockSize());
        std::ostringstream oss;
        oss << "A_" << coarseLevel.GetLevelID();
        rebalancedAc->setObjectLabel(oss.str());
      }
      Set(coarseLevel, matrixName, rebalancedAc);
    }
    if (!rebalancedAc.is_null() && IsPrint(Statistics2)) {
      int oldRank = SetProcRankVerbose(rebalancedAc->getRowMap()->getComm()->getRank());
 
      RCP<ParameterList> params = rcp(new ParameterList());
      params->set("printLoadBalancingInfo", true);
      params->set("printCommInfo", true);
      GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedAc, coarseMatrixName + " (rebalanced)", params);
 
      SetProcRankVerbose(oldRank);
    }
 
  } else {
    // Ac already built by the load balancing process and no load balancing needed
    GetOStream(Runtime1) << "No rebalancing" << std::endl;
    Set(coarseLevel, matrixName, originalAc);
  }
 
  if (rebalanceFacts_.begin() != rebalanceFacts_.end()) {
    SubFactoryMonitor m2(*this, "Rebalance additional data", coarseLevel);
 
    // call Build of all user-given transfer factories
    for (std::vector<RCP<const FactoryBase> >::const_iterator it = rebalanceFacts_.begin(); it != rebalanceFacts_.end(); ++it) {
      GetOStream(Runtime0) << "RebalanceAc: call rebalance factory " << (*it).get() << ": " << (*it)->description() << std::endl;
      (*it)->CallBuild(coarseLevel);
    }
  }
}  // Build()
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void RebalanceAcFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddRebalanceFactory(const RCP<const FactoryBase> &factory) {
  /*TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast,
                             "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. "
                             "This is very strange. (Note: you can remove this exception if there's a good reason for)");
  TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_, Exceptions::RuntimeError, "MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");*/
  rebalanceFacts_.push_back(factory);
}  // AddRebalanceFactory()
 
}  // namespace MueLu
 
#endif  // MUELU_REBALANCEACFACTORY_DEF_HPP