Xpetra_MatrixMatrix.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 
48 #ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_
49 #define PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_
50 
51 #include "Xpetra_ConfigDefs.hpp"
52 
53 #include "Xpetra_BlockedCrsMatrix.hpp"
54 #include "Xpetra_CrsMatrixWrap.hpp"
55 #include "Xpetra_MapExtractor.hpp"
56 #include "Xpetra_Map.hpp"
57 #include "Xpetra_MatrixFactory.hpp"
58 #include "Xpetra_Matrix.hpp"
59 #include "Xpetra_StridedMapFactory.hpp"
60 #include "Xpetra_StridedMap.hpp"
61 
62 #ifdef HAVE_XPETRA_EPETRA
63 #include <Xpetra_EpetraCrsMatrix_fwd.hpp>
64 #endif
65 
66 #ifdef HAVE_XPETRA_EPETRAEXT
67 #include <EpetraExt_MatrixMatrix.h>
68 #include <EpetraExt_RowMatrixOut.h>
69 #include <Epetra_RowMatrixTransposer.h>
70 #endif // HAVE_XPETRA_EPETRAEXT
71 
72 #ifdef HAVE_XPETRA_TPETRA
73 #include <TpetraExt_MatrixMatrix.hpp>
74 #include <MatrixMarket_Tpetra.hpp>
75 #include <Xpetra_TpetraCrsMatrix.hpp>
76 #include <Xpetra_TpetraMultiVector.hpp>
77 #include <Xpetra_TpetraVector.hpp>
78 #endif // HAVE_XPETRA_TPETRA
79 
80 namespace Xpetra {
81 
88  template <class Scalar,
89  class LocalOrdinal = int,
90  class GlobalOrdinal = LocalOrdinal,
91  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
92  class Helpers {
93 #include "Xpetra_UseShortNames.hpp"
94 
95  public:
96 
97 #ifdef HAVE_XPETRA_EPETRA
98  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<Matrix> Op) {
99  // Get the underlying Epetra Mtx
100  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
101  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast,
102  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
103 
104  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
105  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
106  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast,
107  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
108 
109  return tmp_ECrsMtx->getEpetra_CrsMatrix();
110  }
111 
112  static RCP<Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) {
113  RCP<Epetra_CrsMatrix> A;
114  // Get the underlying Epetra Mtx
115  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
116  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast,
117  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
118 
119  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
120  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
121  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
122 
123  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
124  }
125 
126  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) {
127  // Get the underlying Epetra Mtx
128  try {
129  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
130  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
131  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
132  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast,
133  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
134 
135  return *tmp_ECrsMtx->getEpetra_CrsMatrix();
136 
137  } catch(...) {
138  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
139  }
140  }
141 
142  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(const Matrix& Op) {
143  RCP<Epetra_CrsMatrix> A;
144  // Get the underlying Epetra Mtx
145  try {
146  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
147  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
148  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
149  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
150 
151  return *Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
152 
153  } catch(...) {
154  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
155  }
156  }
157 #endif // HAVE_XPETRA_EPETRA
158 
159 #ifdef HAVE_XPETRA_TPETRA
160  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<Matrix> Op) {
161  // Get the underlying Tpetra Mtx
162  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
163  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
164 
165  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
166  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
167  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
168 
169  return tmp_ECrsMtx->getTpetra_CrsMatrix();
170  }
171 
172  static RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) {
173  // Get the underlying Tpetra Mtx
174  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
175  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
176  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
177  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
178  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
179 
180  return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
181  }
182 
183  static const Tpetra::CrsMatrix<SC,LO,GO,NO> & Op2TpetraCrs(const Matrix& Op) {
184  // Get the underlying Tpetra Mtx
185  try{
186  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
187  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
188  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
189  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
190 
191  return *tmp_TCrsMtx->getTpetra_CrsMatrix();
192 
193  } catch (...) {
194  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
195  }
196  }
197 
198  static Tpetra::CrsMatrix<SC,LO,GO,NO> & Op2NonConstTpetraCrs(const Matrix& Op) {
199  // Get the underlying Tpetra Mtx
200  try{
201  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap& >(Op);
202  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
203  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
204  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
205 
206  return *Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_TCrsMtx->getTpetra_CrsMatrix());
207 
208  } catch (...) {
209  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
210  }
211  }
212 #endif // HAVE_XPETRA_TPETRA
213 
214  };
215 
216  template <class Scalar,
217  class LocalOrdinal /*= int*/,
218  class GlobalOrdinal /*= LocalOrdinal*/,
219  class Node /*= KokkosClassic::DefaultNode::DefaultNodeType*/>
220  class MatrixMatrix {
221 #undef XPETRA_MATRIXMATRIX_SHORT
222 #include "Xpetra_UseShortNames.hpp"
223 
224  public:
225 
250  static void Multiply(const Matrix& A, bool transposeA,
251  const Matrix& B, bool transposeB,
252  Matrix& C,
253  bool call_FillComplete_on_result = true,
254  bool doOptimizeStorage = true,
255  const std::string & label = std::string(),
256  const RCP<ParameterList>& params = null) {
257 
258  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
259  Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
260  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
261  Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
262 
263  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
264  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
265 
266  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
267 
268  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
269 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
270  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
271 #else
272  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
273 
274 #endif
275  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
276 #ifdef HAVE_XPETRA_TPETRA
277  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
278  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
279  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
280 
281  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
282  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
283  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
284 #else
285  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
286 #endif
287  }
288 
289  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
290  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
291  fillParams->set("Optimize Storage", doOptimizeStorage);
292  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
293  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
294  fillParams);
295  }
296 
297  // transfer striding information
298  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
299  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
300  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
301  } // end Multiply
302 
325  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, RCP<Matrix> C_in,
326  Teuchos::FancyOStream& fos,
327  bool doFillComplete = true,
328  bool doOptimizeStorage = true,
329  const std::string & label = std::string(),
330  const RCP<ParameterList>& params = null) {
331 
332  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
333  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
334 
335  // Default case: Xpetra Multiply
336  RCP<Matrix> C = C_in;
337 
338  if (C == Teuchos::null) {
339  double nnzPerRow = Teuchos::as<double>(0);
340 
341 #if 0
342  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
343  // For now, follow what ML and Epetra do.
344  GO numRowsA = A.getGlobalNumRows();
345  GO numRowsB = B.getGlobalNumRows();
346  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
347  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
348  nnzPerRow *= nnzPerRow;
349  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
350  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
351  if (totalNnz < minNnz)
352  totalNnz = minNnz;
353  nnzPerRow = totalNnz / A.getGlobalNumRows();
354 
355  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
356  }
357 #endif
358  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
359  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
360 
361  } else {
362  C->resumeFill(); // why this is not done inside of Tpetra MxM?
363  fos << "Reuse C pattern" << std::endl;
364  }
365 
366  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
367 
368  return C;
369  }
370 
381  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Teuchos::FancyOStream &fos,
382  bool callFillCompleteOnResult = true, bool doOptimizeStorage = true, const std::string& label = std::string(),
383  const RCP<ParameterList>& params = null) {
384  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
385  }
386 
387 #ifdef HAVE_XPETRA_EPETRAEXT
388  // Michael Gee's MLMultiply
389  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
390  const Epetra_CrsMatrix& epB,
391  Teuchos::FancyOStream& fos) {
392  throw(Xpetra::Exceptions::RuntimeError("MLTwoMatrixMultiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
393  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
394  }
395 #endif //ifdef HAVE_XPETRA_EPETRAEXT
396 
407  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
408  const BlockedCrsMatrix& B, bool transposeB,
409  Teuchos::FancyOStream& fos,
410  bool doFillComplete = true,
411  bool doOptimizeStorage = true) {
412  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
413  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
414 
415  // Preconditions
416  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
417  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
418 
419  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
420  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
421 
422  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
423 
424  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
425  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
426  RCP<Matrix> Cij;
427 
428  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
429  RCP<Matrix> crmat1 = A.getMatrix(i,l);
430  RCP<Matrix> crmat2 = B.getMatrix(l,j);
431 
432  if (crmat1.is_null() || crmat2.is_null())
433  continue;
434 
435  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
436  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
437  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
438  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
439 
440  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
441  if (!crop1.is_null())
442  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
443  if (!crop2.is_null())
444  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
445 
446  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
447  "crmat1 does not have global constants");
448  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
449  "crmat2 does not have global constants");
450 
451  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
452  continue;
453 
454  // temporary matrix containing result of local block multiplication
455  RCP<Matrix> temp = Teuchos::null;
456 
457  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
458  temp = Multiply (*crop1, false, *crop2, false, fos);
459  else {
460  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
461  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
462  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
463  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
464  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)");
465  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)");
466 
467  // recursive multiplication call
468  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
469 
470  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
471  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)");
472  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)");
473  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)");
474  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)");
475  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)");
476  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)");
477  }
478 
479  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
480 
481  RCP<Matrix> addRes = null;
482  if (Cij.is_null ())
483  Cij = temp;
484  else {
485  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
486  Cij = addRes;
487  }
488  }
489 
490  if (!Cij.is_null()) {
491  if (Cij->isFillComplete())
492  Cij->resumeFill();
493  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
494  C->setMatrix(i, j, Cij);
495  } else {
496  C->setMatrix(i, j, Teuchos::null);
497  }
498  }
499  }
500 
501  if (doFillComplete)
502  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
503 
504  return C;
505  } // TwoMatrixMultiplyBlock
506 
519  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
520  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
521  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
522 
523  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
524  throw Exceptions::RuntimeError("TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
525  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
526 #ifdef HAVE_XPETRA_TPETRA
527  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
528  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
529 
530  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
531 #else
532  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
533 #endif
534  }
535  } //MatrixMatrix::TwoMatrixAdd()
536 
537 
550  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
551  const Matrix& B, bool transposeB, const SC& beta,
552  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
553 
554  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
555  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
556  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
557  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
558 
559  // We have to distinguish 4 cases:
560  // both matrices are CrsMatrixWrap based, only one of them is CrsMatrixWrap based, or none.
561 
562  // Both matrices are CrsMatrixWrap
563  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
564 
565  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
566  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
567 
568 
569  if (C == Teuchos::null) {
570  size_t maxNzInA = 0;
571  size_t maxNzInB = 0;
572  size_t numLocalRows = 0;
573  if (A.isFillComplete() && B.isFillComplete()) {
574  maxNzInA = A.getNodeMaxNumRowEntries();
575  maxNzInB = B.getNodeMaxNumRowEntries();
576  numLocalRows = A.getNodeNumRows();
577  }
578  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
579  // first check if either A or B has at most 1 nonzero per row
580  // the case of both having at most 1 nz per row is handled by the ``else''
581  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
582 
583  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
584  for (size_t i = 0; i < numLocalRows; ++i)
585  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
586 
587  } else {
588  for (size_t i = 0; i < numLocalRows; ++i)
589  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
590  }
591 
592  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
593  << ", using static profiling" << std::endl;
594  C = rcp(new CrsMatrixWrap(A.getRowMap(), exactNnzPerRow, Xpetra::StaticProfile));
595  } else {
596  // general case
597  double nnzPerRowInA = Teuchos::as<double>(A.getNodeNumEntries()) / A.getNodeNumRows();
598  double nnzPerRowInB = Teuchos::as<double>(B.getNodeNumEntries()) / B.getNodeNumRows();
599  LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
600 
601  LO maxPossible = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
602  //Use static profiling (more efficient) if the estimate is at least as big as the max
603  //possible nnz's in any single row of the result.
604  Xpetra::ProfileType pft = (maxPossible) > nnzToAllocate ? Xpetra::DynamicProfile : Xpetra::StaticProfile;
605 
606  fos << "nnzPerRowInA = " << nnzPerRowInA << ", nnzPerRowInB = " << nnzPerRowInB << std::endl;
607  fos << "MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
608  << ", max possible nnz per row in sum = " << maxPossible
609  << ", using " << (pft == Xpetra::DynamicProfile ? "dynamic" : "static" ) << " profiling"
610  << std::endl;
611  C = rcp(new CrsMatrixWrap(A.getRowMap(), nnzToAllocate, pft));
612  }
613  if (transposeB)
614  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
615  }
616 
617  if (C->getRowMap()->lib() == Xpetra::UseEpetra) {
618  throw Exceptions::RuntimeError("MatrixMatrix::Add for Epetra only available with Scalar = double, LO = GO = int.");
619  } else if (C->getRowMap()->lib() == Xpetra::UseTpetra) {
620  #ifdef HAVE_XPETRA_TPETRA
621  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
622  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
623  RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
624 
625  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
626  #else
627  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
628  #endif
629  }
631  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
632  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
634  }
635  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
636  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
637  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
638  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
639 
640  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
641  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
642 
643  size_t i = 0;
644  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
645  RCP<Matrix> Cij = Teuchos::null;
646  if(rcpA != Teuchos::null &&
647  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
648  // recursive call
649  TwoMatrixAdd(*rcpA, transposeA, alpha,
650  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
651  Cij, fos, AHasFixedNnzPerRow);
652  } else if (rcpA == Teuchos::null &&
653  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
654  Cij = rcpBopB->getMatrix(i,j);
655  } else if (rcpA != Teuchos::null &&
656  rcpBopB->getMatrix(i,j)==Teuchos::null) {
657  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
658  } else {
659  Cij = Teuchos::null;
660  }
661 
662  if (!Cij.is_null()) {
663  if (Cij->isFillComplete())
664  Cij->resumeFill();
665  Cij->fillComplete();
666  bC->setMatrix(i, j, Cij);
667  } else {
668  bC->setMatrix(i, j, Teuchos::null);
669  }
670  } // loop over columns j
671  }
672  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
673  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
674  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
675  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
676 
677  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
678  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
679  size_t j = 0;
680  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
681  RCP<Matrix> Cij = Teuchos::null;
682  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
683  rcpB!=Teuchos::null) {
684  // recursive call
685  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
686  *rcpB, transposeB, beta,
687  Cij, fos, AHasFixedNnzPerRow);
688  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
689  rcpB!=Teuchos::null) {
690  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
691  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
692  rcpB==Teuchos::null) {
693  Cij = rcpBopA->getMatrix(i,j);
694  } else {
695  Cij = Teuchos::null;
696  }
697 
698  if (!Cij.is_null()) {
699  if (Cij->isFillComplete())
700  Cij->resumeFill();
701  Cij->fillComplete();
702  bC->setMatrix(i, j, Cij);
703  } else {
704  bC->setMatrix(i, j, Teuchos::null);
705  }
706  } // loop over rows i
707  }
708  else {
709  // This is the version for blocked matrices
710 
711  // check the compatibility of the blocked operators
712  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
713  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
714  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)");
715  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)");
716  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)");
717  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)");
718  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)");
719  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)");
720 
721  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
722  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
723 
724  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
725  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
726  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
727  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
728 
729  RCP<Matrix> Cij = Teuchos::null;
730  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
731  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
732  // recursive call
733  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
734  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
735  Cij, fos, AHasFixedNnzPerRow);
736  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
737  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
738  Cij = rcpBopB->getMatrix(i,j);
739  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
740  rcpBopB->getMatrix(i,j)==Teuchos::null) {
741  Cij = rcpBopA->getMatrix(i,j);
742  } else {
743  Cij = Teuchos::null;
744  }
745 
746  if (!Cij.is_null()) {
747  if (Cij->isFillComplete())
748  Cij->resumeFill();
749  Cij->fillComplete();
750  bC->setMatrix(i, j, Cij);
751  } else {
752  bC->setMatrix(i, j, Teuchos::null);
753  }
754  } // loop over columns j
755  } // loop over rows i
756 
757  } // end blocked recursive algorithm
758  } //MatrixMatrix::TwoMatrixAdd()
759 
760 
761  }; // class MatrixMatrix
762 
763 
764 #ifdef HAVE_XPETRA_EPETRA
765  // specialization MatrixMatrix for SC=double, LO=GO=int
766  template<>
767  class MatrixMatrix<double,int,int,EpetraNode> {
768  typedef double Scalar;
769  typedef int LocalOrdinal;
770  typedef int GlobalOrdinal;
771  typedef EpetraNode Node;
772 #include "Xpetra_UseShortNames.hpp"
773 
774  public:
775 
800  static void Multiply(const Matrix& A, bool transposeA,
801  const Matrix& B, bool transposeB,
802  Matrix& C,
803  bool call_FillComplete_on_result = true,
804  bool doOptimizeStorage = true,
805  const std::string & label = std::string(),
806  const RCP<ParameterList>& params = null) {
807  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
808  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
809  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
810  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
811 
812  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
813  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
814 
815  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
816 
817  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
818 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
819  Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(A);
820  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
821  Epetra_CrsMatrix& epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
822 
823  int i = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, epC, haveMultiplyDoFillComplete);
824  if (haveMultiplyDoFillComplete) {
825  // Due to Epetra wrapper intricacies, we need to explicitly call
826  // fillComplete on Xpetra matrix here. Specifically, EpetraCrsMatrix
827  // only keeps an internal variable to check whether we are in resumed
828  // state or not, but never touches the underlying Epetra object. As
829  // such, we need to explicitly update the state of Xpetra matrix to
830  // that of Epetra one afterwords
831  C.fillComplete();
832  }
833 
834  if (i != 0) {
835  std::ostringstream buf;
836  buf << i;
837  std::string msg = "EpetraExt::MatrixMatrix::Multiply return value of " + buf.str();
838  throw(Exceptions::RuntimeError(msg));
839  }
840 
841 #else
842  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
843 #endif
844  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
845 #ifdef HAVE_XPETRA_TPETRA
846 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
847  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
848  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
849 # else
850  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
851  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
852  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
853 
854  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
855  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
856  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
857 # endif
858 #else
859  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
860 #endif
861  }
862 
863  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
864  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
865  fillParams->set("Optimize Storage",doOptimizeStorage);
866  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
867  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
868  fillParams);
869  }
870 
871  // transfer striding information
872  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
873  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
874  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
875  } // end Multiply
876 
899  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
900  const Matrix& B, bool transposeB,
901  RCP<Matrix> C_in,
902  Teuchos::FancyOStream& fos,
903  bool doFillComplete = true,
904  bool doOptimizeStorage = true,
905  const std::string & label = std::string(),
906  const RCP<ParameterList>& params = null) {
907 
908  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
909  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
910 
911  // Optimization using ML Multiply when available and requested
912  // This feature is currently not supported. We would have to introduce the HAVE_XPETRA_ML_MMM flag
913 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) && defined(HAVE_XPETRA_ML_MMM)
914  if (B.getDomainMap()->lib() == Xpetra::UseEpetra && !transposeA && !transposeB) {
915  RCP<const Epetra_CrsMatrix> epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(rcpFromRef(A));
916  RCP<const Epetra_CrsMatrix> epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(rcpFromRef(B));
917  RCP<Epetra_CrsMatrix> epC = MLTwoMatrixMultiply(*epA, *epB, fos);
918 
919  RCP<Matrix> C = Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO> (epC);
920  if (doFillComplete) {
921  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
922  fillParams->set("Optimize Storage", doOptimizeStorage);
923  C->fillComplete(B.getDomainMap(), A.getRangeMap(), fillParams);
924  }
925 
926  // Fill strided maps information
927  // This is necessary since the ML matrix matrix multiplication routine has no handling for this
928  // TODO: move this call to MLMultiply...
929  C->CreateView("stridedMaps", rcpFromRef(A), transposeA, rcpFromRef(B), transposeB);
930 
931  return C;
932  }
933 #endif // EPETRA + EPETRAEXT + ML
934 
935  // Default case: Xpetra Multiply
936  RCP<Matrix> C = C_in;
937 
938  if (C == Teuchos::null) {
939  double nnzPerRow = Teuchos::as<double>(0);
940 
941 #if 0
942  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
943  // For now, follow what ML and Epetra do.
944  GO numRowsA = A.getGlobalNumRows();
945  GO numRowsB = B.getGlobalNumRows();
946  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
947  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
948  nnzPerRow *= nnzPerRow;
949  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
950  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
951  if (totalNnz < minNnz)
952  totalNnz = minNnz;
953  nnzPerRow = totalNnz / A.getGlobalNumRows();
954 
955  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
956  }
957 #endif
958 
959  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
960  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
961 
962  } else {
963  C->resumeFill(); // why this is not done inside of Tpetra MxM?
964  fos << "Reuse C pattern" << std::endl;
965  }
966 
967  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
968 
969  return C;
970  }
971 
982  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
983  const Matrix& B, bool transposeB,
984  Teuchos::FancyOStream &fos,
985  bool callFillCompleteOnResult = true,
986  bool doOptimizeStorage = true,
987  const std::string & label = std::string(),
988  const RCP<ParameterList>& params = null) {
989  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
990  }
991 
992 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
993  // Michael Gee's MLMultiply
994  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
995  const Epetra_CrsMatrix& epB,
996  Teuchos::FancyOStream& fos) {
997 #if defined(HAVE_XPETRA_ML_MMM) // Note: this is currently not supported
998  ML_Comm* comm;
999  ML_Comm_Create(&comm);
1000  fos << "****** USING ML's MATRIX MATRIX MULTIPLY ******" << std::endl;
1001 #ifdef HAVE_MPI
1002  // ML_Comm uses MPI_COMM_WORLD, so try to use the same communicator as epA.
1003  const Epetra_MpiComm * Mcomm = dynamic_cast<const Epetra_MpiComm*>(&(epA.Comm()));
1004  if (Mcomm)
1005  ML_Comm_Set_UsrComm(comm,Mcomm->GetMpiComm());
1006 # endif
1007  //in order to use ML, there must be no indices missing from the matrix column maps.
1008  EpetraExt::CrsMatrix_SolverMap Atransform;
1009  EpetraExt::CrsMatrix_SolverMap Btransform;
1010  const Epetra_CrsMatrix& A = Atransform(const_cast<Epetra_CrsMatrix&>(epA));
1011  const Epetra_CrsMatrix& B = Btransform(const_cast<Epetra_CrsMatrix&>(epB));
1012 
1013  if (!A.Filled()) throw Exceptions::RuntimeError("A has to be FillCompleted");
1014  if (!B.Filled()) throw Exceptions::RuntimeError("B has to be FillCompleted");
1015 
1016  // create ML operators from EpetraCrsMatrix
1017  ML_Operator* ml_As = ML_Operator_Create(comm);
1018  ML_Operator* ml_Bs = ML_Operator_Create(comm);
1019  ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&A),ml_As); // Should we test if the lightweight wrapper is actually used or if WrapEpetraCrsMatrix fall backs to the heavy one?
1020  ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&B),ml_Bs);
1021  ML_Operator* ml_AtimesB = ML_Operator_Create(comm);
1022  {
1023  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("ML_2matmult kernel"));
1024  ML_2matmult(ml_As,ml_Bs,ml_AtimesB,ML_CSR_MATRIX); // do NOT use ML_EpetraCRS_MATRIX!!!
1025  }
1026  ML_Operator_Destroy(&ml_As);
1027  ML_Operator_Destroy(&ml_Bs);
1028 
1029  // For ml_AtimesB we have to reconstruct the column map in global indexing,
1030  // The following is going down to the salt-mines of ML ...
1031  // note: we use integers, since ML only works for Epetra...
1032  int N_local = ml_AtimesB->invec_leng;
1033  ML_CommInfoOP* getrow_comm = ml_AtimesB->getrow->pre_comm;
1034  if (!getrow_comm) throw(Exceptions::RuntimeError("ML_Operator does not have a CommInfo"));
1035  ML_Comm* comm_AB = ml_AtimesB->comm; // get comm object
1036  if (N_local != B.DomainMap().NumMyElements())
1037  throw(Exceptions::RuntimeError("Mismatch in local row dimension between ML and Epetra"));
1038  int N_rcvd = 0;
1039  int N_send = 0;
1040  int flag = 0;
1041  for (int i = 0; i < getrow_comm->N_neighbors; i++) {
1042  N_rcvd += (getrow_comm->neighbors)[i].N_rcv;
1043  N_send += (getrow_comm->neighbors)[i].N_send;
1044  if ( ((getrow_comm->neighbors)[i].N_rcv != 0) &&
1045  ((getrow_comm->neighbors)[i].rcv_list != NULL) ) flag = 1;
1046  }
1047  // For some unknown reason, ML likes to have stuff one larger than
1048  // neccessary...
1049  std::vector<double> dtemp(N_local + N_rcvd + 1); // "double" vector for comm function
1050  std::vector<int> cmap (N_local + N_rcvd + 1); // vector for gids
1051  for (int i = 0; i < N_local; ++i) {
1052  cmap[i] = B.DomainMap().GID(i);
1053  dtemp[i] = (double) cmap[i];
1054  }
1055  ML_cheap_exchange_bdry(&dtemp[0],getrow_comm,N_local,N_send,comm_AB); // do communication
1056  if (flag) { // process received data
1057  int count = N_local;
1058  const int neighbors = getrow_comm->N_neighbors;
1059  for (int i = 0; i < neighbors; i++) {
1060  const int nrcv = getrow_comm->neighbors[i].N_rcv;
1061  for (int j = 0; j < nrcv; j++)
1062  cmap[getrow_comm->neighbors[i].rcv_list[j]] = (int) dtemp[count++];
1063  }
1064  } else {
1065  for (int i = 0; i < N_local+N_rcvd; ++i)
1066  cmap[i] = (int)dtemp[i];
1067  }
1068  dtemp.clear(); // free double array
1069 
1070  // we can now determine a matching column map for the result
1071  Epetra_Map gcmap(-1, N_local+N_rcvd, &cmap[0], B.ColMap().IndexBase(), A.Comm());
1072 
1073  int allocated = 0;
1074  int rowlength;
1075  double* val = NULL;
1076  int* bindx = NULL;
1077 
1078  const int myrowlength = A.RowMap().NumMyElements();
1079  const Epetra_Map& rowmap = A.RowMap();
1080 
1081  // Determine the maximum bandwith for the result matrix.
1082  // replaces the old, very(!) memory-consuming guess:
1083  // int guessnpr = A.MaxNumEntries()*B.MaxNumEntries();
1084  int educatedguess = 0;
1085  for (int i = 0; i < myrowlength; ++i) {
1086  // get local row
1087  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1088  if (rowlength>educatedguess)
1089  educatedguess = rowlength;
1090  }
1091 
1092  // allocate our result matrix and fill it
1093  RCP<Epetra_CrsMatrix> result = rcp(new Epetra_CrsMatrix(::Copy, A.RangeMap(), gcmap, educatedguess, false));
1094 
1095  std::vector<int> gcid(educatedguess);
1096  for (int i = 0; i < myrowlength; ++i) {
1097  const int grid = rowmap.GID(i);
1098  // get local row
1099  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1100  if (!rowlength) continue;
1101  if ((int)gcid.size() < rowlength) gcid.resize(rowlength);
1102  for (int j = 0; j < rowlength; ++j) {
1103  gcid[j] = gcmap.GID(bindx[j]);
1104  if (gcid[j] < 0)
1105  throw Exceptions::RuntimeError("Error: cannot find gcid!");
1106  }
1107  int err = result->InsertGlobalValues(grid, rowlength, val, &gcid[0]);
1108  if (err != 0 && err != 1) {
1109  std::ostringstream errStr;
1110  errStr << "Epetra_CrsMatrix::InsertGlobalValues returned err=" << err;
1111  throw Exceptions::RuntimeError(errStr.str());
1112  }
1113  }
1114  // free memory
1115  if (bindx) ML_free(bindx);
1116  if (val) ML_free(val);
1117  ML_Operator_Destroy(&ml_AtimesB);
1118  ML_Comm_Destroy(&comm);
1119 
1120  return result;
1121 #else // no MUELU_ML
1122  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1123  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1124  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1125 #endif
1126  }
1127 #endif //ifdef HAVE_XPETRA_EPETRAEXT
1128 
1139  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
1140  const BlockedCrsMatrix& B, bool transposeB,
1141  Teuchos::FancyOStream& fos,
1142  bool doFillComplete = true,
1143  bool doOptimizeStorage = true) {
1144  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1145  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1146 
1147  // Preconditions
1148  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1149  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1150 
1151  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1152  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
1153 
1154  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1155 
1156  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1157  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1158  RCP<Matrix> Cij;
1159 
1160  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1161  RCP<Matrix> crmat1 = A.getMatrix(i,l);
1162  RCP<Matrix> crmat2 = B.getMatrix(l,j);
1163 
1164  if (crmat1.is_null() || crmat2.is_null())
1165  continue;
1166 
1167  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1168  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1169  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
1170  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1171 
1172  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
1173  if (!crop1.is_null())
1174  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
1175  if (!crop2.is_null())
1176  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
1177 
1178  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
1179  "crmat1 does not have global constants");
1180  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
1181  "crmat2 does not have global constants");
1182 
1183  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1184  continue;
1185 
1186 
1187  // temporary matrix containing result of local block multiplication
1188  RCP<Matrix> temp = Teuchos::null;
1189 
1190  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1191  temp = Multiply (*crop1, false, *crop2, false, fos);
1192  else {
1193  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1194  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1195  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1196  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1197  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)");
1198  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)");
1199 
1200  // recursive multiplication call
1201  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1202 
1203  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1204  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)");
1205  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)");
1206  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)");
1207  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)");
1208  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)");
1209  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)");
1210  }
1211 
1212  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1213 
1214  RCP<Matrix> addRes = null;
1215  if (Cij.is_null ())
1216  Cij = temp;
1217  else {
1218  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1219  Cij = addRes;
1220  }
1221  }
1222 
1223  if (!Cij.is_null()) {
1224  if (Cij->isFillComplete())
1225  Cij->resumeFill();
1226  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1227  C->setMatrix(i, j, Cij);
1228  } else {
1229  C->setMatrix(i, j, Teuchos::null);
1230  }
1231  }
1232  }
1233 
1234  if (doFillComplete)
1235  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1236 
1237  return C;
1238  } // TwoMatrixMultiplyBlock
1239 
1252  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1253  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*B.getRowMap())), Exceptions::Incompatible,
1254  "TwoMatrixAdd: matrix row maps are not the same.");
1255 
1256  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1257 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1258  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1259  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1260 
1261  //FIXME is there a bug if beta=0?
1262  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1263 
1264  if (rv != 0)
1265  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1266  std::ostringstream buf;
1267 #else
1268  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1269 #endif
1270  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1271 #ifdef HAVE_XPETRA_TPETRA
1272 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1273  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1274  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1275 # else
1276  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1277  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1278 
1279  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1280 # endif
1281 #else
1282  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1283 #endif
1284  }
1285  } //MatrixMatrix::TwoMatrixAdd()
1286 
1287 
1300  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1301  const Matrix& B, bool transposeB, const SC& beta,
1302  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1303  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1304  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1305  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1306  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1307 
1308  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1309 
1310 
1311  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
1312  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
1313 
1314  if (C == Teuchos::null) {
1315  size_t maxNzInA = 0;
1316  size_t maxNzInB = 0;
1317  size_t numLocalRows = 0;
1318  if (A.isFillComplete() && B.isFillComplete()) {
1319 
1320  maxNzInA = A.getNodeMaxNumRowEntries();
1321  maxNzInB = B.getNodeMaxNumRowEntries();
1322  numLocalRows = A.getNodeNumRows();
1323  }
1324 
1325  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1326  // first check if either A or B has at most 1 nonzero per row
1327  // the case of both having at most 1 nz per row is handled by the ``else''
1328  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1329 
1330  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1331  for (size_t i = 0; i < numLocalRows; ++i)
1332  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1333 
1334  } else {
1335  for (size_t i = 0; i < numLocalRows; ++i)
1336  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1337  }
1338 
1339  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1340  << ", using static profiling" << std::endl;
1341  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow, Xpetra::StaticProfile));
1342 
1343  } else {
1344  // general case
1345  double nnzPerRowInA = Teuchos::as<double>(A.getNodeNumEntries()) / A.getNodeNumRows();
1346  double nnzPerRowInB = Teuchos::as<double>(B.getNodeNumEntries()) / B.getNodeNumRows();
1347  LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
1348 
1349  LO maxPossible = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
1350  //Use static profiling (more efficient) if the estimate is at least as big as the max
1351  //possible nnz's in any single row of the result.
1352  Xpetra::ProfileType pft = (maxPossible) > nnzToAllocate ? Xpetra::DynamicProfile : Xpetra::StaticProfile;
1353 
1354  fos << "nnzPerRowInA = " << nnzPerRowInA << ", nnzPerRowInB = " << nnzPerRowInB << std::endl;
1355  fos << "MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
1356  << ", max possible nnz per row in sum = " << maxPossible
1357  << ", using " << (pft == Xpetra::DynamicProfile ? "dynamic" : "static" ) << " profiling"
1358  << std::endl;
1359  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), nnzToAllocate, pft));
1360  }
1361  if (transposeB)
1362  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1363  }
1364 
1365  if (C->getRowMap()->lib() == Xpetra::UseEpetra) {
1366  #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1367  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1368  const Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(B);
1369  RCP<Epetra_CrsMatrix> epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
1370  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1371 
1372  //FIXME is there a bug if beta=0?
1373  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1374 
1375  if (rv != 0)
1376  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1377  #else
1378  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1379  #endif
1380  } else if (C->getRowMap()->lib() == Xpetra::UseTpetra) {
1381  #ifdef HAVE_XPETRA_TPETRA
1382  # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1383  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1384  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1385  # else
1386  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1387  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
1388  RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
1389 
1390  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
1391  # endif
1392  #else
1393  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
1394  #endif
1395  }
1396 
1398  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
1399  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
1401  }
1402  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
1403  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1404  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1405  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1406 
1407  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1408  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1409 
1410  size_t i = 0;
1411  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1412  RCP<Matrix> Cij = Teuchos::null;
1413  if(rcpA != Teuchos::null &&
1414  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1415  // recursive call
1416  TwoMatrixAdd(*rcpA, transposeA, alpha,
1417  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1418  Cij, fos, AHasFixedNnzPerRow);
1419  } else if (rcpA == Teuchos::null &&
1420  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1421  Cij = rcpBopB->getMatrix(i,j);
1422  } else if (rcpA != Teuchos::null &&
1423  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1424  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
1425  } else {
1426  Cij = Teuchos::null;
1427  }
1428 
1429  if (!Cij.is_null()) {
1430  if (Cij->isFillComplete())
1431  Cij->resumeFill();
1432  Cij->fillComplete();
1433  bC->setMatrix(i, j, Cij);
1434  } else {
1435  bC->setMatrix(i, j, Teuchos::null);
1436  }
1437  } // loop over columns j
1438  }
1439  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
1440  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
1441  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1442  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
1443 
1444  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1445  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1446 
1447  size_t j = 0;
1448  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1449  RCP<Matrix> Cij = Teuchos::null;
1450  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1451  rcpB!=Teuchos::null) {
1452  // recursive call
1453  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1454  *rcpB, transposeB, beta,
1455  Cij, fos, AHasFixedNnzPerRow);
1456  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1457  rcpB!=Teuchos::null) {
1458  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
1459  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1460  rcpB==Teuchos::null) {
1461  Cij = rcpBopA->getMatrix(i,j);
1462  } else {
1463  Cij = Teuchos::null;
1464  }
1465 
1466  if (!Cij.is_null()) {
1467  if (Cij->isFillComplete())
1468  Cij->resumeFill();
1469  Cij->fillComplete();
1470  bC->setMatrix(i, j, Cij);
1471  } else {
1472  bC->setMatrix(i, j, Teuchos::null);
1473  }
1474  } // loop over rows i
1475  }
1476  else {
1477  // This is the version for blocked matrices
1478 
1479  // check the compatibility of the blocked operators
1480  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1481  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1482  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)");
1483  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)");
1484  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)");
1485  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)");
1486  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)");
1487  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)");
1488 
1489  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1490  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1491 
1492  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1493  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1494 
1495  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1496  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1497 
1498  RCP<Matrix> Cij = Teuchos::null;
1499  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1500  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1501  // recursive call
1502 
1503  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1504  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1505  Cij, fos, AHasFixedNnzPerRow);
1506  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1507  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1508  Cij = rcpBopB->getMatrix(i,j);
1509  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1510  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1511  Cij = rcpBopA->getMatrix(i,j);
1512  } else {
1513  Cij = Teuchos::null;
1514  }
1515 
1516  if (!Cij.is_null()) {
1517  if (Cij->isFillComplete())
1518  Cij->resumeFill();
1519  //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
1520  Cij->fillComplete();
1521  bC->setMatrix(i, j, Cij);
1522  } else {
1523  bC->setMatrix(i, j, Teuchos::null);
1524  }
1525  } // loop over columns j
1526  } // loop over rows i
1527 
1528  } // end blocked recursive algorithm
1529  } //MatrixMatrix::TwoMatrixAdd()
1530  }; // end specialization on SC=double, GO=int and NO=EpetraNode
1531 
1532  // specialization MatrixMatrix for SC=double, GO=long long and NO=EptraNode
1533  template<>
1534  class MatrixMatrix<double,int,long long,EpetraNode> {
1535  typedef double Scalar;
1536  typedef int LocalOrdinal;
1537  typedef long long GlobalOrdinal;
1538  typedef EpetraNode Node;
1539 #include "Xpetra_UseShortNames.hpp"
1540 
1541  public:
1542 
1567  static void Multiply(const Matrix& A, bool transposeA,
1568  const Matrix& B, bool transposeB,
1569  Matrix& C,
1570  bool call_FillComplete_on_result = true,
1571  bool doOptimizeStorage = true,
1572  const std::string & label = std::string(),
1573  const RCP<ParameterList>& params = null) {
1574  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
1575  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
1576  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
1577  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
1578 
1579 
1580  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
1581  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
1582 
1583  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1584 
1585  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
1586 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1587  Epetra_CrsMatrix & epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(A);
1588  Epetra_CrsMatrix & epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1589  Epetra_CrsMatrix & epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
1590 
1591  int i = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, epC, haveMultiplyDoFillComplete);
1592  if (haveMultiplyDoFillComplete) {
1593  // Due to Epetra wrapper intricacies, we need to explicitly call
1594  // fillComplete on Xpetra matrix here. Specifically, EpetraCrsMatrix
1595  // only keeps an internal variable to check whether we are in resumed
1596  // state or not, but never touches the underlying Epetra object. As
1597  // such, we need to explicitly update the state of Xpetra matrix to
1598  // that of Epetra one afterwords
1599  C.fillComplete();
1600  }
1601 
1602  if (i != 0) {
1603  std::ostringstream buf;
1604  buf << i;
1605  std::string msg = "EpetraExt::MatrixMatrix::Multiply return value of " + buf.str();
1606  throw(Exceptions::RuntimeError(msg));
1607  }
1608 
1609 #else
1610  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
1611 #endif
1612  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
1613 #ifdef HAVE_XPETRA_TPETRA
1614 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1615  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1616  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long, EpetraNode> ETI enabled."));
1617 # else
1618  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1619  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
1620  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
1621 
1622  //18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
1623  //Previously, Tpetra's matrix matrix multiply did not support fillComplete.
1624  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1625 # endif
1626 #else
1627  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
1628 #endif
1629  }
1630 
1631  if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1632  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
1633  fillParams->set("Optimize Storage",doOptimizeStorage);
1634  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
1635  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
1636  fillParams);
1637  }
1638 
1639  // transfer striding information
1640  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1641  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1642  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
1643  } // end Multiply
1644 
1667  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1668  const Matrix& B, bool transposeB,
1669  RCP<Matrix> C_in,
1670  Teuchos::FancyOStream& fos,
1671  bool doFillComplete = true,
1672  bool doOptimizeStorage = true,
1673  const std::string & label = std::string(),
1674  const RCP<ParameterList>& params = null) {
1675 
1676  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1677  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1678 
1679  // Default case: Xpetra Multiply
1680  RCP<Matrix> C = C_in;
1681 
1682  if (C == Teuchos::null) {
1683  double nnzPerRow = Teuchos::as<double>(0);
1684 
1685 #if 0
1686  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
1687  // For now, follow what ML and Epetra do.
1688  GO numRowsA = A.getGlobalNumRows();
1689  GO numRowsB = B.getGlobalNumRows();
1690  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
1691  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
1692  nnzPerRow *= nnzPerRow;
1693  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1694  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1695  if (totalNnz < minNnz)
1696  totalNnz = minNnz;
1697  nnzPerRow = totalNnz / A.getGlobalNumRows();
1698 
1699  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1700  }
1701 #endif
1702  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1703  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1704 
1705  } else {
1706  C->resumeFill(); // why this is not done inside of Tpetra MxM?
1707  fos << "Reuse C pattern" << std::endl;
1708  }
1709 
1710  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1711 
1712  return C;
1713  }
1714 
1725  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1726  const Matrix& B, bool transposeB,
1727  Teuchos::FancyOStream &fos,
1728  bool callFillCompleteOnResult = true,
1729  bool doOptimizeStorage = true,
1730  const std::string & label = std::string(),
1731  const RCP<ParameterList>& params = null) {
1732  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1733  }
1734 
1735 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1736  // Michael Gee's MLMultiply
1737  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
1738  const Epetra_CrsMatrix& epB,
1739  Teuchos::FancyOStream& fos) {
1740  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1741  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1742  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1743  }
1744 #endif //ifdef HAVE_XPETRA_EPETRAEXT
1745 
1756  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
1757  const BlockedCrsMatrix& B, bool transposeB,
1758  Teuchos::FancyOStream& fos,
1759  bool doFillComplete = true,
1760  bool doOptimizeStorage = true) {
1761  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1762  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1763 
1764  // Preconditions
1765  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1766  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1767 
1768  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1769  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
1770 
1771  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1772 
1773  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1774  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1775  RCP<Matrix> Cij;
1776 
1777  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1778  RCP<Matrix> crmat1 = A.getMatrix(i,l);
1779  RCP<Matrix> crmat2 = B.getMatrix(l,j);
1780 
1781  if (crmat1.is_null() || crmat2.is_null())
1782  continue;
1783 
1784  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1785  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1786  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
1787  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1788 
1789  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
1790  if (!crop1.is_null())
1791  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
1792  if (!crop2.is_null())
1793  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
1794 
1795  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
1796  "crmat1 does not have global constants");
1797  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
1798  "crmat2 does not have global constants");
1799 
1800  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1801  continue;
1802 
1803  // temporary matrix containing result of local block multiplication
1804  RCP<Matrix> temp = Teuchos::null;
1805 
1806  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1807  temp = Multiply (*crop1, false, *crop2, false, fos);
1808  else {
1809  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1810  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1811  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1812  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1813  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)");
1814  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)");
1815 
1816  // recursive multiplication call
1817  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1818 
1819  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1820  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)");
1821  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)");
1822  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)");
1823  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)");
1824  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)");
1825  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)");
1826  }
1827 
1828  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1829 
1830  RCP<Matrix> addRes = null;
1831  if (Cij.is_null ())
1832  Cij = temp;
1833  else {
1834  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1835  Cij = addRes;
1836  }
1837  }
1838 
1839  if (!Cij.is_null()) {
1840  if (Cij->isFillComplete())
1841  Cij->resumeFill();
1842  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1843  C->setMatrix(i, j, Cij);
1844  } else {
1845  C->setMatrix(i, j, Teuchos::null);
1846  }
1847  }
1848  }
1849 
1850  if (doFillComplete)
1851  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1852 
1853  return C;
1854  } // TwoMatrixMultiplyBlock
1855 
1868  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1869  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
1870  "TwoMatrixAdd: matrix row maps are not the same.");
1871 
1872  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1873 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1874  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1875  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1876 
1877  //FIXME is there a bug if beta=0?
1878  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1879 
1880  if (rv != 0)
1881  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1882  std::ostringstream buf;
1883 #else
1884  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1885 #endif
1886  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1887 #ifdef HAVE_XPETRA_TPETRA
1888 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1889  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1890  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
1891 # else
1892  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1893  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1894 
1895  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1896 # endif
1897 #else
1898  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1899 #endif
1900  }
1901  } //MatrixMatrix::TwoMatrixAdd()
1902 
1903 
1916  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1917  const Matrix& B, bool transposeB, const SC& beta,
1918  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1919  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1920  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1921  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1922  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1923 
1924  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1925  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
1926  "TwoMatrixAdd: matrix row maps are not the same.");
1927 
1928  if (C == Teuchos::null) {
1929  size_t maxNzInA = 0;
1930  size_t maxNzInB = 0;
1931  size_t numLocalRows = 0;
1932  if (A.isFillComplete() && B.isFillComplete()) {
1933  maxNzInA = A.getNodeMaxNumRowEntries();
1934  maxNzInB = B.getNodeMaxNumRowEntries();
1935  numLocalRows = A.getNodeNumRows();
1936  }
1937 
1938  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1939  // first check if either A or B has at most 1 nonzero per row
1940  // the case of both having at most 1 nz per row is handled by the ``else''
1941  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1942 
1943  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1944  for (size_t i = 0; i < numLocalRows; ++i)
1945  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1946 
1947  } else {
1948  for (size_t i = 0; i < numLocalRows; ++i)
1949  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1950  }
1951 
1952  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1953  << ", using static profiling" << std::endl;
1954  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow, Xpetra::StaticProfile));
1955 
1956  } else {
1957  // general case
1958  double nnzPerRowInA = Teuchos::as<double>(A.getNodeNumEntries()) / A.getNodeNumRows();
1959  double nnzPerRowInB = Teuchos::as<double>(B.getNodeNumEntries()) / B.getNodeNumRows();
1960  LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
1961 
1962  LO maxPossible = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
1963  //Use static profiling (more efficient) if the estimate is at least as big as the max
1964  //possible nnz's in any single row of the result.
1965  Xpetra::ProfileType pft = (maxPossible) > nnzToAllocate ? Xpetra::DynamicProfile : Xpetra::StaticProfile;
1966 
1967  fos << "nnzPerRowInA = " << nnzPerRowInA << ", nnzPerRowInB = " << nnzPerRowInB << std::endl;
1968  fos << "MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
1969  << ", max possible nnz per row in sum = " << maxPossible
1970  << ", using " << (pft == Xpetra::DynamicProfile ? "dynamic" : "static" ) << " profiling"
1971  << std::endl;
1972  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), nnzToAllocate, pft));
1973  }
1974  if (transposeB)
1975  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1976  }
1977 
1978  if (C->getRowMap()->lib() == Xpetra::UseEpetra) {
1979  #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1980  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1981  const Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(B);
1982  RCP<Epetra_CrsMatrix> epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
1983  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1984 
1985  //FIXME is there a bug if beta=0?
1986  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1987 
1988  if (rv != 0)
1989  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1990  #else
1991  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1992  #endif
1993  } else if (C->getRowMap()->lib() == Xpetra::UseTpetra) {
1994  #ifdef HAVE_XPETRA_TPETRA
1995  # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1996  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1997  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
1998  # else
1999  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
2000  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
2001  RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
2002 
2003  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
2004  # endif
2005  #else
2006  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
2007  #endif
2008  }
2009 
2011  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
2012  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
2014  }
2015  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
2016  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
2017  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
2018  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2019 
2020  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2021  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2022 
2023  size_t i = 0;
2024  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2025  RCP<Matrix> Cij = Teuchos::null;
2026  if(rcpA != Teuchos::null &&
2027  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2028  // recursive call
2029  TwoMatrixAdd(*rcpA, transposeA, alpha,
2030  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2031  Cij, fos, AHasFixedNnzPerRow);
2032  } else if (rcpA == Teuchos::null &&
2033  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2034  Cij = rcpBopB->getMatrix(i,j);
2035  } else if (rcpA != Teuchos::null &&
2036  rcpBopB->getMatrix(i,j)==Teuchos::null) {
2037  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
2038  } else {
2039  Cij = Teuchos::null;
2040  }
2041 
2042  if (!Cij.is_null()) {
2043  if (Cij->isFillComplete())
2044  Cij->resumeFill();
2045  Cij->fillComplete();
2046  bC->setMatrix(i, j, Cij);
2047  } else {
2048  bC->setMatrix(i, j, Teuchos::null);
2049  }
2050  } // loop over columns j
2051  }
2052  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
2053  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
2054  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2055  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
2056 
2057  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2058  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2059 
2060  size_t j = 0;
2061  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2062  RCP<Matrix> Cij = Teuchos::null;
2063  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2064  rcpB!=Teuchos::null) {
2065  // recursive call
2066  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2067  *rcpB, transposeB, beta,
2068  Cij, fos, AHasFixedNnzPerRow);
2069  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2070  rcpB!=Teuchos::null) {
2071  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
2072  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2073  rcpB==Teuchos::null) {
2074  Cij = rcpBopA->getMatrix(i,j);
2075  } else {
2076  Cij = Teuchos::null;
2077  }
2078 
2079  if (!Cij.is_null()) {
2080  if (Cij->isFillComplete())
2081  Cij->resumeFill();
2082  Cij->fillComplete();
2083  bC->setMatrix(i, j, Cij);
2084  } else {
2085  bC->setMatrix(i, j, Teuchos::null);
2086  }
2087  } // loop over rows i
2088  }
2089  else {
2090  // This is the version for blocked matrices
2091 
2092  // check the compatibility of the blocked operators
2093  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2094  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2095  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)");
2096  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)");
2097  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)");
2098  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)");
2099  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)");
2100  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)");
2101 
2102  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2103  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2104 
2105  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2106  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2107 
2108  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2109  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2110 
2111  RCP<Matrix> Cij = Teuchos::null;
2112  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2113  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2114  // recursive call
2115 
2116  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2117  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2118  Cij, fos, AHasFixedNnzPerRow);
2119  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2120  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2121  Cij = rcpBopB->getMatrix(i,j);
2122  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2123  rcpBopB->getMatrix(i,j)==Teuchos::null) {
2124  Cij = rcpBopA->getMatrix(i,j);
2125  } else {
2126  Cij = Teuchos::null;
2127  }
2128 
2129  if (!Cij.is_null()) {
2130  if (Cij->isFillComplete())
2131  Cij->resumeFill();
2132  //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
2133  Cij->fillComplete();
2134  bC->setMatrix(i, j, Cij);
2135  } else {
2136  bC->setMatrix(i, j, Teuchos::null);
2137  }
2138  } // loop over columns j
2139  } // loop over rows i
2140 
2141  } // end blocked recursive algorithm
2142  } //MatrixMatrix::TwoMatrixAdd()
2143  }; // end specialization on GO=long long and NO=EpetraNode
2144 
2145 #endif // HAVE_XPETRA_EPETRA
2146 
2147 } // end namespace Xpetra
2148 
2149 #define XPETRA_MATRIXMATRIX_SHORT
2150 
2151 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_ */
Epetra_CrsMatrix::RowMap
const Epetra_Map & RowMap() const
Xpetra::Matrix::getNodeMaxNumRowEntries
virtual size_t getNodeMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
Xpetra::Matrix::getNodeNumRows
virtual size_t getNodeNumRows() const =0
Returns the number of matrix rows owned on the calling node.
Xpetra::toString
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.
Definition: Xpetra_Utils.cpp:51
Epetra_Map
Xpetra::MatrixMatrix::TwoMatrixAdd
static void TwoMatrixAdd(const Matrix &A, bool transposeA, const SC &alpha, const Matrix &B, bool transposeB, const SC &beta, RCP< Matrix > &C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow=false)
Helper function to calculate C = alpha*A + beta*B.
Definition: Xpetra_MatrixMatrix.hpp:550
Xpetra::Matrix::getNodeNumEntries
virtual size_t getNodeNumEntries() const =0
Returns the local number of entries in this matrix.
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::TwoMatrixAdd
static void TwoMatrixAdd(const Matrix &A, bool transposeA, const SC &alpha, const Matrix &B, bool transposeB, const SC &beta, RCP< Matrix > &C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow=false)
Helper function to calculate C = alpha*A + beta*B.
Definition: Xpetra_MatrixMatrix.hpp:1300
Xpetra::Matrix::getGlobalNumRows
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
Copy
Copy
Xpetra::Matrix::getRowMap
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Xpetra_Matrix.hpp:320
Xpetra::BlockedCrsMatrix::isFillComplete
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Definition: Xpetra_BlockedCrsMatrix.hpp:701
Xpetra::MatrixMatrix< double, int, long long, EpetraNode >::Multiply
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, RCP< Matrix > C_in, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Helper function to do matrix-matrix multiply.
Definition: Xpetra_MatrixMatrix.hpp:1667
GO
GlobalOrdinal GO
Definition: Xpetra_UseShortNamesOrdinal.hpp:139
Xpetra::Helpers::Op2NonConstEpetraCrs
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:142
Xpetra_UseShortNames.hpp
Xpetra::Helpers
Xpetra utility class containing transformation routines between Xpetra::Matrix and Epetra/Tpetra obje...
Definition: Xpetra_MatrixMatrix.hpp:92
Xpetra
Xpetra namespace
Definition: Xpetra_BlockedCrsMatrix.hpp:86
Xpetra_EpetraCrsMatrix_fwd.hpp
Xpetra::Matrix::resumeFill
virtual void resumeFill(const RCP< ParameterList > &params=null)=0
Epetra_BlockMap::IndexBase
int IndexBase() const
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:101
LO
LocalOrdinal LO
Definition: Xpetra_UseShortNamesOrdinal.hpp:138
Xpetra::MatrixMatrix< double, int, long long, EpetraNode >::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.hpp:1756
Epetra_CrsMatrix::Comm
const Epetra_Comm & Comm() const
Epetra_MpiComm
Xpetra_Matrix.hpp
Epetra_CrsMatrix::RangeMap
const Epetra_Map & RangeMap() const
Xpetra_StridedMapFactory.hpp
Xpetra::Helpers::Op2EpetraCrs
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:126
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1073
Xpetra::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: Xpetra_Exceptions.hpp:86
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::MatrixMatrix< double, int, long long, EpetraNode >::TwoMatrixAdd
static void TwoMatrixAdd(const Matrix &A, bool transposeA, const SC &alpha, const Matrix &B, bool transposeB, const SC &beta, RCP< Matrix > &C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow=false)
Helper function to calculate C = alpha*A + beta*B.
Definition: Xpetra_MatrixMatrix.hpp:1916
Xpetra::MatrixMatrix::Multiply
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool callFillCompleteOnResult=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Helper function to do matrix-matrix multiply.
Definition: Xpetra_MatrixMatrix.hpp:381
Xpetra_TpetraMultiVector.hpp
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Xpetra_StridedMap.hpp
Xpetra::Matrix::CreateView
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
Definition: Xpetra_Matrix.hpp:128
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::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.hpp:800
Epetra_MpiComm::GetMpiComm
MPI_Comm GetMpiComm() const
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition: Xpetra_CrsMatrixWrap.hpp:621
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::Multiply
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, RCP< Matrix > C_in, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Helper function to do matrix-matrix multiply.
Definition: Xpetra_MatrixMatrix.hpp:899
Xpetra::Helpers::Op2NonConstTpetraCrs
static Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:198
Xpetra::MatrixMatrix::Multiply
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, RCP< Matrix > C_in, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Helper function to do matrix-matrix multiply.
Definition: Xpetra_MatrixMatrix.hpp:325
Xpetra::Matrix::getNumEntriesInLocalRow
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
Returns the current number of entries on this node in the specified local row.
Epetra_CrsMatrix::Filled
bool Filled() const
Epetra_CrsMatrix::ColMap
const Epetra_Map & ColMap() const
Xpetra::MatrixMatrix< double, int, long long, EpetraNode >::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.hpp:1567
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.hpp:407
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.hpp:519
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
Xpetra::Helpers::Op2EpetraCrs
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:98
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.hpp:250
Xpetra::Helpers::Op2NonConstEpetraCrs
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:112
Epetra_CrsMatrix
Xpetra_ConfigDefs.hpp
Xpetra::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: Xpetra_Exceptions.hpp:108
SC
Scalar SC
Definition: Xpetra_UseShortNamesScalar.hpp:180
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::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.hpp:1252
Xpetra::BlockedCrsMatrix::getRangeMap
RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1064
Xpetra_MapExtractor.hpp
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_fwd.hpp:51
Xpetra::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition: Xpetra_BlockedCrsMatrix.hpp:1308
Xpetra::Operator::getRangeMap
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
Xpetra::Matrix::IsView
bool IsView(const viewLabel_t viewLabel) const
Definition: Xpetra_Matrix.hpp:207
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
Definition: Xpetra_MatrixMatrix.hpp:994
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1076
Xpetra_Map.hpp
Xpetra_CrsMatrixWrap.hpp
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor specifying the number of non-zeros for all rows.
Definition: Xpetra_MatrixFactory.hpp:235
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::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.hpp:1139
Xpetra::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1272
Xpetra::Helpers::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:172
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1275
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::Multiply
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool callFillCompleteOnResult=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Helper function to do matrix-matrix multiply.
Definition: Xpetra_MatrixMatrix.hpp:982
Xpetra::BlockedCrsMatrix::getDomainMap
RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1049
Xpetra_TpetraVector.hpp
Xpetra::MatrixMatrix::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
Definition: Xpetra_MatrixMatrix.hpp:389
Xpetra::MatrixMatrix< double, int, long long, EpetraNode >::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
Definition: Xpetra_MatrixMatrix.hpp:1737
Epetra_BlockMap::GID
int GID(int LID) const
Xpetra::MatrixMatrix< double, int, long long, EpetraNode >::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.hpp:1868
Xpetra::Helpers::Op2TpetraCrs
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:183
Xpetra::Operator::getDomainMap
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_fwd.hpp:51
Xpetra::Helpers::Op2TpetraCrs
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:160
Epetra_CrsMatrix::DomainMap
const Epetra_Map & DomainMap() const
Xpetra::MatrixMatrix< double, int, long long, EpetraNode >::Multiply
static RCP< Matrix > Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool callFillCompleteOnResult=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Helper function to do matrix-matrix multiply.
Definition: Xpetra_MatrixMatrix.hpp:1725
Xpetra_TpetraCrsMatrix.hpp
Generated by   doxygen 1.8.14