MueLu_IsorropiaInterface_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
 
/*
 * MueLu_IsorropiaInterface_def.hpp
 *
 *  Created on: Jun 10, 2013
 *      Author: tobias
 */
 
#ifndef MUELU_ISORROPIAINTERFACE_DEF_HPP_
#define MUELU_ISORROPIAINTERFACE_DEF_HPP_
 
#include "MueLu_IsorropiaInterface_decl.hpp"
 
#include <Teuchos_Utils.hpp>
//#include <Teuchos_DefaultMpiComm.hpp> //TODO: fwd decl.
//#include <Teuchos_OpaqueWrapper.hpp>  //TODO: fwd decl.
 
#include <Xpetra_MapFactory.hpp>
#include <Xpetra_CrsGraphFactory.hpp>
#include <Xpetra_BlockedMap.hpp>
#include <Xpetra_BlockedCrsMatrix.hpp>
 
#include "MueLu_Level.hpp"
#include "MueLu_Exceptions.hpp"
#include "MueLu_Monitor.hpp"
 
#include "MueLu_AmalgamationInfo.hpp"
 
namespace MueLu {
 
template <class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const ParameterList> IsorropiaInterface<LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
  RCP<ParameterList> validParamList = rcp(new ParameterList());
 
  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Factory of the matrix A");
  validParamList->set<RCP<const FactoryBase> >("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
  validParamList->set<RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Generating factory of UnAmalgamationInfo");
 
  return validParamList;
}
 
template <class LocalOrdinal, class GlobalOrdinal, class Node>
void IsorropiaInterface<LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
  Input(currentLevel, "A");
  Input(currentLevel, "number of partitions");
  Input(currentLevel, "UnAmalgamationInfo");
}
 
template <class LocalOrdinal, class GlobalOrdinal, class Node>
void IsorropiaInterface<LocalOrdinal, GlobalOrdinal, Node>::Build(Level& level) const {
  FactoryMonitor m(*this, "Build", level);
  typedef Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node> BlockMap;
 
  RCP<Matrix> A                  = Get<RCP<Matrix> >(level, "A");
  RCP<AmalgamationInfo> amalInfo = Get<RCP<AmalgamationInfo> >(level, "UnAmalgamationInfo");
  GO numParts                    = Get<int>(level, "number of partitions");
 
  RCP<const Map> rowMap = A->getRowMap();
  RCP<const Map> colMap = A->getColMap();
 
  if (numParts == 1 || numParts == -1) {
    // Running on one processor, so decomposition is the trivial one, all zeros.
    RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
    Set(level, "AmalgamatedPartition", decomposition);
    return;
  }
 
  // ok, reconstruct graph information of matrix A
  // Note, this is the non-rebalanced matrix A and we need the graph
  // of the non-rebalanced matrix A. We cannot make use of the "Graph"
  // that is/will be built for the aggregates later for several reasons
  // 1) the "Graph" for the aggregates is meant to be based on the rebalanced matrix A
  // 2) we cannot call a CoalesceDropFactory::Build here since this is very complicated and
  //    completely messes up the whole factory chain
  // 3) CoalesceDropFactory is meant to provide some minimal Graph information for the aggregation
  //    (LWGraph), but here we need the full CrsGraph for Isorropia as input
 
  // That is, why this code is somewhat redundant to CoalesceDropFactory
 
  LO blockdim       = 1;                       // block dim for fixed size blocks
  GO indexBase      = rowMap->getIndexBase();  // index base of maps
  GO offset         = 0;
  LO blockid        = -1;  // block id in strided map
  LO nStridedOffset = 0;   // DOF offset for strided block id "blockid" (default = 0)
  // LO stridedblocksize = blockdim; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
 
  // 1) check for blocking/striding information
  //    fill above variables
  if (A->IsView("stridedMaps") &&
      Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
    Xpetra::viewLabel_t oldView  = A->SwitchToView("stridedMaps");  // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
    RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
    TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null, Exceptions::BadCast, "MueLu::IsorropiaInterface::Build: cast to strided row map failed.");
    blockdim = strMap->getFixedBlockSize();
    offset   = strMap->getOffset();
    blockid  = strMap->getStridedBlockId();
    if (blockid > -1) {
      std::vector<size_t> stridingInfo = strMap->getStridingData();
      for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
        nStridedOffset += stridingInfo[j];
      // stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
 
    }  // else {
    //  stridedblocksize = blockdim;
    //}
    oldView = A->SwitchToView(oldView);
    // GetOStream(Statistics0) << "IsorropiaInterface::Build():" << " found blockdim=" << blockdim << " from strided maps (blockid=" << blockid << ", strided block size=" << stridedblocksize << "). offset=" << offset << std::endl;
  } else
    GetOStream(Statistics0) << "IsorropiaInterface::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
 
  // 2) get row map for amalgamated matrix (graph of A)
  //    with same distribution over all procs as row map of A
  RCP<const Map> nodeMap         = amalInfo->getNodeRowMap();
  RCP<const BlockedMap> bnodeMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(nodeMap);
  if (!bnodeMap.is_null()) nodeMap = bnodeMap->getMap();
 
  GetOStream(Statistics0) << "IsorropiaInterface:Build(): nodeMap " << nodeMap->getLocalNumElements() << "/" << nodeMap->getGlobalNumElements() << " elements" << std::endl;
 
  // 3) create graph of amalgamated matrix
  RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getLocalMaxNumRowEntries() * blockdim);
 
  // 4) do amalgamation. generate graph of amalgamated matrix
  for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); row++) {
    // get global DOF id
    GO grid = rowMap->getGlobalElement(row);
 
    // translate grid to nodeid
    // JHU 2019-20-May this is identical to AmalgamationFactory::DOFGid2NodeId(), and is done
    // to break a circular dependence when libraries are built statically
    GO nodeId = (grid - offset - indexBase) / blockdim + indexBase;
 
    size_t nnz = A->getNumEntriesInLocalRow(row);
    Teuchos::ArrayView<const LO> indices;
    Teuchos::ArrayView<const SC> vals;
    A->getLocalRowView(row, indices, vals);
 
    RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(new std::vector<GO>);  // global column block ids
    LO realnnz                     = 0;
    for (LO col = 0; col < Teuchos::as<LO>(nnz); col++) {
      GO gcid = colMap->getGlobalElement(indices[col]);  // global column id
 
      if (vals[col] != 0.0) {
        // JHU 2019-20-May this is identical to AmalgamationFactory::DOFGid2NodeId(), and is done
        // to break a circular dependence when libraries are built statically
        GO cnodeId = (gcid - offset - indexBase) / blockdim + indexBase;
        cnodeIds->push_back(cnodeId);
        realnnz++;  // increment number of nnz in matrix row
      }
    }
 
    Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp(cnodeIds);
 
    if (arr_cnodeIds.size() > 0)
      crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
  }
  // fill matrix graph
  crsGraph->fillComplete(nodeMap, nodeMap);
 
#ifdef HAVE_MPI
#else  // if we don't have MPI
 
  // Running on one processor, so decomposition is the trivial one, all zeros.
  RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
  Set(level, "AmalgamatedPartition", decomposition);
 
#endif  // HAVE_MPI
  // throw a more helpful error message if something failed
  // TEUCHOS_TEST_FOR_EXCEPTION(!level.IsAvailable("AmalgamatedPartition"), Exceptions::RuntimeError, "IsorropiaInterface::Build : no \'Partition\' vector available on level. Isorropia failed to build a partition of the non-repartitioned graph of A. Please make sure, that Isorropia is correctly compiled (Epetra/Tpetra).");
 
}  // Build()
 
}  // namespace MueLu
 
//#endif //if defined(HAVE_MUELU_ISORROPIA) && defined(HAVE_MPI)
 
#endif /* MUELU_ISORROPIAINTERFACE_DEF_HPP_ */