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