Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_ReorderedBlockedCrsMatrix.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Xpetra: A linear algebra interface package
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef XPETRA_REORDEREDBLOCKEDCRSMATRIX_HPP
47#define XPETRA_REORDEREDBLOCKEDCRSMATRIX_HPP
48
49#include <KokkosCompat_DefaultNode.hpp>
50
51#include "Xpetra_ConfigDefs.hpp"
52#include "Xpetra_Exceptions.hpp"
53
54#include "Xpetra_MapUtils.hpp"
55
56#include "Xpetra_BlockedMultiVector.hpp"
58#include "Xpetra_CrsMatrixWrap.hpp"
60
61
66namespace Xpetra {
67
68 typedef std::string viewLabel_t;
69
70 template <class Scalar,
71 class LocalOrdinal,
72 class GlobalOrdinal,
75 public BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
76 public:
77 typedef Scalar scalar_type;
78 typedef LocalOrdinal local_ordinal_type;
79 typedef GlobalOrdinal global_ordinal_type;
80 typedef Node node_type;
81
82 private:
83#undef XPETRA_REORDEREDBLOCKEDCRSMATRIX_SHORT
85
86 public:
87
89
90
92
102 size_t npr,
105 : Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rangeMaps, domainMaps, npr) {
106 brm_ = brm;
107 fullOp_ = bmat;
108 }
109
110 //protected:
111
114
116
117 private:
119 RCP<const MapExtractor> fullRangeMapExtractor = fullOp_->getRangeMapExtractor();
120
121 // number of sub blocks
122 size_t numBlocks = brm->GetNumBlocks();
123
124 Teuchos::RCP<const Map> map = Teuchos::null;
125
126 if(numBlocks == 0) {
127 // it is a leaf node
128 Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);
129
130 // never extract Thyra style maps (since we have to merge them)
131 map = fullRangeMapExtractor->getMap(Teuchos::as<size_t>(leaf->GetIndex()), false);
132 } else {
133 // initialize vector for sub maps
134 std::vector<Teuchos::RCP<const Map> > subMaps (numBlocks, Teuchos::null);
135
136 for(size_t i = 0; i < numBlocks; i++) {
137 Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
138 subMaps[i] = mergeSubBlockMaps(blkMgr);
139 TEUCHOS_ASSERT(subMaps[i].is_null()==false);
140 }
141
142 map = MapUtils::concatenateMaps(subMaps);
143 }
144 TEUCHOS_ASSERT(map.is_null()==false);
145 return map;
146 }
147
148 public:
150
151
153 virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode, Scalar alpha, Scalar beta, bool sumInterfaceValues,
154 const RCP<Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> >& regionInterfaceImporter,
155 const Teuchos::ArrayRCP<LocalOrdinal>& regionInterfaceLIDs) const
156 { }
158
161 virtual void apply(const MultiVector& X, MultiVector& Y,
163 Scalar alpha = ScalarTraits<Scalar>::one(),
164 Scalar beta = ScalarTraits<Scalar>::zero()) const
165 {
166
167 // Nested sub-operators should just use the provided X and B vectors
168 if(fullOp_->getLocalNumRows() != this->getLocalNumRows()) {
170 return;
171 }
172
173 // Special handling for the top level of the nested operator
174
175 // check whether input parameters are blocked or not
176 RCP<const MultiVector> refX = rcpFromRef(X);
177 RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
178 RCP<MultiVector> tmpY = rcpFromRef(Y);
179 RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
180
181
182 bool bCopyResultX = false;
183 bool bCopyResultY = false;
184
185 // TODO create a nested ReorderedBlockedMultiVector out of a nested BlockedMultiVector!
186 // probably not necessary, is it?
187 // check whether X and B are blocked but not ReorderedBlocked operators
188 /*if (refbX != Teuchos::null && fullOp_->getLocalNumRows() == this->getLocalNumRows()) {
189 RCP<const ReorderedBlockedMultiVector> rbCheck = Teuchos::rcp_dynamic_cast<const ReorderedBlockedMultiVector>(refbX);
190 if(rbCheck == Teuchos::null) {
191 RCP<const BlockedMultiVector> bX =
192 Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(Xpetra::buildReorderedBlockedMultiVector(brm_, refbX));
193 TEUCHOS_ASSERT(bX.is_null()==false);
194 refbX.swap(bX);
195 }
196 }
197 if (tmpbY != Teuchos::null && fullOp_->getLocalNumRows() == this->getLocalNumRows()) {
198 RCP<ReorderedBlockedMultiVector> rbCheck = Teuchos::rcp_dynamic_cast<ReorderedBlockedMultiVector>(tmpbY);
199 if(rbCheck == Teuchos::null) {
200 RCP<BlockedMultiVector> bY =
201 Teuchos::rcp_dynamic_cast<BlockedMultiVector>(Xpetra::buildReorderedBlockedMultiVector(brm_, tmpbY));
202 TEUCHOS_ASSERT(bY.is_null()==false);
203 tmpbY.swap(bY);
204 }
205 }*/
206
207 // if X and B are not blocked, create a blocked version here and use the blocked vectors
208 // for the internal (nested) apply call.
209
210 // Check whether "this" operator is the reordered variant of the underlying fullOp_.
211 // Note, that nested ReorderedBlockedCrsMatrices always have the same full operator "fullOp_"
212 // stored underneath for being able to "translate" the block ids.
213 if (refbX == Teuchos::null && fullOp_->getLocalNumRows() == this->getLocalNumRows())
214 {
215 // create a new (non-nested) blocked multi vector (using the blocked range map of fullOp_)
216 RCP<const BlockedMap> blkRgMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(fullOp_->getRangeMap());
217 TEUCHOS_ASSERT(blkRgMap.is_null()==false);
219 TEUCHOS_ASSERT(bXtemp.is_null()==false);
221 Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(Xpetra::buildReorderedBlockedMultiVector(brm_, bXtemp));
222 TEUCHOS_ASSERT(bX.is_null()==false);
223 refbX.swap(bX);
224 bCopyResultX = true;
225 }
226
227 if (tmpbY == Teuchos::null && fullOp_->getLocalNumRows() == this->getLocalNumRows()) {
228 // create a new (non-nested) blocked multi vector (using the blocked range map of fullOp_)
229 RCP<const BlockedMap> blkRgMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(fullOp_->getRangeMap());
230 TEUCHOS_ASSERT(blkRgMap.is_null()==false);
231 RCP<BlockedMultiVector> tmpbYtemp = Teuchos::rcp(new BlockedMultiVector(blkRgMap, tmpY));
232 TEUCHOS_ASSERT(tmpbYtemp.is_null()==false);
234 Teuchos::rcp_dynamic_cast<BlockedMultiVector>(Xpetra::buildReorderedBlockedMultiVector(brm_, tmpbYtemp));
235 TEUCHOS_ASSERT(bY.is_null()==false);
236 tmpbY.swap(bY);
237 bCopyResultY = true;
238 }
239
240 TEUCHOS_ASSERT(refbX.is_null()==false);
241 TEUCHOS_ASSERT(tmpbY.is_null()==false);
242
244
245 if (bCopyResultX == true) {
246 RCP<const MultiVector> Xmerged = refbX->Merge();
247 RCP<MultiVector> nonconstX = Teuchos::rcp_const_cast<MultiVector>(refX);
249 }
250 if (bCopyResultY == true) {
251 RCP< MultiVector> Ymerged = tmpbY->Merge();
253 }
254
255 }
256
257 // @}
258
259
260
262
263
266
267 // @}
268
270
271
273 std::string description() const { return "ReorderedBlockedCrsMatrix"; }
274
277
278 out << "Xpetra::ReorderedBlockedCrsMatrix: " << BlockedCrsMatrix::Rows() << " x " << BlockedCrsMatrix::Cols() << std::endl;
279
281 out << "ReorderedBlockMatrix is fillComplete" << std::endl;
282
283 out << "fullRowMap" << std::endl;
284 BlockedCrsMatrix::getRangeMap(0,false)->describe(out,verbLevel);
285
286 //out << "fullColMap" << std::endl;
287 //fullcolmap_->describe(out,verbLevel);
288
289 } else {
290 out << "Xpetra::ReorderedBlockedCrsMatrix is NOT fillComplete" << std::endl;
291 }
292
293 for (size_t r = 0; r < BlockedCrsMatrix::Rows(); ++r)
294 for (size_t c = 0; c < BlockedCrsMatrix::Cols(); ++c) {
295 out << "Block(" << r << "," << c << ")" << std::endl;
296 BlockedCrsMatrix::getMatrix(r,c)->describe(out,verbLevel);
297 }
298 }
299
301
302 private:
305
306
307};
308
309template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
312
313 // TODO distinguish between range and domain map extractor! provide MapExtractor as parameter!
314 RCP<const Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > fullRangeMapExtractor = bmat->getRangeMapExtractor();
315
316 // number of sub blocks
317 size_t numBlocks = brm->GetNumBlocks();
318
320
321 if(numBlocks == 0) {
322 // it is a leaf node
323 Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);
324
325 map = fullRangeMapExtractor->getMap(Teuchos::as<size_t>(leaf->GetIndex()), bThyraMode);
326 } else {
327 // initialize vector for sub maps
328 std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subMaps (numBlocks, Teuchos::null);
329
330 for(size_t i = 0; i < numBlocks; i++) {
331 Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
332 subMaps[i] = mergeSubBlockMaps(blkMgr,bmat,bThyraMode);
333 TEUCHOS_ASSERT(subMaps[i].is_null()==false);
334 }
335
336#if 1
337 // concatenate submaps
338 // for Thyra mode this map isn't important
340
341 // create new BlockedMap (either in Thyra Mode or Xpetra mode)
342 map = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(fullMap, subMaps, bThyraMode));
343#else
344 // TAW: 11/27/16 we just concatenate the submaps to one monolithic Map object.
345 // Alternatively, we could create a new BlockedMap using the concatenated map and the submaps
346 // However, the block smoothers only need the concatenated map for creating MultiVectors...
347 // But for the Thyra mode version concatenating would not be ok for the whole map
348 map = MapUtils::concatenateMaps(subMaps);
349#endif
350 }
351 TEUCHOS_ASSERT(map.is_null()==false);
352 return map;
353}
354
355template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
357
364
365 // number of sub blocks
366 size_t rowSz = rowMgr->GetNumBlocks();
367 size_t colSz = colMgr->GetNumBlocks();
368
369 Teuchos::RCP<BlockedCrsMatrix> rbmat = Teuchos::null;
370
371 if(rowSz == 0 && colSz == 0) {
372 // it is a leaf node
373 Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
374 Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
375
376 // extract leaf node
377 Teuchos::RCP<Matrix> mat = bmat->getMatrix(rowleaf->GetIndex(), colleaf->GetIndex());
378
379 if (mat == Teuchos::null) return Teuchos::null;
380
381 // check, whether leaf node is of type Xpetra::CrsMatrixWrap
382 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> > matwrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
383 if(matwrap != Teuchos::null) {
384 // If the leaf node is of type Xpetra::CrsMatrixWrap wrap it into a 1x1 ReorderedBlockMatrix
385 // with the corresponding MapExtractors for translating Thyra to Xpetra GIDs if necessary
386 RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
387 Teuchos::RCP<const Map> submap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
388 std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
389 Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::rcp(new MapExtractor(submap, rowSubMaps, false));
390
391 RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
392 Teuchos::RCP<const Map> submap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
393 std::vector<Teuchos::RCP<const Map> > colSubMaps (1, submap2);
394 Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::rcp(new MapExtractor(submap2, colSubMaps, false));
395
396 rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
397 rbmat->setMatrix(0,0,mat);
398 } else {
399 // If leaf node is already wrapped into a blocked matrix do not wrap it again.
400 rbmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
401 TEUCHOS_ASSERT(rbmat != Teuchos::null);
402 }
403 TEUCHOS_ASSERT(mat->getLocalNumEntries() == rbmat->getLocalNumEntries());
404 } else {
405 // create the map extractors
406 // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
407 Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::null;
408 if(rowSz > 0) {
409 std::vector<Teuchos::RCP<const Map> > rowSubMaps (rowSz, Teuchos::null);
410 for(size_t i = 0; i < rowSz; i++) {
411 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
412 rowSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,false /*xpetra*/);
413 TEUCHOS_ASSERT(rowSubMaps[i].is_null()==false);
414 }
415 Teuchos::RCP<const Map> rgMergedSubMaps = MapUtils::concatenateMaps(rowSubMaps);
416 rgMapExtractor = Teuchos::rcp(new MapExtractor(rgMergedSubMaps, rowSubMaps, false));
417 } else {
418 Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
419 RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
420 // TODO think about Thyra style maps: we cannot use thyra style maps when recombining several blocks!!!
421 // The GIDs might not start with 0 and may not be consecutive!
422 Teuchos::RCP<const Map> submap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
423 std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
424 rgMapExtractor = Teuchos::rcp(new MapExtractor(submap, rowSubMaps, false));
425 }
426
427 Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::null;
428 if(colSz > 0) {
429 std::vector<Teuchos::RCP<const Map> > colSubMaps (colSz, Teuchos::null);
430 for(size_t j = 0; j < colSz; j++) {
431 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
432 colSubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,false/*xpetra*/);
433 TEUCHOS_ASSERT(colSubMaps[j].is_null()==false);
434 }
435 Teuchos::RCP<const Map> doMergedSubMaps = MapUtils::concatenateMaps(colSubMaps);
436 doMapExtractor = Teuchos::rcp(new MapExtractor(doMergedSubMaps, colSubMaps, false));
437 } else {
438 Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
439 RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
440 // TODO think about Thyra style maps: we cannot use thyra style maps when recombining several blocks!!!
441 // The GIDs might not start with 0 and may not be consecutive!
442 Teuchos::RCP<const Map> submap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
443 std::vector<Teuchos::RCP<const Map> > colSubMaps (1, submap);
444 doMapExtractor = Teuchos::rcp(new MapExtractor(submap, colSubMaps, false));
445 }
446
447 rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
448
449 size_t cntNNZ = 0;
450
451 if (rowSz == 0 && colSz > 0) {
452 for(size_t j = 0; j < colSz; j++) {
453 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
454 Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowMgr, colSubMgr, bmat);
455 rbmat->setMatrix(0,j,Teuchos::rcp_const_cast<Matrix>(submat));
456 if(submat != Teuchos::null) cntNNZ += submat->getLocalNumEntries();
457 }
458 } else if (rowSz > 0 && colSz == 0) {
459 for(size_t i = 0; i < rowSz; i++) {
460 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
461 Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowSubMgr, colMgr, bmat);
462 rbmat->setMatrix(i,0,Teuchos::rcp_const_cast<Matrix>(submat));
463 if(submat != Teuchos::null) cntNNZ += submat->getLocalNumEntries();
464 }
465 } else {
466 for(size_t i = 0; i < rowSz; i++) {
467 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
468 for(size_t j = 0; j < colSz; j++) {
469 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
470 Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowSubMgr, colSubMgr, bmat);
471 rbmat->setMatrix(i,j,Teuchos::rcp_const_cast<Matrix>(submat));
472 if(submat != Teuchos::null) cntNNZ += submat->getLocalNumEntries();
473 }
474 }
475 }
476 TEUCHOS_ASSERT(rbmat->getLocalNumEntries() == cntNNZ);
477 }
478 rbmat->fillComplete();
479 return rbmat;
480}
481
482 //MapExtractor(const std::vector<RCP<const Map> >& maps, const std::vector<RCP<const Map> >& thyramaps);
483
484template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
486
492
493 TEUCHOS_ASSERT(bmat->getRangeMapExtractor()->getThyraMode() == true);
494 TEUCHOS_ASSERT(bmat->getDomainMapExtractor()->getThyraMode() == true);
495
496 // number of sub blocks
497 size_t rowSz = rowMgr->GetNumBlocks();
498 size_t colSz = colMgr->GetNumBlocks();
499
500 Teuchos::RCP<BlockedCrsMatrix> rbmat = Teuchos::null;
501
502 if(rowSz == 0 && colSz == 0) {
503 // it is a leaf node
504 Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
505 Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
506
507 // this matrix uses Thyra style GIDs as global row, range, domain and column indices
508 Teuchos::RCP<Matrix> mat = bmat->getMatrix(rowleaf->GetIndex(), colleaf->GetIndex());
509
510 if(mat == Teuchos::null) return Teuchos::null; //std::cout << "Block " << rowleaf->GetIndex() << "," << colleaf->GetIndex() << " is zero" << std::endl;
511
512 // check, whether leaf node is of type Xpetra::CrsMatrixWrap
513 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> > matwrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
514 if(matwrap != Teuchos::null) {
516 // build map extractors
517 RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
518 // extract Xpetra and Thyra based GIDs
519 Teuchos::RCP<const Map> xpsubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
520 Teuchos::RCP<const Map> thysubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),true);
521 std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
522 std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
523 // use expert constructor
524 Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
525
526 RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
527 // extract Xpetra and Thyra based GIDs
528 Teuchos::RCP<const Map> xpsubmap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
529 Teuchos::RCP<const Map> tysubmap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),true);
530 std::vector<Teuchos::RCP<const Map> > colXpSubMaps (1, xpsubmap2);
531 std::vector<Teuchos::RCP<const Map> > colTySubMaps (1, tysubmap2);
532 // use expert constructor
533 Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps, colTySubMaps));
534
536 // build reordered block operator
537 rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
538 rbmat->setMatrix(0,0,mat);
539 } else {
540 // If leaf node is already wrapped into a blocked matrix do not wrap it again.
541 rbmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
542 TEUCHOS_ASSERT(rbmat != Teuchos::null);
543 }
544 TEUCHOS_ASSERT(mat->getLocalNumEntries() == rbmat->getLocalNumEntries());
545 } else {
546 // create the map extractors
547 // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
548 Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::null;
549 if(rowSz > 0) {
550 std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (rowSz, Teuchos::null);
551 std::vector<Teuchos::RCP<const Map> > rowTySubMaps (rowSz, Teuchos::null);
552 for(size_t i = 0; i < rowSz; i++) {
553 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
554 // extract Xpetra and Thyra based merged GIDs
555 rowXpSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,false);
556 rowTySubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,true);
557 TEUCHOS_ASSERT(rowXpSubMaps[i].is_null()==false);
558 TEUCHOS_ASSERT(rowTySubMaps[i].is_null()==false);
559 }
560 // use expert constructor
561 rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
562 } else {
563 Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
564 RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
565 // extract Xpetra and Thyra based GIDs
566 Teuchos::RCP<const Map> xpsubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
567 Teuchos::RCP<const Map> thysubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),true);
568 std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
569 std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
570 // use expert constructor
571 rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
572 }
573
574 Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::null;
575 if(colSz > 0) {
576 std::vector<Teuchos::RCP<const Map> > colXpSubMaps (colSz, Teuchos::null);
577 std::vector<Teuchos::RCP<const Map> > colTySubMaps (colSz, Teuchos::null);
578 for(size_t j = 0; j < colSz; j++) {
579 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
580 // extract Xpetra and Thyra based merged GIDs
581 colXpSubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,false);
582 colTySubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,true);
583 TEUCHOS_ASSERT(colXpSubMaps[j].is_null()==false);
584 TEUCHOS_ASSERT(colTySubMaps[j].is_null()==false);
585 }
586 // use expert constructor
587 doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps,colTySubMaps));
588 } else {
589 Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
590 RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
591 // extract Xpetra and Thyra based GIDs
592 Teuchos::RCP<const Map> xpsubmap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
593 Teuchos::RCP<const Map> tysubmap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),true);
594 std::vector<Teuchos::RCP<const Map> > colXpSubMaps (1, xpsubmap);
595 std::vector<Teuchos::RCP<const Map> > colTySubMaps (1, tysubmap);
596 // use expert constructor
597 doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps, colTySubMaps));
598 }
599
600 // TODO matrix should have both rowMgr and colMgr??
601 rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
602
603 size_t cntNNZ = 0;
604
605 if (rowSz == 0 && colSz > 0) {
606 for(size_t j = 0; j < colSz; j++) {
607 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
608 Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowMgr, colSubMgr, bmat);
609 rbmat->setMatrix(0,j,Teuchos::rcp_const_cast<Matrix>(submat));
610 if(submat != Teuchos::null) cntNNZ += submat->getLocalNumEntries();
611 }
612 } else if (rowSz > 0 && colSz == 0) {
613 for(size_t i = 0; i < rowSz; i++) {
614 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
615 Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowSubMgr, colMgr, bmat);
616 rbmat->setMatrix(i,0,Teuchos::rcp_const_cast<Matrix>(submat));
617 if(submat != Teuchos::null) cntNNZ += submat->getLocalNumEntries();
618 }
619 } else {
620 for(size_t i = 0; i < rowSz; i++) {
621 Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
622 for(size_t j = 0; j < colSz; j++) {
623 Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
624 Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowSubMgr, colSubMgr, bmat);
625 rbmat->setMatrix(i,j,Teuchos::rcp_const_cast<Matrix>(submat));
626 if(submat != Teuchos::null) cntNNZ += submat->getLocalNumEntries();
627 }
628 }
629 }
630 TEUCHOS_ASSERT(rbmat->getLocalNumEntries() == cntNNZ);
631 }
632
633 rbmat->fillComplete();
634 return rbmat;
635}
636
637template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
639 TEUCHOS_ASSERT(bmat->getRangeMapExtractor()->getThyraMode() == bmat->getDomainMapExtractor()->getThyraMode());
641 if(bmat->getRangeMapExtractor()->getThyraMode() == false) {
642 rbmat = mergeSubBlocks(brm, brm, bmat);
643 } else {
644 rbmat = mergeSubBlocksThyra(brm, brm, bmat);
645 }
646
647 // TAW, 6/7/2016: rbmat might be Teuchos::null for empty blocks!
648 return rbmat;
649}
650
651} //namespace Xpetra
652
653#define XPETRA_REORDEREDBLOCKEDCRSMATRIX_SHORT
654#endif /* XPETRA_REORDEREDBLOCKEDCRSMATRIX_HPP */
static const EVerbosityLevel verbLevel_default
void swap(RCP< T > &r_ptr)
bool is_null() const
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual size_t Rows() const
number of row blocks
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode.
RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode, Scalar alpha, Scalar beta, bool sumInterfaceValues, const RCP< Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > > &regionInterfaceImporter, const Teuchos::ArrayRCP< LocalOrdinal > &regionInterfaceLIDs) const
sparse matrix-multivector multiplication for the region layout matrices (currently no blocked impleme...
virtual size_t Cols() const
number of column blocks
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.
Xpetra-specific matrix class.
virtual void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta)=0
Update multi-vector values with scaled values of A, this = beta*this + alpha*A.
virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode, Scalar alpha, Scalar beta, bool sumInterfaceValues, const RCP< Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > > &regionInterfaceImporter, const Teuchos::ArrayRCP< LocalOrdinal > &regionInterfaceLIDs) const
sparse matrix-multivector multiplication for the region layout matrices (currently no blocked impleme...
std::string description() const
Return a simple one-line description of this object.
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.
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.
Teuchos::RCP< const Xpetra::BlockReorderManager > brm_
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)
Constructor.
Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlockMaps(Teuchos::RCP< const Xpetra::BlockReorderManager > brm)
Teuchos::RCP< const Xpetra::BlockReorderManager > getBlockReorderManager()
Returns internal BlockReorderManager object.
Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > fullOp_
#define TEUCHOS_ASSERT(assertion_test)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace
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)
std::string viewLabel_t
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)
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)
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::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > buildReorderedBlockedMultiVector(Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bvec)