MueLu_CoarseMapFactory_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
/*
 * MueLu_CoarseMapFactory_def.hpp
 *
 *  Created on: Oct 12, 2012
 *      Author: wiesner
 */
 
#ifndef MUELU_COARSEMAPFACTORY_DEF_HPP_
#define MUELU_COARSEMAPFACTORY_DEF_HPP_
 
#include <Teuchos_Array.hpp>
 
#include <Xpetra_MultiVector.hpp>
#include <Xpetra_StridedMapFactory.hpp>
 
#include "MueLu_CoarseMapFactory_decl.hpp"
#include "MueLu_Level.hpp"
#include "MueLu_Aggregates.hpp"
#include "MueLu_Monitor.hpp"
 
namespace MueLu {
 
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  RCP<const ParameterList> CoarseMapFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const
  {
    RCP<ParameterList> validParamList = rcp(new ParameterList());
 
    validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Generating factory for aggregates.");
    validParamList->set< RCP<const FactoryBase> >("Nullspace",  Teuchos::null, "Generating factory for null space.");
 
    validParamList->set< std::string  >("Striding info", "{}", "Striding information");
    validParamList->set< LocalOrdinal >("Strided block id", -1, "Strided block id");
 
    // Domain GID offsets: list of offset values for domain gids (coarse gids) of tentative prolongator  (default = 0).
    // For each multigrid level we can define a different offset value, ie for the prolongator between
    // level 0 and level 1 we use the offset entry domainGidOffsets_[0], for the prolongator between
    // level 1 and level 2 we use domainGidOffsets_[1]...
    // If the vector domainGidOffsets_ is too short and does not contain a sufficient number of gid offset
    // values for all levels, a gid offset of 0 is used for all the remaining levels!
    // The GIDs for the domain dofs of Ptent start with domainGidOffset, are contiguous and distributed
    // equally over the procs (unless some reordering is done).
    validParamList->set< std::string > ("Domain GID offsets", "{0}", "vector with offsets for GIDs for each level. If no offset GID value is given for the level we use 0 as default.");
 
    return validParamList;
  }
 
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void CoarseMapFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &currentLevel) const
  {
    Input(currentLevel, "Aggregates");
    Input(currentLevel, "Nullspace");
  }
 
  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
  void CoarseMapFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setStridingData(std::vector<size_t> stridingInfo)
  {
    // store striding map in internal variable
    stridingInfo_ = stridingInfo;
 
    // try to remove string "Striding info" from parameter list to make sure,
    // that stridingInfo_ is not replaced in the Build call.
    std::string strStridingInfo; strStridingInfo.clear();
    SetParameter("Striding info", ParameterEntry(strStridingInfo));
  }
 
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void CoarseMapFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &currentLevel) const
  {
    FactoryMonitor m(*this, "Build", currentLevel);
 
    GlobalOrdinal domainGIDOffset = GetDomainGIDOffset(currentLevel);
    BuildCoarseMap(currentLevel, domainGIDOffset);
  }
 
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void CoarseMapFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCoarseMap(
    Level& currentLevel, const GlobalOrdinal domainGIDOffset) const
  {
 
    RCP<Aggregates> aggregates = Get< RCP<Aggregates> >(currentLevel, "Aggregates");
    GlobalOrdinal numAggs = aggregates->GetNumAggregates();
    RCP<const Map> aggMap = aggregates->GetMap();
 
    RCP<MultiVector> nullspace = Get< RCP<MultiVector> >(currentLevel, "Nullspace");
 
    const size_t NSDim = nullspace->getNumVectors();
    RCP<const Teuchos::Comm<int> > comm = aggMap->getComm();
    const ParameterList & pL = GetParameterList();
 
    LocalOrdinal stridedBlockId = pL.get<LocalOrdinal>("Strided block id");
 
    // read in stridingInfo from parameter list and fill the internal member variable
    // read the data only if the parameter "Striding info" exists and is non-empty
    if(pL.isParameter("Striding info")) {
      std::string strStridingInfo = pL.get<std::string>("Striding info");
      if(strStridingInfo.empty() == false) {
        Teuchos::Array<size_t> arrayVal = Teuchos::fromStringToArray<size_t>(strStridingInfo);
        stridingInfo_ = Teuchos::createVector(arrayVal);
      }
    }
 
    CheckForConsistentStridingInformation(stridedBlockId, NSDim);
 
    GetOStream(Statistics2) << "domainGIDOffset: " << domainGIDOffset << " block size: " << getFixedBlockSize() << " stridedBlockId: " << stridedBlockId << std::endl;
 
    // number of coarse level dofs (fixed by number of aggregates and blocksize data)
    GlobalOrdinal nCoarseDofs = numAggs * getFixedBlockSize();
    GlobalOrdinal indexBase = aggMap->getIndexBase();
 
    RCP<const Map> coarseMap = StridedMapFactory::Build(aggMap->lib(),
        Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
        nCoarseDofs,
        indexBase,
        stridingInfo_,
        comm,
        stridedBlockId,
        domainGIDOffset);
 
    Set(currentLevel, "CoarseMap", coarseMap);
  }
 
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  GlobalOrdinal CoarseMapFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetDomainGIDOffset(
    Level& currentLevel) const
  {
    GlobalOrdinal domainGidOffset = Teuchos::ScalarTraits<GlobalOrdinal>::zero();
 
    std::vector<GlobalOrdinal> domainGidOffsets;
    domainGidOffsets.clear();
    const ParameterList & pL = GetParameterList();
    if(pL.isParameter("Domain GID offsets")) {
      std::string strDomainGIDs = pL.get<std::string>("Domain GID offsets");
      if(strDomainGIDs.empty() == false) {
        Teuchos::Array<GlobalOrdinal> arrayVal = Teuchos::fromStringToArray<GlobalOrdinal>(strDomainGIDs);
        domainGidOffsets = Teuchos::createVector(arrayVal);
        if(currentLevel.GetLevelID() < Teuchos::as<int>(domainGidOffsets.size()) ) {
          domainGidOffset = domainGidOffsets[currentLevel.GetLevelID()];
        }
      }
    }
 
    return domainGidOffset;
  }
 
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void CoarseMapFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CheckForConsistentStridingInformation(
    const LocalOrdinal stridedBlockId, const size_t nullspaceDimension) const
  {
    // check for consistency of striding information with NSDim and nCoarseDofs
    if (stridedBlockId == -1) {
      // this means we have no real strided map but only a block map with constant blockSize "nullspaceDimension"
      TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo_.size() > 1, Exceptions::RuntimeError, "MueLu::CoarseMapFactory::Build(): stridingInfo_.size() but must be one");
      stridingInfo_.clear();
      stridingInfo_.push_back(nullspaceDimension);
      TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo_.size() != 1, Exceptions::RuntimeError, "MueLu::CoarseMapFactory::Build(): stridingInfo_.size() but must be one");
 
    } else {
      // stridedBlockId > -1, set by user
      TEUCHOS_TEST_FOR_EXCEPTION(stridedBlockId > Teuchos::as<LO>(stridingInfo_.size() - 1), Exceptions::RuntimeError, "MueLu::CoarseMapFactory::Build(): it is stridingInfo_.size() <= stridedBlockId_. error.");
      size_t stridedBlockSize = stridingInfo_[stridedBlockId];
      TEUCHOS_TEST_FOR_EXCEPTION(stridedBlockSize != nullspaceDimension , Exceptions::RuntimeError, "MueLu::CoarseMapFactory::Build(): dimension of strided block != nullspaceDimension. error.");
    }
  }
 
} //namespace MueLu
 
#endif /* MUELU_COARSEMAPFACTORY_DEF_HPP_ */