46 #ifndef MUELU_COMBINEPFACTORY_DEF_HPP
47 #define MUELU_COMBINEPFACTORY_DEF_HPP
57 #include <Xpetra_CrsMatrixWrap.hpp>
60 #include <Xpetra_MapFactory.hpp>
61 #include <Xpetra_MultiVectorFactory.hpp>
74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 return validParamList;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
91 Level& coarseLevel)
const {
92 return BuildP(fineLevel, coarseLevel);
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 Level& coarseLevel)
const {
101 const LO nBlks = as<LO>(pL.
get<
int>(
"combine: numBlks"));
110 size_t nComboRowMap = 0, nnzCombo = 0, nComboColMap = 0, nComboDomMap = 0;
112 LO nTotalNumberLocalColMapEntries = 0;
119 for (
int j = 0; j < nBlks; j++) {
123 nComboRowMap += Teuchos::as<size_t>((arrayOfMatrices[j])->getRowMap()->getLocalNumElements());
124 DomMapSizePerBlk[j] = Teuchos::as<size_t>((arrayOfMatrices[j])->getDomainMap()->getLocalNumElements());
125 ColMapSizePerBlk[j] = Teuchos::as<size_t>((arrayOfMatrices[j])->getColMap()->getLocalNumElements());
126 nComboDomMap += DomMapSizePerBlk[j];
127 nComboColMap += ColMapSizePerBlk[j];
128 nnzCombo += Teuchos::as<size_t>((arrayOfMatrices[j])->getLocalNumEntries());
133 for (
int i = 0; i < (int)DomMapSizePerBlk[j]; i++) {
135 if ((arrayOfMatrices[j])->getDomainMap()->getGlobalElement(i) == (arrayOfMatrices[j])->getColMap()->getGlobalElement(tempii)) tempii++;
137 nTotalNumberLocalColMapEntries += tempii;
138 ColMapLocalSizePerBlk[j] = tempii;
139 ColMapRemoteSizePerBlk[j] = ColMapSizePerBlk[j] - ColMapLocalSizePerBlk[j];
141 arrayOfMatrices[j] = Teuchos::null;
142 ColMapLocalSizePerBlk[j] = 0;
143 ColMapRemoteSizePerBlk[j] = 0;
145 ColMapLocalCumulativePerBlk[j + 1] = ColMapLocalSizePerBlk[j] + ColMapLocalCumulativePerBlk[j];
146 ColMapRemoteCumulativePerBlk[j + 1] = ColMapRemoteSizePerBlk[j] + ColMapRemoteCumulativePerBlk[j];
155 size_t nnzCnt = 0, nrowCntFromPrevBlks = 0, ncolCntFromPrevBlks = 0;
157 for (
int j = 0; j < nBlks; j++) {
159 if (arrayOfMatrices[j] != Teuchos::null) {
163 if ((
int)(arrayOfMatrices[j])->getLocalMaxNumRowEntries() > maxNzPerRow) maxNzPerRow = (int)(arrayOfMatrices[j])->getLocalMaxNumRowEntries();
166 subblockcrs->getAllValues(subblockRowPtr, subblockCols, subblockVals);
170 for (decltype(subblockRowPtr.
size()) i = 0; i < subblockRowPtr.
size() - 1; i++) {
171 size_t rowLength = subblockRowPtr[i + 1] - subblockRowPtr[i];
172 comboPRowPtr[nrowCntFromPrevBlks + i] = nnzCnt;
173 for (
size_t k = 0; k < rowLength; k++) {
174 if ((
int)subblockCols[k + subblockRowPtr[i]] < (int)ColMapLocalSizePerBlk[j]) {
175 comboPCols[nnzCnt] = subblockCols[k + subblockRowPtr[i]] + ColMapLocalCumulativePerBlk[j];
176 if ((
int)comboPCols[nnzCnt] >= (int)ColMapLocalCumulativePerBlk[nBlks]) {
184 comboPCols[nnzCnt] = subblockCols[k + subblockRowPtr[i]] - ColMapLocalSizePerBlk[j] + ColMapLocalCumulativePerBlk[nBlks] + ColMapRemoteCumulativePerBlk[j];
185 if ((
int)comboPCols[nnzCnt] >= (
int)(ColMapLocalCumulativePerBlk[nBlks] + ColMapRemoteCumulativePerBlk[nBlks])) {
190 comboPVals[nnzCnt++] = subblockVals[k + subblockRowPtr[i]];
194 nrowCntFromPrevBlks += Teuchos::as<size_t>((arrayOfMatrices[j])->getRowMap()->getLocalNumElements());
195 ncolCntFromPrevBlks += DomMapSizePerBlk[j];
198 comboPRowPtr[nComboRowMap] = nnzCnt;
207 LO nTotalNumberRemoteColMapEntries = 0;
209 size_t domainMapIndex = 0;
210 int nComboColIndex = 0;
211 for (
int j = 0; j < nBlks; j++) {
212 int nThisBlkColIndex = 0;
215 if (arrayOfMatrices[j] != Teuchos::null) tempMax = (arrayOfMatrices[j])->getDomainMap()->getMaxGlobalIndex();
216 Teuchos::reduceAll(*(A->getDomainMap()->getComm()),
Teuchos::REDUCE_MAX, tempMax, Teuchos::ptr(&maxGID));
218 if (arrayOfMatrices[j] != Teuchos::null) {
222 for (
size_t c = 0; c < DomMapSizePerBlk[j]; ++c) {
225 priorDomGID = (arrayOfMatrices[j])->getDomainMap()->getGlobalElement(c);
226 comboDomainMapGIDs[domainMapIndex++] = offset + priorDomGID;
227 if (priorDomGID == (arrayOfMatrices[j])->getColMap()->getGlobalElement(nThisBlkColIndex)) {
228 comboColMapGIDs[nComboColIndex++] = offset + priorDomGID;
233 for (
size_t cc = nThisBlkColIndex; cc < ColMapSizePerBlk[j]; ++cc) {
234 comboColMapGIDs[nTotalNumberLocalColMapEntries + nTotalNumberRemoteColMapEntries++] = offset + (arrayOfMatrices[j])->getColMap()->getGlobalElement(cc);
237 offset += maxGID + 1;
255 for (
size_t i = 0; i < nComboRowMap; i++) {
256 comboPCrs->insertLocalValues(i, comboPCols.
view(comboPRowPtr[i], comboPRowPtr[i + 1] - comboPRowPtr[i]),
257 comboPVals.
view(comboPRowPtr[i], comboPRowPtr[i + 1] - comboPRowPtr[i]));
259 comboPCrs->fillComplete(coarseDomainMap, A->getRowMap());
263 Set(coarseLevel,
"P", comboP);
268 #define MUELU_COMBINEPFACTORY_SHORT
269 #endif // MUELU_COMBINEPFACTORY_DEF_HPP
ParameterList & setEntry(const std::string &name, U &&entry)
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static const NoFactory * get()
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception throws to report errors in the internal logical of the program.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
std::string toString(const T &t)
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
ArrayView< T > view(size_type lowerOffset, size_type size) const