MueLu_BlockedPFactory_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_BLOCKEDPFACTORY_DEF_HPP_
#define MUELU_BLOCKEDPFACTORY_DEF_HPP_

#include <Xpetra_MultiVectorFactory.hpp>
#include <Xpetra_VectorFactory.hpp>
#include <Xpetra_ImportFactory.hpp>
#include <Xpetra_ExportFactory.hpp>
#include <Xpetra_CrsMatrixWrap.hpp>

#include <Xpetra_BlockedCrsMatrix.hpp>
#include <Xpetra_Map.hpp>
#include <Xpetra_MapFactory.hpp>
#include <Xpetra_MapExtractor.hpp>
#include <Xpetra_MapExtractorFactory.hpp>

#include "MueLu_BlockedPFactory_decl.hpp"
#include "MueLu_TentativePFactory.hpp"
#include "MueLu_FactoryBase.hpp"
#include "MueLu_SmootherFactory.hpp"
#include "MueLu_FactoryManager.hpp"
#include "MueLu_Utilities.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_HierarchyUtils.hpp"

namespace MueLu {

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

    validParamList->set< RCP<const FactoryBase> >("A",          null, "Generating factory of the matrix A (block matrix)");
    validParamList->set< bool >                  ("backwards", false, "Forward/backward order");

    return validParamList;
  }

  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
  void BlockedPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &coarseLevel) const {
    Input(fineLevel, "A");

    const ParameterList& pL = GetParameterList();
    const bool backwards = pL.get<bool>("backwards");

    const int numFactManagers = FactManager_.size();
    for (int k = 0; k < numFactManagers; k++) {
      int i = (backwards ? numFactManagers-1 - k : k);
      const RCP<const FactoryManagerBase>& factManager = FactManager_[i];

      SetFactoryManager fineSFM  (rcpFromRef(fineLevel),   factManager);
      SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), factManager);

      if (!restrictionMode_)
        coarseLevel.DeclareInput("P", factManager->GetFactory("P").get(), this);
      else
        coarseLevel.DeclareInput("R", factManager->GetFactory("R").get(), this);
    }
  }

  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
  void BlockedPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level& fineLevel, Level& coarseLevel) const
  { }

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

    RCP<Matrix> Ain = Get< RCP<Matrix> >(fineLevel, "A");

    RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
    TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::BadCast, "Input matrix A is not a BlockedCrsMatrix.");

    const int numFactManagers = FactManager_.size();

    // Plausibility check
    TEUCHOS_TEST_FOR_EXCEPTION(A->Rows() != as<size_t>(numFactManagers), Exceptions::RuntimeError,
                               "Number of block rows [" << A->Rows() << "] does not match the number of SubFactorManagers [" << numFactManagers << "]");
    TEUCHOS_TEST_FOR_EXCEPTION(A->Cols() != as<size_t>(numFactManagers), Exceptions::RuntimeError,
                               "Number of block cols [" << A->Cols() << "] does not match the number of SubFactorManagers [" << numFactManagers << "]");

    // Build blocked prolongator
    std::vector<RCP<Matrix> >     subBlockP          (numFactManagers);
    std::vector<RCP<const Map> >  subBlockPRangeMaps (numFactManagers);
    std::vector<RCP<const Map> >  subBlockPDomainMaps(numFactManagers);

    std::vector<GO> fullRangeMapVector;
    std::vector<GO> fullDomainMapVector;

    RCP<const MapExtractor> rangeAMapExtractor = A->getRangeMapExtractor();
    RCP<const MapExtractor> domainAMapExtractor = A->getDomainMapExtractor();

    const ParameterList& pL = GetParameterList();
    const bool backwards = pL.get<bool>("backwards");

    // Build and store the subblocks and the corresponding range and domain
    // maps.  Since we put together the full range and domain map from the
    // submaps, we do not have to use the maps from blocked A
    for (int k = 0; k < numFactManagers; k++) {
      int i = (backwards ? numFactManagers-1 - k : k);
      const RCP<const FactoryManagerBase>& factManager = FactManager_[i];

      SetFactoryManager fineSFM  (rcpFromRef(fineLevel),   factManager);
      SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), factManager);

      if (!restrictionMode_) subBlockP[i] = coarseLevel.Get<RCP<Matrix> >("P", factManager->GetFactory("P").get());
      else                   subBlockP[i] = coarseLevel.Get<RCP<Matrix> >("R", factManager->GetFactory("R").get());

      // Check if prolongator/restrictor operators have strided maps
      TEUCHOS_TEST_FOR_EXCEPTION(subBlockP[i]->IsView("stridedMaps") == false, Exceptions::BadCast,
                                 "subBlock P operator [" << i << "] has no strided map information.");

      // Append strided row map (= range map) to list of range maps.
      // TAW the row map GIDs extracted here represent the actual row map GIDs.
      // No support for nested operators! (at least if Thyra style gids are used)
      RCP<const StridedMap> strPartialMap = Teuchos::rcp_dynamic_cast<const StridedMap>(subBlockP[i]->getRowMap("stridedMaps"));
      std::vector<size_t> stridedRgData = strPartialMap->getStridingData();
      subBlockPRangeMaps[i] = StridedMapFactory::Build(
          subBlockP[i]->getRangeMap(), // actual global IDs (Thyra or Xpetra)
          stridedRgData,
          strPartialMap->getStridedBlockId(),
          strPartialMap->getOffset());
      //subBlockPRangeMaps[i] = subBlockP[i]->getRowMap("stridedMaps");

      // Use plain range map to determine the DOF ids
      ArrayView<const GO> nodeRangeMap = subBlockPRangeMaps[i]->getNodeElementList();
      fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());

      // Append strided col map (= domain map) to list of range maps.
      RCP<const StridedMap> strPartialMap2 = Teuchos::rcp_dynamic_cast<const StridedMap>(subBlockP[i]->getColMap("stridedMaps"));
      std::vector<size_t> stridedRgData2 = strPartialMap2->getStridingData();
      subBlockPDomainMaps[i] = StridedMapFactory::Build(
          subBlockP[i]->getDomainMap(), // actual global IDs (Thyra or Xpetra)
          stridedRgData2,
          strPartialMap2->getStridedBlockId(),
          strPartialMap2->getOffset());
      //subBlockPDomainMaps[i] = subBlockP[i]->getColMap("stridedMaps");

      // Use plain domain map to determine the DOF ids
      ArrayView<const GO> nodeDomainMap = subBlockPDomainMaps[i]->getNodeElementList();
      fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
    }

    // check if GIDs for full maps have to be sorted:
    // For the Thyra mode ordering they do not have to be sorted since the GIDs are
    // numbered as 0...n1,0...,n2 (starting with zero for each subblock). The MapExtractor
    // generates unique GIDs during the construction.
    // For Xpetra style, the GIDs have to be reordered. Such that one obtains a ordered
    // list of GIDs in an increasing ordering. In Xpetra, the GIDs are all unique through
    // out all submaps.
    bool bDoSortRangeGIDs  = rangeAMapExtractor->getThyraMode()  ? false : true;
    bool bDoSortDomainGIDs = domainAMapExtractor->getThyraMode() ? false : true;

    if (bDoSortRangeGIDs) {
      sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
      fullRangeMapVector.erase(std::unique(fullRangeMapVector.begin(), fullRangeMapVector.end()), fullRangeMapVector.end());
    }
    if (bDoSortDomainGIDs) {
      sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
      fullDomainMapVector.erase(std::unique(fullDomainMapVector.begin(), fullDomainMapVector.end()), fullDomainMapVector.end());
    }

    // extract map index base from maps of blocked A
    GO rangeIndexBase  = 0;
    GO domainIndexBase = 0;
    if (!restrictionMode_) {
      // Prolongation mode: just use index base of range and domain map of A
      rangeIndexBase  = A->getRangeMap() ->getIndexBase();
      domainIndexBase = A->getDomainMap()->getIndexBase();

    } else {
      // Restriction mode: switch range and domain map for blocked restriction operator
      rangeIndexBase  = A->getDomainMap()->getIndexBase();
      domainIndexBase = A->getRangeMap()->getIndexBase();
    }

    // Build full range map.
    // If original range map has striding information, then transfer it to the
    // new range map
    RCP<const StridedMap>   stridedRgFullMap = rcp_dynamic_cast<const StridedMap>(rangeAMapExtractor->getFullMap());
    RCP<const Map >         fullRangeMap = Teuchos::null;

    ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0, fullRangeMapVector.size());
    if (stridedRgFullMap != Teuchos::null) {
      std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
      fullRangeMap = StridedMapFactory::Build(
                                   A->getRangeMap()->lib(),
                                   Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
                                   fullRangeMapGIDs,
                                   rangeIndexBase,
                                   stridedData,
                                   A->getRangeMap()->getComm(),
                                   -1, /* the full map vector should always have strided block id -1! */
                                   stridedRgFullMap->getOffset());
    } else {
      fullRangeMap = MapFactory::Build(
                            A->getRangeMap()->lib(),
                            Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
                            fullRangeMapGIDs,
                            rangeIndexBase,
                            A->getRangeMap()->getComm());
    }

    ArrayView<GO>         fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
    RCP<const StridedMap> stridedDoFullMap = rcp_dynamic_cast<const StridedMap>(domainAMapExtractor->getFullMap());
    RCP<const Map >       fullDomainMap    = null;
    if (stridedDoFullMap != null) {
      TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap == Teuchos::null, Exceptions::BadCast,
        "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");

      std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
      fullDomainMap = StridedMapFactory::Build(
                                   A->getDomainMap()->lib(),
                                   Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
                                   fullDomainMapGIDs,
                                   domainIndexBase,
                                   stridedData2,
                                   A->getDomainMap()->getComm(),
                                   -1, /* the full map vector should always have strided block id -1! */
                                   stridedDoFullMap->getOffset());
    } else {
      fullDomainMap = MapFactory::Build(
                            A->getDomainMap()->lib(),
                            Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
                            fullDomainMapGIDs,
                            domainIndexBase,
                            A->getDomainMap()->getComm());
    }

    // Build map extractors
    RCP<const MapExtractor> rangeMapExtractor  = MapExtractorFactory::Build(fullRangeMap,  subBlockPRangeMaps, rangeAMapExtractor->getThyraMode());
    RCP<const MapExtractor> domainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, domainAMapExtractor->getThyraMode());

    RCP<BlockedCrsMatrix> P = rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10));
    for (size_t i = 0; i < subBlockPRangeMaps.size(); i++)
      for (size_t j = 0; j < subBlockPRangeMaps.size(); j++)
        if (i == j) {
          RCP<CrsMatrixWrap> crsOpii  = rcp_dynamic_cast<CrsMatrixWrap>(subBlockP[i]);
          TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null, Xpetra::Exceptions::BadCast,
                                     "Block [" << i << ","<< j << "] must be of type CrsMatrixWrap.");
          P->setMatrix(i, i, crsOpii);
        } else {
          P->setMatrix(i, j, Teuchos::null);
        }

    P->fillComplete();

    // Level Set
    if (!restrictionMode_) {
      // Prolongation mode
      coarseLevel.Set("P", rcp_dynamic_cast<Matrix>(P), this);

    } else {
      // Restriction mode
      // We do not have to transpose the blocked R operator since the subblocks
      // on the diagonal are already valid R subblocks
      coarseLevel.Set("R", Teuchos::rcp_dynamic_cast<Matrix>(P), this);
    }

  }

} // namespace MueLu

#endif /* MUELU_BLOCKEDPFACTORY_DEF_HPP_ */