MueLu_MapTransferFactory_def.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @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_ */