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"
62 #ifdef HAVE_XPETRA_TPETRA
63 #include <TpetraExt_TripleMatrixMultiply.hpp>
64 #include <Xpetra_TpetraCrsMatrix.hpp>
67 #endif // HAVE_XPETRA_TPETRA
71 template <
class Scalar,
76 #undef XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
106 const Matrix& A,
bool transposeA,
107 const Matrix& P,
bool transposeP,
109 bool call_FillComplete_on_result =
true,
110 bool doOptimizeStorage =
true,
111 const std::string & label = std::string(),
112 const RCP<ParameterList>& params = null) {
114 TEUCHOS_TEST_FOR_EXCEPTION(transposeR ==
false && Ac.
getRowMap()->isSameAs(*R.
getRowMap()) ==
false,
115 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
116 TEUCHOS_TEST_FOR_EXCEPTION(transposeR ==
true && Ac.
getRowMap()->isSameAs(*R.
getDomainMap()) ==
false,
117 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
123 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
128 #ifdef HAVE_XPETRA_TPETRA
136 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
142 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
143 RCP<Teuchos::ParameterList> fillParams = rcp(
new Teuchos::ParameterList());
144 fillParams->set(
"Optimize Storage", doOptimizeStorage);
151 RCP<const Map> domainMap = Teuchos::null;
152 RCP<const Map> rangeMap = Teuchos::null;
154 const std::string stridedViewLabel(
"stridedMaps");
155 const size_t blkSize = 1;
156 std::vector<size_t> stridingInfo(1, blkSize);
157 LocalOrdinal stridedBlockId = -1;
159 if (R.
IsView(stridedViewLabel)) {
160 rangeMap = transposeR ? R.
getColMap(stridedViewLabel) : R.
getRowMap(stridedViewLabel);
166 if (P.
IsView(stridedViewLabel)) {
167 domainMap = transposeP ? P.
getRowMap(stridedViewLabel) : P.
getColMap(stridedViewLabel);
172 Ac.
CreateView(stridedViewLabel, rangeMap, domainMap);
178 #ifdef HAVE_XPETRA_EPETRA
191 const Matrix& A,
bool transposeA,
192 const Matrix& P,
bool transposeP,
194 bool call_FillComplete_on_result =
true,
195 bool doOptimizeStorage =
true,
196 const std::string & label = std::string(),
197 const RCP<ParameterList>& params = null) {
199 TEUCHOS_TEST_FOR_EXCEPTION(transposeR ==
false && Ac.
getRowMap()->isSameAs(*R.
getRowMap()) ==
false,
200 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
201 TEUCHOS_TEST_FOR_EXCEPTION(transposeR ==
true && Ac.
getRowMap()->isSameAs(*R.
getDomainMap()) ==
false,
202 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
208 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
213 #ifdef HAVE_XPETRA_TPETRA
214 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
215 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
225 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
230 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
231 RCP<Teuchos::ParameterList> fillParams = rcp(
new Teuchos::ParameterList());
232 fillParams->set(
"Optimize Storage", doOptimizeStorage);
239 RCP<const Map> domainMap = Teuchos::null;
240 RCP<const Map> rangeMap = Teuchos::null;
242 const std::string stridedViewLabel(
"stridedMaps");
243 const size_t blkSize = 1;
244 std::vector<size_t> stridingInfo(1, blkSize);
247 if (R.
IsView(stridedViewLabel)) {
248 rangeMap = transposeR ? R.
getColMap(stridedViewLabel) : R.
getRowMap(stridedViewLabel);
254 if (P.
IsView(stridedViewLabel)) {
255 domainMap = transposeP ? P.
getRowMap(stridedViewLabel) : P.
getColMap(stridedViewLabel);
260 Ac.
CreateView(stridedViewLabel, rangeMap, domainMap);
279 const Matrix& A,
bool transposeA,
280 const Matrix& P,
bool transposeP,
282 bool call_FillComplete_on_result =
true,
283 bool doOptimizeStorage =
true,
284 const std::string & label = std::string(),
285 const RCP<ParameterList>& params = null) {
287 TEUCHOS_TEST_FOR_EXCEPTION(transposeR ==
false && Ac.
getRowMap()->isSameAs(*R.
getRowMap()) ==
false,
288 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
289 TEUCHOS_TEST_FOR_EXCEPTION(transposeR ==
true && Ac.
getRowMap()->isSameAs(*R.
getDomainMap()) ==
false,
290 Exceptions::RuntimeError,
"XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
296 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
301 #ifdef HAVE_XPETRA_TPETRA
302 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
303 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
313 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
318 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
319 RCP<Teuchos::ParameterList> fillParams = rcp(
new Teuchos::ParameterList());
320 fillParams->set(
"Optimize Storage", doOptimizeStorage);
327 RCP<const Map> domainMap = Teuchos::null;
328 RCP<const Map> rangeMap = Teuchos::null;
330 const std::string stridedViewLabel(
"stridedMaps");
331 const size_t blkSize = 1;
332 std::vector<size_t> stridingInfo(1, blkSize);
335 if (R.
IsView(stridedViewLabel)) {
336 rangeMap = transposeR ? R.
getColMap(stridedViewLabel) : R.
getRowMap(stridedViewLabel);
342 if (P.
IsView(stridedViewLabel)) {
343 domainMap = transposeP ? P.
getRowMap(stridedViewLabel) : P.
getColMap(stridedViewLabel);
348 Ac.
CreateView(stridedViewLabel, rangeMap, domainMap);
358 #define XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
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)
Exception throws to report errors in the internal logical of the program.
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 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)
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
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.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< 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)