48#ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_
49#define PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_
54#include "Xpetra_CrsMatrixWrap.hpp"
55#include "Xpetra_MapExtractor.hpp"
56#include "Xpetra_Map.hpp"
59#include "Xpetra_StridedMapFactory.hpp"
60#include "Xpetra_StridedMap.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>
72#ifdef HAVE_XPETRA_TPETRA
73#include <TpetraExt_MatrixMatrix.hpp>
74#include <Tpetra_RowMatrixTransposer.hpp>
75#include <MatrixMarket_Tpetra.hpp>
76#include <Xpetra_TpetraCrsMatrix.hpp>
77#include <Xpetra_TpetraBlockCrsMatrix.hpp>
78#include <Tpetra_BlockCrsMatrix_Helpers.hpp>
79#include <Xpetra_TpetraMultiVector.hpp>
80#include <Xpetra_TpetraVector.hpp>
91 template <
class Scalar,
92 class LocalOrdinal = int,
93 class GlobalOrdinal = LocalOrdinal,
100#ifdef HAVE_XPETRA_EPETRA
105 "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
110 "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
112 return tmp_ECrsMtx->getEpetra_CrsMatrix();
120 "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
126 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
136 "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
138 return *tmp_ECrsMtx->getEpetra_CrsMatrix();
154 return *Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
162#ifdef HAVE_XPETRA_TPETRA
172 return tmp_ECrsMtx->getTpetra_CrsMatrix();
183 return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
194 return *tmp_TCrsMtx->getTpetra_CrsMatrix();
209 return *Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_TCrsMtx->getTpetra_CrsMatrix());
218 if(crsOp == Teuchos::null)
return false;
221 if(tmp_ECrsMtx == Teuchos::null)
return false;
231 if(tmp_ECrsMtx == Teuchos::null)
return false;
245 return tmp_BlockCrs->getTpetra_BlockCrsMatrix();
255 return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
264 return *tmp_BlockCrs->getTpetra_BlockCrsMatrix();
276 return *tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
284 if(crsOp == Teuchos::null)
return false;
287 if(tmp_BlockCrs == Teuchos::null)
return false;
296 if(tmp_BlockCrs == Teuchos::null)
return false;
313#ifdef HAVE_XPETRA_TPETRA
316 const tcrs_matrix_type& A,
bool transposeA,
const typename tcrs_matrix_type::scalar_type alpha,
317 const tcrs_matrix_type& B,
bool transposeB,
const typename tcrs_matrix_type::scalar_type beta)
322 using Teuchos::rcp_implicit_cast;
323 using Teuchos::rcpFromRef;
328 using transposer_type = Tpetra::RowMatrixTransposer<SC,LO,GO,NO>;
329 using import_type = Tpetra::Import<LO,GO,NO>;
332 Aprime = transposer_type(Aprime).createTranspose();
334 if(A.isFillComplete() && B.isFillComplete())
338 addParams->set(
"Call fillComplete",
false);
340 Tpetra::MatrixMatrix::add<SC,LO,GO,NO>(beta, transposeB, B, alpha,
false, *Aprime, *C, Teuchos::null, Teuchos::null, addParams);
341 return rcp_implicit_cast<Matrix>(
rcp(
new CrsWrap(rcp_implicit_cast<CrsType>(
rcp(
new XTCrsType(C))))));
351 Bprime = transposer_type(Bprime).createTranspose();
353 if(!(Aprime->getRowMap()->isSameAs(*(Bprime->getRowMap()))))
355 auto import =
rcp(
new import_type(Bprime->getRowMap(), Aprime->getRowMap()));
356 Bprime = Tpetra::importAndFillCompleteCrsMatrix<tcrs_matrix_type>(Bprime, *
import, Aprime->getDomainMap(), Aprime->getRangeMap());
359 LO numLocalRows = Aprime->getLocalNumRows();
362 for(
LO i = 0; i < numLocalRows; i++)
364 allocPerRow[i] = Aprime->getNumEntriesInLocalRow(i) + Bprime->getNumEntriesInLocalRow(i);
369 Tpetra::MatrixMatrix::Add<SC,LO,GO,NO>(
370 *Aprime,
false, alpha,
371 *Bprime,
false, beta,
373 return rcp(
new CrsWrap(rcp_implicit_cast<CrsType>(
rcp(
new XTCrsType(C)))));
378#ifdef HAVE_XPETRA_EPETRAEXT
389 if(fillCompleteResult) {
390 errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp,
true);
397 errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp,
false);
400 long long* globalElementList =
nullptr;
401 Crowmap.MyGlobalElementsPtr(globalElementList);
403 for(
int i = 0; i < numLocalRows; i++)
409 for(
int i = 0; i < numLocalRows; i++)
411 int gid = globalElementList[i];
421 std::ostringstream buf;
423 std::string msg =
"EpetraExt::MatrixMatrix::Multiply returned nonzero error code " + buf.str();
430 template <
class Scalar,
432 class GlobalOrdinal ,
435#undef XPETRA_MATRIXMATRIX_SHORT
465 const Matrix& B,
bool transposeB,
467 bool call_FillComplete_on_result =
true,
468 bool doOptimizeStorage =
true,
469 const std::string & label = std::string(),
480 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
483#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
484 throw(
Xpetra::Exceptions::RuntimeError(
"Xpetra::MatrixMatrix::Multiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
490#ifdef HAVE_XPETRA_TPETRA
492 if(helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
494 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
495 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
496 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
500 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
502 else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
507 std::cout<<
"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
511 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
517 if(!params.is_null()) {
519 new_params->set(
"compute global constants",
true);
524 Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
533 Ac_w->replaceCrsMatrix(Ac_p);
544 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
546 fillParams->set(
"Optimize Storage", doOptimizeStorage);
553 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
554 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
555 C.
CreateView(
"stridedMaps", rcpA, transposeA, rcpB, transposeB);
582 bool doFillComplete =
true,
583 bool doOptimizeStorage =
true,
584 const std::string & label = std::string(),
593 if (C == Teuchos::null) {
594 double nnzPerRow = Teuchos::as<double>(0);
603 nnzPerRow *= nnzPerRow;
606 if (totalNnz < minNnz)
610 fos <<
"Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
618 fos <<
"Reuse C pattern" << std::endl;
621 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params);
637 bool callFillCompleteOnResult =
true,
bool doOptimizeStorage =
true,
const std::string& label = std::string(),
639 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
642#ifdef HAVE_XPETRA_EPETRAEXT
647 throw(
Xpetra::Exceptions::RuntimeError(
"MLTwoMatrixMultiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
665 bool doFillComplete =
true,
666 bool doOptimizeStorage =
true) {
668 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
679 for (
size_t i = 0; i < A.
Rows(); ++i) {
680 for (
size_t j = 0; j < B.
Cols(); ++j) {
683 for (
size_t l = 0; l < B.
Rows(); ++l) {
692 auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
693 auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
695 if(unwrap1.is_null() != unwrap2.is_null()) {
696 if(unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
697 crmat1 = unwrap1->getCrsMatrix();
698 if(unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
699 crmat2 = unwrap2->getCrsMatrix();
706 "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
710 Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
712 Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
715 "crmat1 does not have global constants");
717 "crmat2 does not have global constants");
719 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
725 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
726 temp =
Multiply (*crop1,
false, *crop2,
false, fos);
736 temp =
TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
759 if (Cij->isFillComplete())
762 C->setMatrix(i, j, Cij);
764 C->setMatrix(i, j, Teuchos::null);
792 throw Exceptions::RuntimeError(
"TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
794#ifdef HAVE_XPETRA_TPETRA
798 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
821 const Matrix& B,
bool transposeB,
const SC& beta,
835 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
841 #ifdef HAVE_XPETRA_TPETRA
842 using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
844 const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
845 const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
846 C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
852 if (A.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(A));
853 if (B.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(B));
857 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
865 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
867 if(rcpA != Teuchos::null &&
868 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
871 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
872 Cij, fos, AHasFixedNnzPerRow);
873 }
else if (rcpA == Teuchos::null &&
874 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
875 Cij = rcpBopB->getMatrix(i,j);
876 }
else if (rcpA != Teuchos::null &&
877 rcpBopB->getMatrix(i,j)==Teuchos::null) {
878 Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
884 if (Cij->isFillComplete())
887 bC->setMatrix(i, j, Cij);
889 bC->setMatrix(i, j, Teuchos::null);
894 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
901 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
903 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
904 rcpB!=Teuchos::null) {
906 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
907 *rcpB, transposeB, beta,
908 Cij, fos, AHasFixedNnzPerRow);
909 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
910 rcpB!=Teuchos::null) {
911 Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
912 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
913 rcpB==Teuchos::null) {
914 Cij = rcpBopA->getMatrix(i,j);
920 if (Cij->isFillComplete())
923 bC->setMatrix(i, j, Cij);
925 bC->setMatrix(i, j, Teuchos::null);
947 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
948 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
951 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
952 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
954 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
955 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
956 Cij, fos, AHasFixedNnzPerRow);
957 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
958 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
959 Cij = rcpBopB->getMatrix(i,j);
960 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
961 rcpBopB->getMatrix(i,j)==Teuchos::null) {
962 Cij = rcpBopA->getMatrix(i,j);
968 if (Cij->isFillComplete())
971 bC->setMatrix(i, j, Cij);
973 bC->setMatrix(i, j, Teuchos::null);
985#ifdef HAVE_XPETRA_EPETRA
1022 const Matrix& B,
bool transposeB,
1024 bool call_FillComplete_on_result =
true,
1025 bool doOptimizeStorage =
true,
1026 const std::string & label = std::string(),
1036 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1041#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1042 helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
1047#ifdef HAVE_XPETRA_TPETRA
1048# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1049 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1052 if(helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
1054 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
1055 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
1056 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
1060 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1062 else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
1067 if(!A.
getRowMap()->getComm()->getRank())
1068 std::cout<<
"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
1072 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
1078 if(!params.is_null()) {
1080 new_params->set(
"compute global constants",
true);
1085 Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
1094 Ac_w->replaceCrsMatrix(Ac_p);
1106 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1108 fillParams->set(
"Optimize Storage",doOptimizeStorage);
1115 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1116 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1117 C.
CreateView(
"stridedMaps", rcpA, transposeA, rcpB, transposeB);
1143 const Matrix& B,
bool transposeB,
1146 bool doFillComplete =
true,
1147 bool doOptimizeStorage =
true,
1148 const std::string & label = std::string(),
1156#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) && defined(HAVE_XPETRA_ML_MMM)
1162 RCP<Matrix> C = Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO> (epC);
1163 if (doFillComplete) {
1165 fillParams->set(
"Optimize Storage", doOptimizeStorage);
1172 C->CreateView(
"stridedMaps", rcpFromRef(A), transposeA, rcpFromRef(B), transposeB);
1181 if (C == Teuchos::null) {
1182 double nnzPerRow = Teuchos::as<double>(0);
1191 nnzPerRow *= nnzPerRow;
1194 if (totalNnz < minNnz)
1198 fos <<
"Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1207 fos <<
"Reuse C pattern" << std::endl;
1210 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params);
1226 const Matrix& B,
bool transposeB,
1228 bool callFillCompleteOnResult =
true,
1229 bool doOptimizeStorage =
true,
1230 const std::string & label = std::string(),
1232 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1235#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1240#if defined(HAVE_XPETRA_ML_MMM)
1242 ML_Comm_Create(&comm);
1243 fos <<
"****** USING ML's MATRIX MATRIX MULTIPLY ******" << std::endl;
1248 ML_Comm_Set_UsrComm(comm,Mcomm->
GetMpiComm());
1251 EpetraExt::CrsMatrix_SolverMap Atransform;
1252 EpetraExt::CrsMatrix_SolverMap Btransform;
1260 ML_Operator* ml_As = ML_Operator_Create(comm);
1261 ML_Operator* ml_Bs = ML_Operator_Create(comm);
1264 ML_Operator* ml_AtimesB = ML_Operator_Create(comm);
1267 ML_2matmult(ml_As,ml_Bs,ml_AtimesB,ML_CSR_MATRIX);
1269 ML_Operator_Destroy(&ml_As);
1270 ML_Operator_Destroy(&ml_Bs);
1275 int N_local = ml_AtimesB->invec_leng;
1276 ML_CommInfoOP* getrow_comm = ml_AtimesB->getrow->pre_comm;
1278 ML_Comm* comm_AB = ml_AtimesB->comm;
1284 for (
int i = 0; i < getrow_comm->N_neighbors; i++) {
1285 N_rcvd += (getrow_comm->neighbors)[i].N_rcv;
1286 N_send += (getrow_comm->neighbors)[i].N_send;
1287 if ( ((getrow_comm->neighbors)[i].N_rcv != 0) &&
1288 ((getrow_comm->neighbors)[i].rcv_list != NULL) ) flag = 1;
1292 std::vector<double> dtemp(N_local + N_rcvd + 1);
1293 std::vector<int> cmap (N_local + N_rcvd + 1);
1294 for (
int i = 0; i < N_local; ++i) {
1296 dtemp[i] = (double) cmap[i];
1298 ML_cheap_exchange_bdry(&dtemp[0],getrow_comm,N_local,N_send,comm_AB);
1300 int count = N_local;
1301 const int neighbors = getrow_comm->N_neighbors;
1302 for (
int i = 0; i < neighbors; i++) {
1303 const int nrcv = getrow_comm->neighbors[i].N_rcv;
1304 for (
int j = 0; j < nrcv; j++)
1305 cmap[getrow_comm->neighbors[i].rcv_list[j]] = (
int) dtemp[count++];
1308 for (
int i = 0; i < N_local+N_rcvd; ++i)
1309 cmap[i] = (
int)dtemp[i];
1327 int educatedguess = 0;
1328 for (
int i = 0; i < myrowlength; ++i) {
1330 ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1331 if (rowlength>educatedguess)
1332 educatedguess = rowlength;
1338 std::vector<int> gcid(educatedguess);
1339 for (
int i = 0; i < myrowlength; ++i) {
1340 const int grid = rowmap.
GID(i);
1342 ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1343 if (!rowlength)
continue;
1344 if ((
int)gcid.size() < rowlength) gcid.resize(rowlength);
1345 for (
int j = 0; j < rowlength; ++j) {
1346 gcid[j] = gcmap.
GID(bindx[j]);
1350 int err = result->InsertGlobalValues(grid, rowlength, val, &gcid[0]);
1351 if (err != 0 && err != 1) {
1352 std::ostringstream errStr;
1353 errStr <<
"Epetra_CrsMatrix::InsertGlobalValues returned err=" << err;
1358 if (bindx) ML_free(bindx);
1359 if (val) ML_free(val);
1360 ML_Operator_Destroy(&ml_AtimesB);
1361 ML_Comm_Destroy(&comm);
1369 "No ML multiplication available. This feature is currently not supported by Xpetra.");
1388 bool doFillComplete =
true,
1389 bool doOptimizeStorage =
true) {
1391 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1402 for (
size_t i = 0; i < A.
Rows(); ++i) {
1403 for (
size_t j = 0; j < B.
Cols(); ++j) {
1406 for (
size_t l = 0; l < B.
Rows(); ++l) {
1415 auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1416 auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1418 if(unwrap1.is_null() != unwrap2.is_null()) {
1419 if(unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
1420 crmat1 = unwrap1->getCrsMatrix();
1421 if(unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
1422 crmat2 = unwrap2->getCrsMatrix();
1429 "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1433 Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
1435 Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
1438 "crmat1 does not have global constants");
1440 "crmat2 does not have global constants");
1442 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1449 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1450 temp =
Multiply (*crop1,
false, *crop2,
false, fos);
1460 temp =
TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1483 if (Cij->isFillComplete())
1486 C->setMatrix(i, j, Cij);
1488 C->setMatrix(i, j, Teuchos::null);
1513 "TwoMatrixAdd: matrix row maps are not the same.");
1516#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1521 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1525 std::ostringstream buf;
1530#ifdef HAVE_XPETRA_TPETRA
1531# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1532 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1538 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1562 const Matrix& B,
bool transposeB,
const SC& beta,
1571 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1579 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1584 size_t maxNzInA = 0;
1585 size_t maxNzInB = 0;
1586 size_t numLocalRows = 0;
1594 if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1599 if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1600 for (
size_t i = 0; i < numLocalRows; ++i)
1604 for (
size_t i = 0; i < numLocalRows; ++i)
1608 fos <<
"MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1609 <<
", using static profiling" << std::endl;
1618 fos <<
"MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1624 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1632 #ifdef HAVE_XPETRA_TPETRA
1633 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1634 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1637 using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
1638 const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
1639 const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
1640 C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
1648 if (A.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(A));
1649 if (B.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(B));
1653 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1661 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
1663 if(rcpA != Teuchos::null &&
1664 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1667 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1668 Cij, fos, AHasFixedNnzPerRow);
1669 }
else if (rcpA == Teuchos::null &&
1670 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1671 Cij = rcpBopB->getMatrix(i,j);
1672 }
else if (rcpA != Teuchos::null &&
1673 rcpBopB->getMatrix(i,j)==Teuchos::null) {
1674 Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
1676 Cij = Teuchos::null;
1680 if (Cij->isFillComplete())
1682 Cij->fillComplete();
1683 bC->setMatrix(i, j, Cij);
1685 bC->setMatrix(i, j, Teuchos::null);
1690 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
1698 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
1700 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1701 rcpB!=Teuchos::null) {
1703 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1704 *rcpB, transposeB, beta,
1705 Cij, fos, AHasFixedNnzPerRow);
1706 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1707 rcpB!=Teuchos::null) {
1708 Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
1709 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1710 rcpB==Teuchos::null) {
1711 Cij = rcpBopA->getMatrix(i,j);
1713 Cij = Teuchos::null;
1717 if (Cij->isFillComplete())
1719 Cij->fillComplete();
1720 bC->setMatrix(i, j, Cij);
1722 bC->setMatrix(i, j, Teuchos::null);
1745 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
1746 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
1749 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1750 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1753 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1754 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1755 Cij, fos, AHasFixedNnzPerRow);
1756 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1757 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1758 Cij = rcpBopB->getMatrix(i,j);
1759 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1760 rcpBopB->getMatrix(i,j)==Teuchos::null) {
1761 Cij = rcpBopA->getMatrix(i,j);
1763 Cij = Teuchos::null;
1767 if (Cij->isFillComplete())
1770 Cij->fillComplete();
1771 bC->setMatrix(i, j, Cij);
1773 bC->setMatrix(i, j, Teuchos::null);
1818 const Matrix& B,
bool transposeB,
1820 bool call_FillComplete_on_result =
true,
1821 bool doOptimizeStorage =
true,
1822 const std::string & label = std::string(),
1833 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1836#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1837 helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
1842#ifdef HAVE_XPETRA_TPETRA
1843# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1844 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1847 if(helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
1849 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
1850 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
1851 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
1855 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1857 else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
1862 if(!A.
getRowMap()->getComm()->getRank())
1863 std::cout<<
"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
1867 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
1873 if(!params.is_null()) {
1875 new_params->set(
"compute global constants",
true);
1880 Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
1889 Ac_w->replaceCrsMatrix(Ac_p);
1902 if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1904 fillParams->set(
"Optimize Storage",doOptimizeStorage);
1911 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1912 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1913 C.
CreateView(
"stridedMaps", rcpA, transposeA, rcpB, transposeB);
1939 const Matrix& B,
bool transposeB,
1942 bool doFillComplete =
true,
1943 bool doOptimizeStorage =
true,
1944 const std::string & label = std::string(),
1953 if (C == Teuchos::null) {
1954 double nnzPerRow = Teuchos::as<double>(0);
1963 nnzPerRow *= nnzPerRow;
1966 if (totalNnz < minNnz)
1970 fos <<
"Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1978 fos <<
"Reuse C pattern" << std::endl;
1981 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params);
1997 const Matrix& B,
bool transposeB,
1999 bool callFillCompleteOnResult =
true,
2000 bool doOptimizeStorage =
true,
2001 const std::string & label = std::string(),
2003 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
2006#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
2012 "No ML multiplication available. This feature is currently not supported by Xpetra.");
2030 bool doFillComplete =
true,
2031 bool doOptimizeStorage =
true) {
2033 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
2044 for (
size_t i = 0; i < A.
Rows(); ++i) {
2045 for (
size_t j = 0; j < B.
Cols(); ++j) {
2048 for (
size_t l = 0; l < B.
Rows(); ++l) {
2057 auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
2058 auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
2060 if(unwrap1.is_null() != unwrap2.is_null()) {
2061 if(unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
2062 crmat1 = unwrap1->getCrsMatrix();
2063 if(unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
2064 crmat2 = unwrap2->getCrsMatrix();
2071 "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
2075 Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
2077 Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
2080 "crmat1 does not have global constants");
2082 "crmat2 does not have global constants");
2084 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
2090 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
2091 temp =
Multiply (*crop1,
false, *crop2,
false, fos);
2101 temp =
TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
2124 if (Cij->isFillComplete())
2127 C->setMatrix(i, j, Cij);
2129 C->setMatrix(i, j, Teuchos::null);
2154 "TwoMatrixAdd: matrix row maps are not the same.");
2157#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
2162 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
2166 std::ostringstream buf;
2171#ifdef HAVE_XPETRA_TPETRA
2172# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
2173 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
2179 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
2203 const Matrix& B,
bool transposeB,
const SC& beta,
2210 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
2212 "TwoMatrixAdd: matrix row maps are not the same.");
2215 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
2220 size_t maxNzInA = 0;
2221 size_t maxNzInB = 0;
2222 size_t numLocalRows = 0;
2230 if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
2235 if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
2236 for (
size_t i = 0; i < numLocalRows; ++i)
2240 for (
size_t i = 0; i < numLocalRows; ++i)
2244 fos <<
"MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
2245 <<
", using static profiling" << std::endl;
2254 fos <<
"MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
2260 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
2268 #ifdef HAVE_XPETRA_TPETRA
2269 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
2270 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
2274 using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
2277 C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
2285 if (A.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(A));
2286 if (B.
IsView(
"stridedMaps")) C->CreateView(
"stridedMaps", rcpFromRef(B));
2290 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
2298 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
2300 if(rcpA != Teuchos::null &&
2301 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2304 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2305 Cij, fos, AHasFixedNnzPerRow);
2306 }
else if (rcpA == Teuchos::null &&
2307 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2308 Cij = rcpBopB->getMatrix(i,j);
2309 }
else if (rcpA != Teuchos::null &&
2310 rcpBopB->getMatrix(i,j)==Teuchos::null) {
2311 Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
2313 Cij = Teuchos::null;
2317 if (Cij->isFillComplete())
2319 Cij->fillComplete();
2320 bC->setMatrix(i, j, Cij);
2322 bC->setMatrix(i, j, Teuchos::null);
2327 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
2335 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
2337 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2338 rcpB!=Teuchos::null) {
2340 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2341 *rcpB, transposeB, beta,
2342 Cij, fos, AHasFixedNnzPerRow);
2343 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2344 rcpB!=Teuchos::null) {
2345 Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
2346 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2347 rcpB==Teuchos::null) {
2348 Cij = rcpBopA->getMatrix(i,j);
2350 Cij = Teuchos::null;
2354 if (Cij->isFillComplete())
2356 Cij->fillComplete();
2357 bC->setMatrix(i, j, Cij);
2359 bC->setMatrix(i, j, Teuchos::null);
2382 for (
size_t i = 0; i < rcpBopA->Rows(); ++i) {
2383 for (
size_t j = 0; j < rcpBopB->Cols(); ++j) {
2386 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2387 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2390 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2391 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2392 Cij, fos, AHasFixedNnzPerRow);
2393 }
else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2394 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2395 Cij = rcpBopB->getMatrix(i,j);
2396 }
else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2397 rcpBopB->getMatrix(i,j)==Teuchos::null) {
2398 Cij = rcpBopA->getMatrix(i,j);
2400 Cij = Teuchos::null;
2404 if (Cij->isFillComplete())
2407 Cij->fillComplete();
2408 bC->setMatrix(i, j, Cij);
2410 bC->setMatrix(i, j, Teuchos::null);
2423#define XPETRA_MATRIXMATRIX_SHORT
int NumMyElements() const
int ExtractGlobalRowView(int GlobalRow, int &NumEntries, double *&Values, int *&Indices) const
int NumGlobalEntries(long long Row) const
const Epetra_Map & RowMap() const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
const Epetra_Comm & Comm() const
const Epetra_Map & RangeMap() const
const Epetra_Map & ColMap() const
const Epetra_Map & DomainMap() const
MPI_Comm GetMpiComm() const
static RCP< Time > getNewTimer(const std::string &name)
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual size_t Rows() const
number of row blocks
RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode.
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.
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
virtual size_t Cols() const
number of column blocks
Concrete implementation of Xpetra::Matrix.
RCP< CrsMatrix > getCrsMatrix() const
Exception indicating invalid cast attempted.
Exception throws to report incompatible objects (like maps).
Exception throws to report errors in the internal logical of the program.
Tpetra::CrsMatrix< SC, LO, GO, NO > tcrs_matrix_type
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
static Tpetra::BlockCrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraBlockCrs(const Matrix &Op)
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
static bool isTpetraCrs(const Matrix &Op)
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Matrix &Op)
static Teuchos::RCP< Matrix > tpetraAdd(const tcrs_matrix_type &A, bool transposeA, const typename tcrs_matrix_type::scalar_type alpha, const tcrs_matrix_type &B, bool transposeB, const typename tcrs_matrix_type::scalar_type beta)
static bool isTpetraBlockCrs(RCP< Matrix > Op)
static bool isTpetraBlockCrs(const Matrix &Op)
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Matrix > Op)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
static bool isTpetraCrs(RCP< Matrix > Op)
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
static Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraCrs(const Matrix &Op)
static RCP< Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraBlockCrs(RCP< Matrix > Op)
static const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > & Op2TpetraBlockCrs(const Matrix &Op)
static void epetraExtMult(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool fillCompleteResult)
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
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, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
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 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 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.
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 > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
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.
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< 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 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< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &, const Epetra_CrsMatrix &, Teuchos::FancyOStream &)
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< 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 RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
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.
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, 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 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)
Xpetra-specific matrix class.
virtual global_size_t getGlobalNumEntries() const =0
Returns the global number of entries in this matrix.
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode.
virtual size_t getLocalNumRows() const =0
Returns the number of matrix rows owned on the calling node.
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.
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
Returns the current number of entries on this node in the specified local row.
bool IsView(const viewLabel_t viewLabel) const
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
virtual size_t getLocalMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
virtual LocalOrdinal GetStorageBlockSize() const =0
Returns the block size of the storage mechanism, which is usually 1, except for Tpetra::BlockCrsMatri...
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y....
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X....
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
std::string toString(const T &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)