119 for(
size_t i = 0; i<input.
getColMap()->getLocalNumElements(); i++) {
120 GlobalOrdinal gcid = input.
getColMap()->getGlobalElement(i);
128 std::vector<int> myUnknownDofGIDs(comm->getSize(),0);
129 std::vector<int> numUnknownDofGIDs(comm->getSize(),0);
130 myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.
size();
131 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myUnknownDofGIDs[0],&numUnknownDofGIDs[0]);
134 size_t cntUnknownDofGIDs = 0;
135 for(
int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
136 std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs,-1);
137 std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs,-1);
139 size_t cntUnknownOffset = 0;
140 for(
int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
141 for(
size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.
size()); k++) {
142 lUnknownDofGIDs[k+cntUnknownOffset] = ovlUnknownStatusGids[k];
144 if(cntUnknownDofGIDs > 0)
145 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lUnknownDofGIDs[0],&gUnknownDofGIDs[0]);
146 std::sort(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end());
147 gUnknownDofGIDs.erase(std::unique(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end()), gUnknownDofGIDs.end());
150 for(
size_t k=0; k < gUnknownDofGIDs.size(); k++) {
151 GlobalOrdinal curgid = gUnknownDofGIDs[k];
157 std::vector<int> myFoundDofGIDs(comm->getSize(),0);
158 std::vector<int> numFoundDofGIDs(comm->getSize(),0);
159 myFoundDofGIDs[comm->getRank()] = ovlFoundStatusGids.
size();
160 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myFoundDofGIDs[0],&numFoundDofGIDs[0]);
163 size_t cntFoundDofGIDs = 0;
164 for(
int p = 0; p < comm->getSize(); p++) cntFoundDofGIDs += numFoundDofGIDs[p];
165 std::vector<GlobalOrdinal> lFoundDofGIDs(cntFoundDofGIDs,-1);
166 std::vector<GlobalOrdinal> gFoundDofGIDs(cntFoundDofGIDs,-1);
168 size_t cntFoundOffset = 0;
169 for(
int p = 0; p < comm->getRank(); p++) cntFoundOffset += numFoundDofGIDs[p];
170 for(
size_t k=0; k < Teuchos::as<size_t>(ovlFoundStatusGids.
size()); k++) {
171 lFoundDofGIDs[k+cntFoundOffset] = ovlFoundStatusGids[k];
173 if(cntFoundDofGIDs > 0)
174 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntFoundDofGIDs),&lFoundDofGIDs[0],&gFoundDofGIDs[0]);
177 for(
size_t i = 0; i<input.
getColMap()->getLocalNumElements(); i++) {
178 GlobalOrdinal gcid = input.
getColMap()->getGlobalElement(i);
180 std::find(gFoundDofGIDs.begin(), gFoundDofGIDs.end(), gcid) != gFoundDofGIDs.end()) {
203 bool bThyraMode =
false) {
205 size_t numRows = rangeMapExtractor->NumMaps();
206 size_t numCols = domainMapExtractor->NumMaps();
208 TEUCHOS_TEST_FOR_EXCEPTION(rangeMapExtractor->getThyraMode() ==
true,
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor must not use Thyra style numbering of GIDs. The MapExtractor must contain all GIDs of the full range map in order to define a proper splitting.")
209 TEUCHOS_TEST_FOR_EXCEPTION(domainMapExtractor->getThyraMode() ==
true,
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor must not use Thyra style numbering of GIDs. The MapExtractor must contain all GIDs of the full domain map in order to define a proper splitting.")
221 TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getMaxAllGlobalIndex() != input.
getColMap()->getMaxAllGlobalIndex(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix. fullDomainMap->getMaxAllGlobalIndex() = " << fullDomainMap->getMaxAllGlobalIndex() <<
" vs. input.getColMap()->getMaxAllGlobalIndex() = " << input.
getColMap()->getMaxAllGlobalIndex())
228 if(columnMapExtractor == Teuchos::null) {
231 std::vector<Teuchos::RCP<const Map> > ovlxmaps(numCols, Teuchos::null);
232 for(
size_t c = 0; c < numCols; c++) {
235 ovlxmaps[c] = colMap;
241 myColumnMapExtractor = columnMapExtractor;
248 if(bThyraMode ==
true) {
250 std::vector<Teuchos::RCP<const Map> > thyRgMapExtractorMaps(numRows, Teuchos::null);
251 for (
size_t r = 0; r < numRows; r++) {
255 if(strRangeMap != Teuchos::null) {
256 std::vector<size_t> strInfo = strRangeMap->getStridingData();
257 GlobalOrdinal offset = strRangeMap->getOffset();
258 LocalOrdinal stridedBlockId = strRangeMap->getStridedBlockId();
260 thyRgMapExtractorMaps[r] = strShrinkedMap;
262 thyRgMapExtractorMaps[r] = shrinkedMap;
268 std::vector<Teuchos::RCP<const Map> > thyDoMapExtractorMaps (numCols, Teuchos::null);
269 std::vector<Teuchos::RCP<const Map> > thyColMapExtractorMaps(numCols, Teuchos::null);
270 for (
size_t c = 0; c < numCols; c++) {
275 if(strDomainMap != Teuchos::null) {
276 std::vector<size_t> strInfo = strDomainMap->getStridingData();
277 GlobalOrdinal offset = strDomainMap->getOffset();
278 LocalOrdinal stridedBlockId = strDomainMap->getStridedBlockId();
280 thyDoMapExtractorMaps[c] = strShrinkedDomainMap;
282 thyDoMapExtractorMaps[c] = shrinkedDomainMap;
287 if(strColMap != Teuchos::null) {
288 std::vector<size_t> strInfo = strColMap->getStridingData();
289 GlobalOrdinal offset = strColMap->getOffset();
290 LocalOrdinal stridedBlockId = strColMap->getStridedBlockId();
292 thyColMapExtractorMaps[c] = strShrinkedColMap;
294 thyColMapExtractorMaps[c] = shrinkedColMap;
306 std::vector<Teuchos::RCP<Matrix> > subMatrices(numRows*numCols, Teuchos::null);
307 for (
size_t r = 0; r < numRows; r++) {
308 for (
size_t c = 0; c < numCols; c++) {
312 if(bThyraMode ==
true)
332 doCheck->putScalar(1.0);
338 doCheck->putScalar(-1.0);
339 coCheck->putScalar(-1.0);
342 for (
size_t rrr = 0; rrr < input.
getDomainMap()->getLocalNumElements(); rrr++) {
344 GlobalOrdinal
id = input.
getDomainMap()->getGlobalElement(rrr);
347 size_t BlockId = domainMapExtractor->getMapIndexForGID(
id);
349 doCheckData[rrr] = Teuchos::as<Scalar>(BlockId);
357 for (
size_t rr = 0; rr < input.
getRowMap()->getLocalNumElements(); rr++) {
360 GlobalOrdinal growid = input.
getRowMap()->getGlobalElement(rr);
369 size_t rowBlockId = rangeMapExtractor->getMapIndexForGID(growid);
372 GlobalOrdinal subblock_growid = growid;
373 if(bThyraMode ==
true) {
375 LocalOrdinal lrowid = rangeMapExtractor->getMap(rowBlockId)->getLocalElement(growid);
377 subblock_growid = thyRangeMapExtractor->getMap(rowBlockId,
true)->getGlobalElement(lrowid);
390 for(
size_t i=0; i<(size_t)indices.
size(); i++) {
392 GlobalOrdinal gcolid = input.
getColMap()->getGlobalElement(indices[i]);
394 size_t colBlockId = myColumnMapExtractor->getMapIndexForGID(gcolid);
398 GlobalOrdinal subblock_gcolid = gcolid;
399 if(bThyraMode ==
true) {
401 LocalOrdinal lcolid = myColumnMapExtractor->getMap(colBlockId)->getLocalElement(gcolid);
403 subblock_gcolid = thyColMapExtractor->getMap(colBlockId,
true)->getGlobalElement(lcolid);
405 blockColIdx [colBlockId].push_back(subblock_gcolid);
406 blockColVals[colBlockId].push_back(vals[i]);
409 for (
size_t c = 0; c < numCols; c++) {
410 subMatrices[rowBlockId*numCols+c]->insertGlobalValues(subblock_growid,blockColIdx[c].view(0,blockColIdx[c].size()),blockColVals[c].view(0,blockColVals[c].size()));
417 if(bThyraMode ==
true) {
418 for (
size_t r = 0; r < numRows; r++) {
419 for (
size_t c = 0; c < numCols; c++) {
420 subMatrices[r*numCols+c]->fillComplete(thyDomainMapExtractor->getMap(c,
true), thyRangeMapExtractor->getMap(r,
true));
425 for (
size_t r = 0; r < numRows; r++) {
426 for (
size_t c = 0; c < numCols; c++) {
427 subMatrices[r*numCols+c]->fillComplete(domainMapExtractor->getMap(c), rangeMapExtractor->getMap(r));
433 for (
size_t r = 0; r < numRows; r++) {
434 for (
size_t c = 0; c < numCols; c++) {
435 bA->setMatrix(r,c,subMatrices[r*numCols+c]);
449 using Teuchos::rcp_dynamic_cast;
451 GlobalOrdinal gZeroDiags;
452 bool usedEfficientPath =
false;
454#ifdef HAVE_MUELU_TPETRA
458 tpCrsAc = rcp_dynamic_cast<TpetraCrsMatrix>(crsWrapAc->getCrsMatrix());
461 auto tpCrsGraph = tpCrsAc->getTpetra_CrsMatrix()->getCrsGraph();
462 size_t numRows = Ac->getRowMap()->getLocalNumElements();
463 typedef typename Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::offset_device_view_type offset_type;
464 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
465 auto offsets = offset_type(Kokkos::ViewAllocateWithoutInitializing(
"offsets"), numRows);
466 tpCrsGraph->getLocalDiagOffsets(offsets);
469 Tpetra::Details::OrdinalTraits<typename offset_type::value_type>::invalid ();
471 if (repairZeroDiagonals) {
475 LO numMissingDiagonalEntries = 0;
477 Kokkos::parallel_reduce(
"countMissingDiagonalEntries",
478 range_type(0, numRows),
479 KOKKOS_LAMBDA (
const LO i, LO& missing) {
480 missing += (offsets(i) == STINV);
481 }, numMissingDiagonalEntries);
483 GlobalOrdinal gNumMissingDiagonalEntries;
484 Teuchos::reduceAll(*(Ac->getRowMap()->getComm()), Teuchos::REDUCE_SUM, Teuchos::as<GlobalOrdinal>(numMissingDiagonalEntries),
485 Teuchos::outArg(gNumMissingDiagonalEntries));
487 if (gNumMissingDiagonalEntries == 0) {
490 auto lclA = tpCrsAc->getTpetra_CrsMatrix()->getLocalMatrixDevice();
492 using ATS = Kokkos::ArithTraits<Scalar>;
493 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
496 typename ATS::val_type impl_replacementValue = replacementValue;
498 Kokkos::parallel_reduce(
"fixSmallDiagonalEntries",
499 range_type(0, numRows),
500 KOKKOS_LAMBDA (
const LO i, LO& fixed) {
501 const auto offset = offsets(i);
502 auto curRow = lclA.row (i);
503 if (impl_ATS::magnitude(curRow.value(offset)) <= threshold) {
504 curRow.value(offset) = impl_replacementValue;
509 Teuchos::reduceAll(*(Ac->getRowMap()->getComm()), Teuchos::REDUCE_SUM, Teuchos::as<GlobalOrdinal>(lZeroDiags),
510 Teuchos::outArg(gZeroDiags));
512 usedEfficientPath =
true;
518 auto lclA = tpCrsAc->getTpetra_CrsMatrix()->getLocalMatrixDevice();
520 using ATS = Kokkos::ArithTraits<Scalar>;
521 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
525 Kokkos::parallel_reduce(
"detectSmallDiagonalEntries",
526 range_type(0, numRows),
527 KOKKOS_LAMBDA (
const LO i, LO& small) {
528 const auto offset = offsets(i);
532 auto curRow = lclA.row (i);
533 if (impl_ATS::magnitude(curRow.value(offset)) <= threshold) {
539 Teuchos::reduceAll(*(Ac->getRowMap()->getComm()), Teuchos::REDUCE_SUM, Teuchos::as<GlobalOrdinal>(lZeroDiags),
540 Teuchos::outArg(gZeroDiags));
542 usedEfficientPath =
true;
548 if (!usedEfficientPath) {
551 Ac->getLocalDiagCopy(*diagVec);
553 LocalOrdinal lZeroDiags = 0;
556 for (
size_t i = 0; i < rowMap->getLocalNumElements(); i++) {
557 if (TST::magnitude(diagVal[i]) <= threshold) {
562 Teuchos::reduceAll(*(rowMap->getComm()), Teuchos::REDUCE_SUM, Teuchos::as<GlobalOrdinal>(lZeroDiags),
563 Teuchos::outArg(gZeroDiags));
565 if (repairZeroDiagonals && gZeroDiags > 0) {
583 for (
size_t r = 0; r < rowMap->getLocalNumElements(); r++) {
584 if (TST::magnitude(diagVal[r]) <= threshold) {
585 GlobalOrdinal grid = rowMap->getGlobalElement(r);
587 valout[0] = replacementValue;
588 fixDiagMatrix->insertGlobalValues(grid,indout(), valout());
593 fixDiagMatrix->fillComplete(Ac->getDomainMap(),Ac->getRangeMap());
597 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*Ac,
false, 1.0, *fixDiagMatrix,
false, 1.0, newAc, fos);
598 if (Ac->IsView(
"stridedMaps"))
599 newAc->CreateView(
"stridedMaps", Ac);
602 fixDiagMatrix = Teuchos::null;
608 p->set(
"DoOptimizeStorage",
true);
617 fos <<
"CheckRepairMainDiagonal: " << (repairZeroDiagonals ?
"repaired " :
"found ")
618 << gZeroDiags <<
" too small entries (threshold = " << threshold <<
") on main diagonal of Ac." << std::endl;
620#ifdef HAVE_XPETRA_DEBUG
626 Ac->getLocalDiagCopy(*diagVec);
627 diagVal = diagVec->getData(0);
628 for (
size_t r = 0; r < Ac->getRowMap()->getLocalNumElements(); r++) {
629 if (TST::magnitude(diagVal[r]) <= threshold) {
630 fos <<
"Error: there are too small entries left on diagonal after repair..." << std::endl;