Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_BlockedCrsMatrix.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_BLOCKEDCRSMATRIX_HPP
47#define XPETRA_BLOCKEDCRSMATRIX_HPP
48
49#include <KokkosCompat_DefaultNode.hpp>
50
52#include <Teuchos_Hashtable.hpp>
53
54#include "Xpetra_ConfigDefs.hpp"
55#include "Xpetra_Exceptions.hpp"
56
57#include "Xpetra_MapFactory.hpp"
58#include "Xpetra_MultiVector.hpp"
59#include "Xpetra_BlockedMultiVector.hpp"
60#include "Xpetra_MultiVectorFactory.hpp"
61#include "Xpetra_BlockedVector.hpp"
62#include "Xpetra_CrsGraph.hpp"
63#include "Xpetra_CrsMatrix.hpp"
65
66#include "Xpetra_MapExtractor.hpp"
68
69#include "Xpetra_Matrix.hpp"
71#include "Xpetra_CrsMatrixWrap.hpp"
72
73#ifdef HAVE_XPETRA_THYRA
74#include <Thyra_ProductVectorSpaceBase.hpp>
75#include <Thyra_VectorSpaceBase.hpp>
76#include <Thyra_LinearOpBase.hpp>
77#include <Thyra_BlockedLinearOpBase.hpp>
78#include <Thyra_PhysicallyBlockedLinearOpBase.hpp>
79#include "Xpetra_ThyraUtils.hpp"
80#endif
81
82#include "Xpetra_VectorFactory.hpp"
83
84
89namespace Xpetra {
90
91#ifdef HAVE_XPETRA_THYRA
92 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> class ThyraUtils;
93#endif
94
95 typedef std::string viewLabel_t;
96
97 template <class Scalar,
98 class LocalOrdinal,
99 class GlobalOrdinal,
102 public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
103 public:
104 typedef Scalar scalar_type;
105 typedef LocalOrdinal local_ordinal_type;
106 typedef GlobalOrdinal global_ordinal_type;
107 typedef Node node_type;
108
109 private:
110#undef XPETRA_BLOCKEDCRSMATRIX_SHORT
112
113 public:
114
116
117
119
125 const Teuchos::RCP<const BlockedMap>& domainMaps,
126 size_t numEntriesPerRow)
127 : is_diagonal_(true)
128 {
131 bRangeThyraMode_ = rangeMaps->getThyraMode();
132 bDomainThyraMode_ = domainMaps->getThyraMode();
133
134 blocks_.reserve(Rows()*Cols());
135
136 // add CrsMatrix objects in row,column order
137 for (size_t r = 0; r < Rows(); ++r)
138 for (size_t c = 0; c < Cols(); ++c)
139 blocks_.push_back(MatrixFactory::Build(getRangeMap(r, bRangeThyraMode_), numEntriesPerRow));
140
141 // Default view
143 }
144
146
154 Teuchos::RCP<const MapExtractor>& domainMapExtractor,
155 size_t numEntriesPerRow)
156 : is_diagonal_(true), domainmaps_(domainMapExtractor), rangemaps_(rangeMapExtractor)
157 {
158 bRangeThyraMode_ = rangeMapExtractor->getThyraMode();
159 bDomainThyraMode_ = domainMapExtractor->getThyraMode();
160
161 blocks_.reserve(Rows()*Cols());
162
163 // add CrsMatrix objects in row,column order
164 for (size_t r = 0; r < Rows(); ++r)
165 for (size_t c = 0; c < Cols(); ++c)
166 blocks_.push_back(MatrixFactory::Build(getRangeMap(r, bRangeThyraMode_), numEntriesPerRow));
167
168 // Default view
170 }
171
172#ifdef HAVE_XPETRA_THYRA
174
179 BlockedCrsMatrix(const Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> >& thyraOp, const Teuchos::RCP<const Teuchos::Comm<int> >& /* comm */)
180 : is_diagonal_(true), thyraOp_(thyraOp)
181 {
182 // extract information from Thyra blocked operator and rebuilt information
183 const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productRangeSpace = thyraOp->productRange();
184 const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productDomainSpace = thyraOp->productDomain();
185
186 int numRangeBlocks = productRangeSpace->numBlocks();
187 int numDomainBlocks = productDomainSpace->numBlocks();
188
189 // build range map extractor from Thyra::BlockedLinearOpBase object
190 std::vector<Teuchos::RCP<const Map> > subRangeMaps(numRangeBlocks);
191 for (size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
192 for (size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
193 if (thyraOp->blockExists(r,c)) {
194 // we only need at least one block in each block row to extract the range map
195 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
197 Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(const_op);
198 subRangeMaps[r] = xop->getRangeMap();
199 if(r!=c) is_diagonal_ = false;
200 break;
201 }
202 }
203 }
204 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullRangeMap = mergeMaps(subRangeMaps);
205
206 // check whether the underlying Thyra operator uses Thyra-style numbering (default for most applications) or
207 // Xpetra style numbering
208 bool bRangeUseThyraStyleNumbering = false;
209 size_t numAllElements = 0;
210 for(size_t v = 0; v < subRangeMaps.size(); ++v) {
211 numAllElements += subRangeMaps[v]->getGlobalNumElements();
212 }
213 if ( fullRangeMap->getGlobalNumElements() != numAllElements) bRangeUseThyraStyleNumbering = true;
214 rangemaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullRangeMap, subRangeMaps, bRangeUseThyraStyleNumbering));
215
216 // build domain map extractor from Thyra::BlockedLinearOpBase object
217 std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subDomainMaps(numDomainBlocks);
218 for (size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
219 for (size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
220 if (thyraOp->blockExists(r,c)) {
221 // we only need at least one block in each block row to extract the range map
222 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
224 Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toXpetra(const_op);
225 subDomainMaps[c] = xop->getDomainMap();
226 break;
227 }
228 }
229 }
230 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullDomainMap = mergeMaps(subDomainMaps);
231 // plausibility check for numbering style (Xpetra or Thyra)
232 bool bDomainUseThyraStyleNumbering = false;
233 numAllElements = 0;
234 for(size_t v = 0; v < subDomainMaps.size(); ++v) {
235 numAllElements += subDomainMaps[v]->getGlobalNumElements();
236 }
237 if (fullDomainMap->getGlobalNumElements() != numAllElements) bDomainUseThyraStyleNumbering = true;
238 domainmaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullDomainMap, subDomainMaps, bDomainUseThyraStyleNumbering));
239
240 // store numbering mode
241 bRangeThyraMode_ = bRangeUseThyraStyleNumbering;
242 bDomainThyraMode_ = bDomainUseThyraStyleNumbering;
243
244 // add CrsMatrix objects in row,column order
245 blocks_.reserve(Rows()*Cols());
246 for (size_t r = 0; r < Rows(); ++r) {
247 for (size_t c = 0; c < Cols(); ++c) {
248 if(thyraOp->blockExists(r,c)) {
249 // TODO we do not support nested Thyra operators here!
250 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
251 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar> >(const_op); // cast away const
253 Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toXpetra(op);
255 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LO,GO,Node> >(xop, true);
256 blocks_.push_back(xwrap);
257 } else {
258 // add empty block
260 }
261 }
262 }
263 // Default view
265 }
266
267 private:
269
277 // TODO merging for Thyra mode is missing (similar to what we do in constructor of MapExtractor
278
279 // merge submaps to global map
280 std::vector<GlobalOrdinal> gids;
281 for(size_t tt = 0; tt<subMaps.size(); ++tt) {
283#if 1
284 Teuchos::ArrayView< const GlobalOrdinal > subMapGids = subMap->getLocalElementList();
285 gids.insert(gids.end(), subMapGids.begin(), subMapGids.end());
286#else
287 size_t myNumElements = subMap->getLocalNumElements();
288 for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(myNumElements); ++l) {
289 GlobalOrdinal gid = subMap->getGlobalElement(l);
290 gids.push_back(gid);
291 }
292#endif
293 }
294
295 // we have to sort the matrix entries and get rid of the double entries
296 // since we use this to detect Thyra-style numbering or Xpetra-style
297 // numbering. In Thyra-style numbering mode, the Xpetra::MapExtractor builds
298 // the correct row maps.
300 std::sort(gids.begin(), gids.end());
301 gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
302 Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
303 Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullMap = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(subMaps[0]->lib(), INVALID, gidsView, subMaps[0]->getIndexBase(), subMaps[0]->getComm());
304 return fullMap;
305 }
306
307 public:
308#endif
309
311 virtual ~BlockedCrsMatrix() {}
312
314
315
317
318
320
345 void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
346 XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertGlobalValues");
347 if (Rows() == 1 && Cols () == 1) {
348 getMatrix(0,0)->insertGlobalValues(globalRow, cols, vals);
349 return;
350 }
351 throw Xpetra::Exceptions::RuntimeError("insertGlobalValues not supported by BlockedCrsMatrix");
352 }
353
355
365 void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
366 XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertLocalValues");
367 if (Rows() == 1 && Cols () == 1) {
368 getMatrix(0,0)->insertLocalValues(localRow, cols, vals);
369 return;
370 }
371 throw Xpetra::Exceptions::RuntimeError("insertLocalValues not supported by BlockedCrsMatrix");
372 }
373
375 XPETRA_MONITOR("XpetraBlockedCrsMatrix::removeEmptyProcessesInPlace");
376 if (Rows() == 1 && Cols () == 1) {
377 getMatrix(0,0)->removeEmptyProcessesInPlace(newMap);
378 return;
379 }
380 throw Xpetra::Exceptions::RuntimeError("removeEmptyProcesses not supported by BlockedCrsMatrix");
381 }
382
384
392 void replaceGlobalValues(GlobalOrdinal globalRow,
394 const ArrayView<const Scalar> &vals) {
395 XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceGlobalValues");
396 if (Rows() == 1 && Cols () == 1) {
397 getMatrix(0,0)->replaceGlobalValues(globalRow,cols,vals);
398 return;
399 }
400 throw Xpetra::Exceptions::RuntimeError("replaceGlobalValues not supported by BlockedCrsMatrix");
401 }
402
404
408 void replaceLocalValues(LocalOrdinal localRow,
410 const ArrayView<const Scalar> &vals) {
411 XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceLocalValues");
412 if (Rows() == 1 && Cols () == 1) {
413 getMatrix(0,0)->replaceLocalValues(localRow,cols,vals);
414 return;
415 }
416 throw Xpetra::Exceptions::RuntimeError("replaceLocalValues not supported by BlockedCrsMatrix");
417 }
418
420 virtual void setAllToScalar(const Scalar& alpha) {
421 XPETRA_MONITOR("XpetraBlockedCrsMatrix::setAllToScalar");
422 for (size_t row = 0; row < Rows(); row++) {
423 for (size_t col = 0; col < Cols(); col++) {
424 if (!getMatrix(row,col).is_null()) {
425 getMatrix(row,col)->setAllToScalar(alpha);
426 }
427 }
428 }
429 }
430
432 void scale(const Scalar& alpha) {
433 XPETRA_MONITOR("XpetraBlockedCrsMatrix::scale");
434 for (size_t row = 0; row < Rows(); row++) {
435 for (size_t col = 0; col < Cols(); col++) {
436 if (!getMatrix(row,col).is_null()) {
437 getMatrix(row,col)->scale(alpha);
438 }
439 }
440 }
441 }
442
444
446
447
455 void resumeFill(const RCP< ParameterList >& params = null) {
456 XPETRA_MONITOR("XpetraBlockedCrsMatrix::resumeFill");
457 for (size_t row = 0; row < Rows(); row++) {
458 for (size_t col = 0; col < Cols(); col++) {
459 if (!getMatrix(row,col).is_null()) {
460 getMatrix(row,col)->resumeFill(params);
461 }
462 }
463 }
464 }
465
471 void fillComplete(const RCP<const Map>& domainMap, const RCP<const Map>& rangeMap, const RCP<ParameterList>& params = null) {
472 XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
473 if (Rows() == 1 && Cols () == 1) {
474 getMatrix(0,0)->fillComplete(domainMap, rangeMap, params);
475 return;
476 }
477 fillComplete(params);
478 }
479
494 void fillComplete(const RCP<ParameterList>& params = null) {
495 XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
496 TEUCHOS_TEST_FOR_EXCEPTION(rangemaps_==Teuchos::null, Xpetra::Exceptions::RuntimeError,"BlockedCrsMatrix::fillComplete: rangemaps_ is not set. Error.");
497
498 for (size_t r = 0; r < Rows(); ++r)
499 for (size_t c = 0; c < Cols(); ++c) {
500 if(getMatrix(r,c) != Teuchos::null) {
501 Teuchos::RCP<Matrix> Ablock = getMatrix(r,c);
502 if(r!=c) is_diagonal_ = false;
503 if (Ablock != Teuchos::null && !Ablock->isFillComplete())
504 Ablock->fillComplete(getDomainMap(c, bDomainThyraMode_), getRangeMap(r, bRangeThyraMode_), params);
505 }
506 }
507
508#if 0
509 // get full row map
510 RCP<const Map> rangeMap = rangemaps_->getFullMap();
511 fullrowmap_ = MapFactory::Build(rangeMap()->lib(), rangeMap()->getGlobalNumElements(), rangeMap()->getLocalElementList(), rangeMap()->getIndexBase(), rangeMap()->getComm());
512
513 // build full col map
514 fullcolmap_ = Teuchos::null; // delete old full column map
515
516 std::vector<GO> colmapentries;
517 for (size_t c = 0; c < Cols(); ++c) {
518 // copy all local column lids of all block rows to colset
519 std::set<GO> colset;
520 for (size_t r = 0; r < Rows(); ++r) {
522
523 if (Ablock != Teuchos::null) {
524 Teuchos::ArrayView<const GO> colElements = Ablock->getColMap()->getLocalElementList();
525 Teuchos::RCP<const Map> colmap = Ablock->getColMap();
526 copy(colElements.getRawPtr(), colElements.getRawPtr() + colElements.size(), inserter(colset, colset.begin()));
527 }
528 }
529
530 // remove duplicates (entries which are in column maps of more than one block row)
531 colmapentries.reserve(colmapentries.size() + colset.size());
532 copy(colset.begin(), colset.end(), back_inserter(colmapentries));
533 sort(colmapentries.begin(), colmapentries.end());
534 typename std::vector<GO>::iterator gendLocation;
535 gendLocation = std::unique(colmapentries.begin(), colmapentries.end());
536 colmapentries.erase(gendLocation,colmapentries.end());
537 }
538
539 // sum up number of local elements
540 size_t numGlobalElements = 0;
541 Teuchos::reduceAll(*(rangeMap->getComm()), Teuchos::REDUCE_SUM, colmapentries.size(), Teuchos::outArg(numGlobalElements));
542
543 // store global full column map
545 fullcolmap_ = MapFactory::Build(rangeMap->lib(), numGlobalElements, aView, 0, rangeMap->getComm());
546#endif
547 }
548
550
552
555 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumRows");
556 global_size_t globalNumRows = 0;
557
558 for (size_t row = 0; row < Rows(); row++)
559 for (size_t col = 0; col < Cols(); col++)
560 if (!getMatrix(row,col).is_null()) {
561 globalNumRows += getMatrix(row,col)->getGlobalNumRows();
562 break; // we need only one non-null matrix in a row
563 }
564
565 return globalNumRows;
566 }
567
569
572 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumCols");
573 global_size_t globalNumCols = 0;
574
575 for (size_t col = 0; col < Cols(); col++)
576 for (size_t row = 0; row < Rows(); row++)
577 if (!getMatrix(row,col).is_null()) {
578 globalNumCols += getMatrix(row,col)->getGlobalNumCols();
579 break; // we need only one non-null matrix in a col
580 }
581
582 return globalNumCols;
583 }
584
586 size_t getLocalNumRows() const {
587 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalNumRows");
588 global_size_t nodeNumRows = 0;
589
590 for (size_t row = 0; row < Rows(); ++row)
591 for (size_t col = 0; col < Cols(); col++)
592 if (!getMatrix(row,col).is_null()) {
593 nodeNumRows += getMatrix(row,col)->getLocalNumRows();
594 break; // we need only one non-null matrix in a row
595 }
596
597 return nodeNumRows;
598 }
599
602 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumEntries");
603 global_size_t globalNumEntries = 0;
604
605 for (size_t row = 0; row < Rows(); ++row)
606 for (size_t col = 0; col < Cols(); ++col)
607 if (!getMatrix(row,col).is_null())
608 globalNumEntries += getMatrix(row,col)->getGlobalNumEntries();
609
610 return globalNumEntries;
611 }
612
614 size_t getLocalNumEntries() const {
615 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalNumEntries");
616 global_size_t nodeNumEntries = 0;
617
618 for (size_t row = 0; row < Rows(); ++row)
619 for (size_t col = 0; col < Cols(); ++col)
620 if (!getMatrix(row,col).is_null())
621 nodeNumEntries += getMatrix(row,col)->getLocalNumEntries();
622
623 return nodeNumEntries;
624 }
625
627
628 size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
629 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNumEntriesInLocalRow");
630 GlobalOrdinal gid = this->getRowMap()->getGlobalElement(localRow);
631 size_t row = getBlockedRangeMap()->getMapIndexForGID(gid);
632 LocalOrdinal lid = getBlockedRangeMap()->getMap(row)->getLocalElement(gid);
633 size_t numEntriesInLocalRow = 0;
634 for (size_t col = 0; col < Cols(); ++col)
635 if (!getMatrix(row,col).is_null())
636 numEntriesInLocalRow += getMatrix(row,col)->getNumEntriesInLocalRow(lid);
637 return numEntriesInLocalRow;
638 }
639
641
642 size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const {
643 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNumEntriesInGlobalRow");
644 size_t row = getBlockedRangeMap()->getMapIndexForGID(globalRow);
645 size_t numEntriesInGlobalRow = 0;
646 for (size_t col = 0; col < Cols(); ++col)
647 if (!getMatrix(row,col).is_null())
648 numEntriesInGlobalRow += getMatrix(row,col)->getNumEntriesInGlobalRow(globalRow);
649 return numEntriesInGlobalRow;
650 }
651
653
656 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalMaxNumRowEntries");
657 global_size_t globalMaxEntries = 0;
658
659 for (size_t row = 0; row < Rows(); row++) {
660 global_size_t globalMaxEntriesBlockRows = 0;
661 for (size_t col = 0; col < Cols(); col++) {
662 if (!getMatrix(row,col).is_null()) {
663 globalMaxEntriesBlockRows += getMatrix(row,col)->getGlobalMaxNumRowEntries();
664 }
665 }
666 if(globalMaxEntriesBlockRows > globalMaxEntries)
667 globalMaxEntries = globalMaxEntriesBlockRows;
668 }
669 return globalMaxEntries;
670 }
671
673
676 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalMaxNumRowEntries");
677 size_t localMaxEntries = 0;
678
679 for (size_t row = 0; row < Rows(); row++) {
680 size_t localMaxEntriesBlockRows = 0;
681 for (size_t col = 0; col < Cols(); col++) {
682 if (!getMatrix(row,col).is_null()) {
683 localMaxEntriesBlockRows += getMatrix(row,col)->getLocalMaxNumRowEntries();
684 }
685 }
686 if(localMaxEntriesBlockRows > localMaxEntries)
687 localMaxEntries = localMaxEntriesBlockRows;
688 }
689 return localMaxEntries;
690 }
691
693
696 bool isLocallyIndexed() const {
697 XPETRA_MONITOR("XpetraBlockedCrsMatrix::isLocallyIndexed");
698 for (size_t i = 0; i < blocks_.size(); ++i)
699 if (blocks_[i] != Teuchos::null && !blocks_[i]->isLocallyIndexed())
700 return false;
701 return true;
702 }
703
705
708 bool isGloballyIndexed() const {
709 XPETRA_MONITOR("XpetraBlockedCrsMatrix::isGloballyIndexed");
710 for (size_t i = 0; i < blocks_.size(); i++)
711 if (blocks_[i] != Teuchos::null && !blocks_[i]->isGloballyIndexed())
712 return false;
713 return true;
714 }
715
717 bool isFillComplete() const {
718 XPETRA_MONITOR("XpetraBlockedCrsMatrix::isFillComplete");
719 for (size_t i = 0; i < blocks_.size(); i++)
720 if (blocks_[i] != Teuchos::null && !blocks_[i]->isFillComplete())
721 return false;
722 return true;
723 }
724
726
740 virtual void getLocalRowCopy(LocalOrdinal LocalRow,
741 const ArrayView<LocalOrdinal>& Indices,
742 const ArrayView<Scalar>& Values,
743 size_t &NumEntries) const {
744 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowCopy");
745 if (Rows() == 1 && Cols () == 1) {
746 getMatrix(0,0)->getLocalRowCopy(LocalRow, Indices, Values, NumEntries);
747 return;
748 }
749 throw Xpetra::Exceptions::RuntimeError("getLocalRowCopy not supported by BlockedCrsMatrix" );
750 }
751
753
762 void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values) const {
763 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalRowView");
764 if (Rows() == 1 && Cols () == 1) {
765 getMatrix(0,0)->getGlobalRowView(GlobalRow, indices, values);
766 return;
767 }
768 throw Xpetra::Exceptions::RuntimeError("getGlobalRowView not supported by BlockedCrsMatrix");
769 }
770
772
781 void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values) const {
782 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowView");
783 if (Rows() == 1 && Cols () == 1) {
784 getMatrix(0,0)->getLocalRowView(LocalRow, indices, values);
785 return;
786 }
787 else if(is_diagonal_) {
788 GlobalOrdinal gid = this->getRowMap()->getGlobalElement(LocalRow);
789 size_t row = getBlockedRangeMap()->getMapIndexForGID(gid);
790 getMatrix(row,row)->getLocalRowView(getMatrix(row,row)->getRowMap()->getLocalElement(gid),indices,values);
791 return;
792 }
793 throw Xpetra::Exceptions::RuntimeError("getLocalRowView not supported by BlockedCrsMatrix");
794 }
795
797
800 void getLocalDiagCopy(Vector& diag) const {
801 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalDiagCopy");
802
803 //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
804
805 Teuchos::RCP<Vector> rcpdiag = Teuchos::rcpFromRef(diag);
806 Teuchos::RCP<BlockedVector> bdiag = Teuchos::rcp_dynamic_cast<BlockedVector>(rcpdiag);
807
808 // special treatment for 1x1 block matrices
809 // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
810 // BlockedVectors have Vector objects as Leaf objects.
811 if(Rows() == 1 && Cols() == 1 && bdiag.is_null() == true) {
813 rm->getLocalDiagCopy(diag);
814 return;
815 }
816
817 TEUCHOS_TEST_FOR_EXCEPTION(bdiag.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::getLocalDiagCopy: diag must be a Blocked(Multi)Vector.");
818 TEUCHOS_TEST_FOR_EXCEPTION(bdiag->getNumVectors() != 1, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::getLocalDiagCopy: diag must be a Blocked(Multi)Vector with exactly one vector. However, the number of stored vectors is " << bdiag->getNumVectors());
819 TEUCHOS_TEST_FOR_EXCEPTION(bdiag->getBlockedMap()->getNumMaps() != this->Rows(), Xpetra::Exceptions::RuntimeError,
820 "BlockedCrsMatrix::getLocalDiagCopy(): the number of blocks in diag differ from the number of blocks in this operator." );
821 //XPETRA_TEST_FOR_EXCEPTION(bdiag->getMap()->isSameAs(*(getMap())) == false, Xpetra::Exceptions::RuntimeError,
822 // "BlockedCrsMatrix::getLocalDiagCopy(): the map of the vector diag is not compatible with the map of the blocked operator." );
823
824 for (size_t row = 0; row < Rows(); row++) {
826 if (!rm.is_null()) {
827 Teuchos::RCP<Vector> rv = VectorFactory::Build(bdiag->getBlockedMap()->getMap(row,bdiag->getBlockedMap()->getThyraMode()));
828 rm->getLocalDiagCopy(*rv);
829 bdiag->setMultiVector(row,rv,bdiag->getBlockedMap()->getThyraMode());
830 }
831 }
832 }
833
835 void leftScale (const Vector& x) {
836 XPETRA_MONITOR("XpetraBlockedCrsMatrix::leftScale");
837
838 Teuchos::RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
839 Teuchos::RCP<const BlockedVector> bx = Teuchos::rcp_dynamic_cast<const BlockedVector>(rcpx);
840
841 //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
842
843 // special treatment for 1xn block matrices
844 // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
845 // BlockedVectors have Vector objects as Leaf objects.
846 if(Rows() == 1 && bx.is_null() == true) {
847 for (size_t col = 0; col < Cols(); ++col) {
849 rm->leftScale(x);
850 }
851 return;
852 }
853
854 TEUCHOS_TEST_FOR_EXCEPTION(bx.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Blocked(Multi)Vector.");
855 TEUCHOS_TEST_FOR_EXCEPTION(bx->getNumVectors() != 1, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Blocked(Multi)Vector with exactly one vector. However, the number of stored vectors is " << bx->getNumVectors());
856 TEUCHOS_TEST_FOR_EXCEPTION(bx->getBlockedMap()->getNumMaps() != this->Rows(), Xpetra::Exceptions::RuntimeError,
857 "BlockedCrsMatrix::leftScale(): the number of blocks in diag differ from the number of blocks in this operator." );
858
859 for (size_t row = 0; row < Rows(); row++) {
860 Teuchos::RCP<const MultiVector> rmv = bx->getMultiVector(row);
861 Teuchos::RCP<const Vector> rscale = rmv->getVector(0);
862 XPETRA_TEST_FOR_EXCEPTION(rscale.is_null()==true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Vector.");
863 for (size_t col = 0; col < Cols(); ++col) {
864 Teuchos::RCP<Matrix> rm = getMatrix(row,col);
865 if (!rm.is_null()) {
866 rm->leftScale(*rscale);
867 }
868 }
869 }
870 }
871
873 void rightScale (const Vector& x) {
874 XPETRA_MONITOR("XpetraBlockedCrsMatrix::rightScale");
875
876 Teuchos::RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
877 Teuchos::RCP<const BlockedVector> bx = Teuchos::rcp_dynamic_cast<const BlockedVector>(rcpx);
878
879 //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
880
881 // special treatment for nx1 block matrices
882 // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
883 // BlockedVectors have Vector objects as Leaf objects.
884 if(Cols() == 1 && bx.is_null() == true) {
885 for (size_t row = 0; row < Rows(); ++row) {
887 rm->rightScale(x);
888 }
889 return;
890 }
891
892 TEUCHOS_TEST_FOR_EXCEPTION(bx.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::rightScale: x must be a Blocked(Multi)Vector.");
893 TEUCHOS_TEST_FOR_EXCEPTION(bx->getNumVectors() != 1, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::rightScale: x must be a Blocked(Multi)Vector with exactly one vector. However, the number of stored vectors is " << bx->getNumVectors());
894 TEUCHOS_TEST_FOR_EXCEPTION(bx->getBlockedMap()->getNumMaps() != this->Cols(), Xpetra::Exceptions::RuntimeError,
895 "BlockedCrsMatrix::rightScale(): the number of blocks in diag differ from the number of blocks in this operator." );
896
897 for (size_t col = 0; col < Cols(); ++col) {
898 Teuchos::RCP<const MultiVector> rmv = bx->getMultiVector(col);
899 Teuchos::RCP<const Vector> rscale = rmv->getVector(0);
900 XPETRA_TEST_FOR_EXCEPTION(rscale.is_null()==true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Vector.");
901 for (size_t row = 0; row < Rows(); row++) {
902 Teuchos::RCP<Matrix> rm = getMatrix(row,col);
903 if (!rm.is_null()) {
904 rm->rightScale(*rscale);
905 }
906 }
907 }
908 }
909
910
913 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getFrobeniusNorm");
915 for (size_t col = 0; col < Cols(); ++col) {
916 for (size_t row = 0; row < Rows(); ++row) {
917 if(getMatrix(row,col)!=Teuchos::null) {
918 typename ScalarTraits<Scalar>::magnitudeType n = getMatrix(row,col)->getFrobeniusNorm();
919 ret += n * n;
920 }
921 }
922 }
924 }
925
926
928 virtual bool haveGlobalConstants() const {return true;}
929
930
932
934
935
937
957
959
960
962 virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode, Scalar alpha, Scalar beta, bool sumInterfaceValues,
963 const RCP<Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> >& regionInterfaceImporter,
964 const Teuchos::ArrayRCP<LocalOrdinal>& regionInterfaceLIDs) const
965 { }
966
968
971 virtual void apply(const MultiVector& X, MultiVector& Y,
973 Scalar alpha = ScalarTraits<Scalar>::one(),
974 Scalar beta = ScalarTraits<Scalar>::zero()) const
975 {
976 XPETRA_MONITOR("XpetraBlockedCrsMatrix::apply");
977 //using Teuchos::RCP;
978
980 "apply() only supports the following modes: NO_TRANS and TRANS." );
981
982 // check whether input parameters are blocked or not
983 RCP<const MultiVector> refX = rcpFromRef(X);
984 RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
985 //RCP<MultiVector> tmpY = rcpFromRef(Y);
986 //RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
987
988 // TODO get rid of me: adapt MapExtractor
989 bool bBlockedX = (refbX != Teuchos::null) ? true : false;
990
991 // create (temporary) vectors for output
992 // In the end we call Y.update(alpha, *tmpY, beta). Therefore we need a new vector storing the temporary results
994
995 //RCP<Teuchos::FancyOStream> out = rcp(new Teuchos::FancyOStream(rcp(&std::cout,false)));
996
997 SC one = ScalarTraits<SC>::one();
998
999 if (mode == Teuchos::NO_TRANS) {
1000
1001 for (size_t row = 0; row < Rows(); row++) {
1002 RCP<MultiVector> Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, true);
1003 for (size_t col = 0; col < Cols(); col++) {
1004
1005 // extract matrix block
1006 RCP<Matrix> Ablock = getMatrix(row, col);
1007
1008 if (Ablock.is_null())
1009 continue;
1010
1011 // check whether Ablock is itself a blocked operator
1012 // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
1013 bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
1014
1015 // input/output vectors for local block operation
1016 RCP<const MultiVector> Xblock = Teuchos::null; // subpart of X vector to be applied to subblock of A
1017
1018 // extract sub part of X using Xpetra or Thyra GIDs
1019 // if submatrix is again blocked, we extract it using Xpetra style gids. If it is a single
1020 // block matrix we use the Thyra or Xpetra style GIDs that are used to store the matrix
1021 if(bBlockedX) Xblock = domainmaps_->ExtractVector(refbX, col, bDomainThyraMode_);
1022 else Xblock = domainmaps_->ExtractVector(refX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
1023
1024 RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, false); // subpart of Y vector containing part of solution of Xblock applied to Ablock
1025 Ablock->apply(*Xblock, *tmpYblock);
1026
1027 Yblock->update(one, *tmpYblock, one);
1028 }
1029 rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
1030 }
1031
1032 } else if (mode == Teuchos::TRANS) {
1033 // TODO: test me!
1034 for (size_t col = 0; col < Cols(); col++) {
1035 RCP<MultiVector> Yblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_, true);
1036
1037 for (size_t row = 0; row<Rows(); row++) {
1038 RCP<Matrix> Ablock = getMatrix(row, col);
1039
1040 if (Ablock.is_null())
1041 continue;
1042
1043 // check whether Ablock is itself a blocked operator
1044 // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
1045 bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
1046
1047 RCP<const MultiVector> Xblock = Teuchos::null;
1048
1049 // extract sub part of X using Xpetra or Thyra GIDs
1050 if(bBlockedX) Xblock = rangemaps_->ExtractVector(refbX, row, bRangeThyraMode_);
1051 else Xblock = rangemaps_->ExtractVector(refX, row, bBlockedSubMatrix == true ? false : bRangeThyraMode_);
1052 RCP<MultiVector> tmpYblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_, false);
1053 Ablock->apply(*Xblock, *tmpYblock, Teuchos::TRANS);
1054
1055 Yblock->update(one, *tmpYblock, one);
1056 }
1057 domainmaps_->InsertVector(Yblock, col, tmpY, bDomainThyraMode_);
1058 }
1059 }
1060 Y.update(alpha, *tmpY, beta);
1061 }
1062
1064 RCP<const Map > getFullDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getFullDomainMap()"); return domainmaps_->getFullMap(); }
1065
1067 RCP<const BlockedMap > getBlockedDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getBlockedDomainMap()"); return domainmaps_->getBlockedMap(); }
1068
1070 RCP<const Map > getDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap()"); return domainmaps_->getMap(); /*domainmaps_->getFullMap();*/ }
1071
1073 RCP<const Map > getDomainMap(size_t i) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t)"); return domainmaps_->getMap(i, bDomainThyraMode_); }
1074
1076 RCP<const Map > getDomainMap(size_t i, bool bThyraMode) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t,bool)"); return domainmaps_->getMap(i, bThyraMode); }
1077
1079 RCP<const Map > getFullRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap()"); return rangemaps_->getFullMap(); }
1080
1082 RCP<const BlockedMap > getBlockedRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getBlockedRangeMap()"); return rangemaps_->getBlockedMap(); }
1083
1085 RCP<const Map > getRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap()"); return rangemaps_->getMap(); /*rangemaps_->getFullMap();*/ }
1086
1088 RCP<const Map > getRangeMap(size_t i) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t)"); return rangemaps_->getMap(i, bRangeThyraMode_); }
1089
1091 RCP<const Map > getRangeMap(size_t i, bool bThyraMode) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t,bool)"); return rangemaps_->getMap(i, bThyraMode); }
1092
1094 RCP<const MapExtractor> getRangeMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMapExtractor()"); return rangemaps_; }
1095
1097 RCP<const MapExtractor> getDomainMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMapExtractor()"); return domainmaps_; }
1098
1100
1102 //{@
1103
1112 virtual void bgs_apply(
1113 const MultiVector& X,
1114 MultiVector& Y,
1115 size_t row,
1117 Scalar alpha = ScalarTraits<Scalar>::one(),
1118 Scalar beta = ScalarTraits<Scalar>::zero()
1119 ) const
1120 {
1121 XPETRA_MONITOR("XpetraBlockedCrsMatrix::bgs_apply");
1122 //using Teuchos::RCP;
1123
1125 "apply() only supports the following modes: NO_TRANS and TRANS." );
1126
1127 // check whether input parameters are blocked or not
1128 RCP<const MultiVector> refX = rcpFromRef(X);
1129 RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
1130 //RCP<MultiVector> tmpY = rcpFromRef(Y);
1131 //RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
1132
1133 bool bBlockedX = (refbX != Teuchos::null) ? true : false;
1134
1135 // create (temporary) vectors for output
1136 // In the end we call Y.update(alpha, *tmpY, beta). Therefore we need a new vector storing the temporary results
1138
1139 SC one = ScalarTraits<SC>::one();
1140
1141 if (mode == Teuchos::NO_TRANS) {
1142 RCP<MultiVector> Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, true);
1143 for (size_t col = 0; col < Cols(); col++) {
1144
1145 // extract matrix block
1146 RCP<Matrix> Ablock = getMatrix(row, col);
1147
1148 if (Ablock.is_null())
1149 continue;
1150
1151 // check whether Ablock is itself a blocked operator
1152 // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
1153 bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
1154
1155 // input/output vectors for local block operation
1156 RCP<const MultiVector> Xblock = Teuchos::null; // subpart of X vector to be applied to subblock of A
1157
1158 // extract sub part of X using Xpetra or Thyra GIDs
1159 // if submatrix is again blocked, we extract it using Xpetra style gids. If it is a single
1160 // block matrix we use the Thyra or Xpetra style GIDs that are used to store the matrix
1161 if(bBlockedX) Xblock = domainmaps_->ExtractVector(refbX, col, bDomainThyraMode_);
1162 else Xblock = domainmaps_->ExtractVector(refX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
1163
1164 RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, false); // subpart of Y vector containing part of solution of Xblock applied to Ablock
1165 Ablock->apply(*Xblock, *tmpYblock);
1166
1167 Yblock->update(one, *tmpYblock, one);
1168 }
1169 rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
1170 } else {
1171 TEUCHOS_TEST_FOR_EXCEPTION(true,Xpetra::Exceptions::NotImplemented,"Xpetar::BlockedCrsMatrix::bgs_apply: not implemented for transpose case.");
1172 }
1173 Y.update(alpha, *tmpY, beta);
1174 }
1175
1176
1178
1179
1181 //{@
1182
1185 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMap");
1186 if (Rows() == 1 && Cols () == 1) {
1187 return getMatrix(0,0)->getMap();
1188 }
1189 throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getMap(): operation not supported.");
1190 }
1191
1193 void doImport(const Matrix &source, const Import& importer, CombineMode CM) {
1194 XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
1195 if (Rows() == 1 && Cols () == 1) {
1196 getMatrix(0,0)->doImport(source, importer, CM);
1197 return;
1198 }
1199 throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
1200 }
1201
1203 void doExport(const Matrix& dest, const Import& importer, CombineMode CM) {
1204 XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
1205 if (Rows() == 1 && Cols () == 1) {
1206 getMatrix(0,0)->doExport(dest, importer, CM);
1207 return;
1208 }
1209 throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
1210 }
1211
1213 void doImport(const Matrix& source, const Export& exporter, CombineMode CM) {
1214 XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
1215 if (Rows() == 1 && Cols () == 1) {
1216 getMatrix(0,0)->doImport(source, exporter, CM);
1217 return;
1218 }
1219 throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
1220 }
1221
1223 void doExport(const Matrix& dest, const Export& exporter, CombineMode CM) {
1224 XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
1225 if (Rows() == 1 && Cols () == 1) {
1226 getMatrix(0,0)->doExport(dest, exporter, CM);
1227 return;
1228 }
1229 throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
1230 }
1231
1232 // @}
1233
1235
1236
1238 std::string description() const { return "Xpetra_BlockedCrsMatrix.description()"; }
1239
1242 out << "Xpetra::BlockedCrsMatrix: " << Rows() << " x " << Cols() << std::endl;
1243
1244 if (isFillComplete()) {
1245 out << "BlockMatrix is fillComplete" << std::endl;
1246
1247 /*if(fullrowmap_ != Teuchos::null) {
1248 out << "fullRowMap" << std::endl;
1249 fullrowmap_->describe(out,verbLevel);
1250 } else {
1251 out << "fullRowMap not set. Check whether block matrix is properly fillCompleted!" << std::endl;
1252 }*/
1253
1254 //out << "fullColMap" << std::endl;
1255 //fullcolmap_->describe(out,verbLevel);
1256
1257 } else {
1258 out << "BlockMatrix is NOT fillComplete" << std::endl;
1259 }
1260
1261 for (size_t r = 0; r < Rows(); ++r)
1262 for (size_t c = 0; c < Cols(); ++c) {
1263 if(getMatrix(r,c)!=Teuchos::null) {
1264 out << "Block(" << r << "," << c << ")" << std::endl;
1265 getMatrix(r,c)->describe(out,verbLevel);
1266 } else out << "Block(" << r << "," << c << ") = null" << std::endl;
1267 }
1268 }
1269
1271
1272 void setObjectLabel( const std::string &objectLabel ) {
1273 XPETRA_MONITOR("TpetraBlockedCrsMatrix::setObjectLabel");
1274 for (size_t r = 0; r < Rows(); ++r)
1275 for (size_t c = 0; c < Cols(); ++c) {
1276 if(getMatrix(r,c)!=Teuchos::null) {
1277 std::ostringstream oss; oss<< objectLabel << "(" << r << "," << c << ")";
1278 getMatrix(r,c)->setObjectLabel(oss.str());
1279 }
1280 }
1281 }
1283
1284
1286 bool hasCrsGraph() const {
1287 if (Rows() == 1 && Cols () == 1) return true;
1288 else return false;
1289 }
1290
1293 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsGraph");
1294 if (Rows() == 1 && Cols () == 1) {
1295 return getMatrix(0,0)->getCrsGraph();
1296 }
1297 throw Xpetra::Exceptions::RuntimeError("getCrsGraph() not supported by BlockedCrsMatrix");
1298 }
1299
1301
1303
1304
1305 virtual bool isDiagonal() const {return is_diagonal_;}
1306
1308 virtual size_t Rows() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Rows"); return rangemaps_->NumMaps(); }
1309
1311 virtual size_t Cols() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Cols"); return domainmaps_->NumMaps(); }
1312
1315 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsMatrix");
1316 TEUCHOS_TEST_FOR_EXCEPTION(Rows()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Rows() << " block rows, though.");
1317 TEUCHOS_TEST_FOR_EXCEPTION(Cols()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Cols() << " block columns, though.");
1318
1319 RCP<Matrix> mat = getMatrix(0,0);
1320 RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
1321 if (bmat == Teuchos::null) return mat;
1322 return bmat->getCrsMatrix();
1323 }
1324
1327 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getInnermostCrsMatrix");
1328 size_t row = Rows()+1, col = Cols()+1;
1329 for (size_t r = 0; r < Rows(); ++r)
1330 for(size_t c = 0; c < Cols(); ++c)
1331 if (getMatrix(r,c) != Teuchos::null) {
1332 row = r;
1333 col = c;
1334 break;
1335 }
1336 TEUCHOS_TEST_FOR_EXCEPTION(row == Rows()+1 || col == Cols()+1, Xpetra::Exceptions::Incompatible, "Xpetra::BlockedCrsMatrix::getInnermostCrsMatrix: Could not find a non-zero sub-block in blocked operator.")
1337 RCP<Matrix> mm = getMatrix(row,col);
1338 RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mm);
1339 if (bmat == Teuchos::null) return mm;
1340 return bmat->getInnermostCrsMatrix();
1341 }
1342
1344 Teuchos::RCP<Matrix> getMatrix(size_t r, size_t c) const {
1345 XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMatrix");
1346 TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
1347 TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
1348
1349 // transfer strided/blocked map information
1350 /* if (blocks_[r*Cols()+c] != Teuchos::null &&
1351 blocks_[r*Cols()+c]->IsView("stridedMaps") == false)
1352 blocks_[r*Cols()+c]->CreateView("stridedMaps", getRangeMap(r,bRangeThyraMode_), getDomainMap(c,bDomainThyraMode_));*/
1353 return blocks_[r*Cols()+c];
1354 }
1355
1357 //void setMatrix(size_t r, size_t c, Teuchos::RCP<CrsMatrix> mat) {
1358 void setMatrix(size_t r, size_t c, Teuchos::RCP<Matrix> mat) {
1359 XPETRA_MONITOR("XpetraBlockedCrsMatrix::setMatrix");
1360 // TODO: if filled -> return error
1361
1362 TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
1363 TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
1364 if(!mat.is_null() && r!=c) is_diagonal_ = false;
1365 // set matrix
1366 blocks_[r*Cols() + c] = mat;
1367 }
1368
1370 // NOTE: This is a rather expensive operation, since all blocks are copied
1371 // into a new big CrsMatrix
1373 XPETRA_MONITOR("XpetraBlockedCrsMatrix::Merge");
1374 using Teuchos::RCP;
1375 using Teuchos::rcp_dynamic_cast;
1376 Scalar one = ScalarTraits<SC>::one();
1377
1379 "BlockedCrsMatrix::Merge: only implemented for Xpetra-style or Thyra-style numbering. No mixup allowed!" );
1380
1382 "BlockedCrsMatrix::Merge: BlockMatrix must be fill-completed." );
1383
1384 LocalOrdinal lclNumRows = getFullRangeMap()->getLocalNumElements();
1385 Teuchos::ArrayRCP<size_t> numEntPerRow (lclNumRows);
1386 for (LocalOrdinal lclRow = 0; lclRow < lclNumRows; ++lclRow)
1387 numEntPerRow[lclRow] = getNumEntriesInLocalRow(lclRow);
1388
1389 RCP<Matrix> sparse = MatrixFactory::Build(getFullRangeMap(), numEntPerRow);
1390
1391 if(bRangeThyraMode_ == false) {
1392 // Xpetra mode
1393 for (size_t i = 0; i < Rows(); i++) {
1394 for (size_t j = 0; j < Cols(); j++) {
1395 if (getMatrix(i,j) != Teuchos::null) {
1396 RCP<const Matrix> mat = getMatrix(i,j);
1397
1398 // recursively call Merge routine
1399 RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1400 if (bMat != Teuchos::null) mat = bMat->Merge();
1401
1402 bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1404 "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1405
1406 // jump over empty blocks
1407 if(mat->getLocalNumEntries() == 0) continue;
1408
1409 this->Add(*mat, one, *sparse, one);
1410 }
1411 }
1412 }
1413 } else {
1414 // Thyra mode
1415 for (size_t i = 0; i < Rows(); i++) {
1416 for (size_t j = 0; j < Cols(); j++) {
1417 if (getMatrix(i,j) != Teuchos::null) {
1418 RCP<const Matrix> mat = getMatrix(i,j);
1419 // recursively call Merge routine
1420 RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1421 if (bMat != Teuchos::null) mat = bMat->Merge();
1422
1423 bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1425 "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1426
1427 // check whether we have a CrsMatrix block (no blocked operator)
1428 RCP<const CrsMatrixWrap> crsMat = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(mat);
1429 TEUCHOS_ASSERT(crsMat != Teuchos::null);
1430
1431 // these are the thyra style maps of the matrix
1432 RCP<const Map> trowMap = mat->getRowMap();
1433 RCP<const Map> tcolMap = mat->getColMap();
1434 RCP<const Map> tdomMap = mat->getDomainMap();
1435
1436 // get Xpetra maps
1437 RCP<const Map> xrowMap = getRangeMapExtractor()->getMap(i,false);
1438 RCP<const Map> xdomMap = getDomainMapExtractor()->getMap(j,false);
1439
1440 // generate column map with Xpetra GIDs
1441 // We have to do this separately for each block since the column
1442 // map of each block might be different in the same block column
1444 *tcolMap,
1445 *tdomMap,
1446 *xdomMap);
1447
1448 // jump over empty blocks
1449 if(mat->getLocalNumEntries() == 0) continue;
1450
1451 size_t maxNumEntries = mat->getLocalMaxNumRowEntries();
1452
1453 size_t numEntries;
1454 Array<GO> inds (maxNumEntries);
1455 Array<GO> inds2(maxNumEntries);
1456 Array<SC> vals (maxNumEntries);
1457
1458 // loop over all rows and add entries
1459 for (size_t k = 0; k < mat->getLocalNumRows(); k++) {
1460 GlobalOrdinal rowTGID = trowMap->getGlobalElement(k);
1461 crsMat->getCrsMatrix()->getGlobalRowCopy(rowTGID, inds(), vals(), numEntries);
1462
1463 // create new indices array
1464 for (size_t l = 0; l < numEntries; ++l) {
1465 LocalOrdinal lid = tcolMap->getLocalElement(inds[l]);
1466 inds2[l] = xcolMap->getGlobalElement(lid);
1467 }
1468
1469 GlobalOrdinal rowXGID = xrowMap->getGlobalElement(k);
1470 sparse->insertGlobalValues(
1471 rowXGID, inds2(0, numEntries),
1472 vals(0, numEntries));
1473 }
1474 }
1475 }
1476 }
1477 }
1478
1479 sparse->fillComplete(getFullDomainMap(), getFullRangeMap());
1480
1482 "BlockedCrsMatrix::Merge: Local number of entries of merged matrix does not coincide with local number of entries of blocked operator." );
1483
1485 "BlockedCrsMatrix::Merge: Global number of entries of merged matrix does not coincide with global number of entries of blocked operator." );
1486
1487 return sparse;
1488 }
1490
1494 if (Rows() == 1 && Cols () == 1) {
1495 return getMatrix(0,0)->getLocalMatrixDevice();
1496 }
1497 throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getLocalMatrix(): operation not supported.");
1498 }
1500 typename local_matrix_type::HostMirror getLocalMatrixHost () const {
1501 if (Rows() == 1 && Cols () == 1) {
1502 return getMatrix(0,0)->getLocalMatrixHost();
1503 }
1504 throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getLocalMatrix(): operation not supported.");
1505 }
1506
1507
1508#ifdef HAVE_XPETRA_THYRA
1511 Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(Teuchos::rcpFromRef(*this));
1513
1515 Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<Scalar> >(thOp);
1517 return thbOp;
1518 }
1519#endif
1521 LocalOrdinal GetStorageBlockSize() const {return 1;}
1522
1523
1525 void residual(const MultiVector & X,
1526 const MultiVector & B,
1527 MultiVector & R) const {
1529 R.update(STS::one(),B,STS::zero());
1530 this->apply (X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
1531 }
1532
1533 private:
1534
1537
1539
1549 void Add(const Matrix& A, const Scalar scalarA, Matrix& B, const Scalar scalarB) const {
1550 XPETRA_MONITOR("XpetraBlockedCrsMatrix::Add");
1552 "Matrix A is not completed");
1553 using Teuchos::Array;
1554 using Teuchos::ArrayView;
1555
1556 B.scale(scalarB);
1557
1558 Scalar one = ScalarTraits<SC>::one();
1559 Scalar zero = ScalarTraits<SC>::zero();
1560
1561 if (scalarA == zero)
1562 return;
1563
1564 Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1565 Teuchos::RCP<const CrsMatrixWrap> rcpAwrap = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(rcpA);
1567 "BlockedCrsMatrix::Add: matrix A must be of type CrsMatrixWrap.");
1568 Teuchos::RCP<const CrsMatrix> crsA = rcpAwrap->getCrsMatrix();
1569
1570 size_t maxNumEntries = crsA->getLocalMaxNumRowEntries();
1571
1572 size_t numEntries;
1573 Array<GO> inds(maxNumEntries);
1574 Array<SC> vals(maxNumEntries);
1575
1576 RCP<const Map> rowMap = crsA->getRowMap();
1577 RCP<const Map> colMap = crsA->getColMap();
1578
1579 ArrayView<const GO> rowGIDs = crsA->getRowMap()->getLocalElementList();
1580 for (size_t i = 0; i < crsA->getLocalNumRows(); i++) {
1581 GO row = rowGIDs[i];
1582 crsA->getGlobalRowCopy(row, inds(), vals(), numEntries);
1583
1584 if (scalarA != one)
1585 for (size_t j = 0; j < numEntries; ++j)
1586 vals[j] *= scalarA;
1587
1588 B.insertGlobalValues(row, inds(0, numEntries), vals(0, numEntries)); // insert should be ok, since blocks in BlockedCrsOpeartor do not overlap!
1589 }
1590 }
1591
1593
1594 // Default view is created after fillComplete()
1595 // Because ColMap might not be available before fillComplete().
1597 XPETRA_MONITOR("XpetraBlockedCrsMatrix::CreateDefaultView");
1598
1599 // Create default view
1600 this->defaultViewLabel_ = "point";
1602
1603 // Set current view
1604 this->currentViewLabel_ = this->GetDefaultViewLabel();
1605 }
1606
1607 private:
1611
1612 std::vector<Teuchos::RCP<Matrix> > blocks_;
1613#ifdef HAVE_XPETRA_THYRA
1615#endif
1618
1619};
1620
1621} //namespace Xpetra
1622
1623#define XPETRA_BLOCKEDCRSMATRIX_SHORT
1624#endif /* XPETRA_BLOCKEDCRSMATRIX_HPP */
#define XPETRA_MONITOR(funcName)
#define XPETRA_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
iterator end() const
iterator begin() const
size_type size() const
T * getRawPtr() const
static const EVerbosityLevel verbLevel_default
bool is_null() const
Teuchos::RCP< const MapExtractor > domainmaps_
full domain map together with all partial domain maps
void doImport(const Matrix &source, const Export &exporter, CombineMode CM)
Import (using an Exporter).
void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView< const GlobalOrdinal > &indices, ArrayView< const Scalar > &values) const
Extract a const, non-persisting view of global indices in a specified row of the matrix.
virtual ScalarTraits< Scalar >::magnitudeType getFrobeniusNorm() const
Get Frobenius norm of the matrix.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Teuchos::RCP< Matrix > getInnermostCrsMatrix()
helper routine recursively returns the first inner-most non-null matrix block from a (nested) blocked...
bool isLocallyIndexed() const
If matrix indices of all matrix blocks are in the local range, this function returns true....
Teuchos::RCP< const MapExtractor > rangemaps_
full range map together with all partial domain maps
RCP< const BlockedMap > getBlockedRangeMap() const
Returns the BlockedMap associated with the range of this operator.
void setMatrix(size_t r, size_t c, Teuchos::RCP< Matrix > mat)
set matrix block
RCP< const BlockedMap > getBlockedDomainMap() const
Returns the BlockedMap associated with the domain of this operator.
virtual size_t Rows() const
number of row blocks
void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using global IDs.
size_t getLocalNumRows() const
Returns the number of matrix rows owned on the calling node.
Teuchos::RCP< Matrix > Merge() const
merge BlockedCrsMatrix blocks in a CrsMatrix
BlockedCrsMatrix(const Teuchos::RCP< const BlockedMap > &rangeMaps, const Teuchos::RCP< const BlockedMap > &domainMaps, size_t numEntriesPerRow)
Constructor.
void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > &params=null)
Signal that data entry is complete.
bool bDomainThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...
size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
size_t getLocalNumEntries() const
Returns the local number of entries in this matrix.
RCP< const Map > getRangeMap(size_t i) const
Returns the Map associated with the i'th block range of this operator.
virtual void bgs_apply(const MultiVector &X, MultiVector &Y, size_t row, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=ScalarTraits< Scalar >::one(), Scalar beta=ScalarTraits< Scalar >::zero()) const
Special multiplication routine (for BGS/Jacobi smoother)
void doExport(const Matrix &dest, const Export &exporter, CombineMode CM)
Export (using an Importer).
RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
bool hasCrsGraph() const
Supports the getCrsGraph() call.
std::vector< Teuchos::RCP< Matrix > > blocks_
row major matrix block storage
virtual void getLocalRowCopy(LocalOrdinal LocalRow, const ArrayView< LocalOrdinal > &Indices, const ArrayView< Scalar > &Values, size_t &NumEntries) const
Extract a list of entries in a specified local row of the matrix. Put into storage allocated by calli...
RCP< const CrsGraph > getCrsGraph() const
Returns the CrsGraph associated with this matrix.
void replaceGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Replace matrix entries, using global IDs.
void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const
Extract a const, non-persisting view of local indices in a specified row of the matrix.
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode.
Teuchos::RCP< Matrix > getCrsMatrix() const
return unwrap 1x1 blocked operators
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
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...
RCP< const Map > getFullRangeMap() const
Returns the Map associated with the full range of this operator.
void leftScale(const Vector &x)
Left scale matrix using the given vector entries.
virtual void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to scalar.
RCP< const Map > getRangeMap(size_t i, bool bThyraMode) const
Returns the Map associated with the i'th block range of this operator.
void replaceLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Replace matrix entries, using local IDs.
void doImport(const Matrix &source, const Import &importer, CombineMode CM)
Import.
global_size_t getGlobalNumCols() const
Returns the number of global columns in the matrix.
size_t getLocalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
void residual(const MultiVector &X, const MultiVector &B, MultiVector &R) const
Compute a residual R = B - (*this) * X.
bool is_diagonal_
If we're diagonal, a bunch of the extraction stuff should work.
const Teuchos::RCP< const Map > getMap() const
Implements DistObject interface.
void Add(const Matrix &A, const Scalar scalarA, Matrix &B, const Scalar scalarB) const
Add a Xpetra::CrsMatrix to another: B = B*scalarB + A*scalarA.
bool bRangeThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...
BlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMapExtractor, Teuchos::RCP< const MapExtractor > &domainMapExtractor, size_t numEntriesPerRow)
Constructor.
void scale(const Scalar &alpha)
Scale the current values of a matrix, this = alpha*this.
local_matrix_type getLocalMatrixDevice() const
Access the underlying local Kokkos::CrsMatrix object.
virtual bool haveGlobalConstants() const
Returns true if globalConstants have been computed; false otherwise.
RCP< const Map > getDomainMap(size_t i) const
Returns the Map associated with the i'th block domain of this operator.
LocalOrdinal GetStorageBlockSize() const
Returns the block size of the storage mechanism.
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.
size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
Returns the current number of entries in the specified (locally owned) global row.
CrsMatrix::local_matrix_type local_matrix_type
void insertLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using local IDs.
void rightScale(const Vector &x)
Right scale matrix using the given vector entries.
void setObjectLabel(const std::string &objectLabel)
local_matrix_type::HostMirror getLocalMatrixHost() const
Access the underlying local Kokkos::CrsMatrix object.
void fillComplete(const RCP< ParameterList > &params=null)
Signal that data entry is complete.
void resumeFill(const RCP< ParameterList > &params=null)
size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
void getLocalDiagCopy(Vector &diag) const
Get a copy of the diagonal entries owned by this node, with local row indices.
bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
global_size_t getGlobalNumRows() const
Returns the number of global rows.
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
virtual size_t Cols() const
number of column blocks
RCP< const Map > getDomainMap(size_t i, bool bThyraMode) const
Returns the Map associated with the i'th block domain of this operator.
RCP< const Map > getFullDomainMap() const
Returns the Map associated with the full domain of this operator.
global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
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.
void removeEmptyProcessesInPlace(const Teuchos::RCP< const Map > &newMap)
void doExport(const Matrix &dest, const Import &importer, CombineMode CM)
Export.
std::string description() const
Return a simple one-line description of this object.
KokkosSparse::CrsMatrix< impl_scalar_type, LocalOrdinal, execution_space, void, typename local_graph_type::size_type > local_matrix_type
The specialization of Kokkos::CrsMatrix that represents the part of the sparse matrix on each MPI pro...
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
Exception indicating invalid cast attempted.
Exception throws to report incompatible objects (like maps).
Exception throws when you call an unimplemented method of Xpetra.
Exception throws to report errors in the internal logical of the program.
static Teuchos::RCP< Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map > &fullmap, const std::vector< Teuchos::RCP< const Map > > &maps, bool bThyraMode=false)
Build MapExtrasctor from a given set of partial maps.
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.
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > transformThyra2XpetraGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlInput, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlReferenceInput)
replace set of global ids by new global ids
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Xpetra-specific matrix class.
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
viewLabel_t currentViewLabel_
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode.
viewLabel_t defaultViewLabel_
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
const viewLabel_t & GetDefaultViewLabel() const
virtual void scale(const Scalar &alpha)=0
Scale the current values of a matrix, this = alpha*this.
virtual void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)=0
Insert matrix entries, using global IDs.
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.
virtual size_t getNumVectors() const =0
Number of columns in the multivector.
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.
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.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool is_null(const std::shared_ptr< T > &p)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace
std::string viewLabel_t
size_t global_size_t
Global size_t object.
CombineMode
Xpetra::Combine Mode enumerable type.
static magnitudeType magnitude(T a)