48#ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_HPP_
49#define PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_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"
63#ifdef HAVE_XPETRA_TPETRA
64#include <TpetraExt_TripleMatrixMultiply.hpp>
65#include <Xpetra_TpetraCrsMatrix.hpp>
66#include <Tpetra_BlockCrsMatrix.hpp>
67#include <Tpetra_BlockCrsMatrix_Helpers.hpp>
74 template <
class Scalar,
79#undef XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
109 const Matrix& A,
bool transposeA,
110 const Matrix& P,
bool transposeP,
112 bool call_FillComplete_on_result =
true,
113 bool doOptimizeStorage =
true,
114 const std::string & label = std::string(),
118 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
120 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
126 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
131#ifdef HAVE_XPETRA_TPETRA
133 if(helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P)) {
142 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
144 else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
150 std::cout<<
"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
157 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
165 const bool do_fill_complete=
true;
166 Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
175 Ac_w->replaceCrsMatrix(Ac_p);
186 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
188 fillParams->set(
"Optimize Storage", doOptimizeStorage);
198 const std::string stridedViewLabel(
"stridedMaps");
199 const size_t blkSize = 1;
200 std::vector<size_t> stridingInfo(1, blkSize);
201 LocalOrdinal stridedBlockId = -1;
203 if (R.
IsView(stridedViewLabel)) {
204 rangeMap = transposeR ? R.
getColMap(stridedViewLabel) : R.
getRowMap(stridedViewLabel);
210 if (P.
IsView(stridedViewLabel)) {
211 domainMap = transposeP ? P.
getRowMap(stridedViewLabel) : P.
getColMap(stridedViewLabel);
216 Ac.
CreateView(stridedViewLabel, rangeMap, domainMap);
222#ifdef HAVE_XPETRA_EPETRA
235 const Matrix& A,
bool transposeA,
236 const Matrix& P,
bool transposeP,
238 bool call_FillComplete_on_result =
true,
239 bool doOptimizeStorage =
true,
240 const std::string & label = std::string(),
244 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
246 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
252 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
257#ifdef HAVE_XPETRA_TPETRA
258# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
259 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
263 if(helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P) && helpers::isTpetraCrs(Ac)) {
272 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
274 else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
279 std::cout<<
"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
286 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
294 const bool do_fill_complete=
true;
295 Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
304 Ac_w->replaceCrsMatrix(Ac_p);
315 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
317 fillParams->set(
"Optimize Storage", doOptimizeStorage);
327 const std::string stridedViewLabel(
"stridedMaps");
328 const size_t blkSize = 1;
329 std::vector<size_t> stridingInfo(1, blkSize);
332 if (R.
IsView(stridedViewLabel)) {
333 rangeMap = transposeR ? R.
getColMap(stridedViewLabel) : R.
getRowMap(stridedViewLabel);
339 if (P.
IsView(stridedViewLabel)) {
340 domainMap = transposeP ? P.
getRowMap(stridedViewLabel) : P.
getColMap(stridedViewLabel);
345 Ac.
CreateView(stridedViewLabel, rangeMap, domainMap);
364 const Matrix& A,
bool transposeA,
365 const Matrix& P,
bool transposeP,
367 bool call_FillComplete_on_result =
true,
368 bool doOptimizeStorage =
true,
369 const std::string & label = std::string(),
373 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
375 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
381 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
386#ifdef HAVE_XPETRA_TPETRA
387# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
388 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
392 if(helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P)) {
401 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
403 else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
408 std::cout<<
"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
415 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
423 const bool do_fill_complete=
true;
424 Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
433 Ac_w->replaceCrsMatrix(Ac_p);
444 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
446 fillParams->set(
"Optimize Storage", doOptimizeStorage);
456 const std::string stridedViewLabel(
"stridedMaps");
457 const size_t blkSize = 1;
458 std::vector<size_t> stridingInfo(1, blkSize);
461 if (R.
IsView(stridedViewLabel)) {
462 rangeMap = transposeR ? R.
getColMap(stridedViewLabel) : R.
getRowMap(stridedViewLabel);
468 if (P.
IsView(stridedViewLabel)) {
469 domainMap = transposeP ? P.
getRowMap(stridedViewLabel) : P.
getColMap(stridedViewLabel);
474 Ac.
CreateView(stridedViewLabel, rangeMap, domainMap);
484#define XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
Exception throws to report errors in the internal logical of the program.
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Xpetra-specific matrix class.
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode.
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
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.
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 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....
static RCP< Xpetra::StridedMap< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
static void MultiplyRAP(const Matrix &R, bool transposeR, const Matrix &A, bool transposeA, const Matrix &P, bool transposeP, Matrix &Ac, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
static void MultiplyRAP(const Matrix &R, bool transposeR, const Matrix &A, bool transposeA, const Matrix &P, bool transposeP, Matrix &Ac, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
static void MultiplyRAP(const Matrix &R, bool transposeR, const Matrix &A, bool transposeA, const Matrix &P, bool transposeP, Matrix &Ac, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)