MueLu_MapTransferFactory_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_MAPTRANSFERFACTORY_DEF_HPP_
#define MUELU_MAPTRANSFERFACTORY_DEF_HPP_

#include "MueLu_MapTransferFactory_decl.hpp"

#include <Xpetra_Map.hpp>
#include <Xpetra_MapFactory.hpp>
#include <Xpetra_Matrix.hpp>

#include "MueLu_Level.hpp"
#include "MueLu_FactoryManagerBase.hpp"
#include "MueLu_Monitor.hpp"

namespace MueLu {

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const ParameterList> MapTransferFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
  RCP<ParameterList> validParamList = rcp(new ParameterList());

  validParamList->setEntry("map: name", Teuchos::ParameterEntry(std::string("")));
  validParamList->setEntry("map: factory", Teuchos::ParameterEntry(std::string("null")));

  validParamList->set<RCP<const FactoryBase>>("P", Teuchos::null, "Tentative prolongator factory");
  validParamList->set<std::string>("nullspace vectors: limit to", "all", "Limit the number of nullspace vectors to be used for the map transfer (especially to exclude rotational vectors).");

  return validParamList;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MapTransferFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& fineLevel, Level& coarseLevel) const {
  const ParameterList& pL       = GetParameterList();
  const std::string mapFactName = pL.get<std::string>("map: factory");
  const std::string mapName     = pL.get<std::string>("map: name");

  if (fineLevel.GetLevelID() == 0) {
    // Not needed, if the map is provided as user data
    fineLevel.DeclareInput(mapName, NoFactory::get(), this);
  } else {
    // check whether user has provided a specific name for the MapFactory
    if (mapFactName == "" || mapFactName == "NoFactory")
      mapFact_ = MueLu::NoFactory::getRCP();
    else if (mapFactName != "null")
      mapFact_ = coarseLevel.GetFactoryManager()->GetFactory(mapFactName);

    // request map generated by mapFact_
    fineLevel.DeclareInput(mapName, mapFact_.get(), this);
  }

  // request Ptent
  // note that "P" provided by the user (through XML file) is supposed to be of type TentativePFactory
  Teuchos::RCP<const FactoryBase> tentPFact = GetFactory("P");
  if (tentPFact == Teuchos::null)
    tentPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent");
  coarseLevel.DeclareInput("P", tentPFact.get(), this);
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MapTransferFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& fineLevel, Level& coarseLevel) const {
  Monitor m(*this, "MapTransferFactory");

  const ParameterList& pL     = GetParameterList();
  const std::string mapName   = pL.get<std::string>("map: name");
  const int maxNumProlongCols = GetLimitOfProlongatorColumns(pL);

  // fetch map from level
  RCP<const Map> transferMap = Teuchos::null;
  if (fineLevel.GetLevelID() == 0) {
    transferMap = fineLevel.Get<RCP<const Map>>(mapName, NoFactory::get());
  } else {
    if (fineLevel.IsAvailable(mapName, mapFact_.get()) == false)
      GetOStream(Runtime0) << "MapTransferFactory::Build: User provided map \"" << mapName << "\" not found in Level class on level " << fineLevel.GetLevelID() << "." << std::endl;
    transferMap = fineLevel.Get<RCP<const Map>>(mapName, mapFact_.get());
  }

  // Get default tentative prolongator factory
  // Getting it that way ensures that the same factory instance will be used for both SaPFactory and NullspaceFactory.
  RCP<const FactoryBase> tentPFact = GetFactory("P");
  if (tentPFact == Teuchos::null)
    tentPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent");
  TEUCHOS_TEST_FOR_EXCEPTION(!coarseLevel.IsAvailable("P", tentPFact.get()), Exceptions::RuntimeError,
                             "MueLu::MapTransferFactory::Build(): P (generated by TentativePFactory) not available.");
  RCP<Matrix> Ptent = coarseLevel.Get<RCP<Matrix>>("P", tentPFact.get());

  // loop over local rows of Ptent and figure out the corresponding coarse GIDs
  Array<GO> coarseMapGids;
  RCP<const Map> prolongColMap = Ptent->getColMap();
  GO gRowID                    = -1;
  int numColEntries            = 0;
  for (size_t row = 0; row < Ptent->getLocalNumRows(); ++row) {
    gRowID = Ptent->getRowMap()->getGlobalElement(row);

    if (transferMap->isNodeGlobalElement(gRowID)) {
      Teuchos::ArrayView<const LO> indices;
      Teuchos::ArrayView<const SC> vals;
      Ptent->getLocalRowView(row, indices, vals);

      numColEntries = as<int>(indices.size());
      if (maxNumProlongCols > 0)
        numColEntries = std::min(numColEntries, maxNumProlongCols);

      for (size_t col = 0; col < as<size_t>(numColEntries); ++col) {
        // mark all (selected) columns in Ptent(gRowID,*) to be coarse Dofs of next level transferMap
        GO gcid = prolongColMap->getGlobalElement(indices[col]);
        coarseMapGids.push_back(gcid);
      }
    }
  }

  // build coarse version of the input map
  const GO INVALID = Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
  std::sort(coarseMapGids.begin(), coarseMapGids.end());
  coarseMapGids.erase(std::unique(coarseMapGids.begin(), coarseMapGids.end()), coarseMapGids.end());
  RCP<const Map> coarseTransferMap = MapFactory::Build(prolongColMap->lib(), INVALID, coarseMapGids(),
                                                       prolongColMap->getIndexBase(), prolongColMap->getComm());

  // store map in coarse level
  if (fineLevel.GetLevelID() == 0) {
    const std::string mapFactName  = pL.get<std::string>("map: factory");
    RCP<const FactoryBase> mapFact = coarseLevel.GetFactoryManager()->GetFactory(mapFactName);
    coarseLevel.Set(mapName, coarseTransferMap, mapFact.get());
  } else
    coarseLevel.Set(mapName, coarseTransferMap, mapFact_.get());
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
int MapTransferFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetLimitOfProlongatorColumns(const ParameterList& pL) const {
  const std::string useTheseNspVectors = pL.get<std::string>("nullspace vectors: limit to");

  // Leave right away, if no limit is prescribed by the user
  if (useTheseNspVectors == "all" || useTheseNspVectors == "")
    return -1;

  // Simplify? Maybe replace by boolean flag "nullspace: exclude rotations"
  int maxNumProlongCols = -1;
  if (useTheseNspVectors == "translations")
    maxNumProlongCols = 1;
  else
    TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::InvalidArgument, "Unknown subset of nullspace vectors to be used, when performing a map transfer.")

  return maxNumProlongCols;
}

}  // namespace MueLu

#endif /* MUELU_MAPTRANSFERFACTORY_DEF_HPP_ */