47 #ifndef PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
48 #define PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
52 #include "Xpetra_Map.hpp"
54 #include "Xpetra_StridedMap.hpp"
55 #include "Xpetra_MapFactory.hpp"
56 #include "Xpetra_MapExtractor.hpp"
65 #ifdef HAVE_XPETRA_TPETRA
66 #include "Xpetra_TpetraMultiVector.hpp"
67 #include <Tpetra_RowMatrixTransposer.hpp>
68 #include <Tpetra_Details_extractBlockDiagonal.hpp>
69 #include <Tpetra_Details_scaleBlockDiagonal.hpp>
83 template <
class Scalar,
88 #undef XPETRA_MATRIXUTILS_SHORT
97 Teuchos::ArrayRCP<const Scalar> data = input.
getData(c);
99 ret->replaceLocalValue(Teuchos::as<LocalOrdinal>(r), c, data[r]);
109 RCP<const Teuchos::Comm<int>> comm = input.
getRowMap()->getComm();
112 Teuchos::Array<GlobalOrdinal> ovlUnknownStatusGids;
113 Teuchos::Array<GlobalOrdinal> ovlFoundStatusGids;
116 for (
size_t i = 0; i < input.
getColMap()->getLocalNumElements(); i++) {
117 GlobalOrdinal gcid = input.
getColMap()->getGlobalElement(i);
119 ovlUnknownStatusGids.push_back(gcid);
125 std::vector<int> myUnknownDofGIDs(comm->getSize(), 0);
126 std::vector<int> numUnknownDofGIDs(comm->getSize(), 0);
127 myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.size();
128 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, comm->getSize(), &myUnknownDofGIDs[0], &numUnknownDofGIDs[0]);
131 size_t cntUnknownDofGIDs = 0;
132 for (
int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
133 std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs, -1);
134 std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs, -1);
136 size_t cntUnknownOffset = 0;
137 for (
int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
138 for (
size_t k = 0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
139 lUnknownDofGIDs[k + cntUnknownOffset] = ovlUnknownStatusGids[k];
141 if (cntUnknownDofGIDs > 0)
142 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(cntUnknownDofGIDs), &lUnknownDofGIDs[0], &gUnknownDofGIDs[0]);
143 std::sort(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end());
144 gUnknownDofGIDs.erase(std::unique(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end()), gUnknownDofGIDs.end());
147 for (
size_t k = 0; k < gUnknownDofGIDs.size(); k++) {
148 GlobalOrdinal curgid = gUnknownDofGIDs[k];
150 ovlFoundStatusGids.push_back(curgid);
154 std::vector<int> myFoundDofGIDs(comm->getSize(), 0);
155 std::vector<int> numFoundDofGIDs(comm->getSize(), 0);
156 myFoundDofGIDs[comm->getRank()] = ovlFoundStatusGids.size();
157 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, comm->getSize(), &myFoundDofGIDs[0], &numFoundDofGIDs[0]);
160 size_t cntFoundDofGIDs = 0;
161 for (
int p = 0; p < comm->getSize(); p++) cntFoundDofGIDs += numFoundDofGIDs[p];
162 std::vector<GlobalOrdinal> lFoundDofGIDs(cntFoundDofGIDs, -1);
163 std::vector<GlobalOrdinal> gFoundDofGIDs(cntFoundDofGIDs, -1);
165 size_t cntFoundOffset = 0;
166 for (
int p = 0; p < comm->getRank(); p++) cntFoundOffset += numFoundDofGIDs[p];
167 for (
size_t k = 0; k < Teuchos::as<size_t>(ovlFoundStatusGids.size()); k++) {
168 lFoundDofGIDs[k + cntFoundOffset] = ovlFoundStatusGids[k];
170 if (cntFoundDofGIDs > 0)
171 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(cntFoundDofGIDs), &lFoundDofGIDs[0], &gFoundDofGIDs[0]);
173 Teuchos::Array<GlobalOrdinal> ovlDomainMapArray;
174 for (
size_t i = 0; i < input.
getColMap()->getLocalNumElements(); i++) {
175 GlobalOrdinal gcid = input.
getColMap()->getGlobalElement(i);
177 std::find(gFoundDofGIDs.begin(), gFoundDofGIDs.end(), gcid) != gFoundDofGIDs.end()) {
178 ovlDomainMapArray.push_back(gcid);
181 RCP<Map> ovlDomainMap =
MapFactory::Build(domainMap.
lib(), Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(), ovlDomainMapArray(), 0, comm);
195 static Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
SplitMatrix(
200 bool bThyraMode =
false) {
201 size_t numRows = rangeMapExtractor->NumMaps();
202 size_t numCols = domainMapExtractor->NumMaps();
204 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.")
205 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.")
207 RCP<const Map> fullRangeMap = rangeMapExtractor->getFullMap();
208 RCP<const Map> fullDomainMap = domainMapExtractor->getFullMap();
210 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getMaxAllGlobalIndex() != input.
getRowMap()->getMaxAllGlobalIndex(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
211 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getGlobalNumElements() != input.
getRowMap()->getGlobalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
212 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getLocalNumElements() != input.
getRowMap()->getLocalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
213 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getMaxAllGlobalIndex() != input.
getRangeMap()->getMaxAllGlobalIndex(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
214 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getGlobalNumElements() != input.
getRangeMap()->getGlobalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
215 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getLocalNumElements() != input.
getRangeMap()->getLocalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
217 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())
218 TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getMaxAllGlobalIndex() != input.
getDomainMap()->getMaxAllGlobalIndex(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
219 TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getGlobalNumElements() != input.
getDomainMap()->getGlobalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
220 TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getLocalNumElements() != input.
getDomainMap()->getLocalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
223 Teuchos::RCP<const MapExtractor> myColumnMapExtractor = Teuchos::null;
224 if (columnMapExtractor == Teuchos::null) {
225 TEUCHOS_TEST_FOR_EXCEPTION(domainMapExtractor->getThyraMode() ==
true,
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: Auto generation of column map extractor not supported for Thyra style numbering.");
227 std::vector<Teuchos::RCP<const Map>> ovlxmaps(numCols, Teuchos::null);
228 for (
size_t c = 0; c < numCols; c++) {
231 ovlxmaps[c] = colMap;
237 myColumnMapExtractor = columnMapExtractor;
241 Teuchos::RCP<const MapExtractor> thyRangeMapExtractor = Teuchos::null;
242 Teuchos::RCP<const MapExtractor> thyDomainMapExtractor = Teuchos::null;
243 Teuchos::RCP<const MapExtractor> thyColMapExtractor = Teuchos::null;
244 if (bThyraMode ==
true) {
246 std::vector<Teuchos::RCP<const Map>> thyRgMapExtractorMaps(numRows, Teuchos::null);
247 for (
size_t r = 0; r < numRows; r++) {
248 RCP<const Map> rMap = rangeMapExtractor->getMap(r);
250 RCP<const StridedMap> strRangeMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rMap);
251 if (strRangeMap != Teuchos::null) {
253 GlobalOrdinal offset = strRangeMap->getOffset();
254 LocalOrdinal stridedBlockId = strRangeMap->getStridedBlockId();
255 RCP<const StridedMap> strShrinkedMap = Teuchos::rcp(
new StridedMap(shrinkedMap, strInfo, shrinkedMap->getIndexBase(), stridedBlockId, offset));
256 thyRgMapExtractorMaps[r] = strShrinkedMap;
258 thyRgMapExtractorMaps[r] = shrinkedMap;
260 TEUCHOS_TEST_FOR_EXCEPTION(thyRgMapExtractorMaps[r]->getLocalNumElements() != rMap->getLocalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: Thyra-style range map extractor contains faulty data.")
264 std::vector<Teuchos::RCP<const Map>> thyDoMapExtractorMaps(numCols, Teuchos::null);
265 std::vector<Teuchos::RCP<const Map>> thyColMapExtractorMaps(numCols, Teuchos::null);
266 for (
size_t c = 0; c < numCols; c++) {
267 RCP<const Map> cMap = domainMapExtractor->getMap(c);
270 RCP<const StridedMap> strDomainMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(cMap);
271 if (strDomainMap != Teuchos::null) {
273 GlobalOrdinal offset = strDomainMap->getOffset();
274 LocalOrdinal stridedBlockId = strDomainMap->getStridedBlockId();
275 RCP<const StridedMap> strShrinkedDomainMap = Teuchos::rcp(
new StridedMap(shrinkedDomainMap, strInfo, shrinkedDomainMap->getIndexBase(), stridedBlockId, offset));
276 thyDoMapExtractorMaps[c] = strShrinkedDomainMap;
278 thyDoMapExtractorMaps[c] = shrinkedDomainMap;
280 RCP<const Map> colMap = myColumnMapExtractor->getMap(c);
282 RCP<const StridedMap> strColMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(colMap);
283 if (strColMap != Teuchos::null) {
285 GlobalOrdinal offset = strColMap->getOffset();
286 LocalOrdinal stridedBlockId = strColMap->getStridedBlockId();
287 RCP<const StridedMap> strShrinkedColMap = Teuchos::rcp(
new StridedMap(shrinkedColMap, strInfo, shrinkedColMap->getIndexBase(), stridedBlockId, offset));
288 thyColMapExtractorMaps[c] = strShrinkedColMap;
290 thyColMapExtractorMaps[c] = shrinkedColMap;
293 TEUCHOS_TEST_FOR_EXCEPTION(thyColMapExtractorMaps[c]->getLocalNumElements() != colMap->getLocalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: Thyra-style column map extractor contains faulty data.")
294 TEUCHOS_TEST_FOR_EXCEPTION(thyDoMapExtractorMaps[c]->getLocalNumElements() != cMap->getLocalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: Thyra-style domain map extractor contains faulty data.")
302 std::vector<Teuchos::RCP<Matrix>> subMatrices(numRows * numCols, Teuchos::null);
303 for (
size_t r = 0; r < numRows; r++) {
304 for (
size_t c = 0; c < numCols; c++) {
308 if (bThyraMode ==
true)
321 #if 0 // TAW needs to be fixed (does not compile for Scalar=complex)
328 doCheck->putScalar(1.0);
332 TEUCHOS_TEST_FOR_EXCEPTION(coCheck->normInf() != Teuchos::ScalarTraits< Scalar >::magnitude(1.0),
Xpetra::Exceptions::RuntimeError,
"Xpetra::MatrixUtils::SplitMatrix: error when distributing data.");
334 doCheck->putScalar(-1.0);
335 coCheck->putScalar(-1.0);
337 Teuchos::ArrayRCP< Scalar > doCheckData = doCheck->getDataNonConst(0);
338 for (
size_t rrr = 0; rrr < input.
getDomainMap()->getLocalNumElements(); rrr++) {
340 GlobalOrdinal
id = input.
getDomainMap()->getGlobalElement(rrr);
343 size_t BlockId = domainMapExtractor->getMapIndexForGID(
id);
345 doCheckData[rrr] = Teuchos::as<Scalar>(BlockId);
350 Teuchos::ArrayRCP< Scalar > coCheckData = coCheck->getDataNonConst(0);
353 for (
size_t rr = 0; rr < input.
getRowMap()->getLocalNumElements(); rr++) {
355 GlobalOrdinal growid = input.
getRowMap()->getGlobalElement(rr);
364 size_t rowBlockId = rangeMapExtractor->getMapIndexForGID(growid);
367 GlobalOrdinal subblock_growid = growid;
368 if (bThyraMode ==
true) {
370 LocalOrdinal lrowid = rangeMapExtractor->getMap(rowBlockId)->getLocalElement(growid);
372 subblock_growid = thyRangeMapExtractor->getMap(rowBlockId,
true)->getGlobalElement(lrowid);
378 Teuchos::ArrayView<const LocalOrdinal> indices;
379 Teuchos::ArrayView<const Scalar> vals;
382 std::vector<Teuchos::Array<GlobalOrdinal>> blockColIdx(numCols, Teuchos::Array<GlobalOrdinal>());
383 std::vector<Teuchos::Array<Scalar>> blockColVals(numCols, Teuchos::Array<Scalar>());
385 for (
size_t i = 0; i < (size_t)indices.size(); i++) {
387 GlobalOrdinal gcolid = input.
getColMap()->getGlobalElement(indices[i]);
389 size_t colBlockId = myColumnMapExtractor->getMapIndexForGID(gcolid);
393 GlobalOrdinal subblock_gcolid = gcolid;
394 if (bThyraMode ==
true) {
396 LocalOrdinal lcolid = myColumnMapExtractor->getMap(colBlockId)->getLocalElement(gcolid);
398 subblock_gcolid = thyColMapExtractor->getMap(colBlockId,
true)->getGlobalElement(lcolid);
400 blockColIdx[colBlockId].push_back(subblock_gcolid);
401 blockColVals[colBlockId].push_back(vals[i]);
404 for (
size_t c = 0; c < numCols; c++) {
405 subMatrices[rowBlockId * numCols + c]->insertGlobalValues(subblock_growid, blockColIdx[c].view(0, blockColIdx[c].size()), blockColVals[c].view(0, blockColVals[c].size()));
410 RCP<BlockedCrsMatrix> bA = Teuchos::null;
411 if (bThyraMode ==
true) {
412 for (
size_t r = 0; r < numRows; r++) {
413 for (
size_t c = 0; c < numCols; c++) {
414 subMatrices[r * numCols + c]->fillComplete(thyDomainMapExtractor->getMap(c,
true), thyRangeMapExtractor->getMap(r,
true));
417 bA = Teuchos::rcp(
new BlockedCrsMatrix(thyRangeMapExtractor, thyDomainMapExtractor, 10 ));
419 for (
size_t r = 0; r < numRows; r++) {
420 for (
size_t c = 0; c < numCols; c++) {
421 subMatrices[r * numCols + c]->fillComplete(domainMapExtractor->getMap(c), rangeMapExtractor->getMap(r));
424 bA = Teuchos::rcp(
new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10 ));
427 for (
size_t r = 0; r < numRows; r++) {
428 for (
size_t c = 0; c < numCols; c++) {
429 bA->setMatrix(r, c, subMatrices[r * numCols + c]);
438 bool const& repairZeroDiagonals, Teuchos::FancyOStream& fos,
439 const typename Teuchos::ScalarTraits<Scalar>::magnitudeType threshold = Teuchos::ScalarTraits<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType>::zero(),
440 const Scalar replacementValue = Teuchos::ScalarTraits<Scalar>::one()) {
441 using TST =
typename Teuchos::ScalarTraits<Scalar>;
442 using Teuchos::rcp_dynamic_cast;
444 GlobalOrdinal gZeroDiags;
445 bool usedEfficientPath =
false;
447 #ifdef HAVE_MUELU_TPETRA
448 RCP<CrsMatrixWrap> crsWrapAc = rcp_dynamic_cast<
CrsMatrixWrap>(Ac);
449 RCP<TpetraCrsMatrix> tpCrsAc;
450 if (!crsWrapAc.is_null())
451 tpCrsAc = rcp_dynamic_cast<TpetraCrsMatrix>(crsWrapAc->getCrsMatrix());
453 if (!tpCrsAc.is_null()) {
454 auto tpCrsGraph = tpCrsAc->getTpetra_CrsMatrix()->
getCrsGraph();
455 size_t numRows = Ac->getRowMap()->getLocalNumElements();
456 typedef typename Tpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>::offset_device_view_type offset_type;
457 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
458 auto offsets = offset_type(Kokkos::ViewAllocateWithoutInitializing(
"offsets"), numRows);
459 tpCrsGraph->getLocalDiagOffsets(offsets);
462 Tpetra::Details::OrdinalTraits<typename offset_type::value_type>::invalid();
464 if (repairZeroDiagonals) {
468 LO numMissingDiagonalEntries = 0;
470 Kokkos::parallel_reduce(
471 "countMissingDiagonalEntries",
472 range_type(0, numRows),
473 KOKKOS_LAMBDA(
const LO i, LO& missing) {
474 missing += (offsets(i) == STINV);
476 numMissingDiagonalEntries);
478 GlobalOrdinal gNumMissingDiagonalEntries;
479 Teuchos::reduceAll(*(Ac->getRowMap()->getComm()), Teuchos::REDUCE_SUM, Teuchos::as<GlobalOrdinal>(numMissingDiagonalEntries),
480 Teuchos::outArg(gNumMissingDiagonalEntries));
482 if (gNumMissingDiagonalEntries == 0) {
485 auto lclA = tpCrsAc->getTpetra_CrsMatrix()->getLocalMatrixDevice();
487 using ATS = Kokkos::ArithTraits<Scalar>;
488 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
491 typename ATS::val_type impl_replacementValue = replacementValue;
493 Kokkos::parallel_reduce(
494 "fixSmallDiagonalEntries",
495 range_type(0, numRows),
496 KOKKOS_LAMBDA(
const LO i, LO& fixed) {
497 const auto offset = offsets(i);
498 auto curRow = lclA.row(i);
499 if (impl_ATS::magnitude(curRow.value(offset)) <= threshold) {
500 curRow.value(offset) = impl_replacementValue;
506 Teuchos::reduceAll(*(Ac->getRowMap()->getComm()), Teuchos::REDUCE_SUM, Teuchos::as<GlobalOrdinal>(lZeroDiags),
507 Teuchos::outArg(gZeroDiags));
509 usedEfficientPath =
true;
515 auto lclA = tpCrsAc->getTpetra_CrsMatrix()->getLocalMatrixDevice();
517 using ATS = Kokkos::ArithTraits<Scalar>;
518 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
522 Kokkos::parallel_reduce(
523 "detectSmallDiagonalEntries",
524 range_type(0, numRows),
525 KOKKOS_LAMBDA(
const LO i, LO& small) {
526 const auto offset = offsets(i);
530 auto curRow = lclA.row(i);
531 if (impl_ATS::magnitude(curRow.value(offset)) <= threshold) {
538 Teuchos::reduceAll(*(Ac->getRowMap()->getComm()), Teuchos::REDUCE_SUM, Teuchos::as<GlobalOrdinal>(lZeroDiags),
539 Teuchos::outArg(gZeroDiags));
541 usedEfficientPath =
true;
546 if (!usedEfficientPath) {
547 RCP<const Map> rowMap = Ac->getRowMap();
549 Ac->getLocalDiagCopy(*diagVec);
551 LocalOrdinal lZeroDiags = 0;
552 Teuchos::ArrayRCP<const Scalar> diagVal = diagVec->getData(0);
554 for (
size_t i = 0; i < rowMap->getLocalNumElements(); i++) {
555 if (TST::magnitude(diagVal[i]) <= threshold) {
560 Teuchos::reduceAll(*(rowMap->getComm()), Teuchos::REDUCE_SUM, Teuchos::as<GlobalOrdinal>(lZeroDiags),
561 Teuchos::outArg(gZeroDiags));
563 if (repairZeroDiagonals && gZeroDiags > 0) {
579 Teuchos::Array<GlobalOrdinal> indout(1);
580 Teuchos::Array<Scalar> valout(1);
581 for (
size_t r = 0; r < rowMap->getLocalNumElements(); r++) {
582 if (TST::magnitude(diagVal[r]) <= threshold) {
583 GlobalOrdinal grid = rowMap->getGlobalElement(r);
585 valout[0] = replacementValue;
586 fixDiagMatrix->insertGlobalValues(grid, indout(), valout());
590 Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer(
"CheckRepairMainDiagonal: fillComplete1"));
591 fixDiagMatrix->fillComplete(Ac->getDomainMap(), Ac->getRangeMap());
595 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Ac,
false, 1.0, *fixDiagMatrix,
false, 1.0, newAc, fos);
596 if (Ac->IsView(
"stridedMaps"))
597 newAc->CreateView(
"stridedMaps", Ac);
600 fixDiagMatrix = Teuchos::null;
605 Teuchos::RCP<Teuchos::ParameterList> p = Teuchos::rcp(
new Teuchos::ParameterList());
606 p->set(
"DoOptimizeStorage",
true);
608 Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer(
"CheckRepairMainDiagonal: fillComplete2"));
615 fos <<
"CheckRepairMainDiagonal: " << (repairZeroDiagonals ?
"repaired " :
"found ")
616 << gZeroDiags <<
" too small entries (threshold = " << threshold <<
") on main diagonal of Ac." << std::endl;
618 #ifdef HAVE_XPETRA_DEBUG // only for debugging
621 RCP<const Map> rowMap = Ac->getRowMap();
623 Teuchos::ArrayRCP<const Scalar> diagVal;
624 Ac->getLocalDiagCopy(*diagVec);
625 diagVal = diagVec->getData(0);
626 for (
size_t r = 0; r < Ac->getRowMap()->getLocalNumElements(); r++) {
627 if (TST::magnitude(diagVal[r]) <= threshold) {
628 fos <<
"Error: there are too small entries left on diagonal after repair..." << std::endl;
643 const Teuchos::ArrayView<const double>& relativeThreshold, Teuchos::FancyOStream& fos) {
644 Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer(
"RelativeDiagonalBoost"));
646 TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() != relativeThreshold.size() && relativeThreshold.size() != 1,
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::RelativeDiagonal Boost: Either A->GetFixedBlockSize() != relativeThreshold.size() OR relativeThreshold.size() == 1");
648 LocalOrdinal numPDEs = A->GetFixedBlockSize();
649 typedef typename Teuchos::ScalarTraits<Scalar> TST;
650 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
651 Scalar zero = TST::zero();
652 Scalar one = TST::one();
656 A->getLocalDiagCopy(*diag);
657 Teuchos::ArrayRCP<const Scalar> dataVal = diag->getData(0);
658 size_t N = A->getRowMap()->getLocalNumElements();
661 std::vector<MT> l_diagMax(numPDEs), g_diagMax(numPDEs);
662 for (
size_t i = 0; i < N; i++) {
663 int pde = (int)(i % numPDEs);
664 if ((
int)i < numPDEs)
665 l_diagMax[pde] = TST::magnitude(dataVal[i]);
667 l_diagMax[pde] = std::max(l_diagMax[pde], TST::magnitude(dataVal[i]));
669 Teuchos::reduceAll(*A->getRowMap()->getComm(), Teuchos::REDUCE_MAX, numPDEs, l_diagMax.data(), g_diagMax.data());
673 Teuchos::Array<GlobalOrdinal> index(1);
674 Teuchos::Array<Scalar> value(1);
675 for (
size_t i = 0; i < N; i++) {
676 GlobalOrdinal GRID = A->getRowMap()->getGlobalElement(i);
677 int pde = (int)(i % numPDEs);
679 if (TST::magnitude(dataVal[i]) < relativeThreshold[pde] * g_diagMax[pde])
680 value[0] = relativeThreshold[pde] * g_diagMax[pde] - TST::magnitude(dataVal[i]);
683 boostMatrix->insertGlobalValues(GRID, index(), value());
685 boostMatrix->fillComplete(A->getDomainMap(), A->getRangeMap());
689 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*A,
false, one, *boostMatrix,
false, one, newA, fos);
690 if (A->IsView(
"stridedMaps"))
691 newA->CreateView(
"stridedMaps", A);
702 #if defined(HAVE_XPETRA_EPETRA)
704 #endif // HAVE_XPETRA_EPETRA
706 #ifdef HAVE_XPETRA_TPETRA
709 Tpetra::Details::extractBlockDiagonal(At, Dt);
710 #endif // HAVE_XPETRA_TPETRA
721 #if defined(HAVE_XPETRA_EPETRA)
723 #endif // HAVE_XPETRA_EPETRA
725 #ifdef HAVE_XPETRA_TPETRA
728 Tpetra::Details::inverseScaleBlockDiagonal(Dt, doTranspose, St);
729 #endif // HAVE_XPETRA_TPETRA
738 auto tpRowMap = Teuchos::rcp_dynamic_cast<
const TpetraMap>(rowMap,
true)->getTpetra_Map();
739 auto tpColMap = Teuchos::rcp_dynamic_cast<
const TpetraMap>(colMap,
true)->getTpetra_Map();
740 fail = !tpColMap->isLocallyFitted(*tpRowMap);
742 RCP<const Teuchos::Comm<int>> comm = rowMap->getComm();
743 LO numRows = Teuchos::as<LocalOrdinal>(rowMap->getLocalNumElements());
745 for (LO rowLID = 0; rowLID < numRows; rowLID++) {
746 GO rowGID = rowMap->getGlobalElement(rowLID);
747 LO colGID = colMap->getGlobalElement(rowLID);
748 if (rowGID != colGID) {
750 std::cerr <<
"On rank " << comm->getRank() <<
", LID " << rowLID <<
" is GID " << rowGID <<
" in the rowmap, but GID " << colGID <<
" in the column map.\n";
755 "Local parts of row and column map do not match!");
767 std::vector<size_t>& rangeStridingInfo, std::vector<size_t>& domainStridingInfo) {
771 if (matrix->IsView(
"stridedMaps") ==
true) matrix->RemoveView(
"stridedMaps");
772 matrix->CreateView(
"stridedMaps", stridedRowMap, stridedColMap);
778 #define XPETRA_MATRIXUTILS_SHORT
780 #endif // PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > findColumnSubMap(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &domainMap)
static void checkLocalRowMapMatchesColMap(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
static void inverseScaleBlockDiagonal(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &blockDiagonal, bool doTranspose, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &toBeScaled)
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
static Teuchos::RCP< Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > SplitMatrix(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &input, Teuchos::RCP< const Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node >> rangeMapExtractor, Teuchos::RCP< const Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node >> domainMapExtractor, Teuchos::RCP< const Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node >> columnMapExtractor=Teuchos::null, bool bThyraMode=false)
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> &map, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Exception throws to report errors in the internal logical of the program.
static Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > concatenateMaps(const std::vector< Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > > &subMaps)
Helper function to concatenate several maps.
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
virtual bool isNodeGlobalElement(GlobalOrdinal globalIndex) const =0
Whether the given global index is valid for this Map on this process.
static void convertMatrixToStridedMaps(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> matrix, std::vector< size_t > &rangeStridingInfo, std::vector< size_t > &domainStridingInfo)
virtual void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const =0
Extract a const, non-persisting view of local indices in a specified row of the matrix.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
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)
Map constructor with Xpetra-defined contiguous uniform distribution.
virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const =0
Const view of the local values in a particular vector of this multivector.
static void RelativeDiagonalBoost(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const Teuchos::ArrayView< const double > &relativeThreshold, Teuchos::FancyOStream &fos)
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
virtual const Teuchos::RCP< const map_type > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
virtual size_t getLocalLength() const =0
Local number of rows on the calling process.
RCP< const CrsGraph > getCrsGraph() const
Returns the CrsGraph associated with this matrix.
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > shrinkMapGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlInput)
Helper function to shrink the GIDs and generate a standard map whith GIDs starting at 0...
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
static Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > xpetraGidNumbering2ThyraGidNumbering(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &input)
std::vector< size_t > getStridingData() const
Exception throws to report incompatible objects (like maps).
static void extractBlockDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &diagonal)
Concrete implementation of Xpetra::Matrix.
virtual const Teuchos::RCP< const map_type > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
virtual size_t getLocalMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
static RCP< Xpetra::StridedMap< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
virtual size_t getNumVectors() const =0
Number of columns in the multivector.
static void CheckRepairMainDiagonal(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &Ac, bool const &repairZeroDiagonals, Teuchos::FancyOStream &fos, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType threshold=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< Scalar >::magnitudeType >::zero(), const Scalar replacementValue=Teuchos::ScalarTraits< Scalar >::one())
virtual UnderlyingLib lib() const =0
Get the library used by this object (Tpetra or Epetra?)
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)
Constructor specifying the number of non-zeros for all rows.
Xpetra-specific matrix class.
Class that stores a strided map.