47 #ifndef PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_ 48 #define PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_ 74 template <
class Scalar,
79 #undef XPETRA_MATRIXUTILS_SHORT 92 Teuchos::ArrayRCP< const Scalar > data = input.
getData(c);
94 ret->replaceLocalValue(Teuchos::as<LocalOrdinal>(r), c, data[r]);
105 RCP< const Teuchos::Comm<int> > comm = input.
getRowMap()->getComm();
108 Teuchos::Array<GlobalOrdinal> ovlUnknownStatusGids;
109 Teuchos::Array<GlobalOrdinal> ovlFoundStatusGids;
112 for(
size_t i = 0; i<input.
getColMap()->getNodeNumElements(); i++) {
113 GlobalOrdinal gcid = input.
getColMap()->getGlobalElement(i);
115 ovlUnknownStatusGids.push_back(gcid);
121 std::vector<int> myUnknownDofGIDs(comm->getSize(),0);
122 std::vector<int> numUnknownDofGIDs(comm->getSize(),0);
123 myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.size();
124 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myUnknownDofGIDs[0],&numUnknownDofGIDs[0]);
127 size_t cntUnknownDofGIDs = 0;
128 for(
int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
129 std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs,-1);
130 std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs,-1);
132 size_t cntUnknownOffset = 0;
133 for(
int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
134 for(
size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
135 lUnknownDofGIDs[k+cntUnknownOffset] = ovlUnknownStatusGids[k];
137 if(cntUnknownDofGIDs > 0)
138 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lUnknownDofGIDs[0],&gUnknownDofGIDs[0]);
139 std::sort(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end());
140 gUnknownDofGIDs.erase(std::unique(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end()), gUnknownDofGIDs.end());
143 for(
size_t k=0; k < gUnknownDofGIDs.size(); k++) {
144 GlobalOrdinal curgid = gUnknownDofGIDs[k];
146 ovlFoundStatusGids.push_back(curgid);
150 std::vector<int> myFoundDofGIDs(comm->getSize(),0);
151 std::vector<int> numFoundDofGIDs(comm->getSize(),0);
152 myFoundDofGIDs[comm->getRank()] = ovlFoundStatusGids.size();
153 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myFoundDofGIDs[0],&numFoundDofGIDs[0]);
156 size_t cntFoundDofGIDs = 0;
157 for(
int p = 0; p < comm->getSize(); p++) cntFoundDofGIDs += numFoundDofGIDs[p];
158 std::vector<GlobalOrdinal> lFoundDofGIDs(cntFoundDofGIDs,-1);
159 std::vector<GlobalOrdinal> gFoundDofGIDs(cntFoundDofGIDs,-1);
161 size_t cntFoundOffset = 0;
162 for(
int p = 0; p < comm->getRank(); p++) cntFoundOffset += numFoundDofGIDs[p];
163 for(
size_t k=0; k < Teuchos::as<size_t>(ovlFoundStatusGids.size()); k++) {
164 lFoundDofGIDs[k+cntFoundOffset] = ovlFoundStatusGids[k];
166 if(cntFoundDofGIDs > 0)
167 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntFoundDofGIDs),&lFoundDofGIDs[0],&gFoundDofGIDs[0]);
169 Teuchos::Array<GlobalOrdinal> ovlDomainMapArray;
170 for(
size_t i = 0; i<input.
getColMap()->getNodeNumElements(); i++) {
171 GlobalOrdinal gcid = input.
getColMap()->getGlobalElement(i);
173 std::find(gFoundDofGIDs.begin(), gFoundDofGIDs.end(), gcid) != gFoundDofGIDs.end()) {
174 ovlDomainMapArray.push_back(gcid);
177 RCP<Map> ovlDomainMap =
MapFactory::Build (domainMap.
lib(),Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),ovlDomainMapArray(),0,comm);
189 static Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
SplitMatrix(
194 bool bThyraMode =
false) {
200 size_t numRows = rangeMapExtractor->NumMaps();
201 size_t numCols = domainMapExtractor->NumMaps();
203 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.")
204 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.")
206 RCP<const Map> fullRangeMap = rangeMapExtractor->getFullMap();
207 RCP<const Map> fullDomainMap = domainMapExtractor->getFullMap();
209 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getMaxAllGlobalIndex() != input.
getRowMap()->getMaxAllGlobalIndex(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
210 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getGlobalNumElements() != input.
getRowMap()->getGlobalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
211 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getNodeNumElements() != input.
getRowMap()->getNodeNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
212 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getMaxAllGlobalIndex() != input.
getRangeMap()->getMaxAllGlobalIndex(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
213 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getGlobalNumElements() != input.
getRangeMap()->getGlobalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
214 TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getNodeNumElements() != input.
getRangeMap()->getNodeNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
216 TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getMaxAllGlobalIndex() != input.
getColMap()->getMaxAllGlobalIndex(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
217 TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getMaxAllGlobalIndex() != input.
getDomainMap()->getMaxAllGlobalIndex(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
218 TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getGlobalNumElements() != input.
getDomainMap()->getGlobalNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
219 TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getNodeNumElements() != input.
getDomainMap()->getNodeNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
222 Teuchos::RCP<const MapExtractor> myColumnMapExtractor = Teuchos::null;
223 if(columnMapExtractor == Teuchos::null) {
224 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.");
226 std::vector<Teuchos::RCP<const Map> > ovlxmaps(numCols, Teuchos::null);
227 for(
size_t c = 0; c < numCols; c++) {
230 ovlxmaps[c] = colMap;
236 myColumnMapExtractor = columnMapExtractor;
240 Teuchos::RCP<const MapExtractor> thyRangeMapExtractor = Teuchos::null;
241 Teuchos::RCP<const MapExtractor> thyDomainMapExtractor = Teuchos::null;
242 Teuchos::RCP<const MapExtractor> thyColMapExtractor = Teuchos::null;
243 if(bThyraMode ==
true) {
245 std::vector<Teuchos::RCP<const Map> > thyRgMapExtractorMaps(numRows, Teuchos::null);
246 for (
size_t r = 0; r < numRows; r++) {
247 RCP<const Map> rMap = rangeMapExtractor->getMap(r);
249 RCP<const StridedMap> strRangeMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rMap);
250 if(strRangeMap != Teuchos::null) {
252 GlobalOrdinal offset = strRangeMap->getOffset();
253 LocalOrdinal stridedBlockId = strRangeMap->getStridedBlockId();
254 RCP<const StridedMap> strShrinkedMap = Teuchos::rcp(
new StridedMap(shrinkedMap, strInfo, shrinkedMap->getIndexBase(), stridedBlockId, offset));
255 thyRgMapExtractorMaps[r] = strShrinkedMap;
257 thyRgMapExtractorMaps[r] = shrinkedMap;
259 TEUCHOS_TEST_FOR_EXCEPTION(thyRgMapExtractorMaps[r]->getNodeNumElements() != rMap->getNodeNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: Thyra-style range map extractor contains faulty data.")
263 std::vector<Teuchos::RCP<const Map> > thyDoMapExtractorMaps (numCols, Teuchos::null);
264 std::vector<Teuchos::RCP<const Map> > thyColMapExtractorMaps(numCols, Teuchos::null);
265 for (
size_t c = 0; c < numCols; c++) {
266 RCP<const Map> cMap = domainMapExtractor->getMap(c);
269 RCP<const StridedMap> strDomainMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(cMap);
270 if(strDomainMap != Teuchos::null) {
272 GlobalOrdinal offset = strDomainMap->getOffset();
273 LocalOrdinal stridedBlockId = strDomainMap->getStridedBlockId();
274 RCP<const StridedMap> strShrinkedDomainMap = Teuchos::rcp(
new StridedMap(shrinkedDomainMap, strInfo, shrinkedDomainMap->getIndexBase(), stridedBlockId, offset));
275 thyDoMapExtractorMaps[c] = strShrinkedDomainMap;
277 thyDoMapExtractorMaps[c] = shrinkedDomainMap;
279 RCP<const Map> colMap = myColumnMapExtractor->getMap(c);
281 RCP<const StridedMap> strColMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(colMap);
282 if(strColMap != Teuchos::null) {
284 GlobalOrdinal offset = strColMap->getOffset();
285 LocalOrdinal stridedBlockId = strColMap->getStridedBlockId();
286 RCP<const StridedMap> strShrinkedColMap = Teuchos::rcp(
new StridedMap(shrinkedColMap, strInfo, shrinkedColMap->getIndexBase(), stridedBlockId, offset));
287 thyColMapExtractorMaps[c] = strShrinkedColMap;
289 thyColMapExtractorMaps[c] = shrinkedColMap;
292 TEUCHOS_TEST_FOR_EXCEPTION(thyColMapExtractorMaps[c]->getNodeNumElements() != colMap->getNodeNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: Thyra-style column map extractor contains faulty data.")
293 TEUCHOS_TEST_FOR_EXCEPTION(thyDoMapExtractorMaps[c]->getNodeNumElements() != cMap->getNodeNumElements(),
Xpetra::Exceptions::Incompatible,
"Xpetra::MatrixUtils::Split: Thyra-style domain map extractor contains faulty data.")
301 std::vector<Teuchos::RCP<Matrix> > subMatrices(numRows*numCols, Teuchos::null);
302 for (
size_t r = 0; r < numRows; r++) {
303 for (
size_t c = 0; c < numCols; c++) {
307 if(bThyraMode ==
true)
320 #if 0 // TAW needs to be fixed (does not compile for Scalar=complex) 327 doCheck->putScalar(1.0);
331 TEUCHOS_TEST_FOR_EXCEPTION(coCheck->normInf() != Teuchos::ScalarTraits< Scalar >::magnitude(1.0),
Xpetra::Exceptions::RuntimeError,
"Xpetra::MatrixUtils::SplitMatrix: error when distributing data.");
333 doCheck->putScalar(-1.0);
334 coCheck->putScalar(-1.0);
336 Teuchos::ArrayRCP< Scalar > doCheckData = doCheck->getDataNonConst(0);
337 for (
size_t rrr = 0; rrr < input.
getDomainMap()->getNodeNumElements(); rrr++) {
339 GlobalOrdinal
id = input.
getDomainMap()->getGlobalElement(rrr);
342 size_t BlockId = domainMapExtractor->getMapIndexForGID(
id);
344 doCheckData[rrr] = Teuchos::as<Scalar>(BlockId);
349 Teuchos::ArrayRCP< Scalar > coCheckData = coCheck->getDataNonConst(0);
352 for (
size_t rr = 0; rr < input.
getRowMap()->getNodeNumElements(); 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()));
411 RCP<BlockedCrsMatrix> bA = Teuchos::null;
412 if(bThyraMode ==
true) {
413 for (
size_t r = 0; r < numRows; r++) {
414 for (
size_t c = 0; c < numCols; c++) {
415 subMatrices[r*numCols+c]->fillComplete(thyDomainMapExtractor->getMap(c,
true), thyRangeMapExtractor->getMap(r,
true));
418 bA = Teuchos::rcp(
new BlockedCrsMatrix(thyRangeMapExtractor, thyDomainMapExtractor, 10 ));
420 for (
size_t r = 0; r < numRows; r++) {
421 for (
size_t c = 0; c < numCols; c++) {
422 subMatrices[r*numCols+c]->fillComplete(domainMapExtractor->getMap(c), rangeMapExtractor->getMap(r));
425 bA = Teuchos::rcp(
new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10 ));
428 for (
size_t r = 0; r < numRows; r++) {
429 for (
size_t c = 0; c < numCols; c++) {
430 bA->setMatrix(r,c,subMatrices[r*numCols+c]);
439 #define XPETRA_MATRIXUTILS_SHORT 441 #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 RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
Constructor specifying the number of non-zeros for all rows.
std::vector< size_t > getStridingData() const
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, const Teuchos::RCP< Node > &node=defaultArgNode())
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.
Exception throws to report errors in the internal logical of the program.
Xpetra utility class for common map-related routines.
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)
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.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const=0
The Map describing the parallel distribution of this object.
virtual size_t getLocalLength() const =0
Local number of rows on the calling process.
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 Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > xpetraGidNumbering2ThyraGidNumbering(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &input)
Exception throws to report incompatible objects (like maps).
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
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 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 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...
Xpetra-specific matrix class.
Class that stores a strided map.