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 C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
116} // end Multiply
117
118template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
119RCP<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,
120 Teuchos::FancyOStream& fos,
121 bool doFillComplete,
122 bool doOptimizeStorage,
123 const std::string& label,
124 const RCP<ParameterList>& params) {
125 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
126 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
127
128 // Default case: Xpetra Multiply
129 RCP<Matrix> C = C_in;
130
131 if (C == Teuchos::null) {
132 double nnzPerRow = Teuchos::as<double>(0);
133
134#if 0
135 if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
136 // For now, follow what ML and Epetra do.
137 GO numRowsA = A.getGlobalNumRows();
138 GO numRowsB = B.getGlobalNumRows();
139 nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
140 sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
141 nnzPerRow *= nnzPerRow;
142 double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
143 double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
144 if (totalNnz < minNnz)
145 totalNnz = minNnz;
146 nnzPerRow = totalNnz / A.getGlobalNumRows();
147
148 fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
149 }
150#endif
151 if (transposeA)
152 C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
153 else
154 C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
155
156 } else {
157 C->resumeFill(); // why this is not done inside of Tpetra MxM?
158 fos << "Reuse C pattern" << std::endl;
159 }
160
161 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
162
163 return C;
164}
165
166template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
167RCP<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,
168 bool callFillCompleteOnResult, bool doOptimizeStorage, const std::string& label,
169 const RCP<ParameterList>& params) {
170 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
171}
172
173template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
174RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
175 const BlockedCrsMatrix& B, bool transposeB,
176 Teuchos::FancyOStream& fos,
177 bool doFillComplete,
178 bool doOptimizeStorage) {
179 TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
180 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
181
182 // Preconditions
183 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
184 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
185
186 RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
187 RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
188
189 RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
190
191 for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
192 for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
193 RCP<Matrix> Cij;
194
195 for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
196 RCP<Matrix> crmat1 = A.getMatrix(i, l);
197 RCP<Matrix> crmat2 = B.getMatrix(l, j);
198
199 if (crmat1.is_null() || crmat2.is_null())
200 continue;
201
202 // try unwrapping 1x1 blocked matrix
203 {
204 auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
205 auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
206
207 if (unwrap1.is_null() != unwrap2.is_null()) {
208 if (unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
209 crmat1 = unwrap1->getCrsMatrix();
210 if (unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
211 crmat2 = unwrap2->getCrsMatrix();
212 }
213 }
214
215 RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
216 RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
217 TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
218 "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
219
220 // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
221 if (!crop1.is_null())
222 Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
223 if (!crop2.is_null())
224 Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
225
226 TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
227 "crmat1 does not have global constants");
228 TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
229 "crmat2 does not have global constants");
230
231 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
232 continue;
233
234 // temporary matrix containing result of local block multiplication
235 RCP<Matrix> temp = Teuchos::null;
236
237 if (crop1 != Teuchos::null && crop2 != Teuchos::null)
238 temp = Multiply(*crop1, false, *crop2, false, fos);
239 else {
240 RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
241 RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
242 TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
243 TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
244 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)");
245 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)");
246
247 // recursive multiplication call
248 temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
249
250 RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
251 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)");
252 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)");
253 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)");
254 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)");
255 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)");
256 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)");
257 }
258
259 TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
260
261 RCP<Matrix> addRes = null;
262 if (Cij.is_null())
263 Cij = temp;
264 else {
265 MatrixMatrix::TwoMatrixAdd(*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
266 Cij = addRes;
267 }
268 }
269
270 if (!Cij.is_null()) {
271 if (Cij->isFillComplete())
272 Cij->resumeFill();
273 Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
274 C->setMatrix(i, j, Cij);
275 } else {
276 C->setMatrix(i, j, Teuchos::null);
277 }
278 }
279 }
280
281 if (doFillComplete)
282 C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
283
284 return C;
285} // TwoMatrixMultiplyBlock
286
287template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
288void MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
289 if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
290 throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
291
292 if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
293 const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
294 Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(B);
295
296 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
297 }
298} // MatrixMatrix::TwoMatrixAdd()
299
300template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
301void MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
302 const Matrix& B, bool transposeB, const SC& beta,
303 RCP<Matrix>& C, Teuchos::FancyOStream& fos, bool AHasFixedNnzPerRow) {
304 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
305 RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
306 RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
307 RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
308 // We have to distinguish 4 cases:
309 // both matrices are CrsMatrixWrap based, only one of them is CrsMatrixWrap based, or none.
310
311 // C can be null, so just use A to get the lib
312 auto lib = A.getRowMap()->lib();
313
314 // Both matrices are CrsMatrixWrap
315 if (rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
316 if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
317 throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
318 if (lib == Xpetra::UseTpetra) {
319 using tcrs_matrix_type = Tpetra::CrsMatrix<SC, LO, GO, NO>;
320 using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
321 const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
322 const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
323 C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
324 }
326 if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
327 if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
329 }
330 // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
331 else if (rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
332 RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
333 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
334
335 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
336 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
337
338 size_t i = 0;
339 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
340 RCP<Matrix> Cij = Teuchos::null;
341 if (rcpA != Teuchos::null &&
342 rcpBopB->getMatrix(i, j) != Teuchos::null) {
343 // recursive call
344 TwoMatrixAdd(*rcpA, transposeA, alpha,
345 *(rcpBopB->getMatrix(i, j)), transposeB, beta,
346 Cij, fos, AHasFixedNnzPerRow);
347 } else if (rcpA == Teuchos::null &&
348 rcpBopB->getMatrix(i, j) != Teuchos::null) {
349 Cij = rcpBopB->getMatrix(i, j);
350 } else if (rcpA != Teuchos::null &&
351 rcpBopB->getMatrix(i, j) == Teuchos::null) {
352 Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
353 } else {
354 Cij = Teuchos::null;
355 }
356
357 if (!Cij.is_null()) {
358 if (Cij->isFillComplete())
359 Cij->resumeFill();
360 Cij->fillComplete();
361 bC->setMatrix(i, j, Cij);
362 } else {
363 bC->setMatrix(i, j, Teuchos::null);
364 }
365 } // loop over columns j
366 }
367 // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
368 else if (rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
369 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
370 RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
371
372 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
373 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
374 size_t j = 0;
375 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
376 RCP<Matrix> Cij = Teuchos::null;
377 if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
378 rcpB != Teuchos::null) {
379 // recursive call
380 TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
381 *rcpB, transposeB, beta,
382 Cij, fos, AHasFixedNnzPerRow);
383 } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
384 rcpB != Teuchos::null) {
385 Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
386 } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
387 rcpB == Teuchos::null) {
388 Cij = rcpBopA->getMatrix(i, j);
389 } else {
390 Cij = Teuchos::null;
391 }
392
393 if (!Cij.is_null()) {
394 if (Cij->isFillComplete())
395 Cij->resumeFill();
396 Cij->fillComplete();
397 bC->setMatrix(i, j, Cij);
398 } else {
399 bC->setMatrix(i, j, Teuchos::null);
400 }
401 } // loop over rows i
402 } else {
403 // This is the version for blocked matrices
404
405 // check the compatibility of the blocked operators
406 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
407 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
408 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)");
409 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)");
410 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)");
411 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)");
412 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)");
413 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)");
414
415 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
416 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
417
418 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
419 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
420 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
421 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
422
423 RCP<Matrix> Cij = Teuchos::null;
424 if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
425 rcpBopB->getMatrix(i, j) != Teuchos::null) {
426 // recursive call
427 TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
428 *(rcpBopB->getMatrix(i, j)), transposeB, beta,
429 Cij, fos, AHasFixedNnzPerRow);
430 } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
431 rcpBopB->getMatrix(i, j) != Teuchos::null) {
432 Cij = rcpBopB->getMatrix(i, j);
433 } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
434 rcpBopB->getMatrix(i, j) == Teuchos::null) {
435 Cij = rcpBopA->getMatrix(i, j);
436 } else {
437 Cij = Teuchos::null;
438 }
439
440 if (!Cij.is_null()) {
441 if (Cij->isFillComplete())
442 Cij->resumeFill();
443 Cij->fillComplete();
444 bC->setMatrix(i, j, Cij);
445 } else {
446 bC->setMatrix(i, j, Teuchos::null);
447 }
448 } // loop over columns j
449 } // loop over rows i
450
451 } // end blocked recursive algorithm
452} // MatrixMatrix::TwoMatrixAdd()
453
454} // end namespace Xpetra
455
456#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:174
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:288
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)