Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_MatrixMatrix_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Xpetra: A linear algebra interface package
4//
5// Copyright 2012 NTESS and the Xpetra contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_DEF_HPP_
11#define PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_DEF_HPP_
12
13#include "Xpetra_ConfigDefs.hpp"
14
15#include "Xpetra_BlockedCrsMatrix.hpp"
16#include "Xpetra_CrsMatrixWrap.hpp"
17#include "Xpetra_MapExtractor.hpp"
18#include "Xpetra_Map.hpp"
19#include "Xpetra_MatrixFactory.hpp"
20#include "Xpetra_Matrix.hpp"
21#include "Xpetra_StridedMapFactory.hpp"
22#include "Xpetra_StridedMap.hpp"
23
24#ifdef HAVE_XPETRA_EPETRA
25#include <Xpetra_EpetraCrsMatrix_fwd.hpp>
26#endif
27
28#ifdef HAVE_XPETRA_EPETRAEXT
29#include <EpetraExt_MatrixMatrix.h>
30#include <EpetraExt_RowMatrixOut.h>
31#include <Epetra_RowMatrixTransposer.h>
32#endif // HAVE_XPETRA_EPETRAEXT
33
34#ifdef HAVE_XPETRA_TPETRA
35#include <TpetraExt_MatrixMatrix.hpp>
36#include <Tpetra_RowMatrixTransposer.hpp>
37#include <MatrixMarket_Tpetra.hpp>
38#include <Xpetra_TpetraCrsMatrix.hpp>
39#include <Xpetra_TpetraBlockCrsMatrix.hpp>
40#include <Tpetra_BlockCrsMatrix_Helpers.hpp>
41#include <Xpetra_TpetraMultiVector.hpp>
42#include <Xpetra_TpetraVector.hpp>
43#endif // HAVE_XPETRA_TPETRA
44
45#include "Xpetra_MatrixMatrix_decl.hpp"
46
47namespace Xpetra {
48
49template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
50void MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(const Matrix& A, bool transposeA,
51 const Matrix& B, bool transposeB,
52 Matrix& C,
53 bool call_FillComplete_on_result,
54 bool doOptimizeStorage,
55 const std::string& label,
56 const RCP<ParameterList>& params) {
57 TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
58 Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
59 TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
60 Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
61
62 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
63 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
64
65 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
66
67 if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
68#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
69 throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
70#else
71 throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
72
73#endif
74 } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
75#ifdef HAVE_XPETRA_TPETRA
76 using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
77 if (helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
78 // All matrices are Crs
79 const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = helpers::Op2TpetraCrs(A);
80 const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = helpers::Op2TpetraCrs(B);
81 Tpetra::CrsMatrix<SC, LO, GO, NO>& tpC = helpers::Op2NonConstTpetraCrs(C);
82
83 // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
84 // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
85 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
86 } else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
87 // All matrices are BlockCrs (except maybe Ac)
88 // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
89 // we'll need to think about refactoring BlockCrs so we can do something smarter here.
90 if (!A.getRowMap()->getComm()->getRank())
91 std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
92
93 const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
94 const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(B);
95 using CRS = Tpetra::CrsMatrix<SC, LO, GO, NO>;
96 RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
97 RCP<const CRS> Bcrs = Tpetra::convertToCrsMatrix(tpB);
98
99 // We need the global constants to do the copy back to BlockCrs
100 RCP<ParameterList> new_params;
101 if (!params.is_null()) {
102 new_params = rcp(new Teuchos::ParameterList(*params));
103 new_params->set("compute global constants", true);
104 }
105
106 // FIXME: The lines below only works because we're assuming Ac is Point
107 RCP<CRS> tempAc = Teuchos::rcp(new CRS(Acrs->getRowMap(), 0));
108 Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
109
110 // Temporary output matrix
111 RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO>> Ac_t = Tpetra::convertToBlockCrsMatrix(*tempAc, A.GetStorageBlockSize());
112 RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>> Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
113 RCP<Xpetra::CrsMatrix<SC, LO, GO, NO>> Ac_p = Ac_x;
114
115 // We can now cheat and replace the innards of Ac
116 RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO>> Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO>>(Teuchos::rcpFromRef(C));
117 Ac_w->replaceCrsMatrix(Ac_p);
118 } else {
119 // Mix and match
120 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
121 }
122#else
123 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
124#endif
125 }
126
127 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
128 RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
129 fillParams->set("Optimize Storage", doOptimizeStorage);
130 C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
131 (transposeA) ? A.getDomainMap() : A.getRangeMap(),
132 fillParams);
133 }
134
135 // transfer striding information
136 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
137 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
138 C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
139} // end Multiply
140
141template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, RCP<Matrix> C_in,
143 Teuchos::FancyOStream& fos,
144 bool doFillComplete,
145 bool doOptimizeStorage,
146 const std::string& label,
147 const RCP<ParameterList>& params) {
148 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
149 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
150
151 // Default case: Xpetra Multiply
152 RCP<Matrix> C = C_in;
153
154 if (C == Teuchos::null) {
155 double nnzPerRow = Teuchos::as<double>(0);
156
157#if 0
158 if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
159 // For now, follow what ML and Epetra do.
160 GO numRowsA = A.getGlobalNumRows();
161 GO numRowsB = B.getGlobalNumRows();
162 nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
163 sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
164 nnzPerRow *= nnzPerRow;
165 double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
166 double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
167 if (totalNnz < minNnz)
168 totalNnz = minNnz;
169 nnzPerRow = totalNnz / A.getGlobalNumRows();
170
171 fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
172 }
173#endif
174 if (transposeA)
175 C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
176 else
177 C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
178
179 } else {
180 C->resumeFill(); // why this is not done inside of Tpetra MxM?
181 fos << "Reuse C pattern" << std::endl;
182 }
183
184 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
185
186 return C;
187}
188
189template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
190RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Teuchos::FancyOStream& fos,
191 bool callFillCompleteOnResult, bool doOptimizeStorage, const std::string& label,
192 const RCP<ParameterList>& params) {
193 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
194}
195
196#ifdef HAVE_XPETRA_EPETRAEXT
197template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
198RCP<Epetra_CrsMatrix> MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
199 const Epetra_CrsMatrix& epB,
200 Teuchos::FancyOStream& fos) {
201 throw(Xpetra::Exceptions::RuntimeError("MLTwoMatrixMultiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
202 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
203}
204#endif // ifdef HAVE_XPETRA_EPETRAEXT
205
206template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
207RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
208 const BlockedCrsMatrix& B, bool transposeB,
209 Teuchos::FancyOStream& fos,
210 bool doFillComplete,
211 bool doOptimizeStorage) {
212 TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
213 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
214
215 // Preconditions
216 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
217 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
218
219 RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
220 RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
221
222 RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
223
224 for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
225 for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
226 RCP<Matrix> Cij;
227
228 for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
229 RCP<Matrix> crmat1 = A.getMatrix(i, l);
230 RCP<Matrix> crmat2 = B.getMatrix(l, j);
231
232 if (crmat1.is_null() || crmat2.is_null())
233 continue;
234
235 // try unwrapping 1x1 blocked matrix
236 {
237 auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
238 auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
239
240 if (unwrap1.is_null() != unwrap2.is_null()) {
241 if (unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
242 crmat1 = unwrap1->getCrsMatrix();
243 if (unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
244 crmat2 = unwrap2->getCrsMatrix();
245 }
246 }
247
248 RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
249 RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
250 TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
251 "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
252
253 // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
254 if (!crop1.is_null())
255 Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
256 if (!crop2.is_null())
257 Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
258
259 TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
260 "crmat1 does not have global constants");
261 TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
262 "crmat2 does not have global constants");
263
264 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
265 continue;
266
267 // temporary matrix containing result of local block multiplication
268 RCP<Matrix> temp = Teuchos::null;
269
270 if (crop1 != Teuchos::null && crop2 != Teuchos::null)
271 temp = Multiply(*crop1, false, *crop2, false, fos);
272 else {
273 RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
274 RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
275 TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
276 TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
277 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)");
278 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)");
279
280 // recursive multiplication call
281 temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
282
283 RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
284 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)");
285 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)");
286 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)");
287 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)");
288 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)");
289 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)");
290 }
291
292 TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
293
294 RCP<Matrix> addRes = null;
295 if (Cij.is_null())
296 Cij = temp;
297 else {
298 MatrixMatrix::TwoMatrixAdd(*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
299 Cij = addRes;
300 }
301 }
302
303 if (!Cij.is_null()) {
304 if (Cij->isFillComplete())
305 Cij->resumeFill();
306 Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
307 C->setMatrix(i, j, Cij);
308 } else {
309 C->setMatrix(i, j, Teuchos::null);
310 }
311 }
312 }
313
314 if (doFillComplete)
315 C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
316
317 return C;
318} // TwoMatrixMultiplyBlock
319
320template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
321void MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
322 if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
323 throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
324
325 if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
326 throw Exceptions::RuntimeError("TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
327 } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
328#ifdef HAVE_XPETRA_TPETRA
329 const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
330 Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(B);
331
332 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
333#else
334 throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
335#endif
336 }
337} // MatrixMatrix::TwoMatrixAdd()
338
339template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
340void MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
341 const Matrix& B, bool transposeB, const SC& beta,
342 RCP<Matrix>& C, Teuchos::FancyOStream& fos, bool AHasFixedNnzPerRow) {
343 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
344 RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
345 RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
346 RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
347 // We have to distinguish 4 cases:
348 // both matrices are CrsMatrixWrap based, only one of them is CrsMatrixWrap based, or none.
349
350 // C can be null, so just use A to get the lib
351 auto lib = A.getRowMap()->lib();
352
353 // Both matrices are CrsMatrixWrap
354 if (rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
355 if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
356 throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
357 if (lib == Xpetra::UseEpetra) {
358 throw Exceptions::RuntimeError("MatrixMatrix::Add for Epetra only available with Scalar = double, LO = GO = int.");
359 } else if (lib == Xpetra::UseTpetra) {
360#ifdef HAVE_XPETRA_TPETRA
361 using tcrs_matrix_type = Tpetra::CrsMatrix<SC, LO, GO, NO>;
362 using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
363 const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
364 const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
365 C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
366#else
367 throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
368#endif
369 }
371 if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
372 if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
374 }
375 // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
376 else if (rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
377 RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
378 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
379
380 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
381 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
382
383 size_t i = 0;
384 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
385 RCP<Matrix> Cij = Teuchos::null;
386 if (rcpA != Teuchos::null &&
387 rcpBopB->getMatrix(i, j) != Teuchos::null) {
388 // recursive call
389 TwoMatrixAdd(*rcpA, transposeA, alpha,
390 *(rcpBopB->getMatrix(i, j)), transposeB, beta,
391 Cij, fos, AHasFixedNnzPerRow);
392 } else if (rcpA == Teuchos::null &&
393 rcpBopB->getMatrix(i, j) != Teuchos::null) {
394 Cij = rcpBopB->getMatrix(i, j);
395 } else if (rcpA != Teuchos::null &&
396 rcpBopB->getMatrix(i, j) == Teuchos::null) {
397 Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
398 } else {
399 Cij = Teuchos::null;
400 }
401
402 if (!Cij.is_null()) {
403 if (Cij->isFillComplete())
404 Cij->resumeFill();
405 Cij->fillComplete();
406 bC->setMatrix(i, j, Cij);
407 } else {
408 bC->setMatrix(i, j, Teuchos::null);
409 }
410 } // loop over columns j
411 }
412 // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
413 else if (rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
414 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
415 RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
416
417 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
418 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
419 size_t j = 0;
420 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
421 RCP<Matrix> Cij = Teuchos::null;
422 if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
423 rcpB != Teuchos::null) {
424 // recursive call
425 TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
426 *rcpB, transposeB, beta,
427 Cij, fos, AHasFixedNnzPerRow);
428 } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
429 rcpB != Teuchos::null) {
430 Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
431 } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
432 rcpB == Teuchos::null) {
433 Cij = rcpBopA->getMatrix(i, j);
434 } else {
435 Cij = Teuchos::null;
436 }
437
438 if (!Cij.is_null()) {
439 if (Cij->isFillComplete())
440 Cij->resumeFill();
441 Cij->fillComplete();
442 bC->setMatrix(i, j, Cij);
443 } else {
444 bC->setMatrix(i, j, Teuchos::null);
445 }
446 } // loop over rows i
447 } else {
448 // This is the version for blocked matrices
449
450 // check the compatibility of the blocked operators
451 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
452 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
453 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)");
454 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)");
455 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)");
456 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)");
457 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)");
458 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)");
459
460 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
461 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
462
463 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
464 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
465 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
466 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
467
468 RCP<Matrix> Cij = Teuchos::null;
469 if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
470 rcpBopB->getMatrix(i, j) != Teuchos::null) {
471 // recursive call
472 TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
473 *(rcpBopB->getMatrix(i, j)), transposeB, beta,
474 Cij, fos, AHasFixedNnzPerRow);
475 } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
476 rcpBopB->getMatrix(i, j) != Teuchos::null) {
477 Cij = rcpBopB->getMatrix(i, j);
478 } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
479 rcpBopB->getMatrix(i, j) == Teuchos::null) {
480 Cij = rcpBopA->getMatrix(i, j);
481 } else {
482 Cij = Teuchos::null;
483 }
484
485 if (!Cij.is_null()) {
486 if (Cij->isFillComplete())
487 Cij->resumeFill();
488 Cij->fillComplete();
489 bC->setMatrix(i, j, Cij);
490 } else {
491 bC->setMatrix(i, j, Teuchos::null);
492 }
493 } // loop over columns j
494 } // loop over rows i
495
496 } // end blocked recursive algorithm
497} // MatrixMatrix::TwoMatrixAdd()
498
499} // end namespace Xpetra
500
501#endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_DEF_HPP_ */
Xpetra_ConfigDefs.hpp
Xpetra_EpetraCrsMatrix_fwd.hpp
Xpetra_MatrixMatrix_decl.hpp
GO
GlobalOrdinal GO
Definition Xpetra_UseShortNamesOrdinal.hpp:103
SC
Scalar SC
Definition Xpetra_UseShortNamesScalar.hpp:150
Epetra_CrsMatrix
Teuchos::RCP::is_null
bool is_null() const
Teuchos::basic_FancyOStream
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition Xpetra_BlockedCrsMatrix_def.hpp:1067
Xpetra::BlockedCrsMatrix::getRangeMap
const RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
Definition Xpetra_BlockedCrsMatrix_def.hpp:849
Xpetra::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition Xpetra_BlockedCrsMatrix_def.hpp:1103
Xpetra::BlockedCrsMatrix::isFillComplete
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode.
Definition Xpetra_BlockedCrsMatrix_def.hpp:532
Xpetra::BlockedCrsMatrix::getDomainMap
const RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
Definition Xpetra_BlockedCrsMatrix_def.hpp:819
Xpetra::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition Xpetra_BlockedCrsMatrix_def.hpp:1061
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition Xpetra_BlockedCrsMatrix_def.hpp:867
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition Xpetra_BlockedCrsMatrix_def.hpp:873
Xpetra::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition Xpetra_Exceptions.hpp:50
Xpetra::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition Xpetra_Exceptions.hpp:72
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition Xpetra_Exceptions.hpp:65
Xpetra::Helpers::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
Definition Xpetra_Helpers_def.hpp:98
Xpetra::Helpers::Op2TpetraBlockCrs
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
Definition Xpetra_Helpers_def.hpp:169
Xpetra::Helpers::Op2TpetraCrs
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Definition Xpetra_Helpers_def.hpp:85
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition Xpetra_MatrixFactory_def.hpp:22
Xpetra::MatrixMatrix::TwoMatrixMultiplyBlock
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".
Definition Xpetra_MatrixMatrix_def.hpp:207
Xpetra::MatrixMatrix::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
Definition Xpetra_MatrixMatrix_def.hpp:198
Xpetra::MatrixMatrix::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 > &params=null)
Definition Xpetra_MatrixMatrix_def.hpp:50
Xpetra::MatrixMatrix::TwoMatrixAdd
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
Definition Xpetra_MatrixMatrix_def.hpp:321
Xpetra::Matrix
Xpetra-specific matrix class.
Definition Xpetra_Matrix_decl.hpp:60
Xpetra::Matrix::getGlobalNumEntries
virtual global_size_t getGlobalNumEntries() const =0
Returns the global number of entries in this matrix.
Xpetra::Matrix::CreateView
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
Definition Xpetra_Matrix_def.hpp:41
Xpetra::Matrix::getGlobalNumRows
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
Xpetra::Matrix::getRowMap
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition Xpetra_Matrix_def.hpp:134
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode.
Xpetra::Matrix::fillComplete
virtual void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > &params=null)=0
Signal that data entry is complete, specifying domain and range maps.
Xpetra::Matrix::IsView
bool IsView(const viewLabel_t viewLabel) const
Definition Xpetra_Matrix_def.hpp:120
Xpetra::Matrix::GetStorageBlockSize
virtual LocalOrdinal GetStorageBlockSize() const =0
Returns the block size of the storage mechanism, which is usually 1, except for Tpetra::BlockCrsMatri...
Xpetra::Operator::getDomainMap
virtual const Teuchos::RCP< const map_type > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X....
Xpetra::Operator::getRangeMap
virtual const Teuchos::RCP< const map_type > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y....
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_UNREACHABLE_RETURN
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)