46 #ifndef XPETRA_REORDEREDBLOCKEDCRSMATRIX_HPP 47 #define XPETRA_REORDEREDBLOCKEDCRSMATRIX_HPP 68 template <
class Scalar,
72 class ReorderedBlockedCrsMatrix :
73 public BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
81 #undef XPETRA_REORDEREDBLOCKEDCRSMATRIX_SHORT 99 (Teuchos::RCP<const MapExtractor>& rangeMaps,
100 Teuchos::RCP<const MapExtractor>& domainMaps,
102 Teuchos::RCP<const Xpetra::BlockReorderManager> brm,
118 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
mergeSubBlockMaps(Teuchos::RCP<const Xpetra::BlockReorderManager> brm) {
119 RCP<const MapExtractor> fullRangeMapExtractor =
fullOp_->getRangeMapExtractor();
122 size_t numBlocks = brm->GetNumBlocks();
124 Teuchos::RCP<const Map> map = Teuchos::null;
131 map = fullRangeMapExtractor->getMap(Teuchos::as<size_t>(leaf->GetIndex()),
false);
134 std::vector<Teuchos::RCP<const Map> > subMaps (numBlocks, Teuchos::null);
136 for(
size_t i = 0; i < numBlocks; i++) {
137 Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
139 TEUCHOS_ASSERT(subMaps[i].is_null()==
false);
144 TEUCHOS_ASSERT(map.is_null()==
false);
157 Teuchos::ETransp mode = Teuchos::NO_TRANS,
158 Scalar alpha = ScalarTraits<Scalar>::one(),
159 Scalar beta = ScalarTraits<Scalar>::zero())
const 170 std::string
description()
const {
return "ReorderedBlockedCrsMatrix"; }
173 void describe(Teuchos::FancyOStream& out,
const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default)
const {
179 out <<
"ReorderedBlockMatrix is fillComplete" << std::endl;
181 out <<
"fullRowMap" << std::endl;
188 out <<
"Xpetra::ReorderedBlockedCrsMatrix is NOT fillComplete" << std::endl;
193 out <<
"Block(" << r <<
"," << c <<
")" << std::endl;
201 Teuchos::RCP<const Xpetra::BlockReorderManager >
brm_;
202 Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
fullOp_;
207 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
212 RCP<const Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > fullRangeMapExtractor = bmat->getRangeMapExtractor();
215 size_t numBlocks = brm->GetNumBlocks();
217 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> map = Teuchos::null;
223 map = fullRangeMapExtractor->getMap(Teuchos::as<size_t>(leaf->GetIndex()), bThyraMode);
226 std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subMaps (numBlocks, Teuchos::null);
228 for(
size_t i = 0; i < numBlocks; i++) {
229 Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
231 TEUCHOS_ASSERT(subMaps[i].is_null()==
false);
236 TEUCHOS_ASSERT(map.is_null()==
false);
240 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
251 size_t rowSz = rowMgr->GetNumBlocks();
252 size_t colSz = colMgr->GetNumBlocks();
254 Teuchos::RCP<BlockedCrsMatrix> rbmat = Teuchos::null;
256 if(rowSz == 0 && colSz == 0) {
258 Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<
const Xpetra::BlockReorderLeaf>(rowMgr);
259 Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<
const Xpetra::BlockReorderLeaf>(colMgr);
262 Teuchos::RCP<Matrix> mat = bmat->getMatrix(rowleaf->GetIndex(), colleaf->GetIndex());
264 if (mat == Teuchos::null)
return Teuchos::null;
268 if(matwrap != Teuchos::null) {
271 RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
272 Teuchos::RCP<const Map> submap = fullRangeMapExtractor->
getMap(rowleaf->GetIndex(),
false);
273 std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
274 Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::rcp(
new MapExtractor(submap, rowSubMaps,
false));
276 RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
277 Teuchos::RCP<const Map> submap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),
false);
278 std::vector<Teuchos::RCP<const Map> > colSubMaps (1, submap2);
279 Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::rcp(
new MapExtractor(submap2, colSubMaps,
false));
282 rbmat->setMatrix(0,0,mat);
286 TEUCHOS_ASSERT(rbmat != Teuchos::null);
288 TEUCHOS_ASSERT(mat->getNodeNumEntries() == rbmat->getNodeNumEntries());
292 Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::null;
294 std::vector<Teuchos::RCP<const Map> > rowSubMaps (rowSz, Teuchos::null);
295 for(
size_t i = 0; i < rowSz; i++) {
296 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
298 TEUCHOS_ASSERT(rowSubMaps[i].is_null()==
false);
301 rgMapExtractor = Teuchos::rcp(
new MapExtractor(rgMergedSubMaps, rowSubMaps,
false));
303 Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<
const Xpetra::BlockReorderLeaf>(rowMgr);
304 RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
307 Teuchos::RCP<const Map> submap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),
false);
308 std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
309 rgMapExtractor = Teuchos::rcp(
new MapExtractor(submap, rowSubMaps,
false));
312 Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::null;
314 std::vector<Teuchos::RCP<const Map> > colSubMaps (colSz, Teuchos::null);
315 for(
size_t j = 0; j < colSz; j++) {
316 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
318 TEUCHOS_ASSERT(colSubMaps[j].is_null()==
false);
321 doMapExtractor = Teuchos::rcp(
new MapExtractor(doMergedSubMaps, colSubMaps,
false));
323 Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<
const Xpetra::BlockReorderLeaf>(colMgr);
324 RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
327 Teuchos::RCP<const Map> submap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),
false);
328 std::vector<Teuchos::RCP<const Map> > colSubMaps (1, submap);
329 doMapExtractor = Teuchos::rcp(
new MapExtractor(submap, colSubMaps,
false));
336 if (rowSz == 0 && colSz > 0) {
337 for(
size_t j = 0; j < colSz; j++) {
338 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
339 Teuchos::RCP<const Matrix> submat =
mergeSubBlocks(rowMgr, colSubMgr, bmat);
340 rbmat->setMatrix(0,j,Teuchos::rcp_const_cast<Matrix>(submat));
341 if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
343 }
else if (rowSz > 0 && colSz == 0) {
344 for(
size_t i = 0; i < rowSz; i++) {
345 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
346 Teuchos::RCP<const Matrix> submat =
mergeSubBlocks(rowSubMgr, colMgr, bmat);
347 rbmat->setMatrix(i,0,Teuchos::rcp_const_cast<Matrix>(submat));
348 if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
351 for(
size_t i = 0; i < rowSz; i++) {
352 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
353 for(
size_t j = 0; j < colSz; j++) {
354 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
355 Teuchos::RCP<const Matrix> submat =
mergeSubBlocks(rowSubMgr, colSubMgr, bmat);
356 rbmat->setMatrix(i,j,Teuchos::rcp_const_cast<Matrix>(submat));
357 if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
361 TEUCHOS_ASSERT(rbmat->getNodeNumEntries() == cntNNZ);
363 rbmat->fillComplete();
369 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
378 TEUCHOS_ASSERT(bmat->getRangeMapExtractor()->getThyraMode() ==
true);
379 TEUCHOS_ASSERT(bmat->getDomainMapExtractor()->getThyraMode() ==
true);
382 size_t rowSz = rowMgr->GetNumBlocks();
383 size_t colSz = colMgr->GetNumBlocks();
385 Teuchos::RCP<BlockedCrsMatrix> rbmat = Teuchos::null;
387 if(rowSz == 0 && colSz == 0) {
389 Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<
const Xpetra::BlockReorderLeaf>(rowMgr);
390 Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<
const Xpetra::BlockReorderLeaf>(colMgr);
393 Teuchos::RCP<Matrix> mat = bmat->getMatrix(rowleaf->GetIndex(), colleaf->GetIndex());
395 if(mat == Teuchos::null)
return Teuchos::null;
399 if(matwrap != Teuchos::null) {
402 RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
404 Teuchos::RCP<const Map> xpsubmap = fullRangeMapExtractor->
getMap(rowleaf->GetIndex(),
false);
405 Teuchos::RCP<const Map> thysubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),
true);
406 std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
407 std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
409 Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::rcp(
new MapExtractor(rowXpSubMaps, rowTySubMaps));
411 RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
413 Teuchos::RCP<const Map> xpsubmap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),
false);
414 Teuchos::RCP<const Map> tysubmap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),
true);
415 std::vector<Teuchos::RCP<const Map> > colXpSubMaps (1, xpsubmap2);
416 std::vector<Teuchos::RCP<const Map> > colTySubMaps (1, tysubmap2);
418 Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::rcp(
new MapExtractor(colXpSubMaps, colTySubMaps));
423 rbmat->setMatrix(0,0,mat);
427 TEUCHOS_ASSERT(rbmat != Teuchos::null);
429 TEUCHOS_ASSERT(mat->getNodeNumEntries() == rbmat->getNodeNumEntries());
433 Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::null;
435 std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (rowSz, Teuchos::null);
436 std::vector<Teuchos::RCP<const Map> > rowTySubMaps (rowSz, Teuchos::null);
437 for(
size_t i = 0; i < rowSz; i++) {
438 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
442 TEUCHOS_ASSERT(rowXpSubMaps[i].is_null()==
false);
443 TEUCHOS_ASSERT(rowTySubMaps[i].is_null()==
false);
446 rgMapExtractor = Teuchos::rcp(
new MapExtractor(rowXpSubMaps, rowTySubMaps));
448 Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<
const Xpetra::BlockReorderLeaf>(rowMgr);
449 RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
451 Teuchos::RCP<const Map> xpsubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),
false);
452 Teuchos::RCP<const Map> thysubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),
true);
453 std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
454 std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
456 rgMapExtractor = Teuchos::rcp(
new MapExtractor(rowXpSubMaps, rowTySubMaps));
459 Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::null;
461 std::vector<Teuchos::RCP<const Map> > colXpSubMaps (colSz, Teuchos::null);
462 std::vector<Teuchos::RCP<const Map> > colTySubMaps (colSz, Teuchos::null);
463 for(
size_t j = 0; j < colSz; j++) {
464 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
468 TEUCHOS_ASSERT(colXpSubMaps[j].is_null()==
false);
469 TEUCHOS_ASSERT(colTySubMaps[j].is_null()==
false);
472 doMapExtractor = Teuchos::rcp(
new MapExtractor(colXpSubMaps,colTySubMaps));
474 Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<
const Xpetra::BlockReorderLeaf>(colMgr);
475 RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
477 Teuchos::RCP<const Map> xpsubmap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),
false);
478 Teuchos::RCP<const Map> tysubmap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),
true);
479 std::vector<Teuchos::RCP<const Map> > colXpSubMaps (1, xpsubmap);
480 std::vector<Teuchos::RCP<const Map> > colTySubMaps (1, tysubmap);
482 doMapExtractor = Teuchos::rcp(
new MapExtractor(colXpSubMaps, colTySubMaps));
490 if (rowSz == 0 && colSz > 0) {
491 for(
size_t j = 0; j < colSz; j++) {
492 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
494 rbmat->setMatrix(0,j,Teuchos::rcp_const_cast<Matrix>(submat));
495 if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
497 }
else if (rowSz > 0 && colSz == 0) {
498 for(
size_t i = 0; i < rowSz; i++) {
499 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
501 rbmat->setMatrix(i,0,Teuchos::rcp_const_cast<Matrix>(submat));
502 if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
505 for(
size_t i = 0; i < rowSz; i++) {
506 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
507 for(
size_t j = 0; j < colSz; j++) {
508 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
510 rbmat->setMatrix(i,j,Teuchos::rcp_const_cast<Matrix>(submat));
511 if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
515 TEUCHOS_ASSERT(rbmat->getNodeNumEntries() == cntNNZ);
518 rbmat->fillComplete();
522 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
524 TEUCHOS_ASSERT(bmat->getRangeMapExtractor()->getThyraMode() == bmat->getDomainMapExtractor()->getThyraMode());
525 Teuchos::RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbmat = Teuchos::null;
526 if(bmat->getRangeMapExtractor()->getThyraMode() ==
false) {
538 #define XPETRA_REORDEREDBLOCKEDCRSMATRIX_SHORT
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlockMaps(Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat, bool bThyraMode)
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > fullOp_
const Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const
Implements DistObject interface.
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.
GlobalOrdinal global_ordinal_type
Teuchos::RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlocksThyra(Teuchos::RCP< const Xpetra::BlockReorderManager > rowMgr, Teuchos::RCP< const Xpetra::BlockReorderManager > colMgr, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat)
Teuchos::RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlocks(Teuchos::RCP< const Xpetra::BlockReorderManager > rowMgr, Teuchos::RCP< const Xpetra::BlockReorderManager > colMgr, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat)
std::string description() const
Return a simple one-line description of this object.
Teuchos::RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > buildReorderedBlockedCrsMatrix(Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat)
virtual ~ReorderedBlockedCrsMatrix()
Destructor.
Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlockMaps(Teuchos::RCP< const Xpetra::BlockReorderManager > brm)
LocalOrdinal local_ordinal_type
BlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMaps, Teuchos::RCP< const MapExtractor > &domainMaps, size_t npr, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor.
RCP< const Map > getRangeMap() const
Returns the Map associated with the full range of this operator.
Concrete implementation of Xpetra::Matrix.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
ReorderedBlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMaps, Teuchos::RCP< const MapExtractor > &domainMaps, size_t npr, Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor.
virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=ScalarTraits< Scalar >::one(), Scalar beta=ScalarTraits< Scalar >::zero()) const
Computes the sparse matrix-multivector multiplication.
virtual size_t Rows() const
number of row blocks
virtual size_t Cols() const
number of column blocks
Xpetra-specific matrix class.
Teuchos::RCP< const Xpetra::BlockReorderManager > brm_
virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=ScalarTraits< Scalar >::one(), Scalar beta=ScalarTraits< Scalar >::zero()) const
Computes the sparse matrix-multivector multiplication.