Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
47 namespace Xpetra {
48 
49 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
50 void 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 
141 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142 RCP<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 
189 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
190 RCP<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
197 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
198 RCP<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 
206 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
207 RCP<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 
320 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
321 void 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 
339 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
340 void 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::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 &quot;in-place&quot;.
Definition: Xpetra_MatrixMatrix_def.hpp:207
Xpetra::Matrix::getGlobalNumRows
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix_def.hpp:873
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::Helpers
Xpetra utility class containing transformation routines between Xpetra::Matrix and Epetra/Tpetra obje...
Definition: Xpetra_Helpers_decl.hpp:56
Xpetra_EpetraCrsMatrix_fwd.hpp
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::resumeFill
virtual void resumeFill(const RCP< ParameterList > &params=null)=0
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::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:65
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition: Xpetra_BlockedCrsMatrix_def.hpp:1067
Xpetra::Helpers::Op2TpetraCrs
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:85
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::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: Xpetra_Exceptions.hpp:50
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::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
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::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::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix_def.hpp:1061
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::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...
Xpetra::Helpers::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:98
Xpetra_MatrixMatrix_decl.hpp
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix_def.hpp:867
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::CrsMatrixWrap::replaceCrsMatrix
void replaceCrsMatrix(RCP< CrsMatrix > &M)
Expert only.
Definition: Xpetra_CrsMatrixWrap_def.hpp:458
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory_def.hpp:22
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_ConfigDefs.hpp
Xpetra::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: Xpetra_Exceptions.hpp:72
Xpetra::Matrix::getGlobalNumEntries
virtual global_size_t getGlobalNumEntries() const =0
Returns the global number of entries in this matrix.
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_decl.hpp:46
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::Helpers::Op2TpetraBlockCrs
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:169
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::BlockedCrsMatrix::getCrsMatrix
Teuchos::RCP< Matrix > getCrsMatrix() const
return unwrap 1x1 blocked operators
Definition: Xpetra_BlockedCrsMatrix_def.hpp:1073
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_decl.hpp:60