46 #ifndef XPETRA_BLOCKEDCRSMATRIX_HPP 47 #define XPETRA_BLOCKEDCRSMATRIX_HPP 51 #include <Teuchos_SerialDenseMatrix.hpp> 52 #include <Teuchos_Hashtable.hpp> 71 #ifdef HAVE_XPETRA_THYRA 72 #include <Thyra_ProductVectorSpaceBase.hpp> 73 #include <Thyra_VectorSpaceBase.hpp> 74 #include <Thyra_LinearOpBase.hpp> 75 #include <Thyra_BlockedLinearOpBase.hpp> 76 #include <Thyra_PhysicallyBlockedLinearOpBase.hpp> 89 template <
class Scalar,
94 public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
102 #undef XPETRA_BLOCKEDCRSMATRIX_SHORT 118 Teuchos::RCP<const MapExtractor>& domainMaps,
128 for (
size_t r = 0; r <
Rows(); ++r)
129 for (
size_t c = 0; c <
Cols(); ++c)
136 #ifdef HAVE_XPETRA_THYRA 144 BlockedCrsMatrix(
const Teuchos::RCP<
const Thyra::BlockedLinearOpBase<Scalar> >& thyraOp,
const Teuchos::RCP<
const Teuchos::Comm<int> >& comm)
148 const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productRangeSpace = thyraOp->productRange();
149 const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productDomainSpace = thyraOp->productDomain();
151 int numRangeBlocks = productRangeSpace->numBlocks();
152 int numDomainBlocks = productDomainSpace->numBlocks();
155 std::vector<Teuchos::RCP<const Map> > subRangeMaps(numRangeBlocks);
156 for (
size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
157 for (
size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
158 if (thyraOp->blockExists(r,c)) {
160 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c);
161 Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
163 subRangeMaps[r] = xop->getRangeMap();
168 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullRangeMap = mergeMaps(subRangeMaps);
172 bool bRangeUseThyraStyleNumbering =
false;
173 size_t numAllElements = 0;
174 for(
size_t v = 0; v < subRangeMaps.size(); ++v) {
175 numAllElements += subRangeMaps[v]->getGlobalNumElements();
177 if ( fullRangeMap->getGlobalNumElements() != numAllElements) bRangeUseThyraStyleNumbering =
true;
181 std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subDomainMaps(numDomainBlocks);
182 for (
size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
183 for (
size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
184 if (thyraOp->blockExists(r,c)) {
186 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c);
187 Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
189 subDomainMaps[c] = xop->getDomainMap();
194 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullDomainMap = mergeMaps(subDomainMaps);
196 bool bDomainUseThyraStyleNumbering =
false;
198 for(
size_t v = 0; v < subDomainMaps.size(); ++v) {
199 numAllElements += subDomainMaps[v]->getGlobalNumElements();
201 if (fullDomainMap->getGlobalNumElements() != numAllElements) bDomainUseThyraStyleNumbering =
true;
210 for (
size_t r = 0; r <
Rows(); ++r) {
211 for (
size_t c = 0; c <
Cols(); ++c) {
212 if(thyraOp->blockExists(r,c)) {
214 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c);
215 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar> >(const_op);
216 Teuchos::RCP<Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
218 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar,LO,GO,Node> > xwrap =
244 std::vector<GlobalOrdinal> gids;
245 for(
size_t tt = 0; tt<subMaps.size(); ++tt) {
246 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > subMap = subMaps[tt];
247 for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(subMap->getNodeNumElements()); ++l) {
248 GlobalOrdinal gid = subMap->getGlobalElement(l);
257 const GO INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
258 std::sort(gids.begin(), gids.end());
259 gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
260 Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
303 void insertGlobalValues(GlobalOrdinal globalRow,
const ArrayView<const GlobalOrdinal>& cols,
const ArrayView<const Scalar>& vals) {
306 getMatrix(0,0)->insertGlobalValues(globalRow, cols, vals);
323 void insertLocalValues(LocalOrdinal localRow,
const ArrayView<const LocalOrdinal>& cols,
const ArrayView<const Scalar>& vals) {
326 getMatrix(0,0)->insertLocalValues(localRow, cols, vals);
333 XPETRA_MONITOR(
"XpetraBlockedCrsMatrix::removeEmptyProcessesInPlace");
335 getMatrix(0,0)->removeEmptyProcessesInPlace(newMap);
351 const ArrayView<const GlobalOrdinal> &cols,
352 const ArrayView<const Scalar> &vals) {
355 getMatrix(0,0)->replaceGlobalValues(globalRow,cols,vals);
367 const ArrayView<const LocalOrdinal> &cols,
368 const ArrayView<const Scalar> &vals) {
371 getMatrix(0,0)->replaceLocalValues(localRow,cols,vals);
380 for (
size_t row = 0; row <
Rows(); row++) {
381 for (
size_t col = 0; col <
Cols(); col++) {
383 getMatrix(row,col)->setAllToScalar(alpha);
392 for (
size_t row = 0; row <
Rows(); row++) {
393 for (
size_t col = 0; col <
Cols(); col++) {
415 for (
size_t row = 0; row <
Rows(); row++) {
416 for (
size_t col = 0; col <
Cols(); col++) {
429 void fillComplete(
const RCP<const Map>& domainMap,
const RCP<const Map>& rangeMap,
const RCP<ParameterList>& params = null) {
432 getMatrix(0,0)->fillComplete(domainMap, rangeMap, params);
456 for (
size_t r = 0; r <
Rows(); ++r)
457 for (
size_t c = 0; c <
Cols(); ++c) {
459 Teuchos::RCP<Matrix> Ablock =
getMatrix(r,c);
461 if (Ablock != Teuchos::null && !Ablock->isFillComplete())
468 RCP<const Map> rangeMap =
rangemaps_->getFullMap();
469 fullrowmap_ =
MapFactory::Build(rangeMap()->lib(), rangeMap()->getGlobalNumElements(), rangeMap()->getNodeElementList(), rangeMap()->getIndexBase(), rangeMap()->getComm());
472 fullcolmap_ = Teuchos::null;
474 std::vector<GO> colmapentries;
475 for (
size_t c = 0; c <
Cols(); ++c) {
478 for (
size_t r = 0; r <
Rows(); ++r) {
479 Teuchos::RCP<CrsMatrix> Ablock =
getMatrix(r,c);
481 if (Ablock != Teuchos::null) {
482 Teuchos::ArrayView<const GO> colElements = Ablock->getColMap()->getNodeElementList();
483 Teuchos::RCP<const Map> colmap = Ablock->getColMap();
484 copy(colElements.getRawPtr(), colElements.getRawPtr() + colElements.size(), inserter(colset, colset.begin()));
489 colmapentries.reserve(colmapentries.size() + colset.size());
490 copy(colset.begin(), colset.end(), back_inserter(colmapentries));
491 sort(colmapentries.begin(), colmapentries.end());
492 typename std::vector<GO>::iterator gendLocation;
493 gendLocation = std::unique(colmapentries.begin(), colmapentries.end());
494 colmapentries.erase(gendLocation,colmapentries.end());
498 size_t numGlobalElements = 0;
499 Teuchos::reduceAll(*(rangeMap->getComm()), Teuchos::REDUCE_SUM, colmapentries.size(), Teuchos::outArg(numGlobalElements));
502 const Teuchos::ArrayView<const GO> aView = Teuchos::ArrayView<const GO>(colmapentries);
503 fullcolmap_ =
MapFactory::Build(rangeMap->lib(), numGlobalElements, aView, 0, rangeMap->getComm());
516 for (
size_t row = 0; row <
Rows(); row++)
517 for (
size_t col = 0; col <
Cols(); col++)
519 globalNumRows +=
getMatrix(row,col)->getGlobalNumRows();
523 return globalNumRows;
533 for (
size_t col = 0; col <
Cols(); col++)
534 for (
size_t row = 0; row <
Rows(); row++)
536 globalNumCols +=
getMatrix(row,col)->getGlobalNumCols();
540 return globalNumCols;
548 for (
size_t row = 0; row <
Rows(); ++row)
549 for (
size_t col = 0; col <
Cols(); col++)
551 nodeNumRows +=
getMatrix(row,col)->getNodeNumRows();
563 for (
size_t row = 0; row <
Rows(); ++row)
564 for (
size_t col = 0; col <
Cols(); ++col)
566 globalNumEntries +=
getMatrix(row,col)->getGlobalNumEntries();
568 return globalNumEntries;
576 for (
size_t row = 0; row <
Rows(); ++row)
577 for (
size_t col = 0; col <
Cols(); ++col)
579 nodeNumEntries +=
getMatrix(row,col)->getNodeNumEntries();
581 return nodeNumEntries;
587 XPETRA_MONITOR(
"XpetraBlockedCrsMatrix::getNumEntriesInLocalRow");
589 return getMatrix(0,0)->getNumEntriesInLocalRow(localRow);
600 return getMatrix(0,0)->getGlobalNumDiags();
611 return getMatrix(0,0)->getNodeNumDiags();
620 XPETRA_MONITOR(
"XpetraBlockedCrsMatrix::getGlobalMaxNumRowEntries");
623 for (
size_t row = 0; row <
Rows(); row++) {
625 for (
size_t col = 0; col <
Cols(); col++) {
627 globalMaxEntriesBlockRows +=
getMatrix(row,col)->getGlobalMaxNumRowEntries();
630 if(globalMaxEntriesBlockRows > globalMaxEntries)
631 globalMaxEntries = globalMaxEntriesBlockRows;
633 return globalMaxEntries;
640 XPETRA_MONITOR(
"XpetraBlockedCrsMatrix::getNodeMaxNumRowEntries");
641 size_t localMaxEntries = 0;
643 for (
size_t row = 0; row <
Rows(); row++) {
644 size_t localMaxEntriesBlockRows = 0;
645 for (
size_t col = 0; col <
Cols(); col++) {
647 localMaxEntriesBlockRows +=
getMatrix(row,col)->getNodeMaxNumRowEntries();
650 if(localMaxEntriesBlockRows > localMaxEntries)
651 localMaxEntries = localMaxEntriesBlockRows;
653 return localMaxEntries;
662 for (
size_t i = 0; i <
blocks_.size(); ++i)
674 for (
size_t i = 0; i <
blocks_.size(); i++)
683 for (
size_t i = 0; i <
blocks_.size(); i++)
705 const ArrayView<LocalOrdinal>& Indices,
706 const ArrayView<Scalar>& Values,
707 size_t &NumEntries)
const {
710 getMatrix(0,0)->getLocalRowCopy(LocalRow, Indices, Values, NumEntries);
726 void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values)
const {
729 getMatrix(0,0)->getGlobalRowView(GlobalRow, indices, values);
745 void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values)
const {
748 getMatrix(0,0)->getLocalRowView(LocalRow, indices, values);
761 "BlockedCrsMatrix::getLocalDiagCopy(): the map of the vector diag is not compatible with the full map of the blocked operator." );
764 "BlockedCrsMatrix::getLocalDiagCopy(): you cannot extract the diagonal of a "<<
Rows() <<
"x"<<
Cols()<<
" blocked matrix." );
766 for (
size_t row = 0; row <
Rows(); ++row) {
770 RCP<Vector> dd =
rangemaps_->getVector(row, bThyraMode);
771 getMatrix(row,row)->getLocalDiagCopy(*dd);
772 rangemaps_->InsertVector(*dd,row,diag,bThyraMode);
781 "BlockedCrsMatrix::leftScale(): the map of the vector x is not compatible with the full map of the blocked operator." );
783 RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
785 for (
size_t row = 0; row <
Rows(); ++row) {
786 for (
size_t col = 0; col <
Cols(); ++col) {
790 RCP<Vector> xx =
rangemaps_->ExtractVector(rcpx,row,bThyraMode);
801 "BlockedCrsMatrix::rightScale(): the map of the vector x is not compatible with the full map of the blocked operator." );
803 RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
805 for (
size_t col = 0; col <
Cols(); ++col) {
806 for (
size_t row = 0; row <
Rows(); ++row) {
810 RCP<Vector> xx =
domainmaps_->ExtractVector(rcpx,col,bThyraMode);
821 typename ScalarTraits<Scalar>::magnitudeType ret = Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero());
822 for (
size_t col = 0; col <
Cols(); ++col) {
823 for (
size_t row = 0; row <
Rows(); ++row) {
825 typename ScalarTraits<Scalar>::magnitudeType n =
getMatrix(row,col)->getFrobeniusNorm();
830 return Teuchos::ScalarTraits< typename ScalarTraits<Scalar>::magnitudeType >::squareroot(ret);
868 Teuchos::ETransp mode = Teuchos::NO_TRANS,
869 Scalar alpha = ScalarTraits<Scalar>::one(),
870 Scalar beta = ScalarTraits<Scalar>::zero())
const 876 "apply() only supports the following modes: NO_TRANS and TRANS." );
879 RCP<const MultiVector> refX = rcpFromRef(X);
880 RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(refX);
881 RCP<MultiVector> tmpY = rcpFromRef(Y);
884 bool bBlockedX = (refbX != Teuchos::null) ?
true :
false;
885 bool bBlockedY = (tmpbY != Teuchos::null) ?
true :
false;
892 SC one = ScalarTraits<SC>::one();
894 if (mode == Teuchos::NO_TRANS) {
896 for (
size_t row = 0; row <
Rows(); row++) {
898 for (
size_t col = 0; col <
Cols(); col++) {
901 RCP<Matrix> Ablock =
getMatrix(row, col);
903 if (Ablock.is_null())
908 bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(Ablock) == Teuchos::null ?
false :
true;
911 RCP<const MultiVector> Xblock = Teuchos::null;
919 Ablock->apply(*Xblock, *tmpYblock);
924 tmpYblock->replaceMap(
rangemaps_->getMap(row,
true));
926 Yblock->update(one, *tmpYblock, one);
932 }
else if (mode == Teuchos::TRANS) {
934 for (
size_t col = 0; col <
Cols(); col++) {
937 for (
size_t row = 0; row<
Rows(); row++) {
938 RCP<Matrix> Ablock =
getMatrix(row, col);
940 if (Ablock.is_null())
945 bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(Ablock) == Teuchos::null ?
false :
true;
947 RCP<const MultiVector> Xblock = Teuchos::null;
953 Ablock->apply(*Xblock, *tmpYblock, Teuchos::TRANS);
958 tmpYblock->replaceMap(
domainmaps_->getMap(col,
true));
961 Yblock->update(one, *tmpYblock, one);
967 if(bBlockedY) Y.
update(alpha, *tmpbY, beta);
968 else Y.
update(alpha, *tmpY, beta);
1001 const Teuchos::RCP< const Map >
getMap()
const {
1013 getMatrix(0,0)->doImport(source, importer, CM);
1023 getMatrix(0,0)->doExport(dest, importer, CM);
1033 getMatrix(0,0)->doImport(source, exporter, CM);
1043 getMatrix(0,0)->doExport(dest, exporter, CM);
1055 std::string
description()
const {
return "Xpetra_BlockedCrsMatrix.description()"; }
1058 void describe(Teuchos::FancyOStream& out,
const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default)
const {
1059 out <<
"Xpetra::BlockedCrsMatrix: " <<
Rows() <<
" x " <<
Cols() << std::endl;
1062 out <<
"BlockMatrix is fillComplete" << std::endl;
1075 out <<
"BlockMatrix is NOT fillComplete" << std::endl;
1078 for (
size_t r = 0; r <
Rows(); ++r)
1079 for (
size_t c = 0; c <
Cols(); ++c) {
1081 out <<
"Block(" << r <<
"," << c <<
")" << std::endl;
1082 getMatrix(r,c)->describe(out,verbLevel);
1083 }
else out <<
"Block(" << r <<
"," << c <<
") = null" << std::endl;
1110 TEUCHOS_TEST_FOR_EXCEPTION(
Rows()!=1, std::out_of_range,
"Can only unwrap a 1x1 blocked matrix. The matrix has " <<
Rows() <<
" block rows, though.");
1111 TEUCHOS_TEST_FOR_EXCEPTION(
Cols()!=1, std::out_of_range,
"Can only unwrap a 1x1 blocked matrix. The matrix has " <<
Cols() <<
" block columns, though.");
1114 RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(mat);
1115 if (bmat == Teuchos::null)
return mat;
1121 size_t row =
Rows()+1, col =
Cols()+1;
1122 for (
size_t r = 0; r <
Rows(); ++r)
1123 for(
size_t c = 0; c <
Cols(); ++c)
1129 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.")
1131 RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(mm);
1132 if (bmat == Teuchos::null)
return mm;
1139 TEUCHOS_TEST_FOR_EXCEPTION(r >
Rows(), std::out_of_range,
"Error, r = " <<
Rows() <<
" is too big");
1140 TEUCHOS_TEST_FOR_EXCEPTION(c >
Cols(), std::out_of_range,
"Error, c = " <<
Cols() <<
" is too big");
1151 void setMatrix(
size_t r,
size_t c, Teuchos::RCP<Matrix> mat) {
1155 TEUCHOS_TEST_FOR_EXCEPTION(r >
Rows(), std::out_of_range,
"Error, r = " <<
Rows() <<
" is too big");
1156 TEUCHOS_TEST_FOR_EXCEPTION(c >
Cols(), std::out_of_range,
"Error, c = " <<
Cols() <<
" is too big");
1168 using Teuchos::rcp_dynamic_cast;
1169 Scalar one = ScalarTraits<SC>::one();
1172 "BlockedCrsMatrix::Merge: only implemented for Xpetra-style or Thyra-style numbering. No mixup allowed!" );
1175 "BlockedCrsMatrix::Merge: BlockMatrix must be fill-completed." );
1181 for (
size_t i = 0; i <
Rows(); i++) {
1182 for (
size_t j = 0; j <
Cols(); j++) {
1187 RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(mat);
1188 if (bMat != Teuchos::null) mat = bMat->
Merge();
1192 "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1195 if(mat->getNodeNumEntries() == 0)
continue;
1197 this->
Add(*mat, one, *sparse, one);
1203 for (
size_t i = 0; i <
Rows(); i++) {
1204 for (
size_t j = 0; j <
Cols(); j++) {
1208 RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(mat);
1209 if (bMat != Teuchos::null) mat = bMat->
Merge();
1213 "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1216 RCP<const CrsMatrixWrap> crsMat = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(mat);
1217 TEUCHOS_ASSERT(crsMat != Teuchos::null);
1220 RCP<const Map> trowMap = mat->getRowMap();
1221 RCP<const Map> tcolMap = mat->getColMap();
1222 RCP<const Map> tdomMap = mat->getDomainMap();
1237 if(mat->getNodeNumEntries() == 0)
continue;
1239 size_t maxNumEntries = mat->getNodeMaxNumRowEntries();
1242 Array<GO> inds (maxNumEntries);
1243 Array<GO> inds2(maxNumEntries);
1244 Array<SC> vals (maxNumEntries);
1247 for (
size_t k = 0; k < mat->getNodeNumRows(); k++) {
1248 GlobalOrdinal rowTGID = trowMap->getGlobalElement(k);
1249 crsMat->getCrsMatrix()->getGlobalRowCopy(rowTGID, inds(), vals(), numEntries);
1252 for (
size_t l = 0; l < numEntries; ++l) {
1253 LocalOrdinal lid = tcolMap->getLocalElement(inds[l]);
1254 inds2[l] = xcolMap->getGlobalElement(lid);
1257 GlobalOrdinal rowXGID = xrowMap->getGlobalElement(k);
1258 sparse->insertGlobalValues(
1259 rowXGID, inds2(0, numEntries),
1260 vals(0, numEntries));
1270 "BlockedCrsMatrix::Merge: Local number of entries of merged matrix does not coincide with local number of entries of blocked operator." );
1273 "BlockedCrsMatrix::Merge: Global number of entries of merged matrix does not coincide with global number of entries of blocked operator." );
1279 #ifdef HAVE_XPETRA_KOKKOS_REFACTOR 1280 typedef typename CrsMatrix::local_matrix_type local_matrix_type;
1283 local_matrix_type getLocalMatrix ()
const {
1285 return getMatrix(0,0)->getLocalMatrix();
1291 #ifdef HAVE_XPETRA_THYRA 1292 Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > getThyraOperator() {
1293 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thOp =
1294 Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(Teuchos::rcpFromRef(*
this));
1295 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thOp));
1297 Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > thbOp =
1298 Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<Scalar> >(thOp);
1299 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thbOp));
1323 "Matrix A is not completed");
1324 using Teuchos::Array;
1325 using Teuchos::ArrayView;
1329 Scalar one = ScalarTraits<SC>::one();
1330 Scalar zero = ScalarTraits<SC>::zero();
1332 if (scalarA == zero)
1335 Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1336 Teuchos::RCP<const CrsMatrixWrap> rcpAwrap = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(rcpA);
1338 "BlockedCrsMatrix::Add: matrix A must be of type CrsMatrixWrap.");
1339 Teuchos::RCP<const CrsMatrix> crsA = rcpAwrap->getCrsMatrix();
1341 size_t maxNumEntries = crsA->getNodeMaxNumRowEntries();
1344 Array<GO> inds(maxNumEntries);
1345 Array<SC> vals(maxNumEntries);
1347 RCP<const Map> rowMap = crsA->getRowMap();
1348 RCP<const Map> colMap = crsA->getColMap();
1350 ArrayView<const GO> rowGIDs = crsA->getRowMap()->getNodeElementList();
1351 for (
size_t i = 0; i < crsA->getNodeNumRows(); i++) {
1352 GO row = rowGIDs[i];
1353 crsA->getGlobalRowCopy(row, inds(), vals(), numEntries);
1356 for (
size_t j = 0; j < numEntries; ++j)
1383 #ifdef HAVE_XPETRA_THYRA 1384 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > thyraOp_;
1393 #define XPETRA_BLOCKEDCRSMATRIX_SHORT
Teuchos::RCP< Matrix > Merge() const
merge BlockedCrsMatrix blocks in a CrsMatrix
size_t getNodeMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
RCP< const CrsGraph > getCrsGraph() const
Returns the CrsGraph associated with this matrix.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
const Teuchos::RCP< const Map > getMap() const
Implements DistObject interface.
bool isLocallyIndexed() const
If matrix indices of all matrix blocks are in the local range, this function returns true...
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
RCP< const Map > getRangeMap(size_t i, bool bThyraMode) const
Returns the Map associated with the i'th block range of this operator.
void replaceGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Replace matrix entries, using global IDs.
void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > ¶ms=null)
Signal that data entry is complete.
virtual void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to scalar.
RCP< const Map > getDomainMap(size_t i) const
Returns the Map associated with the i'th block domain of this operator.
global_size_t getGlobalNumDiags() const
Returns the number of global diagonal entries, based on global row/column index comparisons.
Teuchos::RCP< const MapExtractor > rangemaps_
void doImport(const Matrix &source, const Import &importer, CombineMode CM)
Import.
Exception throws to report errors in the internal logical of the program.
virtual ScalarTraits< Scalar >::magnitudeType getFrobeniusNorm() const
Get Frobenius norm of the 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.
virtual ~BlockedCrsMatrix()
Destructor.
Teuchos::RCP< Matrix > getCrsMatrix() const
return unwrap 1x1 blocked operators
void doExport(const Matrix &dest, const Import &importer, CombineMode CM)
Export.
LocalOrdinal local_ordinal_type
size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
global_size_t getGlobalNumCols() const
Returns the number of global columns in the matrix.
void resumeFill(const RCP< ParameterList > ¶ms=null)
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.
size_t getNodeNumDiags() const
Returns the number of local diagonal entries, based on global row/column index comparisons.
void rightScale(const Vector &x)
Right scale matrix using the given vector entries.
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Exception indicating invalid cast attempted.
void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using global IDs.
viewLabel_t defaultViewLabel_
virtual void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)=0
Insert matrix entries, using global IDs.
size_t getNodeNumRows() const
Returns the number of matrix rows owned on the calling node.
GlobalOrdinal global_ordinal_type
void doExport(const Matrix &dest, const Export &exporter, CombineMode CM)
Export (using an Importer).
void leftScale(const Vector &x)
Left scale matrix using the given vector entries.
Teuchos::RCP< Matrix > getInnermostCrsMatrix()
void insertLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using local IDs.
bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
RCP< const Map > getDomainMap(size_t i, bool bThyraMode) const
Returns the Map associated with the i'th block domain of this operator.
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
const viewLabel_t & GetDefaultViewLabel() const
void getLocalDiagCopy(Vector &diag) const
Get a copy of the diagonal entries owned by this node, with local row indices.
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.
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
void scale(const Scalar &alpha)
Scale the current values of a matrix, this = alpha*this.
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
viewLabel_t currentViewLabel_
void doImport(const Matrix &source, const Export &exporter, CombineMode CM)
Import (using an Exporter).
RCP< const Map > getRangeMap(size_t i) const
Returns the Map associated with the i'th block range of this operator.
size_t global_size_t
Global size_t object.
global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
std::vector< Teuchos::RCP< Matrix > > blocks_
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...
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.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
#define XPETRA_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setMatrix(size_t r, size_t c, Teuchos::RCP< Matrix > mat)
set matrix block
bool bDomainThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...
BlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMaps, Teuchos::RCP< const MapExtractor > &domainMaps, size_t npr, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor.
Exception throws to report incompatible objects (like maps).
size_t getNodeNumEntries() const
Returns the local number of entries in this matrix.
RCP< const Map > getRangeMap() const
Returns the Map associated with the full range of this operator.
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.
Concrete implementation of Xpetra::Matrix.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
std::string description() const
Return a simple one-line description of this object.
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.
bool IsView(const viewLabel_t viewLabel) const
CombineMode
Xpetra::Combine Mode enumerable type.
global_size_t getGlobalNumRows() const
Returns the number of global rows.
#define XPETRA_MONITOR(funcName)
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor specifying the number of non-zeros for all rows.
virtual void scale(const Scalar &alpha)=0
Scale the current values of a matrix, this = alpha*this.
void fillComplete(const RCP< ParameterList > ¶ms=null)
Signal that data entry is complete.
virtual size_t Rows() const
number of row blocks
virtual size_t Cols() const
number of column blocks
virtual size_t getNumVectors() const =0
Number of columns in the multivector.
RCP< const Map > getDomainMap() const
Returns the Map associated with the full domain of this operator.
void removeEmptyProcessesInPlace(const Teuchos::RCP< const Map > &newMap)
size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
void replaceLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Replace matrix entries, using local IDs.
Xpetra-specific matrix class.
Teuchos::RCP< const MapExtractor > domainmaps_
bool bRangeThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...
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.