48 #ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_ 49 #define PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_ 62 #ifdef HAVE_XPETRA_EPETRA 66 #ifdef HAVE_XPETRA_EPETRAEXT 67 #include <EpetraExt_MatrixMatrix.h> 68 #include <EpetraExt_RowMatrixOut.h> 69 #include <Epetra_RowMatrixTransposer.h> 70 #endif // HAVE_XPETRA_EPETRAEXT 72 #ifdef HAVE_XPETRA_TPETRA 73 #include <TpetraExt_MatrixMatrix.hpp> 74 #include <MatrixMarket_Tpetra.hpp> 78 #endif // HAVE_XPETRA_TPETRA 88 template <
class Scalar,
89 class LocalOrdinal = int,
90 class GlobalOrdinal = LocalOrdinal,
97 #ifdef HAVE_XPETRA_EPETRA 100 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(Op);
102 "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
104 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
105 const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<
const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
107 "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
109 return tmp_ECrsMtx->getEpetra_CrsMatrix();
113 RCP<Epetra_CrsMatrix> A;
115 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(Op);
117 "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
119 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
120 const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<
const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
121 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null,
Exceptions::BadCast,
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
123 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
131 const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<
const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
133 "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
135 return *tmp_ECrsMtx->getEpetra_CrsMatrix();
143 RCP<Epetra_CrsMatrix> A;
148 const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<
const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
149 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null,
Xpetra::Exceptions::BadCast,
"Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
151 return *Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
157 #endif // HAVE_XPETRA_EPETRA 159 #ifdef HAVE_XPETRA_TPETRA 160 static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> >
Op2TpetraCrs(RCP<Matrix> Op) {
162 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(Op);
163 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null,
Xpetra::Exceptions::BadCast,
"Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
165 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
167 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null,
Xpetra::Exceptions::BadCast,
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
169 return tmp_ECrsMtx->getTpetra_CrsMatrix();
174 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(Op);
175 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null,
Xpetra::Exceptions::BadCast,
"Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
176 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
178 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null,
Xpetra::Exceptions::BadCast,
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
180 return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
189 TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null,
Xpetra::Exceptions::BadCast,
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
191 return *tmp_TCrsMtx->getTpetra_CrsMatrix();
204 TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null,
Xpetra::Exceptions::BadCast,
"Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
206 return *Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_TCrsMtx->getTpetra_CrsMatrix());
212 #endif // HAVE_XPETRA_TPETRA 216 template <
class Scalar,
218 class GlobalOrdinal ,
221 #undef XPETRA_MATRIXMATRIX_SHORT 251 const Matrix& B,
bool transposeB,
253 bool call_FillComplete_on_result =
true,
254 bool doOptimizeStorage =
true,
255 const std::string & label = std::string(),
256 const RCP<ParameterList>& params = null) {
258 TEUCHOS_TEST_FOR_EXCEPTION(transposeA ==
false && C.
getRowMap()->isSameAs(*A.
getRowMap()) ==
false,
260 TEUCHOS_TEST_FOR_EXCEPTION(transposeA ==
true && C.
getRowMap()->isSameAs(*A.
getDomainMap()) ==
false,
266 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
269 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 270 throw(
Xpetra::Exceptions::RuntimeError(
"Xpetra::MatrixMatrix::Multiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
276 #ifdef HAVE_XPETRA_TPETRA 283 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
289 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
290 RCP<Teuchos::ParameterList> fillParams = rcp(
new Teuchos::ParameterList());
291 fillParams->set(
"Optimize Storage", doOptimizeStorage);
298 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<
Matrix>(Teuchos::rcpFromRef(A));
299 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<
Matrix>(Teuchos::rcpFromRef(B));
300 C.
CreateView(
"stridedMaps", rcpA, transposeA, rcpB, transposeB);
325 static RCP<Matrix>
Multiply(
const Matrix& A,
bool transposeA,
const Matrix& B,
bool transposeB, RCP<Matrix> C_in,
326 Teuchos::FancyOStream& fos,
327 bool doFillComplete =
true,
328 bool doOptimizeStorage =
true,
329 const std::string & label = std::string(),
330 const RCP<ParameterList>& params = null) {
336 RCP<Matrix> C = C_in;
338 if (C == Teuchos::null) {
339 double nnzPerRow = Teuchos::as<double>(0);
347 nnzPerRow *= nnzPerRow;
350 if (totalNnz < minNnz)
354 fos <<
"Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
362 fos <<
"Reuse C pattern" << std::endl;
365 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params);
380 static RCP<Matrix>
Multiply(
const Matrix& A,
bool transposeA,
const Matrix& B,
bool transposeB, Teuchos::FancyOStream &fos,
381 bool callFillCompleteOnResult =
true,
bool doOptimizeStorage =
true,
const std::string& label = std::string(),
382 const RCP<ParameterList>& params = null) {
383 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
386 #ifdef HAVE_XPETRA_EPETRAEXT 389 const Epetra_CrsMatrix& epB,
390 Teuchos::FancyOStream& fos) {
391 throw(
Xpetra::Exceptions::RuntimeError(
"MLTwoMatrixMultiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
392 return Teuchos::null;
394 #endif //ifdef HAVE_XPETRA_EPETRAEXT 408 Teuchos::FancyOStream& fos,
409 bool doFillComplete =
true,
410 bool doOptimizeStorage =
true) {
412 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
421 RCP<BlockedCrsMatrix> C = rcp(
new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 ));
423 for (
size_t i = 0; i < A.
Rows(); ++i) {
424 for (
size_t j = 0; j < B.
Cols(); ++j) {
427 for (
size_t l = 0; l < B.
Rows(); ++l) {
431 if (crmat1.is_null() || crmat2.is_null())
433 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
436 RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<
CrsMatrixWrap>(crmat1);
437 RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<
CrsMatrixWrap>(crmat2);
438 TEUCHOS_TEST_FOR_EXCEPTION((crop1==Teuchos::null) != (crop2==Teuchos::null),
Xpetra::Exceptions::RuntimeError,
"A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
441 RCP<Matrix> temp = Teuchos::null;
443 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
444 temp =
Multiply (*crop1,
false, *crop2,
false, fos);
446 RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(crmat1);
447 RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(crmat2);
448 TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==
true,
Xpetra::Exceptions::BadCast,
"A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
449 TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==
true,
Xpetra::Exceptions::BadCast,
"B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
450 TEUCHOS_TEST_FOR_EXCEPTION(bop1->Cols()!=bop2->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << bop1->Cols() <<
" columns and B has " << bop2->Rows() <<
" rows. Matrices are not compatible! (TwoMatrixMultiplyBlock)");
451 TEUCHOS_TEST_FOR_EXCEPTION(bop1->getDomainMap()->isSameAs(*(bop2->getRangeMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of A is not the same as range map of B. Matrices are not compatible! (TwoMatrixMultiplyBlock)");
454 temp =
TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
456 RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(temp);
457 TEUCHOS_TEST_FOR_EXCEPTION(btemp->Rows()!=bop1->Rows(),
Xpetra::Exceptions::RuntimeError,
"Number of block rows of local blocked operator is " << btemp->Rows() <<
" but should be " << bop1->Rows() <<
". (TwoMatrixMultiplyBlock)");
458 TEUCHOS_TEST_FOR_EXCEPTION(btemp->Cols()!=bop2->Cols(),
Xpetra::Exceptions::RuntimeError,
"Number of block cols of local blocked operator is " << btemp->Cols() <<
" but should be " << bop2->Cols() <<
". (TwoMatrixMultiplyBlock)");
459 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getRangeMapExtractor()->getFullMap()->isSameAs(*(bop1->getRangeMapExtractor()->getFullMap())) ==
false,
Xpetra::Exceptions::RuntimeError,
"Range map of local blocked operator should be same as first operator. (TwoMatrixMultiplyBlock)");
460 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getDomainMapExtractor()->getFullMap()->isSameAs(*(bop2->getDomainMapExtractor()->getFullMap())) ==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of local blocked operator should be same as second operator. (TwoMatrixMultiplyBlock)");
461 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getRangeMapExtractor()->getThyraMode() != bop1->getRangeMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Thyra mode of local range map extractor incompatible with range map extractor of A (TwoMatrixMultiplyBlock)");
462 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getDomainMapExtractor()->getThyraMode() != bop2->getDomainMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Thyra mode of local domain map extractor incompatible with domain map extractor of B (TwoMatrixMultiplyBlock)");
467 RCP<Matrix> addRes = null;
476 if (!Cij.is_null()) {
477 if (Cij->isFillComplete())
480 C->setMatrix(i, j, Cij);
482 C->setMatrix(i, j, Teuchos::null);
510 throw Exceptions::RuntimeError(
"TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
512 #ifdef HAVE_XPETRA_TPETRA 516 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
537 const Matrix& B,
bool transposeB,
const SC& beta,
538 RCP<Matrix>& C, Teuchos::FancyOStream &fos,
bool AHasFixedNnzPerRow =
false) {
540 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
541 RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
542 RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(rcpA);
543 RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(rcpB);
549 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
555 if (C == Teuchos::null) {
558 size_t numLocalRows = 0;
564 if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
567 Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
569 if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
570 for (
size_t i = 0; i < numLocalRows; ++i)
574 for (
size_t i = 0; i < numLocalRows; ++i)
578 fos <<
"MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)" 579 <<
", using static profiling" << std::endl;
585 LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
592 fos <<
"nnzPerRowInA = " << nnzPerRowInA <<
", nnzPerRowInB = " << nnzPerRowInB << std::endl;
593 fos <<
"MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
594 <<
", max possible nnz per row in sum = " << maxPossible
600 fos <<
"MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
606 #ifdef HAVE_XPETRA_TPETRA 611 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
617 if (A.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(A));
618 if (B.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(B));
622 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
623 RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
624 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
630 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
631 RCP<Matrix> Cij = Teuchos::null;
632 if(rcpA != Teuchos::null &&
633 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
636 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
637 Cij, fos, AHasFixedNnzPerRow);
638 }
else if (rcpA == Teuchos::null &&
639 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
640 Cij = rcpBopB->getMatrix(i,j);
641 }
else if (rcpA != Teuchos::null &&
642 rcpBopB->getMatrix(i,j)==Teuchos::null) {
643 Cij = Teuchos::rcp_const_cast<
Matrix>(rcpA);
648 if (!Cij.is_null()) {
649 if (Cij->isFillComplete())
652 bC->setMatrix(i, j, Cij);
654 bC->setMatrix(i, j, Teuchos::null);
659 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
660 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
661 RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
666 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
667 RCP<Matrix> Cij = Teuchos::null;
668 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
669 rcpB!=Teuchos::null) {
671 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
672 *rcpB, transposeB, beta,
673 Cij, fos, AHasFixedNnzPerRow);
674 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
675 rcpB!=Teuchos::null) {
676 Cij = Teuchos::rcp_const_cast<
Matrix>(rcpB);
677 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
678 rcpB==Teuchos::null) {
679 Cij = rcpBopA->getMatrix(i,j);
684 if (!Cij.is_null()) {
685 if (Cij->isFillComplete())
688 bC->setMatrix(i, j, Cij);
690 bC->setMatrix(i, j, Teuchos::null);
700 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->Rows()!=rcpBopB->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << rcpBopA->Rows() <<
" rows and B has " << rcpBopA->Rows() <<
" rows. Matrices are not compatible! (TwoMatrixAdd)");
701 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->Rows()!=rcpBopB->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << rcpBopA->Cols() <<
" cols and B has " << rcpBopA->Cols() <<
" cols. Matrices are not compatible! (TwoMatrixAdd)");
702 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getRangeMap()->isSameAs(*(rcpBopB->getRangeMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Range map of A is not the same as range map of B. Matrices are not compatible! (TwoMatrixAdd)");
703 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getDomainMap()->isSameAs(*(rcpBopB->getDomainMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of A is not the same as domain map of B. Matrices are not compatible! (TwoMatrixAdd)");
704 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getRangeMapExtractor()->getThyraMode() != rcpBopB->getRangeMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Different Thyra/Xpetra style gids in RangeMapExtractor of A and B. Matrices are not compatible! (TwoMatrixAdd)");
705 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getDomainMapExtractor()->getThyraMode() != rcpBopB->getDomainMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Different Thyra/Xpetra style gids in DomainMapExtractor of A and B. Matrices are not compatible! (TwoMatrixAdd)");
707 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
708 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
712 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
713 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
715 RCP<Matrix> Cij = Teuchos::null;
716 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
717 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
719 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
720 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
721 Cij, fos, AHasFixedNnzPerRow);
722 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
723 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
724 Cij = rcpBopB->getMatrix(i,j);
725 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
726 rcpBopB->getMatrix(i,j)==Teuchos::null) {
727 Cij = rcpBopA->getMatrix(i,j);
732 if (!Cij.is_null()) {
733 if (Cij->isFillComplete())
736 bC->setMatrix(i, j, Cij);
738 bC->setMatrix(i, j, Teuchos::null);
750 #ifdef HAVE_XPETRA_EPETRA 787 const Matrix& B,
bool transposeB,
789 bool call_FillComplete_on_result =
true,
790 bool doOptimizeStorage =
true,
791 const std::string & label = std::string(),
792 const RCP<ParameterList>& params = null) {
793 TEUCHOS_TEST_FOR_EXCEPTION(transposeA ==
false && C.
getRowMap()->isSameAs(*A.
getRowMap()) ==
false,
795 TEUCHOS_TEST_FOR_EXCEPTION(transposeA ==
true && C.
getRowMap()->isSameAs(*A.
getDomainMap()) ==
false,
801 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
804 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 809 int i = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, epC, haveMultiplyDoFillComplete);
810 if (haveMultiplyDoFillComplete) {
821 std::ostringstream buf;
823 std::string msg =
"EpetraExt::MatrixMatrix::Multiply return value of " + buf.str();
831 #ifdef HAVE_XPETRA_TPETRA 832 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \ 833 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT)))) 842 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
849 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
850 RCP<Teuchos::ParameterList> fillParams = rcp(
new Teuchos::ParameterList());
851 fillParams->set(
"Optimize Storage",doOptimizeStorage);
858 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<
Matrix>(Teuchos::rcpFromRef(A));
859 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<
Matrix>(Teuchos::rcpFromRef(B));
860 C.
CreateView(
"stridedMaps", rcpA, transposeA, rcpB, transposeB);
886 const Matrix& B,
bool transposeB,
888 Teuchos::FancyOStream& fos,
889 bool doFillComplete =
true,
890 bool doOptimizeStorage =
true,
891 const std::string & label = std::string(),
892 const RCP<ParameterList>& params = null) {
899 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) && defined(HAVE_XPETRA_ML_MMM) 905 RCP<Matrix> C = Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO> (epC);
906 if (doFillComplete) {
907 RCP<Teuchos::ParameterList> fillParams = rcp(
new Teuchos::ParameterList());
908 fillParams->set(
"Optimize Storage", doOptimizeStorage);
915 C->CreateView(
"stridedMaps", rcpFromRef(A), transposeA, rcpFromRef(B), transposeB);
919 #endif // EPETRA + EPETRAEXT + ML 922 RCP<Matrix> C = C_in;
924 if (C == Teuchos::null) {
925 double nnzPerRow = Teuchos::as<double>(0);
933 nnzPerRow *= nnzPerRow;
936 if (totalNnz < minNnz)
940 fos <<
"Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
948 fos <<
"Reuse C pattern" << std::endl;
951 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params);
967 const Matrix& B,
bool transposeB,
968 Teuchos::FancyOStream &fos,
969 bool callFillCompleteOnResult =
true,
970 bool doOptimizeStorage =
true,
971 const std::string & label = std::string(),
972 const RCP<ParameterList>& params = null) {
973 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
976 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 979 const Epetra_CrsMatrix& epB,
980 Teuchos::FancyOStream& fos) {
981 #if defined(HAVE_XPETRA_ML_MMM) // Note: this is currently not supported 983 ML_Comm_Create(&comm);
984 fos <<
"****** USING ML's MATRIX MATRIX MULTIPLY ******" << std::endl;
987 const Epetra_MpiComm * Mcomm =
dynamic_cast<const Epetra_MpiComm*
>(&(epA.Comm()));
989 ML_Comm_Set_UsrComm(comm,Mcomm->GetMpiComm());
992 EpetraExt::CrsMatrix_SolverMap Atransform;
993 EpetraExt::CrsMatrix_SolverMap Btransform;
994 const Epetra_CrsMatrix& A = Atransform(const_cast<Epetra_CrsMatrix&>(epA));
995 const Epetra_CrsMatrix& B = Btransform(const_cast<Epetra_CrsMatrix&>(epB));
1001 ML_Operator* ml_As = ML_Operator_Create(comm);
1002 ML_Operator* ml_Bs = ML_Operator_Create(comm);
1003 ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&A),ml_As);
1004 ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&B),ml_Bs);
1005 ML_Operator* ml_AtimesB = ML_Operator_Create(comm);
1007 Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer(
"ML_2matmult kernel"));
1008 ML_2matmult(ml_As,ml_Bs,ml_AtimesB,ML_CSR_MATRIX);
1010 ML_Operator_Destroy(&ml_As);
1011 ML_Operator_Destroy(&ml_Bs);
1016 int N_local = ml_AtimesB->invec_leng;
1017 ML_CommInfoOP* getrow_comm = ml_AtimesB->getrow->pre_comm;
1019 ML_Comm* comm_AB = ml_AtimesB->comm;
1020 if (N_local != B.DomainMap().NumMyElements())
1025 for (
int i = 0; i < getrow_comm->N_neighbors; i++) {
1026 N_rcvd += (getrow_comm->neighbors)[i].N_rcv;
1027 N_send += (getrow_comm->neighbors)[i].N_send;
1028 if ( ((getrow_comm->neighbors)[i].N_rcv != 0) &&
1029 ((getrow_comm->neighbors)[i].rcv_list != NULL) ) flag = 1;
1033 std::vector<double> dtemp(N_local + N_rcvd + 1);
1034 std::vector<int> cmap (N_local + N_rcvd + 1);
1035 for (
int i = 0; i < N_local; ++i) {
1036 cmap[i] = B.DomainMap().GID(i);
1037 dtemp[i] = (double) cmap[i];
1039 ML_cheap_exchange_bdry(&dtemp[0],getrow_comm,N_local,N_send,comm_AB);
1041 int count = N_local;
1042 const int neighbors = getrow_comm->N_neighbors;
1043 for (
int i = 0; i < neighbors; i++) {
1044 const int nrcv = getrow_comm->neighbors[i].N_rcv;
1045 for (
int j = 0; j < nrcv; j++)
1046 cmap[getrow_comm->neighbors[i].rcv_list[j]] = (
int) dtemp[count++];
1049 for (
int i = 0; i < N_local+N_rcvd; ++i)
1050 cmap[i] = (
int)dtemp[i];
1055 Epetra_Map gcmap(-1, N_local+N_rcvd, &cmap[0], B.ColMap().IndexBase(), A.Comm());
1062 const int myrowlength = A.RowMap().NumMyElements();
1063 const Epetra_Map& rowmap = A.RowMap();
1068 int educatedguess = 0;
1069 for (
int i = 0; i < myrowlength; ++i) {
1071 ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1072 if (rowlength>educatedguess)
1073 educatedguess = rowlength;
1077 RCP<Epetra_CrsMatrix> result = rcp(
new Epetra_CrsMatrix(::Copy, A.RangeMap(), gcmap, educatedguess,
false));
1079 std::vector<int> gcid(educatedguess);
1080 for (
int i = 0; i < myrowlength; ++i) {
1081 const int grid = rowmap.GID(i);
1083 ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1084 if (!rowlength)
continue;
1085 if ((
int)gcid.size() < rowlength) gcid.resize(rowlength);
1086 for (
int j = 0; j < rowlength; ++j) {
1087 gcid[j] = gcmap.GID(bindx[j]);
1091 int err = result->InsertGlobalValues(grid, rowlength, val, &gcid[0]);
1092 if (err != 0 && err != 1) {
1093 std::ostringstream errStr;
1094 errStr <<
"Epetra_CrsMatrix::InsertGlobalValues returned err=" << err;
1099 if (bindx) ML_free(bindx);
1100 if (val) ML_free(val);
1101 ML_Operator_Destroy(&ml_AtimesB);
1102 ML_Comm_Destroy(&comm);
1105 #else // no MUELU_ML 1107 "No ML multiplication available. This feature is currently not supported by Xpetra.");
1108 return Teuchos::null;
1111 #endif //ifdef HAVE_XPETRA_EPETRAEXT 1125 Teuchos::FancyOStream& fos,
1126 bool doFillComplete =
true,
1127 bool doOptimizeStorage =
true) {
1129 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1138 RCP<BlockedCrsMatrix> C = rcp(
new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 ));
1140 for (
size_t i = 0; i < A.
Rows(); ++i) {
1141 for (
size_t j = 0; j < B.
Cols(); ++j) {
1144 for (
size_t l = 0; l < B.
Rows(); ++l) {
1148 if (crmat1.is_null() || crmat2.is_null())
1150 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1153 RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<
CrsMatrixWrap>(crmat1);
1154 RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<
CrsMatrixWrap>(crmat2);
1155 TEUCHOS_TEST_FOR_EXCEPTION((crop1==Teuchos::null) != (crop2==Teuchos::null),
Xpetra::Exceptions::RuntimeError,
"A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1158 RCP<Matrix> temp = Teuchos::null;
1160 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1161 temp =
Multiply (*crop1,
false, *crop2,
false, fos);
1163 RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(crmat1);
1164 RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(crmat2);
1165 TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==
true,
Xpetra::Exceptions::BadCast,
"A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1166 TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==
true,
Xpetra::Exceptions::BadCast,
"B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1167 TEUCHOS_TEST_FOR_EXCEPTION(bop1->Cols()!=bop2->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << bop1->Cols() <<
" columns and B has " << bop2->Rows() <<
" rows. Matrices are not compatible! (TwoMatrixMultiplyBlock)");
1168 TEUCHOS_TEST_FOR_EXCEPTION(bop1->getDomainMap()->isSameAs(*(bop2->getRangeMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of A is not the same as range map of B. Matrices are not compatible! (TwoMatrixMultiplyBlock)");
1171 temp =
TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1173 RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(temp);
1174 TEUCHOS_TEST_FOR_EXCEPTION(btemp->Rows()!=bop1->Rows(),
Xpetra::Exceptions::RuntimeError,
"Number of block rows of local blocked operator is " << btemp->Rows() <<
" but should be " << bop1->Rows() <<
". (TwoMatrixMultiplyBlock)");
1175 TEUCHOS_TEST_FOR_EXCEPTION(btemp->Cols()!=bop2->Cols(),
Xpetra::Exceptions::RuntimeError,
"Number of block cols of local blocked operator is " << btemp->Cols() <<
" but should be " << bop2->Cols() <<
". (TwoMatrixMultiplyBlock)");
1176 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getRangeMapExtractor()->getFullMap()->isSameAs(*(bop1->getRangeMapExtractor()->getFullMap())) ==
false,
Xpetra::Exceptions::RuntimeError,
"Range map of local blocked operator should be same as first operator. (TwoMatrixMultiplyBlock)");
1177 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getDomainMapExtractor()->getFullMap()->isSameAs(*(bop2->getDomainMapExtractor()->getFullMap())) ==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of local blocked operator should be same as second operator. (TwoMatrixMultiplyBlock)");
1178 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getRangeMapExtractor()->getThyraMode() != bop1->getRangeMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Thyra mode of local range map extractor incompatible with range map extractor of A (TwoMatrixMultiplyBlock)");
1179 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getDomainMapExtractor()->getThyraMode() != bop2->getDomainMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Thyra mode of local domain map extractor incompatible with domain map extractor of B (TwoMatrixMultiplyBlock)");
1184 RCP<Matrix> addRes = null;
1193 if (!Cij.is_null()) {
1194 if (Cij->isFillComplete())
1197 C->setMatrix(i, j, Cij);
1199 C->setMatrix(i, j, Teuchos::null);
1224 "TwoMatrixAdd: matrix row maps are not the same.");
1227 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 1232 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1236 std::ostringstream buf;
1241 #ifdef HAVE_XPETRA_TPETRA 1242 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \ 1243 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT)))) 1249 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1271 const Matrix& B,
bool transposeB,
const SC& beta,
1272 RCP<Matrix>& C, Teuchos::FancyOStream &fos,
bool AHasFixedNnzPerRow =
false) {
1273 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1274 RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1275 RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(rcpA);
1276 RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(rcpB);
1278 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1283 if (C == Teuchos::null) {
1284 size_t maxNzInA = 0;
1285 size_t maxNzInB = 0;
1286 size_t numLocalRows = 0;
1293 if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1296 Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1298 if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1299 for (
size_t i = 0; i < numLocalRows; ++i)
1303 for (
size_t i = 0; i < numLocalRows; ++i)
1307 fos <<
"MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)" 1308 <<
", using static profiling" << std::endl;
1315 LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
1322 fos <<
"nnzPerRowInA = " << nnzPerRowInA <<
", nnzPerRowInB = " << nnzPerRowInB << std::endl;
1323 fos <<
"MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
1324 <<
", max possible nnz per row in sum = " << maxPossible
1330 fos <<
"MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1334 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 1338 Epetra_CrsMatrix* ref2epC = &*epC;
1341 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1349 #ifdef HAVE_XPETRA_TPETRA 1350 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \ 1351 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT)))) 1358 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
1366 if (A.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(A));
1367 if (B.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(B));
1371 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1372 RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1373 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1379 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
1380 RCP<Matrix> Cij = Teuchos::null;
1381 if(rcpA != Teuchos::null &&
1382 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1385 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1386 Cij, fos, AHasFixedNnzPerRow);
1387 }
else if (rcpA == Teuchos::null &&
1388 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1389 Cij = rcpBopB->getMatrix(i,j);
1390 }
else if (rcpA != Teuchos::null &&
1391 rcpBopB->getMatrix(i,j)==Teuchos::null) {
1392 Cij = Teuchos::rcp_const_cast<
Matrix>(rcpA);
1394 Cij = Teuchos::null;
1397 if (!Cij.is_null()) {
1398 if (Cij->isFillComplete())
1400 Cij->fillComplete();
1401 bC->setMatrix(i, j, Cij);
1403 bC->setMatrix(i, j, Teuchos::null);
1408 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
1409 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1410 RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
1416 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
1417 RCP<Matrix> Cij = Teuchos::null;
1418 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1419 rcpB!=Teuchos::null) {
1421 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1422 *rcpB, transposeB, beta,
1423 Cij, fos, AHasFixedNnzPerRow);
1424 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1425 rcpB!=Teuchos::null) {
1426 Cij = Teuchos::rcp_const_cast<
Matrix>(rcpB);
1427 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1428 rcpB==Teuchos::null) {
1429 Cij = rcpBopA->getMatrix(i,j);
1431 Cij = Teuchos::null;
1434 if (!Cij.is_null()) {
1435 if (Cij->isFillComplete())
1437 Cij->fillComplete();
1438 bC->setMatrix(i, j, Cij);
1440 bC->setMatrix(i, j, Teuchos::null);
1450 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->Rows()!=rcpBopB->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << rcpBopA->Rows() <<
" rows and B has " << rcpBopA->Rows() <<
" rows. Matrices are not compatible! (TwoMatrixAdd)");
1451 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->Rows()!=rcpBopB->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << rcpBopA->Cols() <<
" cols and B has " << rcpBopA->Cols() <<
" cols. Matrices are not compatible! (TwoMatrixAdd)");
1452 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getRangeMap()->isSameAs(*(rcpBopB->getRangeMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Range map of A is not the same as range map of B. Matrices are not compatible! (TwoMatrixAdd)");
1453 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getDomainMap()->isSameAs(*(rcpBopB->getDomainMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of A is not the same as domain map of B. Matrices are not compatible! (TwoMatrixAdd)");
1454 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getRangeMapExtractor()->getThyraMode() != rcpBopB->getRangeMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Different Thyra/Xpetra style gids in RangeMapExtractor of A and B. Matrices are not compatible! (TwoMatrixAdd)");
1455 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getDomainMapExtractor()->getThyraMode() != rcpBopB->getDomainMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Different Thyra/Xpetra style gids in DomainMapExtractor of A and B. Matrices are not compatible! (TwoMatrixAdd)");
1457 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1458 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1463 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
1464 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
1466 RCP<Matrix> Cij = Teuchos::null;
1467 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1468 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1471 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1472 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1473 Cij, fos, AHasFixedNnzPerRow);
1474 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1475 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1476 Cij = rcpBopB->getMatrix(i,j);
1477 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1478 rcpBopB->getMatrix(i,j)==Teuchos::null) {
1479 Cij = rcpBopA->getMatrix(i,j);
1481 Cij = Teuchos::null;
1484 if (!Cij.is_null()) {
1485 if (Cij->isFillComplete())
1488 Cij->fillComplete();
1489 bC->setMatrix(i, j, Cij);
1491 bC->setMatrix(i, j, Teuchos::null);
1536 const Matrix& B,
bool transposeB,
1538 bool call_FillComplete_on_result =
true,
1539 bool doOptimizeStorage =
true,
1540 const std::string & label = std::string(),
1541 const RCP<ParameterList>& params = null) {
1542 TEUCHOS_TEST_FOR_EXCEPTION(transposeA ==
false && C.
getRowMap()->isSameAs(*A.
getRowMap()) ==
false,
1544 TEUCHOS_TEST_FOR_EXCEPTION(transposeA ==
true && C.
getRowMap()->isSameAs(*A.
getDomainMap()) ==
false,
1551 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1554 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 1559 int i = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, epC, haveMultiplyDoFillComplete);
1560 if (haveMultiplyDoFillComplete) {
1571 std::ostringstream buf;
1573 std::string msg =
"EpetraExt::MatrixMatrix::Multiply return value of " + buf.str();
1581 #ifdef HAVE_XPETRA_TPETRA 1582 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \ 1583 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG)))) 1592 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1599 if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1600 RCP<Teuchos::ParameterList> fillParams = rcp(
new Teuchos::ParameterList());
1601 fillParams->set(
"Optimize Storage",doOptimizeStorage);
1608 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<
Matrix>(Teuchos::rcpFromRef(A));
1609 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<
Matrix>(Teuchos::rcpFromRef(B));
1610 C.
CreateView(
"stridedMaps", rcpA, transposeA, rcpB, transposeB);
1636 const Matrix& B,
bool transposeB,
1638 Teuchos::FancyOStream& fos,
1639 bool doFillComplete =
true,
1640 bool doOptimizeStorage =
true,
1641 const std::string & label = std::string(),
1642 const RCP<ParameterList>& params = null) {
1648 RCP<Matrix> C = C_in;
1650 if (C == Teuchos::null) {
1651 double nnzPerRow = Teuchos::as<double>(0);
1659 nnzPerRow *= nnzPerRow;
1662 if (totalNnz < minNnz)
1666 fos <<
"Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1674 fos <<
"Reuse C pattern" << std::endl;
1677 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params);
1693 const Matrix& B,
bool transposeB,
1694 Teuchos::FancyOStream &fos,
1695 bool callFillCompleteOnResult =
true,
1696 bool doOptimizeStorage =
true,
1697 const std::string & label = std::string(),
1698 const RCP<ParameterList>& params = null) {
1699 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1702 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 1705 const Epetra_CrsMatrix& epB,
1706 Teuchos::FancyOStream& fos) {
1708 "No ML multiplication available. This feature is currently not supported by Xpetra.");
1709 return Teuchos::null;
1711 #endif //ifdef HAVE_XPETRA_EPETRAEXT 1725 Teuchos::FancyOStream& fos,
1726 bool doFillComplete =
true,
1727 bool doOptimizeStorage =
true) {
1729 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1738 RCP<BlockedCrsMatrix> C = rcp(
new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 ));
1740 for (
size_t i = 0; i < A.
Rows(); ++i) {
1741 for (
size_t j = 0; j < B.
Cols(); ++j) {
1744 for (
size_t l = 0; l < B.
Rows(); ++l) {
1748 if (crmat1.is_null() || crmat2.is_null())
1750 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1753 RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<
CrsMatrixWrap>(crmat1);
1754 RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<
CrsMatrixWrap>(crmat2);
1755 TEUCHOS_TEST_FOR_EXCEPTION((crop1==Teuchos::null) != (crop2==Teuchos::null),
Xpetra::Exceptions::RuntimeError,
"A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1758 RCP<Matrix> temp = Teuchos::null;
1760 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1761 temp =
Multiply (*crop1,
false, *crop2,
false, fos);
1763 RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(crmat1);
1764 RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(crmat2);
1765 TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==
true,
Xpetra::Exceptions::BadCast,
"A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1766 TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==
true,
Xpetra::Exceptions::BadCast,
"B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1767 TEUCHOS_TEST_FOR_EXCEPTION(bop1->Cols()!=bop2->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << bop1->Cols() <<
" columns and B has " << bop2->Rows() <<
" rows. Matrices are not compatible! (TwoMatrixMultiplyBlock)");
1768 TEUCHOS_TEST_FOR_EXCEPTION(bop1->getDomainMap()->isSameAs(*(bop2->getRangeMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of A is not the same as range map of B. Matrices are not compatible! (TwoMatrixMultiplyBlock)");
1771 temp =
TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1773 RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<
BlockedCrsMatrix>(temp);
1774 TEUCHOS_TEST_FOR_EXCEPTION(btemp->Rows()!=bop1->Rows(),
Xpetra::Exceptions::RuntimeError,
"Number of block rows of local blocked operator is " << btemp->Rows() <<
" but should be " << bop1->Rows() <<
". (TwoMatrixMultiplyBlock)");
1775 TEUCHOS_TEST_FOR_EXCEPTION(btemp->Cols()!=bop2->Cols(),
Xpetra::Exceptions::RuntimeError,
"Number of block cols of local blocked operator is " << btemp->Cols() <<
" but should be " << bop2->Cols() <<
". (TwoMatrixMultiplyBlock)");
1776 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getRangeMapExtractor()->getFullMap()->isSameAs(*(bop1->getRangeMapExtractor()->getFullMap())) ==
false,
Xpetra::Exceptions::RuntimeError,
"Range map of local blocked operator should be same as first operator. (TwoMatrixMultiplyBlock)");
1777 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getDomainMapExtractor()->getFullMap()->isSameAs(*(bop2->getDomainMapExtractor()->getFullMap())) ==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of local blocked operator should be same as second operator. (TwoMatrixMultiplyBlock)");
1778 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getRangeMapExtractor()->getThyraMode() != bop1->getRangeMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Thyra mode of local range map extractor incompatible with range map extractor of A (TwoMatrixMultiplyBlock)");
1779 TEUCHOS_TEST_FOR_EXCEPTION(btemp->getDomainMapExtractor()->getThyraMode() != bop2->getDomainMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Thyra mode of local domain map extractor incompatible with domain map extractor of B (TwoMatrixMultiplyBlock)");
1784 RCP<Matrix> addRes = null;
1793 if (!Cij.is_null()) {
1794 if (Cij->isFillComplete())
1797 C->setMatrix(i, j, Cij);
1799 C->setMatrix(i, j, Teuchos::null);
1824 "TwoMatrixAdd: matrix row maps are not the same.");
1827 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 1832 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1836 std::ostringstream buf;
1841 #ifdef HAVE_XPETRA_TPETRA 1842 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \ 1843 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG)))) 1849 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1871 const Matrix& B,
bool transposeB,
const SC& beta,
1872 RCP<Matrix>& C, Teuchos::FancyOStream &fos,
bool AHasFixedNnzPerRow =
false) {
1873 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1874 RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1875 RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(rcpA);
1876 RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(rcpB);
1878 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1880 "TwoMatrixAdd: matrix row maps are not the same.");
1882 if (C == Teuchos::null) {
1883 size_t maxNzInA = 0;
1884 size_t maxNzInB = 0;
1885 size_t numLocalRows = 0;
1892 if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1895 Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1897 if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1898 for (
size_t i = 0; i < numLocalRows; ++i)
1902 for (
size_t i = 0; i < numLocalRows; ++i)
1906 fos <<
"MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)" 1907 <<
", using static profiling" << std::endl;
1914 LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
1921 fos <<
"nnzPerRowInA = " << nnzPerRowInA <<
", nnzPerRowInB = " << nnzPerRowInB << std::endl;
1922 fos <<
"MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
1923 <<
", max possible nnz per row in sum = " << maxPossible
1929 fos <<
"MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1933 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) 1937 Epetra_CrsMatrix* ref2epC = &*epC;
1940 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1948 #ifdef HAVE_XPETRA_TPETRA 1949 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \ 1950 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG)))) 1957 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
1965 if (A.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(A));
1966 if (B.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(B));
1970 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1971 RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1972 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1978 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
1979 RCP<Matrix> Cij = Teuchos::null;
1980 if(rcpA != Teuchos::null &&
1981 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1984 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1985 Cij, fos, AHasFixedNnzPerRow);
1986 }
else if (rcpA == Teuchos::null &&
1987 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1988 Cij = rcpBopB->getMatrix(i,j);
1989 }
else if (rcpA != Teuchos::null &&
1990 rcpBopB->getMatrix(i,j)==Teuchos::null) {
1991 Cij = Teuchos::rcp_const_cast<
Matrix>(rcpA);
1993 Cij = Teuchos::null;
1996 if (!Cij.is_null()) {
1997 if (Cij->isFillComplete())
1999 Cij->fillComplete();
2000 bC->setMatrix(i, j, Cij);
2002 bC->setMatrix(i, j, Teuchos::null);
2007 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
2008 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2009 RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
2015 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
2016 RCP<Matrix> Cij = Teuchos::null;
2017 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2018 rcpB!=Teuchos::null) {
2020 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2021 *rcpB, transposeB, beta,
2022 Cij, fos, AHasFixedNnzPerRow);
2023 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2024 rcpB!=Teuchos::null) {
2025 Cij = Teuchos::rcp_const_cast<
Matrix>(rcpB);
2026 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2027 rcpB==Teuchos::null) {
2028 Cij = rcpBopA->getMatrix(i,j);
2030 Cij = Teuchos::null;
2033 if (!Cij.is_null()) {
2034 if (Cij->isFillComplete())
2036 Cij->fillComplete();
2037 bC->setMatrix(i, j, Cij);
2039 bC->setMatrix(i, j, Teuchos::null);
2049 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->Rows()!=rcpBopB->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << rcpBopA->Rows() <<
" rows and B has " << rcpBopA->Rows() <<
" rows. Matrices are not compatible! (TwoMatrixAdd)");
2050 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->Rows()!=rcpBopB->Rows(),
Xpetra::Exceptions::RuntimeError,
"A has " << rcpBopA->Cols() <<
" cols and B has " << rcpBopA->Cols() <<
" cols. Matrices are not compatible! (TwoMatrixAdd)");
2051 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getRangeMap()->isSameAs(*(rcpBopB->getRangeMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Range map of A is not the same as range map of B. Matrices are not compatible! (TwoMatrixAdd)");
2052 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getDomainMap()->isSameAs(*(rcpBopB->getDomainMap()))==
false,
Xpetra::Exceptions::RuntimeError,
"Domain map of A is not the same as domain map of B. Matrices are not compatible! (TwoMatrixAdd)");
2053 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getRangeMapExtractor()->getThyraMode() != rcpBopB->getRangeMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Different Thyra/Xpetra style gids in RangeMapExtractor of A and B. Matrices are not compatible! (TwoMatrixAdd)");
2054 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA->getDomainMapExtractor()->getThyraMode() != rcpBopB->getDomainMapExtractor()->getThyraMode(),
Xpetra::Exceptions::RuntimeError,
"Different Thyra/Xpetra style gids in DomainMapExtractor of A and B. Matrices are not compatible! (TwoMatrixAdd)");
2056 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2057 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2062 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
2063 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
2065 RCP<Matrix> Cij = Teuchos::null;
2066 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2067 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2070 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2071 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2072 Cij, fos, AHasFixedNnzPerRow);
2073 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2074 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2075 Cij = rcpBopB->getMatrix(i,j);
2076 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2077 rcpBopB->getMatrix(i,j)==Teuchos::null) {
2078 Cij = rcpBopA->getMatrix(i,j);
2080 Cij = Teuchos::null;
2083 if (!Cij.is_null()) {
2084 if (Cij->isFillComplete())
2087 Cij->fillComplete();
2088 bC->setMatrix(i, j, Cij);
2090 bC->setMatrix(i, j, Teuchos::null);
2099 #endif // HAVE_XPETRA_EPETRA 2103 #define XPETRA_MATRIXMATRIX_SHORT virtual size_t getNodeNumRows() const =0
Returns the number of matrix rows owned on the calling node.
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.
static void TwoMatrixAdd(const Matrix &A, bool transposeA, const SC &alpha, const Matrix &B, bool transposeB, const SC &beta, RCP< Matrix > &C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow=false)
Helper function to calculate C = alpha*A + beta*B.
static void TwoMatrixAdd(const Matrix &A, bool transposeA, const SC &alpha, const Matrix &B, bool transposeB, const SC &beta, RCP< Matrix > &C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow=false)
Helper function to calculate C = alpha*A + beta*B.
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, RCP< Matrix > C_in, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
Helper function to do matrix-matrix multiply.
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Matrix &Op)
Xpetra utility class containing transformation routines between Xpetra::Matrix and Epetra/Tpetra obje...
virtual void resumeFill(const RCP< ParameterList > ¶ms=null)=0
Exception throws to report errors in the internal logical of the program.
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
virtual size_t getGlobalMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on all nodes.
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Exception indicating invalid cast attempted.
virtual void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > ¶ms=null)=0
Signal that data entry is complete, specifying domain and range maps.
static void TwoMatrixAdd(const Matrix &A, bool transposeA, const SC &alpha, const Matrix &B, bool transposeB, const SC &beta, RCP< Matrix > &C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow=false)
Helper function to calculate C = alpha*A + beta*B.
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool callFillCompleteOnResult=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
Helper function to do matrix-matrix multiply.
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
RCP< CrsMatrix > getCrsMatrix() const
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, RCP< Matrix > C_in, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
Helper function to do matrix-matrix multiply.
static Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraCrs(const Matrix &Op)
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, RCP< Matrix > C_in, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
Helper function to do matrix-matrix multiply.
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
Returns the current number of entries on this node in the specified local row.
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Matrix > Op)
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
Exception throws to report incompatible objects (like maps).
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
RCP< const Map > getRangeMap() const
Returns the Map associated with the full range of this operator.
virtual global_size_t getGlobalNumEntries() const =0
Returns the global number of entries in this matrix.
Concrete implementation of Xpetra::Matrix.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
bool IsView(const viewLabel_t viewLabel) const
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
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.
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
virtual size_t Rows() const
number of row blocks
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
virtual size_t Cols() const
number of column blocks
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool callFillCompleteOnResult=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
Helper function to do matrix-matrix multiply.
RCP< const Map > getDomainMap() const
Returns the Map associated with the full domain of this operator.
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra-specific matrix class.
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool callFillCompleteOnResult=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
Helper function to do matrix-matrix multiply.