MueLu_CombinePFactory_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_COMBINEPFACTORY_DEF_HPP
#define MUELU_COMBINEPFACTORY_DEF_HPP

#include <stdlib.h>
#include <iomanip>

// #include <Teuchos_LAPACK.hpp>
#include <Teuchos_SerialDenseMatrix.hpp>
#include <Teuchos_SerialDenseVector.hpp>
#include <Teuchos_SerialDenseSolver.hpp>

#include <Xpetra_CrsMatrixWrap.hpp>
#include <Xpetra_ImportFactory.hpp>
#include <Xpetra_Matrix.hpp>
#include <Xpetra_MapFactory.hpp>
#include <Xpetra_MultiVectorFactory.hpp>
#include <Xpetra_VectorFactory.hpp>

#include <Xpetra_IO.hpp>

#include "MueLu_FactoryManagerBase.hpp"
#include "MueLu_CombinePFactory_decl.hpp"

#include "MueLu_MasterList.hpp"
#include "MueLu_Monitor.hpp"

namespace MueLu {

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const ParameterList> CombinePFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
  RCP<ParameterList> validParamList = rcp(new ParameterList());
  validParamList->setEntry("combine: numBlks", ParameterEntry(1));
  validParamList->setEntry("combine: useMaxLevels", ParameterEntry(false));
  validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory of the matrix A");

  return validParamList;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void CombinePFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
    DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {
  //    Input(fineLevel, "subblock");
}

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

namespace {
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
constructIdentityProlongator(const Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& map) {
  using local_matrix_type = typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_type;
  using local_graph_type  = typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_graph_type;
  using row_map_type      = typename local_matrix_type::row_map_type::non_const_type;
  using entries_type      = typename local_graph_type::entries_type::non_const_type;
  using values_type       = typename local_matrix_type::values_type;
  using local_scalar_type = typename values_type::value_type;

  LocalOrdinal numLocal = map->getLocalNumElements();
  row_map_type rowptr("rowptr", numLocal + 1);
  entries_type colind("colind", numLocal);
  values_type values("values", numLocal);

  Kokkos::parallel_for(
      "Setup CRS values for identity prolongator",
      Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>(0, numLocal),
      KOKKOS_LAMBDA(const LocalOrdinal index) {
        rowptr(index) = index;
        colind(index) = index;
        values(index) = local_scalar_type{1.0};
        if (index == (numLocal - 1)) {
          rowptr(numLocal) = numLocal;
        }
      });

  auto eye = Xpetra::CrsMatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(map, map, 0);
  eye->setAllValues(rowptr, colind, values);
  eye->expertStaticFillComplete(map, map);
  return Teuchos::make_rcp<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(eye);
}
}  // namespace

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

  const ParameterList& pL = GetParameterList();
  const LO nBlks          = as<LO>(pL.get<int>("combine: numBlks"));
  const bool useMaxLevels = pL.get<bool>("combine: useMaxLevels");

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

  // Record all matrices that each define a block in block diagonal comboP
  // matrix used for PDE/multiblock interpolation.  Additionally, count
  // total number of local rows, nonzeros, coarseDofs, and colDofs.

  Teuchos::ArrayRCP<RCP<Matrix>> arrayOfMatrices(nBlks);
  size_t nComboRowMap = 0, nnzCombo = 0, nComboColMap = 0, nComboDomMap = 0;

  LO nTotalNumberLocalColMapEntries = 0;
  Teuchos::ArrayRCP<size_t> DomMapSizePerBlk(nBlks);
  Teuchos::ArrayRCP<size_t> ColMapSizePerBlk(nBlks);
  Teuchos::ArrayRCP<size_t> ColMapLocalSizePerBlk(nBlks);
  Teuchos::ArrayRCP<size_t> ColMapRemoteSizePerBlk(nBlks);
  Teuchos::ArrayRCP<size_t> ColMapLocalCumulativePerBlk(nBlks + 1);   // hardwire 0th entry so that it has the value of 0
  Teuchos::ArrayRCP<size_t> ColMapRemoteCumulativePerBlk(nBlks + 1);  // hardwire 0th entry so that it has the value of 0

  bool anyCoarseGridsRemaining = false;

  if (useMaxLevels) {
    for (int j = 0; j < nBlks; j++) {
      std::string blockName = "Psubblock" + Teuchos::toString(j);
      anyCoarseGridsRemaining |= coarseLevel.IsAvailable(blockName, NoFactory::get());
    };

    int localAnyCoarseGridsRemaining  = anyCoarseGridsRemaining;
    int globalAnyCoarseGridsRemaining = localAnyCoarseGridsRemaining;
    Teuchos::reduceAll(*A->getDomainMap()->getComm(), Teuchos::REDUCE_MAX, localAnyCoarseGridsRemaining, Teuchos::ptr(&globalAnyCoarseGridsRemaining));

    anyCoarseGridsRemaining |= globalAnyCoarseGridsRemaining > 0;
  }

  auto setSubblockProlongator = [&](RCP<Matrix> Psubblock, int j) {
    arrayOfMatrices[j] = Psubblock;
    nComboRowMap += Teuchos::as<size_t>((arrayOfMatrices[j])->getRowMap()->getLocalNumElements());
    DomMapSizePerBlk[j] = Teuchos::as<size_t>((arrayOfMatrices[j])->getDomainMap()->getLocalNumElements());
    ColMapSizePerBlk[j] = Teuchos::as<size_t>((arrayOfMatrices[j])->getColMap()->getLocalNumElements());
    nComboDomMap += DomMapSizePerBlk[j];
    nComboColMap += ColMapSizePerBlk[j];
    nnzCombo += Teuchos::as<size_t>((arrayOfMatrices[j])->getLocalNumEntries());
    TEUCHOS_TEST_FOR_EXCEPTION((arrayOfMatrices[j])->getDomainMap()->getIndexBase() != 0, Exceptions::RuntimeError, "interpolation subblocks must use 0 indexbase");

    // figure out how many empty entries in each column map
    int tempii = 0;
    for (int i = 0; i < (int)DomMapSizePerBlk[j]; i++) {
      if ((arrayOfMatrices[j])->getDomainMap()->getGlobalElement(i) == (arrayOfMatrices[j])->getColMap()->getGlobalElement(tempii)) tempii++;
    }
    nTotalNumberLocalColMapEntries += tempii;
    ColMapLocalSizePerBlk[j]  = tempii;
    ColMapRemoteSizePerBlk[j] = ColMapSizePerBlk[j] - ColMapLocalSizePerBlk[j];
  };

  for (int j = 0; j < nBlks; j++) {
    std::string blockName = "Psubblock" + Teuchos::toString(j);
    if (coarseLevel.IsAvailable(blockName, NoFactory::get())) {
      setSubblockProlongator(coarseLevel.Get<RCP<Matrix>>(blockName, NoFactory::get()), j);
    } else {
      std::string subblockOpName = "Operatorsubblock" + Teuchos::toString(j);
      bool hasOperator           = false;
      if (coarseLevel.IsAvailable(subblockOpName)) {
        auto A_blk  = coarseLevel.Get<RCP<Operator>>(subblockOpName);
        hasOperator = A_blk != Teuchos::null;
      }

      if (useMaxLevels && anyCoarseGridsRemaining && hasOperator) {
        // Use Psubblock = I
        auto P_id = constructIdentityProlongator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(coarseLevel.Get<RCP<Operator>>(subblockOpName)->getDomainMap());
        setSubblockProlongator(P_id, j);
      } else {
        arrayOfMatrices[j]        = Teuchos::null;
        ColMapLocalSizePerBlk[j]  = 0;
        ColMapRemoteSizePerBlk[j] = 0;
      }
    }
    ColMapLocalCumulativePerBlk[j + 1]  = ColMapLocalSizePerBlk[j] + ColMapLocalCumulativePerBlk[j];
    ColMapRemoteCumulativePerBlk[j + 1] = ColMapRemoteSizePerBlk[j] + ColMapRemoteCumulativePerBlk[j];
  }

  TEUCHOS_TEST_FOR_EXCEPTION(nComboRowMap != A->getRowMap()->getLocalNumElements(), Exceptions::RuntimeError, "sum of subblock rows != #row's Afine");

  // build up csr arrays for combo block diagonal P
  Teuchos::ArrayRCP<size_t> comboPRowPtr(nComboRowMap + 1);
  Teuchos::ArrayRCP<LocalOrdinal> comboPCols(nnzCombo);
  Teuchos::ArrayRCP<Scalar> comboPVals(nnzCombo);

  size_t nnzCnt = 0, nrowCntFromPrevBlks = 0, ncolCntFromPrevBlks = 0;
  LO maxNzPerRow = 0;
  for (int j = 0; j < nBlks; j++) {
    // grab csr pointers for individual blocks of P
    if (arrayOfMatrices[j] != Teuchos::null) {
      Teuchos::ArrayRCP<const size_t> subblockRowPtr((arrayOfMatrices[j])->getLocalNumRows());
      Teuchos::ArrayRCP<const LocalOrdinal> subblockCols((arrayOfMatrices[j])->getLocalNumEntries());
      Teuchos::ArrayRCP<const Scalar> subblockVals((arrayOfMatrices[j])->getLocalNumEntries());
      if ((int)(arrayOfMatrices[j])->getLocalMaxNumRowEntries() > maxNzPerRow) maxNzPerRow = (int)(arrayOfMatrices[j])->getLocalMaxNumRowEntries();
      Teuchos::RCP<CrsMatrixWrap> subblockwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>((arrayOfMatrices[j]));
      Teuchos::RCP<CrsMatrix> subblockcrs      = subblockwrap->getCrsMatrix();
      subblockcrs->getAllValues(subblockRowPtr, subblockCols, subblockVals);

      // copy jth block into csr arrays of comboP

      for (decltype(subblockRowPtr.size()) i = 0; i < subblockRowPtr.size() - 1; i++) {
        size_t rowLength                      = subblockRowPtr[i + 1] - subblockRowPtr[i];
        comboPRowPtr[nrowCntFromPrevBlks + i] = nnzCnt;
        for (size_t k = 0; k < rowLength; k++) {
          if ((int)subblockCols[k + subblockRowPtr[i]] < (int)ColMapLocalSizePerBlk[j]) {
            comboPCols[nnzCnt] = subblockCols[k + subblockRowPtr[i]] + ColMapLocalCumulativePerBlk[j];
            if ((int)comboPCols[nnzCnt] >= (int)ColMapLocalCumulativePerBlk[nBlks]) {
              printf("ERROR1\n");
              exit(1);
            }
          } else {
            // Here we subtract off the number of local colmap guys ... so this tell us where we are among ghost unknowns. We then want to stick this ghost after
            // all the Local guys ... so we add ColMapLocalCumulativePerBlk[nBlks] .... finally we need to account for the fact that other ghost blocks may have already
            // been handled ... so we then add  + ColMapRemoteCumulativePerBlk[j];
            comboPCols[nnzCnt] = subblockCols[k + subblockRowPtr[i]] - ColMapLocalSizePerBlk[j] + ColMapLocalCumulativePerBlk[nBlks] + ColMapRemoteCumulativePerBlk[j];
            if ((int)comboPCols[nnzCnt] >= (int)(ColMapLocalCumulativePerBlk[nBlks] + ColMapRemoteCumulativePerBlk[nBlks])) {
              printf("ERROR2\n");
              exit(1);
            }
          }
          comboPVals[nnzCnt++] = subblockVals[k + subblockRowPtr[i]];
        }
      }

      nrowCntFromPrevBlks += Teuchos::as<size_t>((arrayOfMatrices[j])->getRowMap()->getLocalNumElements());
      ncolCntFromPrevBlks += DomMapSizePerBlk[j];  // rst: check this
    }
  }
  comboPRowPtr[nComboRowMap] = nnzCnt;

  // Come up with global IDs for the coarse grid maps. We assume that each xxx
  // block has a minimum GID of 0.  Since MueLu is generally creating these
  // GIDS, this assumption is probably correct, but we'll check it.

  Teuchos::Array<GlobalOrdinal> comboDomainMapGIDs(nComboDomMap);
  Teuchos::Array<GlobalOrdinal> comboColMapGIDs(nComboColMap);

  LO nTotalNumberRemoteColMapEntries = 0;
  GlobalOrdinal offset               = 0;
  size_t domainMapIndex              = 0;
  int nComboColIndex                 = 0;
  for (int j = 0; j < nBlks; j++) {
    int nThisBlkColIndex = 0;

    GlobalOrdinal tempMax = 0, maxGID = 0;
    if (arrayOfMatrices[j] != Teuchos::null) tempMax = (arrayOfMatrices[j])->getDomainMap()->getMaxGlobalIndex();
    Teuchos::reduceAll(*(A->getDomainMap()->getComm()), Teuchos::REDUCE_MAX, tempMax, Teuchos::ptr(&maxGID));

    if (arrayOfMatrices[j] != Teuchos::null) {
      TEUCHOS_TEST_FOR_EXCEPTION(arrayOfMatrices[j]->getDomainMap()->getMinAllGlobalIndex() < 0, Exceptions::RuntimeError, "Global ID assumption for domainMap not met within subblock");

      GO priorDomGID = 0;
      for (size_t c = 0; c < DomMapSizePerBlk[j]; ++c) {  // check this
                                                          //  The global ids of jth block are assumed to go between 0 and maxGID_j.  So the 1th blocks DomainGIDs should start at maxGID_0+1. The 2nd
                                                          //  block DomainDIGS starts at maxGID_0+1 + maxGID_1 + 1. We use offset to keep track of these starting points.
        priorDomGID                          = (arrayOfMatrices[j])->getDomainMap()->getGlobalElement(c);
        comboDomainMapGIDs[domainMapIndex++] = offset + priorDomGID;
        if (priorDomGID == (arrayOfMatrices[j])->getColMap()->getGlobalElement(nThisBlkColIndex)) {
          comboColMapGIDs[nComboColIndex++] = offset + priorDomGID;
          nThisBlkColIndex++;
        }
      }

      for (size_t cc = nThisBlkColIndex; cc < ColMapSizePerBlk[j]; ++cc) {
        comboColMapGIDs[nTotalNumberLocalColMapEntries + nTotalNumberRemoteColMapEntries++] = offset + (arrayOfMatrices[j])->getColMap()->getGlobalElement(cc);
      }
    }
    offset += maxGID + 1;
  }
#ifdef out
  RCP<const Map> coarseDomainMap = Xpetra::MapFactory<LO, GO, NO>::Build(A->getDomainMap()->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), comboDomainMapGIDs, 0, A->getDomainMap()->getComm());
  RCP<const Map> coarseColMap    = Xpetra::MapFactory<LO, GO, NO>::Build(A->getDomainMap()->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), comboColMapGIDs, 0, A->getDomainMap()->getComm());
#endif

  RCP<const Map> coarseDomainMap = Xpetra::MapFactory<LO, GO, NO>::Build(A->getDomainMap()->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), comboDomainMapGIDs, 0, A->getRowMap()->getComm());
  RCP<const Map> coarseColMap    = Xpetra::MapFactory<LO, GO, NO>::Build(A->getDomainMap()->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), comboColMapGIDs, 0, A->getRowMap()->getComm());

  Teuchos::RCP<CrsMatrix> comboPCrs = CrsMatrixFactory::Build(A->getRowMap(), coarseColMap, maxNzPerRow);
  // comboPCrs->getCrsGraph(); //.getRowInfo(6122);
  // comboPCrs->getRowInfo(6122);

  //    Teuchos::RCP<CrsMatrix> comboPCrs = CrsMatrixFactory::Build(A->getRowMap(), coarseColMap,nnzCombo+1000);

  //    for (size_t i = 0; i < nComboRowMap; i++) {
  // printf("FIXME\n"); if (nComboRowMap > 6142)  nComboRowMap = 6142;
  for (size_t i = 0; i < nComboRowMap; i++) {
    comboPCrs->insertLocalValues(i, comboPCols.view(comboPRowPtr[i], comboPRowPtr[i + 1] - comboPRowPtr[i]),
                                 comboPVals.view(comboPRowPtr[i], comboPRowPtr[i + 1] - comboPRowPtr[i]));
  }
  comboPCrs->fillComplete(coarseDomainMap, A->getRowMap());

  Teuchos::RCP<Matrix> comboP = Teuchos::rcp(new CrsMatrixWrap(comboPCrs));

  Set(coarseLevel, "P", comboP);
}

}  // namespace MueLu

#define MUELU_COMBINEPFACTORY_SHORT
#endif  // MUELU_COMBINEPFACTORY_DEF_HPP