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  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  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
435  continue;
436 
437  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
438  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
439  TEUCHOS_TEST_FOR_EXCEPTION((crop1==Teuchos::null) != (crop2==Teuchos::null), Xpetra::Exceptions::RuntimeError, "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
440 
441  // temporary matrix containing result of local block multiplication
442  RCP<Matrix> temp = Teuchos::null;
443 
444  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
445  temp = Multiply (*crop1, false, *crop2, false, fos);
446  else {
447  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
448  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
449  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
450  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
451  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)");
452  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)");
453 
454  // recursive multiplication call
455  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
456 
457  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
458  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)");
459  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)");
460  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)");
461  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)");
462  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)");
463  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)");
464  }
465 
466  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
467 
468  RCP<Matrix> addRes = null;
469  if (Cij.is_null ())
470  Cij = temp;
471  else {
472  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
473  Cij = addRes;
474  }
475  }
476 
477  if (!Cij.is_null()) {
478  if (Cij->isFillComplete())
479  Cij->resumeFill();
480  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
481  C->setMatrix(i, j, Cij);
482  } else {
483  C->setMatrix(i, j, Teuchos::null);
484  }
485  }
486  }
487 
488  if (doFillComplete)
489  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
490 
491  return C;
492  } // TwoMatrixMultiplyBlock
493 
506  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
507  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
508  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
509 
510  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
511  throw Exceptions::RuntimeError("TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
512  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
513 #ifdef HAVE_XPETRA_TPETRA
514  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
515  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
516 
517  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
518 #else
519  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
520 #endif
521  }
522  } //MatrixMatrix::TwoMatrixAdd()
523 
524 
537  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
538  const Matrix& B, bool transposeB, const SC& beta,
539  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
540 
541  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
542  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
543  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
544  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
545 
546  // We have to distinguish 4 cases:
547  // both matrices are CrsMatrixWrap based, only one of them is CrsMatrixWrap based, or none.
548 
549  // Both matrices are CrsMatrixWrap
550  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
551 
552  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
553  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
554 
555 
556  if (C == Teuchos::null) {
557  size_t maxNzInA = 0;
558  size_t maxNzInB = 0;
559  size_t numLocalRows = 0;
560  if (A.isFillComplete() && B.isFillComplete()) {
561  maxNzInA = A.getNodeMaxNumRowEntries();
562  maxNzInB = B.getNodeMaxNumRowEntries();
563  numLocalRows = A.getNodeNumRows();
564  }
565  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
566  // first check if either A or B has at most 1 nonzero per row
567  // the case of both having at most 1 nz per row is handled by the ``else''
568  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
569 
570  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
571  for (size_t i = 0; i < numLocalRows; ++i)
572  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
573 
574  } else {
575  for (size_t i = 0; i < numLocalRows; ++i)
576  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
577  }
578 
579  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
580  << ", using static profiling" << std::endl;
581  C = rcp(new CrsMatrixWrap(A.getRowMap(), exactNnzPerRow, Xpetra::StaticProfile));
582  } else {
583  // general case
584  double nnzPerRowInA = Teuchos::as<double>(A.getNodeNumEntries()) / A.getNodeNumRows();
585  double nnzPerRowInB = Teuchos::as<double>(B.getNodeNumEntries()) / B.getNodeNumRows();
586  LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
587 
588  LO maxPossible = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
589  //Use static profiling (more efficient) if the estimate is at least as big as the max
590  //possible nnz's in any single row of the result.
591  Xpetra::ProfileType pft = (maxPossible) > nnzToAllocate ? Xpetra::DynamicProfile : Xpetra::StaticProfile;
592 
593  fos << "nnzPerRowInA = " << nnzPerRowInA << ", nnzPerRowInB = " << nnzPerRowInB << std::endl;
594  fos << "MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
595  << ", max possible nnz per row in sum = " << maxPossible
596  << ", using " << (pft == Xpetra::DynamicProfile ? "dynamic" : "static" ) << " profiling"
597  << std::endl;
598  C = rcp(new CrsMatrixWrap(A.getRowMap(), nnzToAllocate, pft));
599  }
600  if (transposeB)
601  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
602  }
603 
604  if (C->getRowMap()->lib() == Xpetra::UseEpetra) {
605  throw Exceptions::RuntimeError("MatrixMatrix::Add for Epetra only available with Scalar = double, LO = GO = int.");
606  } else if (C->getRowMap()->lib() == Xpetra::UseTpetra) {
607  #ifdef HAVE_XPETRA_TPETRA
608  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
609  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
610  RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
611 
612  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
613  #else
614  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
615  #endif
616  }
618  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
619  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
621  }
622  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
623  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
624  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
625  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
626 
627  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
628  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
629 
630  size_t i = 0;
631  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
632  RCP<Matrix> Cij = Teuchos::null;
633  if(rcpA != Teuchos::null &&
634  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
635  // recursive call
636  TwoMatrixAdd(*rcpA, transposeA, alpha,
637  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
638  Cij, fos, AHasFixedNnzPerRow);
639  } else if (rcpA == Teuchos::null &&
640  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
641  Cij = rcpBopB->getMatrix(i,j);
642  } else if (rcpA != Teuchos::null &&
643  rcpBopB->getMatrix(i,j)==Teuchos::null) {
644  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
645  } else {
646  Cij = Teuchos::null;
647  }
648 
649  if (!Cij.is_null()) {
650  if (Cij->isFillComplete())
651  Cij->resumeFill();
652  Cij->fillComplete();
653  bC->setMatrix(i, j, Cij);
654  } else {
655  bC->setMatrix(i, j, Teuchos::null);
656  }
657  } // loop over columns j
658  }
659  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
660  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
661  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
662  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
663 
664  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
665  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
666  size_t j = 0;
667  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
668  RCP<Matrix> Cij = Teuchos::null;
669  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
670  rcpB!=Teuchos::null) {
671  // recursive call
672  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
673  *rcpB, transposeB, beta,
674  Cij, fos, AHasFixedNnzPerRow);
675  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
676  rcpB!=Teuchos::null) {
677  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
678  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
679  rcpB==Teuchos::null) {
680  Cij = rcpBopA->getMatrix(i,j);
681  } else {
682  Cij = Teuchos::null;
683  }
684 
685  if (!Cij.is_null()) {
686  if (Cij->isFillComplete())
687  Cij->resumeFill();
688  Cij->fillComplete();
689  bC->setMatrix(i, j, Cij);
690  } else {
691  bC->setMatrix(i, j, Teuchos::null);
692  }
693  } // loop over rows i
694  }
695  else {
696  // This is the version for blocked matrices
697 
698  // check the compatibility of the blocked operators
699  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
700  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
701  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)");
702  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)");
703  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)");
704  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)");
705  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)");
706  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)");
707 
708  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
709  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
710 
711  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
712  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
713  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
714  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
715 
716  RCP<Matrix> Cij = Teuchos::null;
717  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
718  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
719  // recursive call
720  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
721  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
722  Cij, fos, AHasFixedNnzPerRow);
723  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
724  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
725  Cij = rcpBopB->getMatrix(i,j);
726  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
727  rcpBopB->getMatrix(i,j)==Teuchos::null) {
728  Cij = rcpBopA->getMatrix(i,j);
729  } else {
730  Cij = Teuchos::null;
731  }
732 
733  if (!Cij.is_null()) {
734  if (Cij->isFillComplete())
735  Cij->resumeFill();
736  Cij->fillComplete();
737  bC->setMatrix(i, j, Cij);
738  } else {
739  bC->setMatrix(i, j, Teuchos::null);
740  }
741  } // loop over columns j
742  } // loop over rows i
743 
744  } // end blocked recursive algorithm
745  } //MatrixMatrix::TwoMatrixAdd()
746 
747 
748  }; // class MatrixMatrix
749 
750 
751 #ifdef HAVE_XPETRA_EPETRA
752  // specialization MatrixMatrix for SC=double, LO=GO=int
753  template<>
754  class MatrixMatrix<double,int,int,EpetraNode> {
755  typedef double Scalar;
756  typedef int LocalOrdinal;
757  typedef int GlobalOrdinal;
758  typedef EpetraNode Node;
759 #include "Xpetra_UseShortNames.hpp"
760 
761  public:
762 
787  static void Multiply(const Matrix& A, bool transposeA,
788  const Matrix& B, bool transposeB,
789  Matrix& C,
790  bool call_FillComplete_on_result = true,
791  bool doOptimizeStorage = true,
792  const std::string & label = std::string(),
793  const RCP<ParameterList>& params = null) {
794  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
795  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
796  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
797  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
798 
799  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
800  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
801 
802  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
803 
804  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
805 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
806  Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(A);
807  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
808  Epetra_CrsMatrix& epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
809 
810  int i = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, epC, haveMultiplyDoFillComplete);
811  if (haveMultiplyDoFillComplete) {
812  // Due to Epetra wrapper intricacies, we need to explicitly call
813  // fillComplete on Xpetra matrix here. Specifically, EpetraCrsMatrix
814  // only keeps an internal variable to check whether we are in resumed
815  // state or not, but never touches the underlying Epetra object. As
816  // such, we need to explicitly update the state of Xpetra matrix to
817  // that of Epetra one afterwords
818  C.fillComplete();
819  }
820 
821  if (i != 0) {
822  std::ostringstream buf;
823  buf << i;
824  std::string msg = "EpetraExt::MatrixMatrix::Multiply return value of " + buf.str();
825  throw(Exceptions::RuntimeError(msg));
826  }
827 
828 #else
829  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
830 #endif
831  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
832 #ifdef HAVE_XPETRA_TPETRA
833 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
834  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
835  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
836 # else
837  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
838  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
839  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
840 
841  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
842  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
843  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
844 # endif
845 #else
846  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
847 #endif
848  }
849 
850  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
851  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
852  fillParams->set("Optimize Storage",doOptimizeStorage);
853  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
854  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
855  fillParams);
856  }
857 
858  // transfer striding information
859  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
860  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
861  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
862  } // end Multiply
863 
886  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
887  const Matrix& B, bool transposeB,
888  RCP<Matrix> C_in,
889  Teuchos::FancyOStream& fos,
890  bool doFillComplete = true,
891  bool doOptimizeStorage = true,
892  const std::string & label = std::string(),
893  const RCP<ParameterList>& params = null) {
894 
895  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
896  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
897 
898  // Optimization using ML Multiply when available and requested
899  // This feature is currently not supported. We would have to introduce the HAVE_XPETRA_ML_MMM flag
900 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) && defined(HAVE_XPETRA_ML_MMM)
901  if (B.getDomainMap()->lib() == Xpetra::UseEpetra && !transposeA && !transposeB) {
902  RCP<const Epetra_CrsMatrix> epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(rcpFromRef(A));
903  RCP<const Epetra_CrsMatrix> epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(rcpFromRef(B));
904  RCP<Epetra_CrsMatrix> epC = MLTwoMatrixMultiply(*epA, *epB, fos);
905 
906  RCP<Matrix> C = Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO> (epC);
907  if (doFillComplete) {
908  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
909  fillParams->set("Optimize Storage", doOptimizeStorage);
910  C->fillComplete(B.getDomainMap(), A.getRangeMap(), fillParams);
911  }
912 
913  // Fill strided maps information
914  // This is necessary since the ML matrix matrix multiplication routine has no handling for this
915  // TODO: move this call to MLMultiply...
916  C->CreateView("stridedMaps", rcpFromRef(A), transposeA, rcpFromRef(B), transposeB);
917 
918  return C;
919  }
920 #endif // EPETRA + EPETRAEXT + ML
921 
922  // Default case: Xpetra Multiply
923  RCP<Matrix> C = C_in;
924 
925  if (C == Teuchos::null) {
926  double nnzPerRow = Teuchos::as<double>(0);
927 
928 #if 0
929  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
930  // For now, follow what ML and Epetra do.
931  GO numRowsA = A.getGlobalNumRows();
932  GO numRowsB = B.getGlobalNumRows();
933  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
934  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
935  nnzPerRow *= nnzPerRow;
936  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
937  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
938  if (totalNnz < minNnz)
939  totalNnz = minNnz;
940  nnzPerRow = totalNnz / A.getGlobalNumRows();
941 
942  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
943  }
944 #endif
945 
946  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
947  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
948 
949  } else {
950  C->resumeFill(); // why this is not done inside of Tpetra MxM?
951  fos << "Reuse C pattern" << std::endl;
952  }
953 
954  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
955 
956  return C;
957  }
958 
969  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
970  const Matrix& B, bool transposeB,
971  Teuchos::FancyOStream &fos,
972  bool callFillCompleteOnResult = true,
973  bool doOptimizeStorage = true,
974  const std::string & label = std::string(),
975  const RCP<ParameterList>& params = null) {
976  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
977  }
978 
979 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
980  // Michael Gee's MLMultiply
981  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
982  const Epetra_CrsMatrix& epB,
983  Teuchos::FancyOStream& fos) {
984 #if defined(HAVE_XPETRA_ML_MMM) // Note: this is currently not supported
985  ML_Comm* comm;
986  ML_Comm_Create(&comm);
987  fos << "****** USING ML's MATRIX MATRIX MULTIPLY ******" << std::endl;
988 #ifdef HAVE_MPI
989  // ML_Comm uses MPI_COMM_WORLD, so try to use the same communicator as epA.
990  const Epetra_MpiComm * Mcomm = dynamic_cast<const Epetra_MpiComm*>(&(epA.Comm()));
991  if (Mcomm)
992  ML_Comm_Set_UsrComm(comm,Mcomm->GetMpiComm());
993 # endif
994  //in order to use ML, there must be no indices missing from the matrix column maps.
995  EpetraExt::CrsMatrix_SolverMap Atransform;
996  EpetraExt::CrsMatrix_SolverMap Btransform;
997  const Epetra_CrsMatrix& A = Atransform(const_cast<Epetra_CrsMatrix&>(epA));
998  const Epetra_CrsMatrix& B = Btransform(const_cast<Epetra_CrsMatrix&>(epB));
999 
1000  if (!A.Filled()) throw Exceptions::RuntimeError("A has to be FillCompleted");
1001  if (!B.Filled()) throw Exceptions::RuntimeError("B has to be FillCompleted");
1002 
1003  // create ML operators from EpetraCrsMatrix
1004  ML_Operator* ml_As = ML_Operator_Create(comm);
1005  ML_Operator* ml_Bs = ML_Operator_Create(comm);
1006  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?
1007  ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&B),ml_Bs);
1008  ML_Operator* ml_AtimesB = ML_Operator_Create(comm);
1009  {
1010  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("ML_2matmult kernel"));
1011  ML_2matmult(ml_As,ml_Bs,ml_AtimesB,ML_CSR_MATRIX); // do NOT use ML_EpetraCRS_MATRIX!!!
1012  }
1013  ML_Operator_Destroy(&ml_As);
1014  ML_Operator_Destroy(&ml_Bs);
1015 
1016  // For ml_AtimesB we have to reconstruct the column map in global indexing,
1017  // The following is going down to the salt-mines of ML ...
1018  // note: we use integers, since ML only works for Epetra...
1019  int N_local = ml_AtimesB->invec_leng;
1020  ML_CommInfoOP* getrow_comm = ml_AtimesB->getrow->pre_comm;
1021  if (!getrow_comm) throw(Exceptions::RuntimeError("ML_Operator does not have a CommInfo"));
1022  ML_Comm* comm_AB = ml_AtimesB->comm; // get comm object
1023  if (N_local != B.DomainMap().NumMyElements())
1024  throw(Exceptions::RuntimeError("Mismatch in local row dimension between ML and Epetra"));
1025  int N_rcvd = 0;
1026  int N_send = 0;
1027  int flag = 0;
1028  for (int i = 0; i < getrow_comm->N_neighbors; i++) {
1029  N_rcvd += (getrow_comm->neighbors)[i].N_rcv;
1030  N_send += (getrow_comm->neighbors)[i].N_send;
1031  if ( ((getrow_comm->neighbors)[i].N_rcv != 0) &&
1032  ((getrow_comm->neighbors)[i].rcv_list != NULL) ) flag = 1;
1033  }
1034  // For some unknown reason, ML likes to have stuff one larger than
1035  // neccessary...
1036  std::vector<double> dtemp(N_local + N_rcvd + 1); // "double" vector for comm function
1037  std::vector<int> cmap (N_local + N_rcvd + 1); // vector for gids
1038  for (int i = 0; i < N_local; ++i) {
1039  cmap[i] = B.DomainMap().GID(i);
1040  dtemp[i] = (double) cmap[i];
1041  }
1042  ML_cheap_exchange_bdry(&dtemp[0],getrow_comm,N_local,N_send,comm_AB); // do communication
1043  if (flag) { // process received data
1044  int count = N_local;
1045  const int neighbors = getrow_comm->N_neighbors;
1046  for (int i = 0; i < neighbors; i++) {
1047  const int nrcv = getrow_comm->neighbors[i].N_rcv;
1048  for (int j = 0; j < nrcv; j++)
1049  cmap[getrow_comm->neighbors[i].rcv_list[j]] = (int) dtemp[count++];
1050  }
1051  } else {
1052  for (int i = 0; i < N_local+N_rcvd; ++i)
1053  cmap[i] = (int)dtemp[i];
1054  }
1055  dtemp.clear(); // free double array
1056 
1057  // we can now determine a matching column map for the result
1058  Epetra_Map gcmap(-1, N_local+N_rcvd, &cmap[0], B.ColMap().IndexBase(), A.Comm());
1059 
1060  int allocated = 0;
1061  int rowlength;
1062  double* val = NULL;
1063  int* bindx = NULL;
1064 
1065  const int myrowlength = A.RowMap().NumMyElements();
1066  const Epetra_Map& rowmap = A.RowMap();
1067 
1068  // Determine the maximum bandwith for the result matrix.
1069  // replaces the old, very(!) memory-consuming guess:
1070  // int guessnpr = A.MaxNumEntries()*B.MaxNumEntries();
1071  int educatedguess = 0;
1072  for (int i = 0; i < myrowlength; ++i) {
1073  // get local row
1074  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1075  if (rowlength>educatedguess)
1076  educatedguess = rowlength;
1077  }
1078 
1079  // allocate our result matrix and fill it
1080  RCP<Epetra_CrsMatrix> result = rcp(new Epetra_CrsMatrix(::Copy, A.RangeMap(), gcmap, educatedguess, false));
1081 
1082  std::vector<int> gcid(educatedguess);
1083  for (int i = 0; i < myrowlength; ++i) {
1084  const int grid = rowmap.GID(i);
1085  // get local row
1086  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1087  if (!rowlength) continue;
1088  if ((int)gcid.size() < rowlength) gcid.resize(rowlength);
1089  for (int j = 0; j < rowlength; ++j) {
1090  gcid[j] = gcmap.GID(bindx[j]);
1091  if (gcid[j] < 0)
1092  throw Exceptions::RuntimeError("Error: cannot find gcid!");
1093  }
1094  int err = result->InsertGlobalValues(grid, rowlength, val, &gcid[0]);
1095  if (err != 0 && err != 1) {
1096  std::ostringstream errStr;
1097  errStr << "Epetra_CrsMatrix::InsertGlobalValues returned err=" << err;
1098  throw Exceptions::RuntimeError(errStr.str());
1099  }
1100  }
1101  // free memory
1102  if (bindx) ML_free(bindx);
1103  if (val) ML_free(val);
1104  ML_Operator_Destroy(&ml_AtimesB);
1105  ML_Comm_Destroy(&comm);
1106 
1107  return result;
1108 #else // no MUELU_ML
1109  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1110  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1111  return Teuchos::null;
1112 #endif
1113  }
1114 #endif //ifdef HAVE_XPETRA_EPETRAEXT
1115 
1126  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
1127  const BlockedCrsMatrix& B, bool transposeB,
1128  Teuchos::FancyOStream& fos,
1129  bool doFillComplete = true,
1130  bool doOptimizeStorage = true) {
1131  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1132  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1133 
1134  // Preconditions
1135  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1136  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1137 
1138  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1139  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
1140 
1141  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1142 
1143  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1144  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1145  RCP<Matrix> Cij;
1146 
1147  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1148  RCP<Matrix> crmat1 = A.getMatrix(i,l);
1149  RCP<Matrix> crmat2 = B.getMatrix(l,j);
1150 
1151  if (crmat1.is_null() || crmat2.is_null())
1152  continue;
1153  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1154  continue;
1155 
1156  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1157  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1158  TEUCHOS_TEST_FOR_EXCEPTION((crop1==Teuchos::null) != (crop2==Teuchos::null), Xpetra::Exceptions::RuntimeError, "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1159 
1160  // temporary matrix containing result of local block multiplication
1161  RCP<Matrix> temp = Teuchos::null;
1162 
1163  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1164  temp = Multiply (*crop1, false, *crop2, false, fos);
1165  else {
1166  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1167  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1168  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1169  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1170  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)");
1171  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)");
1172 
1173  // recursive multiplication call
1174  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1175 
1176  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1177  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)");
1178  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)");
1179  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)");
1180  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)");
1181  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)");
1182  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)");
1183  }
1184 
1185  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1186 
1187  RCP<Matrix> addRes = null;
1188  if (Cij.is_null ())
1189  Cij = temp;
1190  else {
1191  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1192  Cij = addRes;
1193  }
1194  }
1195 
1196  if (!Cij.is_null()) {
1197  if (Cij->isFillComplete())
1198  Cij->resumeFill();
1199  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1200  C->setMatrix(i, j, Cij);
1201  } else {
1202  C->setMatrix(i, j, Teuchos::null);
1203  }
1204  }
1205  }
1206 
1207  if (doFillComplete)
1208  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1209 
1210  return C;
1211  } // TwoMatrixMultiplyBlock
1212 
1225  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1226  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*B.getRowMap())), Exceptions::Incompatible,
1227  "TwoMatrixAdd: matrix row maps are not the same.");
1228 
1229  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1230 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1231  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1232  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1233 
1234  //FIXME is there a bug if beta=0?
1235  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1236 
1237  if (rv != 0)
1238  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1239  std::ostringstream buf;
1240 #else
1241  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1242 #endif
1243  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1244 #ifdef HAVE_XPETRA_TPETRA
1245 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1246  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1247  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1248 # else
1249  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1250  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1251 
1252  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1253 # endif
1254 #else
1255  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1256 #endif
1257  }
1258  } //MatrixMatrix::TwoMatrixAdd()
1259 
1260 
1273  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1274  const Matrix& B, bool transposeB, const SC& beta,
1275  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1276  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1277  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1278  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1279  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1280 
1281  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1282 
1283 
1284  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
1285  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
1286 
1287  if (C == Teuchos::null) {
1288  size_t maxNzInA = 0;
1289  size_t maxNzInB = 0;
1290  size_t numLocalRows = 0;
1291  if (A.isFillComplete() && B.isFillComplete()) {
1292 
1293  maxNzInA = A.getNodeMaxNumRowEntries();
1294  maxNzInB = B.getNodeMaxNumRowEntries();
1295  numLocalRows = A.getNodeNumRows();
1296  }
1297 
1298  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1299  // first check if either A or B has at most 1 nonzero per row
1300  // the case of both having at most 1 nz per row is handled by the ``else''
1301  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1302 
1303  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1304  for (size_t i = 0; i < numLocalRows; ++i)
1305  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1306 
1307  } else {
1308  for (size_t i = 0; i < numLocalRows; ++i)
1309  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1310  }
1311 
1312  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1313  << ", using static profiling" << std::endl;
1314  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow, Xpetra::StaticProfile));
1315 
1316  } else {
1317  // general case
1318  double nnzPerRowInA = Teuchos::as<double>(A.getNodeNumEntries()) / A.getNodeNumRows();
1319  double nnzPerRowInB = Teuchos::as<double>(B.getNodeNumEntries()) / B.getNodeNumRows();
1320  LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
1321 
1322  LO maxPossible = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
1323  //Use static profiling (more efficient) if the estimate is at least as big as the max
1324  //possible nnz's in any single row of the result.
1325  Xpetra::ProfileType pft = (maxPossible) > nnzToAllocate ? Xpetra::DynamicProfile : Xpetra::StaticProfile;
1326 
1327  fos << "nnzPerRowInA = " << nnzPerRowInA << ", nnzPerRowInB = " << nnzPerRowInB << std::endl;
1328  fos << "MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
1329  << ", max possible nnz per row in sum = " << maxPossible
1330  << ", using " << (pft == Xpetra::DynamicProfile ? "dynamic" : "static" ) << " profiling"
1331  << std::endl;
1332  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), nnzToAllocate, pft));
1333  }
1334  if (transposeB)
1335  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1336  }
1337 
1338  if (C->getRowMap()->lib() == Xpetra::UseEpetra) {
1339  #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1340  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1341  const Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(B);
1342  RCP<Epetra_CrsMatrix> epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
1343  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1344 
1345  //FIXME is there a bug if beta=0?
1346  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1347 
1348  if (rv != 0)
1349  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1350  #else
1351  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1352  #endif
1353  } else if (C->getRowMap()->lib() == Xpetra::UseTpetra) {
1354  #ifdef HAVE_XPETRA_TPETRA
1355  # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1356  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1357  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1358  # else
1359  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1360  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
1361  RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
1362 
1363  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
1364  # endif
1365  #else
1366  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
1367  #endif
1368  }
1369 
1371  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
1372  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
1374  }
1375  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
1376  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1377  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1378  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1379 
1380  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1381  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1382 
1383  size_t i = 0;
1384  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1385  RCP<Matrix> Cij = Teuchos::null;
1386  if(rcpA != Teuchos::null &&
1387  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1388  // recursive call
1389  TwoMatrixAdd(*rcpA, transposeA, alpha,
1390  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1391  Cij, fos, AHasFixedNnzPerRow);
1392  } else if (rcpA == Teuchos::null &&
1393  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1394  Cij = rcpBopB->getMatrix(i,j);
1395  } else if (rcpA != Teuchos::null &&
1396  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1397  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
1398  } else {
1399  Cij = Teuchos::null;
1400  }
1401 
1402  if (!Cij.is_null()) {
1403  if (Cij->isFillComplete())
1404  Cij->resumeFill();
1405  Cij->fillComplete();
1406  bC->setMatrix(i, j, Cij);
1407  } else {
1408  bC->setMatrix(i, j, Teuchos::null);
1409  }
1410  } // loop over columns j
1411  }
1412  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
1413  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
1414  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1415  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
1416 
1417  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1418  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1419 
1420  size_t j = 0;
1421  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1422  RCP<Matrix> Cij = Teuchos::null;
1423  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1424  rcpB!=Teuchos::null) {
1425  // recursive call
1426  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1427  *rcpB, transposeB, beta,
1428  Cij, fos, AHasFixedNnzPerRow);
1429  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1430  rcpB!=Teuchos::null) {
1431  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
1432  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1433  rcpB==Teuchos::null) {
1434  Cij = rcpBopA->getMatrix(i,j);
1435  } else {
1436  Cij = Teuchos::null;
1437  }
1438 
1439  if (!Cij.is_null()) {
1440  if (Cij->isFillComplete())
1441  Cij->resumeFill();
1442  Cij->fillComplete();
1443  bC->setMatrix(i, j, Cij);
1444  } else {
1445  bC->setMatrix(i, j, Teuchos::null);
1446  }
1447  } // loop over rows i
1448  }
1449  else {
1450  // This is the version for blocked matrices
1451 
1452  // check the compatibility of the blocked operators
1453  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1454  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1455  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)");
1456  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)");
1457  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)");
1458  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)");
1459  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)");
1460  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)");
1461 
1462  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1463  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1464 
1465  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1466  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1467 
1468  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1469  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1470 
1471  RCP<Matrix> Cij = Teuchos::null;
1472  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1473  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1474  // recursive call
1475 
1476  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1477  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1478  Cij, fos, AHasFixedNnzPerRow);
1479  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1480  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1481  Cij = rcpBopB->getMatrix(i,j);
1482  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1483  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1484  Cij = rcpBopA->getMatrix(i,j);
1485  } else {
1486  Cij = Teuchos::null;
1487  }
1488 
1489  if (!Cij.is_null()) {
1490  if (Cij->isFillComplete())
1491  Cij->resumeFill();
1492  //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
1493  Cij->fillComplete();
1494  bC->setMatrix(i, j, Cij);
1495  } else {
1496  bC->setMatrix(i, j, Teuchos::null);
1497  }
1498  } // loop over columns j
1499  } // loop over rows i
1500 
1501  } // end blocked recursive algorithm
1502  } //MatrixMatrix::TwoMatrixAdd()
1503  }; // end specialization on SC=double, GO=int and NO=EpetraNode
1504 
1505  // specialization MatrixMatrix for SC=double, GO=long long and NO=EptraNode
1506  template<>
1507  class MatrixMatrix<double,int,long long,EpetraNode> {
1508  typedef double Scalar;
1509  typedef int LocalOrdinal;
1510  typedef long long GlobalOrdinal;
1511  typedef EpetraNode Node;
1512 #include "Xpetra_UseShortNames.hpp"
1513 
1514  public:
1515 
1540  static void Multiply(const Matrix& A, bool transposeA,
1541  const Matrix& B, bool transposeB,
1542  Matrix& C,
1543  bool call_FillComplete_on_result = true,
1544  bool doOptimizeStorage = true,
1545  const std::string & label = std::string(),
1546  const RCP<ParameterList>& params = null) {
1547  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
1548  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
1549  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
1550  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
1551 
1552 
1553  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
1554  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
1555 
1556  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1557 
1558  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
1559 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1560  Epetra_CrsMatrix & epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(A);
1561  Epetra_CrsMatrix & epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1562  Epetra_CrsMatrix & epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
1563 
1564  int i = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, epC, haveMultiplyDoFillComplete);
1565  if (haveMultiplyDoFillComplete) {
1566  // Due to Epetra wrapper intricacies, we need to explicitly call
1567  // fillComplete on Xpetra matrix here. Specifically, EpetraCrsMatrix
1568  // only keeps an internal variable to check whether we are in resumed
1569  // state or not, but never touches the underlying Epetra object. As
1570  // such, we need to explicitly update the state of Xpetra matrix to
1571  // that of Epetra one afterwords
1572  C.fillComplete();
1573  }
1574 
1575  if (i != 0) {
1576  std::ostringstream buf;
1577  buf << i;
1578  std::string msg = "EpetraExt::MatrixMatrix::Multiply return value of " + buf.str();
1579  throw(Exceptions::RuntimeError(msg));
1580  }
1581 
1582 #else
1583  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
1584 #endif
1585  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
1586 #ifdef HAVE_XPETRA_TPETRA
1587 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1588  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1589  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long, EpetraNode> ETI enabled."));
1590 # else
1591  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1592  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
1593  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
1594 
1595  //18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
1596  //Previously, Tpetra's matrix matrix multiply did not support fillComplete.
1597  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1598 # endif
1599 #else
1600  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
1601 #endif
1602  }
1603 
1604  if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1605  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
1606  fillParams->set("Optimize Storage",doOptimizeStorage);
1607  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
1608  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
1609  fillParams);
1610  }
1611 
1612  // transfer striding information
1613  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1614  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1615  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
1616  } // end Multiply
1617 
1640  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1641  const Matrix& B, bool transposeB,
1642  RCP<Matrix> C_in,
1643  Teuchos::FancyOStream& fos,
1644  bool doFillComplete = true,
1645  bool doOptimizeStorage = true,
1646  const std::string & label = std::string(),
1647  const RCP<ParameterList>& params = null) {
1648 
1649  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1650  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1651 
1652  // Default case: Xpetra Multiply
1653  RCP<Matrix> C = C_in;
1654 
1655  if (C == Teuchos::null) {
1656  double nnzPerRow = Teuchos::as<double>(0);
1657 
1658 #if 0
1659  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
1660  // For now, follow what ML and Epetra do.
1661  GO numRowsA = A.getGlobalNumRows();
1662  GO numRowsB = B.getGlobalNumRows();
1663  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
1664  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
1665  nnzPerRow *= nnzPerRow;
1666  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1667  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1668  if (totalNnz < minNnz)
1669  totalNnz = minNnz;
1670  nnzPerRow = totalNnz / A.getGlobalNumRows();
1671 
1672  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1673  }
1674 #endif
1675  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1676  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1677 
1678  } else {
1679  C->resumeFill(); // why this is not done inside of Tpetra MxM?
1680  fos << "Reuse C pattern" << std::endl;
1681  }
1682 
1683  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1684 
1685  return C;
1686  }
1687 
1698  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1699  const Matrix& B, bool transposeB,
1700  Teuchos::FancyOStream &fos,
1701  bool callFillCompleteOnResult = true,
1702  bool doOptimizeStorage = true,
1703  const std::string & label = std::string(),
1704  const RCP<ParameterList>& params = null) {
1705  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1706  }
1707 
1708 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1709  // Michael Gee's MLMultiply
1710  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
1711  const Epetra_CrsMatrix& epB,
1712  Teuchos::FancyOStream& fos) {
1713  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1714  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1715  return Teuchos::null;
1716  }
1717 #endif //ifdef HAVE_XPETRA_EPETRAEXT
1718 
1729  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
1730  const BlockedCrsMatrix& B, bool transposeB,
1731  Teuchos::FancyOStream& fos,
1732  bool doFillComplete = true,
1733  bool doOptimizeStorage = true) {
1734  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1735  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1736 
1737  // Preconditions
1738  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1739  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1740 
1741  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1742  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
1743 
1744  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1745 
1746  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1747  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1748  RCP<Matrix> Cij;
1749 
1750  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1751  RCP<Matrix> crmat1 = A.getMatrix(i,l);
1752  RCP<Matrix> crmat2 = B.getMatrix(l,j);
1753 
1754  if (crmat1.is_null() || crmat2.is_null())
1755  continue;
1756  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1757  continue;
1758 
1759  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1760  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1761  TEUCHOS_TEST_FOR_EXCEPTION((crop1==Teuchos::null) != (crop2==Teuchos::null), Xpetra::Exceptions::RuntimeError, "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1762 
1763  // temporary matrix containing result of local block multiplication
1764  RCP<Matrix> temp = Teuchos::null;
1765 
1766  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1767  temp = Multiply (*crop1, false, *crop2, false, fos);
1768  else {
1769  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1770  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1771  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1772  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1773  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)");
1774  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)");
1775 
1776  // recursive multiplication call
1777  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1778 
1779  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1780  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)");
1781  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)");
1782  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)");
1783  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)");
1784  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)");
1785  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)");
1786  }
1787 
1788  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1789 
1790  RCP<Matrix> addRes = null;
1791  if (Cij.is_null ())
1792  Cij = temp;
1793  else {
1794  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1795  Cij = addRes;
1796  }
1797  }
1798 
1799  if (!Cij.is_null()) {
1800  if (Cij->isFillComplete())
1801  Cij->resumeFill();
1802  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1803  C->setMatrix(i, j, Cij);
1804  } else {
1805  C->setMatrix(i, j, Teuchos::null);
1806  }
1807  }
1808  }
1809 
1810  if (doFillComplete)
1811  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1812 
1813  return C;
1814  } // TwoMatrixMultiplyBlock
1815 
1828  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1829  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
1830  "TwoMatrixAdd: matrix row maps are not the same.");
1831 
1832  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1833 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1834  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1835  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1836 
1837  //FIXME is there a bug if beta=0?
1838  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1839 
1840  if (rv != 0)
1841  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1842  std::ostringstream buf;
1843 #else
1844  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1845 #endif
1846  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1847 #ifdef HAVE_XPETRA_TPETRA
1848 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1849  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1850  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
1851 # else
1852  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1853  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1854 
1855  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1856 # endif
1857 #else
1858  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1859 #endif
1860  }
1861  } //MatrixMatrix::TwoMatrixAdd()
1862 
1863 
1876  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1877  const Matrix& B, bool transposeB, const SC& beta,
1878  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1879  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1880  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1881  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1882  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1883 
1884  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1885  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
1886  "TwoMatrixAdd: matrix row maps are not the same.");
1887 
1888  if (C == Teuchos::null) {
1889  size_t maxNzInA = 0;
1890  size_t maxNzInB = 0;
1891  size_t numLocalRows = 0;
1892  if (A.isFillComplete() && B.isFillComplete()) {
1893  maxNzInA = A.getNodeMaxNumRowEntries();
1894  maxNzInB = B.getNodeMaxNumRowEntries();
1895  numLocalRows = A.getNodeNumRows();
1896  }
1897 
1898  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1899  // first check if either A or B has at most 1 nonzero per row
1900  // the case of both having at most 1 nz per row is handled by the ``else''
1901  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1902 
1903  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1904  for (size_t i = 0; i < numLocalRows; ++i)
1905  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1906 
1907  } else {
1908  for (size_t i = 0; i < numLocalRows; ++i)
1909  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1910  }
1911 
1912  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1913  << ", using static profiling" << std::endl;
1914  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow, Xpetra::StaticProfile));
1915 
1916  } else {
1917  // general case
1918  double nnzPerRowInA = Teuchos::as<double>(A.getNodeNumEntries()) / A.getNodeNumRows();
1919  double nnzPerRowInB = Teuchos::as<double>(B.getNodeNumEntries()) / B.getNodeNumRows();
1920  LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
1921 
1922  LO maxPossible = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
1923  //Use static profiling (more efficient) if the estimate is at least as big as the max
1924  //possible nnz's in any single row of the result.
1925  Xpetra::ProfileType pft = (maxPossible) > nnzToAllocate ? Xpetra::DynamicProfile : Xpetra::StaticProfile;
1926 
1927  fos << "nnzPerRowInA = " << nnzPerRowInA << ", nnzPerRowInB = " << nnzPerRowInB << std::endl;
1928  fos << "MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
1929  << ", max possible nnz per row in sum = " << maxPossible
1930  << ", using " << (pft == Xpetra::DynamicProfile ? "dynamic" : "static" ) << " profiling"
1931  << std::endl;
1932  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), nnzToAllocate, pft));
1933  }
1934  if (transposeB)
1935  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1936  }
1937 
1938  if (C->getRowMap()->lib() == Xpetra::UseEpetra) {
1939  #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1940  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1941  const Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(B);
1942  RCP<Epetra_CrsMatrix> epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
1943  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1944 
1945  //FIXME is there a bug if beta=0?
1946  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1947 
1948  if (rv != 0)
1949  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1950  #else
1951  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1952  #endif
1953  } else if (C->getRowMap()->lib() == Xpetra::UseTpetra) {
1954  #ifdef HAVE_XPETRA_TPETRA
1955  # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1956  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1957  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
1958  # else
1959  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1960  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
1961  RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
1962 
1963  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
1964  # endif
1965  #else
1966  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
1967  #endif
1968  }
1969 
1971  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
1972  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
1974  }
1975  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
1976  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1977  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1978  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1979 
1980  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1981  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1982 
1983  size_t i = 0;
1984  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1985  RCP<Matrix> Cij = Teuchos::null;
1986  if(rcpA != Teuchos::null &&
1987  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1988  // recursive call
1989  TwoMatrixAdd(*rcpA, transposeA, alpha,
1990  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1991  Cij, fos, AHasFixedNnzPerRow);
1992  } else if (rcpA == Teuchos::null &&
1993  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1994  Cij = rcpBopB->getMatrix(i,j);
1995  } else if (rcpA != Teuchos::null &&
1996  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1997  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
1998  } else {
1999  Cij = Teuchos::null;
2000  }
2001 
2002  if (!Cij.is_null()) {
2003  if (Cij->isFillComplete())
2004  Cij->resumeFill();
2005  Cij->fillComplete();
2006  bC->setMatrix(i, j, Cij);
2007  } else {
2008  bC->setMatrix(i, j, Teuchos::null);
2009  }
2010  } // loop over columns j
2011  }
2012  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
2013  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
2014  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2015  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
2016 
2017  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2018  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2019 
2020  size_t j = 0;
2021  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2022  RCP<Matrix> Cij = Teuchos::null;
2023  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2024  rcpB!=Teuchos::null) {
2025  // recursive call
2026  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2027  *rcpB, transposeB, beta,
2028  Cij, fos, AHasFixedNnzPerRow);
2029  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2030  rcpB!=Teuchos::null) {
2031  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
2032  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2033  rcpB==Teuchos::null) {
2034  Cij = rcpBopA->getMatrix(i,j);
2035  } else {
2036  Cij = Teuchos::null;
2037  }
2038 
2039  if (!Cij.is_null()) {
2040  if (Cij->isFillComplete())
2041  Cij->resumeFill();
2042  Cij->fillComplete();
2043  bC->setMatrix(i, j, Cij);
2044  } else {
2045  bC->setMatrix(i, j, Teuchos::null);
2046  }
2047  } // loop over rows i
2048  }
2049  else {
2050  // This is the version for blocked matrices
2051 
2052  // check the compatibility of the blocked operators
2053  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2054  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2055  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)");
2056  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)");
2057  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)");
2058  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)");
2059  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)");
2060  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)");
2061 
2062  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2063  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2064 
2065  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2066  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2067 
2068  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2069  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2070 
2071  RCP<Matrix> Cij = Teuchos::null;
2072  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2073  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2074  // recursive call
2075 
2076  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2077  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2078  Cij, fos, AHasFixedNnzPerRow);
2079  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2080  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2081  Cij = rcpBopB->getMatrix(i,j);
2082  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2083  rcpBopB->getMatrix(i,j)==Teuchos::null) {
2084  Cij = rcpBopA->getMatrix(i,j);
2085  } else {
2086  Cij = Teuchos::null;
2087  }
2088 
2089  if (!Cij.is_null()) {
2090  if (Cij->isFillComplete())
2091  Cij->resumeFill();
2092  //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
2093  Cij->fillComplete();
2094  bC->setMatrix(i, j, Cij);
2095  } else {
2096  bC->setMatrix(i, j, Teuchos::null);
2097  }
2098  } // loop over columns j
2099  } // loop over rows i
2100 
2101  } // end blocked recursive algorithm
2102  } //MatrixMatrix::TwoMatrixAdd()
2103  }; // end specialization on GO=long long and NO=EpetraNode
2104 
2105 #endif // HAVE_XPETRA_EPETRA
2106 
2107 } // end namespace Xpetra
2108 
2109 #define XPETRA_MATRIXMATRIX_SHORT
2110 
2111 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_ */
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::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:537
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:1273
Xpetra::Matrix::getGlobalNumRows
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
Xpetra::Matrix::getRowMap
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Xpetra_Matrix.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:716
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:1640
GO
GlobalOrdinal GO
Definition: Xpetra_UseShortNamesOrdinal.hpp:139
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Xpetra::Helpers::Op2NonConstEpetraCrs
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:142
Xpetra_UseShortNames.hpp
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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::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:1729
Xpetra_Matrix.hpp
Xpetra_StridedMapFactory.hpp
Teuchos::RCP::is_null
bool is_null() const
Xpetra::Helpers::Op2EpetraCrs
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:126
Teuchos::TimeMonitor::getNewTimer
static RCP< Time > getNewTimer(const std::string &name)
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1081
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:1876
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
Teuchos::basic_FancyOStream
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
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:787
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:886
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.
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:1540
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:506
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
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:174
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:1225
Copy
Copy
Xpetra::BlockedCrsMatrix::getRangeMap
RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1072
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:1303
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:981
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1084
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:1126
Xpetra::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1267
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:1270
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:969
Xpetra::BlockedCrsMatrix::getDomainMap
RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1057
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:1710
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:1828
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
Teuchos::toString
std::string toString(const T &t)
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:1698
Xpetra_TpetraCrsMatrix.hpp
Generated by   doxygen 1.8.13