47 #ifndef PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
48 #define PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
53 #include "Xpetra_Map.hpp"
55 #include "Xpetra_StridedMap.hpp"
56 #include "Xpetra_MapFactory.hpp"
57 #include "Xpetra_MapExtractor.hpp"
66 #ifdef HAVE_XPETRA_TPETRA
67 #include "Xpetra_TpetraMultiVector.hpp"
68 #include <Tpetra_RowMatrixTransposer.hpp>
69 #include <Tpetra_Details_extractBlockDiagonal.hpp>
70 #include <Tpetra_Details_scaleBlockDiagonal.hpp>
84 template <
class Scalar,
89 #undef XPETRA_MATRIXUTILS_SHORT
104 ret->replaceLocalValue(Teuchos::as<LocalOrdinal>(r), c, data[r]);
122 for(
size_t i = 0; i<input.
getColMap()->getNodeNumElements(); i++) {
123 GlobalOrdinal gcid = input.
getColMap()->getGlobalElement(i);
131 std::vector<int> myUnknownDofGIDs(comm->getSize(),0);
132 std::vector<int> numUnknownDofGIDs(comm->getSize(),0);
133 myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.
size();
137 size_t cntUnknownDofGIDs = 0;
138 for(
int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
139 std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs,-1);
140 std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs,-1);
142 size_t cntUnknownOffset = 0;
143 for(
int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
144 for(
size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.
size()); k++) {
145 lUnknownDofGIDs[k+cntUnknownOffset] = ovlUnknownStatusGids[k];
147 if(cntUnknownDofGIDs > 0)
149 std::sort(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end());
150 gUnknownDofGIDs.erase(std::unique(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end()), gUnknownDofGIDs.end());
153 for(
size_t k=0; k < gUnknownDofGIDs.size(); k++) {
154 GlobalOrdinal curgid = gUnknownDofGIDs[k];
160 std::vector<int> myFoundDofGIDs(comm->getSize(),0);
161 std::vector<int> numFoundDofGIDs(comm->getSize(),0);
162 myFoundDofGIDs[comm->getRank()] = ovlFoundStatusGids.
size();
166 size_t cntFoundDofGIDs = 0;
167 for(
int p = 0; p < comm->getSize(); p++) cntFoundDofGIDs += numFoundDofGIDs[p];
168 std::vector<GlobalOrdinal> lFoundDofGIDs(cntFoundDofGIDs,-1);
169 std::vector<GlobalOrdinal> gFoundDofGIDs(cntFoundDofGIDs,-1);
171 size_t cntFoundOffset = 0;
172 for(
int p = 0; p < comm->getRank(); p++) cntFoundOffset += numFoundDofGIDs[p];
173 for(
size_t k=0; k < Teuchos::as<size_t>(ovlFoundStatusGids.
size()); k++) {
174 lFoundDofGIDs[k+cntFoundOffset] = ovlFoundStatusGids[k];
176 if(cntFoundDofGIDs > 0)
180 for(
size_t i = 0; i<input.
getColMap()->getNodeNumElements(); i++) {
181 GlobalOrdinal gcid = input.
getColMap()->getGlobalElement(i);
183 std::find(gFoundDofGIDs.begin(), gFoundDofGIDs.end(), gcid) != gFoundDofGIDs.end()) {
204 bool bThyraMode =
false) {
210 size_t numRows = rangeMapExtractor->NumMaps();
211 size_t numCols = domainMapExtractor->NumMaps();
213 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.")
214 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.")
226 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())
233 if(columnMapExtractor == Teuchos::null) {
236 std::vector<Teuchos::RCP<const Map> > ovlxmaps(numCols, Teuchos::null);
237 for(
size_t c = 0; c < numCols; c++) {
240 ovlxmaps[c] = colMap;
246 myColumnMapExtractor = columnMapExtractor;
253 if(bThyraMode ==
true) {
255 std::vector<Teuchos::RCP<const Map> > thyRgMapExtractorMaps(numRows, Teuchos::null);
256 for (
size_t r = 0; r < numRows; r++) {
260 if(strRangeMap != Teuchos::null) {
261 std::vector<size_t> strInfo = strRangeMap->getStridingData();
262 GlobalOrdinal offset = strRangeMap->getOffset();
263 LocalOrdinal stridedBlockId = strRangeMap->getStridedBlockId();
265 thyRgMapExtractorMaps[r] = strShrinkedMap;
267 thyRgMapExtractorMaps[r] = shrinkedMap;
273 std::vector<Teuchos::RCP<const Map> > thyDoMapExtractorMaps (numCols, Teuchos::null);
274 std::vector<Teuchos::RCP<const Map> > thyColMapExtractorMaps(numCols, Teuchos::null);
275 for (
size_t c = 0; c < numCols; c++) {
280 if(strDomainMap != Teuchos::null) {
281 std::vector<size_t> strInfo = strDomainMap->getStridingData();
282 GlobalOrdinal offset = strDomainMap->getOffset();
283 LocalOrdinal stridedBlockId = strDomainMap->getStridedBlockId();
285 thyDoMapExtractorMaps[c] = strShrinkedDomainMap;
287 thyDoMapExtractorMaps[c] = shrinkedDomainMap;
292 if(strColMap != Teuchos::null) {
293 std::vector<size_t> strInfo = strColMap->getStridingData();
294 GlobalOrdinal offset = strColMap->getOffset();
295 LocalOrdinal stridedBlockId = strColMap->getStridedBlockId();
297 thyColMapExtractorMaps[c] = strShrinkedColMap;
299 thyColMapExtractorMaps[c] = shrinkedColMap;
311 std::vector<Teuchos::RCP<Matrix> > subMatrices(numRows*numCols, Teuchos::null);
312 for (
size_t r = 0; r < numRows; r++) {
313 for (
size_t c = 0; c < numCols; c++) {
317 if(bThyraMode ==
true)
330 #if 0 // TAW needs to be fixed (does not compile for Scalar=complex)
337 doCheck->putScalar(1.0);
343 doCheck->putScalar(-1.0);
344 coCheck->putScalar(-1.0);
347 for (
size_t rrr = 0; rrr < input.
getDomainMap()->getNodeNumElements(); rrr++) {
349 GlobalOrdinal
id = input.
getDomainMap()->getGlobalElement(rrr);
352 size_t BlockId = domainMapExtractor->getMapIndexForGID(
id);
354 doCheckData[rrr] = Teuchos::as<Scalar>(BlockId);
362 for (
size_t rr = 0; rr < input.
getRowMap()->getNodeNumElements(); rr++) {
365 GlobalOrdinal growid = input.
getRowMap()->getGlobalElement(rr);
374 size_t rowBlockId = rangeMapExtractor->getMapIndexForGID(growid);
377 GlobalOrdinal subblock_growid = growid;
378 if(bThyraMode ==
true) {
380 LocalOrdinal lrowid = rangeMapExtractor->getMap(rowBlockId)->getLocalElement(growid);
382 subblock_growid = thyRangeMapExtractor->getMap(rowBlockId,
true)->getGlobalElement(lrowid);
395 for(
size_t i=0; i<(size_t)indices.
size(); i++) {
397 GlobalOrdinal gcolid = input.
getColMap()->getGlobalElement(indices[i]);
399 size_t colBlockId = myColumnMapExtractor->getMapIndexForGID(gcolid);
403 GlobalOrdinal subblock_gcolid = gcolid;
404 if(bThyraMode ==
true) {
406 LocalOrdinal lcolid = myColumnMapExtractor->getMap(colBlockId)->getLocalElement(gcolid);
408 subblock_gcolid = thyColMapExtractor->getMap(colBlockId,
true)->getGlobalElement(lcolid);
410 blockColIdx [colBlockId].push_back(subblock_gcolid);
411 blockColVals[colBlockId].push_back(vals[i]);
414 for (
size_t c = 0; c < numCols; c++) {
415 subMatrices[rowBlockId*numCols+c]->insertGlobalValues(subblock_growid,blockColIdx[c].view(0,blockColIdx[c].size()),blockColVals[c].view(0,blockColVals[c].size()));
422 if(bThyraMode ==
true) {
423 for (
size_t r = 0; r < numRows; r++) {
424 for (
size_t c = 0; c < numCols; c++) {
425 subMatrices[r*numCols+c]->fillComplete(thyDomainMapExtractor->getMap(c,
true), thyRangeMapExtractor->getMap(r,
true));
430 for (
size_t r = 0; r < numRows; r++) {
431 for (
size_t c = 0; c < numCols; c++) {
432 subMatrices[r*numCols+c]->fillComplete(domainMapExtractor->getMap(c), rangeMapExtractor->getMap(r));
438 for (
size_t r = 0; r < numRows; r++) {
439 for (
size_t c = 0; c < numCols; c++) {
440 bA->setMatrix(r,c,subMatrices[r*numCols+c]);
453 Scalar one = TST::one();
456 p->
set(
"DoOptimizeStorage",
true);
460 Ac->getLocalDiagCopy(*diagVec);
462 LocalOrdinal lZeroDiags = 0;
465 for (
size_t i = 0; i < rowMap->getNodeNumElements(); i++) {
466 if (TST::magnitude(diagVal[i]) <= threshold) {
470 GlobalOrdinal gZeroDiags;
472 Teuchos::outArg(gZeroDiags));
474 if (repairZeroDiagonals && gZeroDiags > 0) {
492 for (
size_t r = 0; r < rowMap->getNodeNumElements(); r++) {
493 if (TST::magnitude(diagVal[r]) <= threshold) {
494 GlobalOrdinal grid = rowMap->getGlobalElement(r);
497 fixDiagMatrix->insertGlobalValues(grid,indout(), valout());
502 fixDiagMatrix->fillComplete(Ac->getDomainMap(),Ac->getRangeMap());
506 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*Ac,
false, 1.0, *fixDiagMatrix,
false, 1.0, newAc, fos);
507 if (Ac->IsView(
"stridedMaps"))
508 newAc->CreateView(
"stridedMaps", Ac);
511 fixDiagMatrix = Teuchos::null;
525 fos <<
"CheckRepairMainDiagonal: " << (repairZeroDiagonals ?
"repaired " :
"found ")
526 << gZeroDiags <<
" too small entries on main diagonal of Ac." << std::endl;
528 #ifdef HAVE_XPETRA_DEBUG // only for debugging
530 Ac->getLocalDiagCopy(*diagVec);
532 for (
size_t r = 0; r < Ac->getRowMap()->getNodeNumElements(); r++) {
533 if (TST::magnitude(diagVal[r]) <= threshold) {
534 fos <<
"Error: there are too small entries left on diagonal after repair..." << std::endl;
556 LocalOrdinal numPDEs = A->GetFixedBlockSize();
559 Scalar zero = TST::zero();
560 Scalar one = TST::one();
564 A->getLocalDiagCopy(*diag);
566 size_t N = A->getRowMap()->getNodeNumElements();
569 std::vector<MT> l_diagMax(numPDEs), g_diagMax(numPDEs);
570 for(
size_t i=0; i<N; i++) {
571 int pde = (int) (i % numPDEs);
573 l_diagMax[pde] = TST::magnitude(dataVal[i]);
575 l_diagMax[pde] = std::max(l_diagMax[pde],TST::magnitude(dataVal[i]));
583 for (
size_t i = 0; i<N; i++) {
584 GlobalOrdinal GRID = A->getRowMap()->getGlobalElement(i);
585 int pde = (int) (i % numPDEs);
587 if (TST::magnitude(dataVal[i]) < relativeThreshold[pde] * g_diagMax[pde])
588 value[0] = relativeThreshold[pde] * g_diagMax[pde] - TST::magnitude(dataVal[i]);
591 boostMatrix->insertGlobalValues(GRID,index(),value());
593 boostMatrix->fillComplete(A->getDomainMap(),A->getRangeMap());
597 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*A,
false,one, *boostMatrix,
false,one,newA,fos);
598 if (A->IsView(
"stridedMaps"))
599 newA->CreateView(
"stridedMaps", A);
611 #if defined(HAVE_XPETRA_EPETRA)
613 #endif // HAVE_XPETRA_EPETRA
615 #ifdef HAVE_XPETRA_TPETRA
618 Tpetra::Details::extractBlockDiagonal(At,Dt);
619 #endif // HAVE_XPETRA_TPETRA
631 #if defined(HAVE_XPETRA_EPETRA)
633 #endif // HAVE_XPETRA_EPETRA
635 #ifdef HAVE_XPETRA_TPETRA
636 const Tpetra::MultiVector<SC,LO,GO,NO> & Dt =
Xpetra::toTpetra(blockDiagonal);
638 Tpetra::Details::inverseScaleBlockDiagonal(Dt,doTranspose,St);
639 #endif // HAVE_XPETRA_TPETRA
648 #define XPETRA_MATRIXUTILS_SHORT
650 #endif // PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
virtual size_t getNodeMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
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 RelativeDiagonalBoost(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayView< const double > &relativeThreshold, Teuchos::FancyOStream &fos)
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.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Exception throws to report errors in the internal logical of the program.
Xpetra utility class for common map-related routines.
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())
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.
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 RCP< Time > getNewTimer(const std::string &name)
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.
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
virtual size_t getLocalLength() const =0
Local number of rows on the calling process.
static void inverseScaleBlockDiagonal(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &blockDiagonal, bool doTranspose, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &toBeScaled)
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 RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
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)
static magnitudeType magnitude(T a)
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
void push_back(const value_type &x)
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)
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor specifying the number of non-zeros for all rows.
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.
virtual UnderlyingLib lib() const =0
Get the library used by this object (Tpetra or Epetra?)
virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const =0
Const view of the local values in a particular vector of this multivector.
Xpetra utility class for common matrix-related routines.
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
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.
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)