Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 #include <execinfo.h>
63 
64 #ifdef HAVE_XPETRA_EPETRA
65 #include <Xpetra_EpetraCrsMatrix_fwd.hpp>
66 #endif
67 
68 #ifdef HAVE_XPETRA_EPETRAEXT
69 #include <EpetraExt_MatrixMatrix.h>
70 #include <EpetraExt_RowMatrixOut.h>
71 #include <Epetra_RowMatrixTransposer.h>
72 #endif // HAVE_XPETRA_EPETRAEXT
73 
74 #ifdef HAVE_XPETRA_TPETRA
75 #include <TpetraExt_MatrixMatrix.hpp>
76 #include <Tpetra_RowMatrixTransposer.hpp>
77 #include <MatrixMarket_Tpetra.hpp>
78 #include <Xpetra_TpetraCrsMatrix.hpp>
79 #include <Xpetra_TpetraMultiVector.hpp>
80 #include <Xpetra_TpetraVector.hpp>
81 #endif // HAVE_XPETRA_TPETRA
82 
83 namespace Xpetra {
84 
91  template <class Scalar,
92  class LocalOrdinal = int,
93  class GlobalOrdinal = LocalOrdinal,
94  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
95  class Helpers {
96 #include "Xpetra_UseShortNames.hpp"
97 
98  public:
99 
100 #ifdef HAVE_XPETRA_EPETRA
101  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<Matrix> Op) {
102  // Get the underlying Epetra Mtx
103  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
104  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast,
105  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
106 
107  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
108  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
109  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast,
110  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
111 
112  return tmp_ECrsMtx->getEpetra_CrsMatrix();
113  }
114 
115  static RCP<Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) {
116  RCP<Epetra_CrsMatrix> A;
117  // Get the underlying Epetra Mtx
118  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
119  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast,
120  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
121 
122  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
123  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
124  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
125 
126  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
127  }
128 
129  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) {
130  // Get the underlying Epetra Mtx
131  try {
132  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
133  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
134  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
135  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast,
136  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
137 
138  return *tmp_ECrsMtx->getEpetra_CrsMatrix();
139 
140  } catch(...) {
141  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
142  }
143  }
144 
145  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(const Matrix& Op) {
146  RCP<Epetra_CrsMatrix> A;
147  // Get the underlying Epetra Mtx
148  try {
149  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
150  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
151  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
152  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
153 
154  return *Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
155 
156  } catch(...) {
157  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
158  }
159  }
160 #endif // HAVE_XPETRA_EPETRA
161 
162 #ifdef HAVE_XPETRA_TPETRA
163  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<Matrix> Op) {
164  // Get the underlying Tpetra Mtx
165  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
166  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
167 
168  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
169  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
170  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
171 
172  return tmp_ECrsMtx->getTpetra_CrsMatrix();
173  }
174 
175  static RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) {
176  // Get the underlying Tpetra Mtx
177  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
178  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
179  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
180  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
181  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
182 
183  return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
184  }
185 
186  static const Tpetra::CrsMatrix<SC,LO,GO,NO> & Op2TpetraCrs(const Matrix& Op) {
187  // Get the underlying Tpetra Mtx
188  try{
189  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
190  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
191  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
192  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
193 
194  return *tmp_TCrsMtx->getTpetra_CrsMatrix();
195 
196  } catch (...) {
197  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
198  }
199  }
200 
201  static Tpetra::CrsMatrix<SC,LO,GO,NO> & Op2NonConstTpetraCrs(const Matrix& Op) {
202  // Get the underlying Tpetra Mtx
203  try{
204  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap& >(Op);
205  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
206  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
207  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
208 
209  return *Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_TCrsMtx->getTpetra_CrsMatrix());
210 
211  } catch (...) {
212  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
213  }
214  }
215 #endif // HAVE_XPETRA_TPETRA
216 
217 #ifdef HAVE_XPETRA_TPETRA
218  using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
219  static Teuchos::RCP<Matrix> tpetraAdd(
220  const tcrs_matrix_type& A, bool transposeA, const typename tcrs_matrix_type::scalar_type alpha,
221  const tcrs_matrix_type& B, bool transposeB, const typename tcrs_matrix_type::scalar_type beta)
222  {
223  using Teuchos::Array;
224  using Teuchos::RCP;
225  using Teuchos::rcp;
226  using Teuchos::rcp_implicit_cast;
227  using Teuchos::rcpFromRef;
228  using Teuchos::ParameterList;
229  using XTCrsType = Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>;
230  using CrsType = Xpetra::CrsMatrix<SC,LO,GO,NO>;
231  using CrsWrap = Xpetra::CrsMatrixWrap<SC,LO,GO,NO>;
232  using transposer_type = Tpetra::RowMatrixTransposer<SC,LO,GO,NO>;
233  using import_type = Tpetra::Import<LO,GO,NO>;
234  RCP<const tcrs_matrix_type> Aprime = rcpFromRef(A);
235  if(transposeA)
236  Aprime = transposer_type(Aprime).createTranspose();
237  //Decide whether the fast code path can be taken.
238  if(A.isFillComplete() && B.isFillComplete())
239  {
240  RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), 0));
241  RCP<ParameterList> addParams = rcp(new ParameterList);
242  addParams->set("Call fillComplete", false);
243  //passing A after B means C will have the same domain/range map as A (or A^T if transposeA)
244  Tpetra::MatrixMatrix::add<SC,LO,GO,NO>(beta, transposeB, B, alpha, false, *Aprime, *C, Teuchos::null, Teuchos::null, addParams);
245  return rcp_implicit_cast<Matrix>(rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C))))));
246  }
247  else
248  {
249  //Slow case - one or both operands are non-fill complete.
250  //TODO: deprecate this.
251  //Need to compute the explicit transpose before add if transposeA and/or transposeB.
252  //This is to count the number of entries to allocate per row in the final sum.
253  RCP<const tcrs_matrix_type> Bprime = rcpFromRef(B);
254  if(transposeB)
255  Bprime = transposer_type(Bprime).createTranspose();
256  //Use Aprime's row map as C's row map.
257  if(!(Aprime->getRowMap()->isSameAs(*(Bprime->getRowMap()))))
258  {
259  auto import = rcp(new import_type(Bprime->getRowMap(), Aprime->getRowMap()));
260  Bprime = Tpetra::importAndFillCompleteCrsMatrix<tcrs_matrix_type>(Bprime, *import, Aprime->getDomainMap(), Aprime->getRangeMap());
261  }
262  //Count the entries to allocate for C in each local row.
263  LO numLocalRows = Aprime->getNodeNumRows();
264  Array<size_t> allocPerRow(numLocalRows);
265  //0 is always the minimum LID
266  for(LO i = 0; i < numLocalRows; i++)
267  {
268  allocPerRow[i] = Aprime->getNumEntriesInLocalRow(i) + Bprime->getNumEntriesInLocalRow(i);
269  }
270  //Construct C
271  RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), allocPerRow(), Tpetra::StaticProfile));
272  //Compute the sum in C (already took care of transposes)
273  Tpetra::MatrixMatrix::Add<SC,LO,GO,NO>(
274  *Aprime, false, alpha,
275  *Bprime, false, beta,
276  C);
277  return rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C)))));
278  }
279  }
280 #endif
281 
282 #ifdef HAVE_XPETRA_EPETRAEXT
283  static void epetraExtMult(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Matrix& C, bool fillCompleteResult)
284  {
285  Epetra_CrsMatrix& epA = Op2NonConstEpetraCrs(A);
286  Epetra_CrsMatrix& epB = Op2NonConstEpetraCrs(B);
287  Epetra_CrsMatrix& epC = Op2NonConstEpetraCrs(C);
288  //Want a static profile (possibly fill complete) matrix as the result.
289  //But, EpetraExt Multiply needs C to be dynamic profile, so compute the product in a temporary DynamicProfile matrix.
290  const Epetra_Map& Crowmap = epC.RowMap();
291  int errCode = 0;
292  Epetra_CrsMatrix Ctemp(::Copy, Crowmap, 0, false);
293  if(fillCompleteResult) {
294  errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, true);
295  if(!errCode) {
296  epC = Ctemp;
297  C.fillComplete();
298  }
299  }
300  else {
301  errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, false);
302  if(!errCode) {
303  int numLocalRows = Crowmap.NumMyElements();
304  long long* globalElementList = nullptr;
305  Crowmap.MyGlobalElementsPtr(globalElementList);
306  Teuchos::Array<int> entriesPerRow(numLocalRows, 0);
307  for(int i = 0; i < numLocalRows; i++)
308  {
309  entriesPerRow[i] = Ctemp.NumGlobalEntries(globalElementList[i]);
310  }
311  //know how many entries to allocate in epC (which must be static profile)
312  epC = Epetra_CrsMatrix(::Copy, Crowmap, entriesPerRow.data(), true);
313  for(int i = 0; i < numLocalRows; i++)
314  {
315  int gid = globalElementList[i];
316  int numEntries;
317  double* values;
318  int* inds;
319  Ctemp.ExtractGlobalRowView(gid, numEntries, values, inds);
320  epC.InsertGlobalValues(gid, numEntries, values, inds);
321  }
322  }
323  }
324  if(errCode) {
325  std::ostringstream buf;
326  buf << errCode;
327  std::string msg = "EpetraExt::MatrixMatrix::Multiply returned nonzero error code " + buf.str();
328  throw(Exceptions::RuntimeError(msg));
329  }
330  }
331 #endif
332  };
333 
334  template <class Scalar,
335  class LocalOrdinal /*= int*/,
336  class GlobalOrdinal /*= LocalOrdinal*/,
337  class Node /*= KokkosClassic::DefaultNode::DefaultNodeType*/>
338  class MatrixMatrix {
339 #undef XPETRA_MATRIXMATRIX_SHORT
340 #include "Xpetra_UseShortNames.hpp"
341 
342  public:
343 
368  static void Multiply(const Matrix& A, bool transposeA,
369  const Matrix& B, bool transposeB,
370  Matrix& C,
371  bool call_FillComplete_on_result = true,
372  bool doOptimizeStorage = true,
373  const std::string & label = std::string(),
374  const RCP<ParameterList>& params = null) {
375 
376  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
377  Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
378  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
379  Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
380 
381  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
382  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
383 
384  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
385 
386  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
387 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
388  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
389 #else
390  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
391 
392 #endif
393  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
394 #ifdef HAVE_XPETRA_TPETRA
395  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
396  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
397  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
398 
399  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
400  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
401  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
402 #else
403  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
404 #endif
405  }
406 
407  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
408  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
409  fillParams->set("Optimize Storage", doOptimizeStorage);
410  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
411  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
412  fillParams);
413  }
414 
415  // transfer striding information
416  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
417  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
418  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
419  } // end Multiply
420 
443  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, RCP<Matrix> C_in,
444  Teuchos::FancyOStream& fos,
445  bool doFillComplete = true,
446  bool doOptimizeStorage = true,
447  const std::string & label = std::string(),
448  const RCP<ParameterList>& params = null) {
449 
450  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
451  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
452 
453  // Default case: Xpetra Multiply
454  RCP<Matrix> C = C_in;
455 
456  if (C == Teuchos::null) {
457  double nnzPerRow = Teuchos::as<double>(0);
458 
459 #if 0
460  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
461  // For now, follow what ML and Epetra do.
462  GO numRowsA = A.getGlobalNumRows();
463  GO numRowsB = B.getGlobalNumRows();
464  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
465  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
466  nnzPerRow *= nnzPerRow;
467  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
468  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
469  if (totalNnz < minNnz)
470  totalNnz = minNnz;
471  nnzPerRow = totalNnz / A.getGlobalNumRows();
472 
473  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
474  }
475 #endif
476  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
477  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
478 
479  } else {
480  C->resumeFill(); // why this is not done inside of Tpetra MxM?
481  fos << "Reuse C pattern" << std::endl;
482  }
483 
484  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
485 
486  return C;
487  }
488 
499  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Teuchos::FancyOStream &fos,
500  bool callFillCompleteOnResult = true, bool doOptimizeStorage = true, const std::string& label = std::string(),
501  const RCP<ParameterList>& params = null) {
502  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
503  }
504 
505 #ifdef HAVE_XPETRA_EPETRAEXT
506  // Michael Gee's MLMultiply
507  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
508  const Epetra_CrsMatrix& epB,
509  Teuchos::FancyOStream& fos) {
510  throw(Xpetra::Exceptions::RuntimeError("MLTwoMatrixMultiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
511  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
512  }
513 #endif //ifdef HAVE_XPETRA_EPETRAEXT
514 
525  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
526  const BlockedCrsMatrix& B, bool transposeB,
527  Teuchos::FancyOStream& fos,
528  bool doFillComplete = true,
529  bool doOptimizeStorage = true) {
530  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
531  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
532 
533  // Preconditions
534  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
535  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
536 
537  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
538  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
539 
540  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
541 
542  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
543  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
544  RCP<Matrix> Cij;
545 
546  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
547  RCP<Matrix> crmat1 = A.getMatrix(i,l);
548  RCP<Matrix> crmat2 = B.getMatrix(l,j);
549 
550  if (crmat1.is_null() || crmat2.is_null())
551  continue;
552 
553  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
554  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
555  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
556  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
557 
558  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
559  if (!crop1.is_null())
560  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
561  if (!crop2.is_null())
562  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
563 
564  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
565  "crmat1 does not have global constants");
566  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
567  "crmat2 does not have global constants");
568 
569  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
570  continue;
571 
572  // temporary matrix containing result of local block multiplication
573  RCP<Matrix> temp = Teuchos::null;
574 
575  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
576  temp = Multiply (*crop1, false, *crop2, false, fos);
577  else {
578  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
579  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
580  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
581  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
582  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)");
583  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)");
584 
585  // recursive multiplication call
586  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
587 
588  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
589  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)");
590  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)");
591  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)");
592  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)");
593  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)");
594  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)");
595  }
596 
597  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
598 
599  RCP<Matrix> addRes = null;
600  if (Cij.is_null ())
601  Cij = temp;
602  else {
603  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
604  Cij = addRes;
605  }
606  }
607 
608  if (!Cij.is_null()) {
609  if (Cij->isFillComplete())
610  Cij->resumeFill();
611  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
612  C->setMatrix(i, j, Cij);
613  } else {
614  C->setMatrix(i, j, Teuchos::null);
615  }
616  }
617  }
618 
619  if (doFillComplete)
620  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
621 
622  return C;
623  } // TwoMatrixMultiplyBlock
624 
637  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
638  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
639  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
640 
641  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
642  throw Exceptions::RuntimeError("TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
643  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
644 #ifdef HAVE_XPETRA_TPETRA
645  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
646  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
647 
648  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
649 #else
650  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
651 #endif
652  }
653  } //MatrixMatrix::TwoMatrixAdd()
654 
655 
668  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
669  const Matrix& B, bool transposeB, const SC& beta,
670  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
671 
672  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
673  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
674  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
675  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
676  // We have to distinguish 4 cases:
677  // both matrices are CrsMatrixWrap based, only one of them is CrsMatrixWrap based, or none.
678 
679  // C can be null, so just use A to get the lib
680  auto lib = A.getRowMap()->lib();
681 
682  // Both matrices are CrsMatrixWrap
683  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
684  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
685  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
686  if (lib == Xpetra::UseEpetra) {
687  throw Exceptions::RuntimeError("MatrixMatrix::Add for Epetra only available with Scalar = double, LO = GO = int.");
688  } else if (lib == Xpetra::UseTpetra) {
689  #ifdef HAVE_XPETRA_TPETRA
690  using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
691  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
692  const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
693  const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
694  C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
695  #else
696  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
697  #endif
698  }
700  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
701  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
703  }
704  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
705  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
706  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
707  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
708 
709  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
710  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
711 
712  size_t i = 0;
713  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
714  RCP<Matrix> Cij = Teuchos::null;
715  if(rcpA != Teuchos::null &&
716  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
717  // recursive call
718  TwoMatrixAdd(*rcpA, transposeA, alpha,
719  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
720  Cij, fos, AHasFixedNnzPerRow);
721  } else if (rcpA == Teuchos::null &&
722  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
723  Cij = rcpBopB->getMatrix(i,j);
724  } else if (rcpA != Teuchos::null &&
725  rcpBopB->getMatrix(i,j)==Teuchos::null) {
726  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
727  } else {
728  Cij = Teuchos::null;
729  }
730 
731  if (!Cij.is_null()) {
732  if (Cij->isFillComplete())
733  Cij->resumeFill();
734  Cij->fillComplete();
735  bC->setMatrix(i, j, Cij);
736  } else {
737  bC->setMatrix(i, j, Teuchos::null);
738  }
739  } // loop over columns j
740  }
741  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
742  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
743  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
744  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
745 
746  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
747  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
748  size_t j = 0;
749  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
750  RCP<Matrix> Cij = Teuchos::null;
751  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
752  rcpB!=Teuchos::null) {
753  // recursive call
754  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
755  *rcpB, transposeB, beta,
756  Cij, fos, AHasFixedNnzPerRow);
757  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
758  rcpB!=Teuchos::null) {
759  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
760  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
761  rcpB==Teuchos::null) {
762  Cij = rcpBopA->getMatrix(i,j);
763  } else {
764  Cij = Teuchos::null;
765  }
766 
767  if (!Cij.is_null()) {
768  if (Cij->isFillComplete())
769  Cij->resumeFill();
770  Cij->fillComplete();
771  bC->setMatrix(i, j, Cij);
772  } else {
773  bC->setMatrix(i, j, Teuchos::null);
774  }
775  } // loop over rows i
776  }
777  else {
778  // This is the version for blocked matrices
779 
780  // check the compatibility of the blocked operators
781  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
782  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
783  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)");
784  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)");
785  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)");
786  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)");
787  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)");
788  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)");
789 
790  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
791  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
792 
793  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
794  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
795  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
796  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
797 
798  RCP<Matrix> Cij = Teuchos::null;
799  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
800  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
801  // recursive call
802  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
803  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
804  Cij, fos, AHasFixedNnzPerRow);
805  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
806  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
807  Cij = rcpBopB->getMatrix(i,j);
808  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
809  rcpBopB->getMatrix(i,j)==Teuchos::null) {
810  Cij = rcpBopA->getMatrix(i,j);
811  } else {
812  Cij = Teuchos::null;
813  }
814 
815  if (!Cij.is_null()) {
816  if (Cij->isFillComplete())
817  Cij->resumeFill();
818  Cij->fillComplete();
819  bC->setMatrix(i, j, Cij);
820  } else {
821  bC->setMatrix(i, j, Teuchos::null);
822  }
823  } // loop over columns j
824  } // loop over rows i
825 
826  } // end blocked recursive algorithm
827  } //MatrixMatrix::TwoMatrixAdd()
828 
829 
830  }; // class MatrixMatrix
831 
832 
833 #ifdef HAVE_XPETRA_EPETRA
834  // specialization MatrixMatrix for SC=double, LO=GO=int
835  template<>
836  class MatrixMatrix<double,int,int,EpetraNode> {
837  typedef double Scalar;
838  typedef int LocalOrdinal;
839  typedef int GlobalOrdinal;
840  typedef EpetraNode Node;
841 #include "Xpetra_UseShortNames.hpp"
842 
843  public:
844 
869  static void Multiply(const Matrix& A, bool transposeA,
870  const Matrix& B, bool transposeB,
871  Matrix& C,
872  bool call_FillComplete_on_result = true,
873  bool doOptimizeStorage = true,
874  const std::string & label = std::string(),
875  const RCP<ParameterList>& params = null) {
876  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
877  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
878  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
879  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
880 
881  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
882  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
883 
884  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
885 
886  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
887 
888  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
889 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
890  helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
891 #else
892  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
893 #endif
894  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
895 #ifdef HAVE_XPETRA_TPETRA
896 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
897  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
898  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
899 # else
900  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
901  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
902  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
903 
904  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
905  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
906  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
907 # endif
908 #else
909  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
910 #endif
911  }
912 
913  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
914  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
915  fillParams->set("Optimize Storage",doOptimizeStorage);
916  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
917  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
918  fillParams);
919  }
920 
921  // transfer striding information
922  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
923  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
924  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
925  } // end Multiply
926 
949  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
950  const Matrix& B, bool transposeB,
951  RCP<Matrix> C_in,
952  Teuchos::FancyOStream& fos,
953  bool doFillComplete = true,
954  bool doOptimizeStorage = true,
955  const std::string & label = std::string(),
956  const RCP<ParameterList>& params = null) {
957 
958  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
959  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
960 
961  // Optimization using ML Multiply when available and requested
962  // This feature is currently not supported. We would have to introduce the HAVE_XPETRA_ML_MMM flag
963 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) && defined(HAVE_XPETRA_ML_MMM)
964  if (B.getDomainMap()->lib() == Xpetra::UseEpetra && !transposeA && !transposeB) {
965  RCP<const Epetra_CrsMatrix> epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(rcpFromRef(A));
966  RCP<const Epetra_CrsMatrix> epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(rcpFromRef(B));
967  RCP<Epetra_CrsMatrix> epC = MLTwoMatrixMultiply(*epA, *epB, fos);
968 
969  RCP<Matrix> C = Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO> (epC);
970  if (doFillComplete) {
971  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
972  fillParams->set("Optimize Storage", doOptimizeStorage);
973  C->fillComplete(B.getDomainMap(), A.getRangeMap(), fillParams);
974  }
975 
976  // Fill strided maps information
977  // This is necessary since the ML matrix matrix multiplication routine has no handling for this
978  // TODO: move this call to MLMultiply...
979  C->CreateView("stridedMaps", rcpFromRef(A), transposeA, rcpFromRef(B), transposeB);
980 
981  return C;
982  }
983 #endif // EPETRA + EPETRAEXT + ML
984 
985  // Default case: Xpetra Multiply
986  RCP<Matrix> C = C_in;
987 
988  if (C == Teuchos::null) {
989  double nnzPerRow = Teuchos::as<double>(0);
990 
991 #if 0
992  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
993  // For now, follow what ML and Epetra do.
994  GO numRowsA = A.getGlobalNumRows();
995  GO numRowsB = B.getGlobalNumRows();
996  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
997  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
998  nnzPerRow *= nnzPerRow;
999  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1000  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1001  if (totalNnz < minNnz)
1002  totalNnz = minNnz;
1003  nnzPerRow = totalNnz / A.getGlobalNumRows();
1004 
1005  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1006  }
1007 #endif
1008 
1009  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1010  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1011 
1012  } else {
1013  C->resumeFill(); // why this is not done inside of Tpetra MxM?
1014  fos << "Reuse C pattern" << std::endl;
1015  }
1016 
1017  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1018 
1019  return C;
1020  }
1021 
1032  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1033  const Matrix& B, bool transposeB,
1034  Teuchos::FancyOStream &fos,
1035  bool callFillCompleteOnResult = true,
1036  bool doOptimizeStorage = true,
1037  const std::string & label = std::string(),
1038  const RCP<ParameterList>& params = null) {
1039  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1040  }
1041 
1042 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1043  // Michael Gee's MLMultiply
1044  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
1045  const Epetra_CrsMatrix& epB,
1046  Teuchos::FancyOStream& fos) {
1047 #if defined(HAVE_XPETRA_ML_MMM) // Note: this is currently not supported
1048  ML_Comm* comm;
1049  ML_Comm_Create(&comm);
1050  fos << "****** USING ML's MATRIX MATRIX MULTIPLY ******" << std::endl;
1051 #ifdef HAVE_MPI
1052  // ML_Comm uses MPI_COMM_WORLD, so try to use the same communicator as epA.
1053  const Epetra_MpiComm * Mcomm = dynamic_cast<const Epetra_MpiComm*>(&(epA.Comm()));
1054  if (Mcomm)
1055  ML_Comm_Set_UsrComm(comm,Mcomm->GetMpiComm());
1056 # endif
1057  //in order to use ML, there must be no indices missing from the matrix column maps.
1058  EpetraExt::CrsMatrix_SolverMap Atransform;
1059  EpetraExt::CrsMatrix_SolverMap Btransform;
1060  const Epetra_CrsMatrix& A = Atransform(const_cast<Epetra_CrsMatrix&>(epA));
1061  const Epetra_CrsMatrix& B = Btransform(const_cast<Epetra_CrsMatrix&>(epB));
1062 
1063  if (!A.Filled()) throw Exceptions::RuntimeError("A has to be FillCompleted");
1064  if (!B.Filled()) throw Exceptions::RuntimeError("B has to be FillCompleted");
1065 
1066  // create ML operators from EpetraCrsMatrix
1067  ML_Operator* ml_As = ML_Operator_Create(comm);
1068  ML_Operator* ml_Bs = ML_Operator_Create(comm);
1069  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?
1070  ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&B),ml_Bs);
1071  ML_Operator* ml_AtimesB = ML_Operator_Create(comm);
1072  {
1073  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("ML_2matmult kernel"));
1074  ML_2matmult(ml_As,ml_Bs,ml_AtimesB,ML_CSR_MATRIX); // do NOT use ML_EpetraCRS_MATRIX!!!
1075  }
1076  ML_Operator_Destroy(&ml_As);
1077  ML_Operator_Destroy(&ml_Bs);
1078 
1079  // For ml_AtimesB we have to reconstruct the column map in global indexing,
1080  // The following is going down to the salt-mines of ML ...
1081  // note: we use integers, since ML only works for Epetra...
1082  int N_local = ml_AtimesB->invec_leng;
1083  ML_CommInfoOP* getrow_comm = ml_AtimesB->getrow->pre_comm;
1084  if (!getrow_comm) throw(Exceptions::RuntimeError("ML_Operator does not have a CommInfo"));
1085  ML_Comm* comm_AB = ml_AtimesB->comm; // get comm object
1086  if (N_local != B.DomainMap().NumMyElements())
1087  throw(Exceptions::RuntimeError("Mismatch in local row dimension between ML and Epetra"));
1088  int N_rcvd = 0;
1089  int N_send = 0;
1090  int flag = 0;
1091  for (int i = 0; i < getrow_comm->N_neighbors; i++) {
1092  N_rcvd += (getrow_comm->neighbors)[i].N_rcv;
1093  N_send += (getrow_comm->neighbors)[i].N_send;
1094  if ( ((getrow_comm->neighbors)[i].N_rcv != 0) &&
1095  ((getrow_comm->neighbors)[i].rcv_list != NULL) ) flag = 1;
1096  }
1097  // For some unknown reason, ML likes to have stuff one larger than
1098  // neccessary...
1099  std::vector<double> dtemp(N_local + N_rcvd + 1); // "double" vector for comm function
1100  std::vector<int> cmap (N_local + N_rcvd + 1); // vector for gids
1101  for (int i = 0; i < N_local; ++i) {
1102  cmap[i] = B.DomainMap().GID(i);
1103  dtemp[i] = (double) cmap[i];
1104  }
1105  ML_cheap_exchange_bdry(&dtemp[0],getrow_comm,N_local,N_send,comm_AB); // do communication
1106  if (flag) { // process received data
1107  int count = N_local;
1108  const int neighbors = getrow_comm->N_neighbors;
1109  for (int i = 0; i < neighbors; i++) {
1110  const int nrcv = getrow_comm->neighbors[i].N_rcv;
1111  for (int j = 0; j < nrcv; j++)
1112  cmap[getrow_comm->neighbors[i].rcv_list[j]] = (int) dtemp[count++];
1113  }
1114  } else {
1115  for (int i = 0; i < N_local+N_rcvd; ++i)
1116  cmap[i] = (int)dtemp[i];
1117  }
1118  dtemp.clear(); // free double array
1119 
1120  // we can now determine a matching column map for the result
1121  Epetra_Map gcmap(-1, N_local+N_rcvd, &cmap[0], B.ColMap().IndexBase(), A.Comm());
1122 
1123  int allocated = 0;
1124  int rowlength;
1125  double* val = NULL;
1126  int* bindx = NULL;
1127 
1128  const int myrowlength = A.RowMap().NumMyElements();
1129  const Epetra_Map& rowmap = A.RowMap();
1130 
1131  // Determine the maximum bandwith for the result matrix.
1132  // replaces the old, very(!) memory-consuming guess:
1133  // int guessnpr = A.MaxNumEntries()*B.MaxNumEntries();
1134  int educatedguess = 0;
1135  for (int i = 0; i < myrowlength; ++i) {
1136  // get local row
1137  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1138  if (rowlength>educatedguess)
1139  educatedguess = rowlength;
1140  }
1141 
1142  // allocate our result matrix and fill it
1143  RCP<Epetra_CrsMatrix> result = rcp(new Epetra_CrsMatrix(::Copy, A.RangeMap(), gcmap, educatedguess, false));
1144 
1145  std::vector<int> gcid(educatedguess);
1146  for (int i = 0; i < myrowlength; ++i) {
1147  const int grid = rowmap.GID(i);
1148  // get local row
1149  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1150  if (!rowlength) continue;
1151  if ((int)gcid.size() < rowlength) gcid.resize(rowlength);
1152  for (int j = 0; j < rowlength; ++j) {
1153  gcid[j] = gcmap.GID(bindx[j]);
1154  if (gcid[j] < 0)
1155  throw Exceptions::RuntimeError("Error: cannot find gcid!");
1156  }
1157  int err = result->InsertGlobalValues(grid, rowlength, val, &gcid[0]);
1158  if (err != 0 && err != 1) {
1159  std::ostringstream errStr;
1160  errStr << "Epetra_CrsMatrix::InsertGlobalValues returned err=" << err;
1161  throw Exceptions::RuntimeError(errStr.str());
1162  }
1163  }
1164  // free memory
1165  if (bindx) ML_free(bindx);
1166  if (val) ML_free(val);
1167  ML_Operator_Destroy(&ml_AtimesB);
1168  ML_Comm_Destroy(&comm);
1169 
1170  return result;
1171 #else // no MUELU_ML
1172  (void)epA;
1173  (void)epB;
1174  (void)fos;
1175  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1176  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1177  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1178 #endif
1179  }
1180 #endif //ifdef HAVE_XPETRA_EPETRAEXT
1181 
1192  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
1193  const BlockedCrsMatrix& B, bool transposeB,
1194  Teuchos::FancyOStream& fos,
1195  bool doFillComplete = true,
1196  bool doOptimizeStorage = true) {
1197  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1198  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1199 
1200  // Preconditions
1201  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1202  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1203 
1204  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1205  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
1206 
1207  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1208 
1209  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1210  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1211  RCP<Matrix> Cij;
1212 
1213  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1214  RCP<Matrix> crmat1 = A.getMatrix(i,l);
1215  RCP<Matrix> crmat2 = B.getMatrix(l,j);
1216 
1217  if (crmat1.is_null() || crmat2.is_null())
1218  continue;
1219 
1220  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1221  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1222  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
1223  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1224 
1225  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
1226  if (!crop1.is_null())
1227  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
1228  if (!crop2.is_null())
1229  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
1230 
1231  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
1232  "crmat1 does not have global constants");
1233  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
1234  "crmat2 does not have global constants");
1235 
1236  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1237  continue;
1238 
1239 
1240  // temporary matrix containing result of local block multiplication
1241  RCP<Matrix> temp = Teuchos::null;
1242 
1243  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1244  temp = Multiply (*crop1, false, *crop2, false, fos);
1245  else {
1246  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1247  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1248  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1249  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1250  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)");
1251  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)");
1252 
1253  // recursive multiplication call
1254  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1255 
1256  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1257  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)");
1258  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)");
1259  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)");
1260  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)");
1261  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)");
1262  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)");
1263  }
1264 
1265  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1266 
1267  RCP<Matrix> addRes = null;
1268  if (Cij.is_null ())
1269  Cij = temp;
1270  else {
1271  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1272  Cij = addRes;
1273  }
1274  }
1275 
1276  if (!Cij.is_null()) {
1277  if (Cij->isFillComplete())
1278  Cij->resumeFill();
1279  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1280  C->setMatrix(i, j, Cij);
1281  } else {
1282  C->setMatrix(i, j, Teuchos::null);
1283  }
1284  }
1285  }
1286 
1287  if (doFillComplete)
1288  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1289 
1290  return C;
1291  } // TwoMatrixMultiplyBlock
1292 
1305  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1306  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*B.getRowMap())), Exceptions::Incompatible,
1307  "TwoMatrixAdd: matrix row maps are not the same.");
1308 
1309  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1310 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1311  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1312  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1313 
1314  //FIXME is there a bug if beta=0?
1315  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1316 
1317  if (rv != 0)
1318  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1319  std::ostringstream buf;
1320 #else
1321  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1322 #endif
1323  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1324 #ifdef HAVE_XPETRA_TPETRA
1325 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1326  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1327  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1328 # else
1329  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1330  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1331 
1332  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1333 # endif
1334 #else
1335  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1336 #endif
1337  }
1338  } //MatrixMatrix::TwoMatrixAdd()
1339 
1340 
1353  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1354  const Matrix& B, bool transposeB, const SC& beta,
1355  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1356  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
1357  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1358  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1359  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1360  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1361 
1362 
1363  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1364 
1365 
1366  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
1367  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
1368 
1369  auto lib = A.getRowMap()->lib();
1370  if (lib == Xpetra::UseEpetra) {
1371  #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1372  const Epetra_CrsMatrix& epA = helpers::Op2EpetraCrs(A);
1373  const Epetra_CrsMatrix& epB = helpers::Op2EpetraCrs(B);
1374  if(C.is_null())
1375  {
1376  size_t maxNzInA = 0;
1377  size_t maxNzInB = 0;
1378  size_t numLocalRows = 0;
1379  if (A.isFillComplete() && B.isFillComplete()) {
1380 
1381  maxNzInA = A.getNodeMaxNumRowEntries();
1382  maxNzInB = B.getNodeMaxNumRowEntries();
1383  numLocalRows = A.getNodeNumRows();
1384  }
1385 
1386  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1387  // first check if either A or B has at most 1 nonzero per row
1388  // the case of both having at most 1 nz per row is handled by the ``else''
1389  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1390 
1391  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1392  for (size_t i = 0; i < numLocalRows; ++i)
1393  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1394 
1395  } else {
1396  for (size_t i = 0; i < numLocalRows; ++i)
1397  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1398  }
1399 
1400  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1401  << ", using static profiling" << std::endl;
1402  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow));
1403 
1404  } else {
1405  // general case
1406  LO maxPossibleNNZ = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
1407  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), maxPossibleNNZ));
1408  }
1409  if (transposeB)
1410  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1411  }
1412  RCP<Epetra_CrsMatrix> epC = helpers::Op2NonConstEpetraCrs(C);
1413  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1414 
1415  //FIXME is there a bug if beta=0?
1416  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1417 
1418  if (rv != 0)
1419  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1420  #else
1421  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1422  #endif
1423  } else if (lib == Xpetra::UseTpetra) {
1424  #ifdef HAVE_XPETRA_TPETRA
1425  # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1426  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1427  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1428  # else
1429  using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
1430  const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
1431  const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
1432  C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
1433  # endif
1434  #else
1435  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
1436  #endif
1437  }
1438 
1440  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
1441  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
1443  }
1444  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
1445  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1446  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1447  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1448 
1449  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1450  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1451 
1452  size_t i = 0;
1453  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1454  RCP<Matrix> Cij = Teuchos::null;
1455  if(rcpA != Teuchos::null &&
1456  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1457  // recursive call
1458  TwoMatrixAdd(*rcpA, transposeA, alpha,
1459  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1460  Cij, fos, AHasFixedNnzPerRow);
1461  } else if (rcpA == Teuchos::null &&
1462  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1463  Cij = rcpBopB->getMatrix(i,j);
1464  } else if (rcpA != Teuchos::null &&
1465  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1466  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
1467  } else {
1468  Cij = Teuchos::null;
1469  }
1470 
1471  if (!Cij.is_null()) {
1472  if (Cij->isFillComplete())
1473  Cij->resumeFill();
1474  Cij->fillComplete();
1475  bC->setMatrix(i, j, Cij);
1476  } else {
1477  bC->setMatrix(i, j, Teuchos::null);
1478  }
1479  } // loop over columns j
1480  }
1481  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
1482  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
1483  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1484  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
1485 
1486  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1487  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1488 
1489  size_t j = 0;
1490  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1491  RCP<Matrix> Cij = Teuchos::null;
1492  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1493  rcpB!=Teuchos::null) {
1494  // recursive call
1495  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1496  *rcpB, transposeB, beta,
1497  Cij, fos, AHasFixedNnzPerRow);
1498  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1499  rcpB!=Teuchos::null) {
1500  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
1501  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1502  rcpB==Teuchos::null) {
1503  Cij = rcpBopA->getMatrix(i,j);
1504  } else {
1505  Cij = Teuchos::null;
1506  }
1507 
1508  if (!Cij.is_null()) {
1509  if (Cij->isFillComplete())
1510  Cij->resumeFill();
1511  Cij->fillComplete();
1512  bC->setMatrix(i, j, Cij);
1513  } else {
1514  bC->setMatrix(i, j, Teuchos::null);
1515  }
1516  } // loop over rows i
1517  }
1518  else {
1519  // This is the version for blocked matrices
1520 
1521  // check the compatibility of the blocked operators
1522  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1523  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1524  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)");
1525  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)");
1526  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)");
1527  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)");
1528  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)");
1529  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)");
1530 
1531  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1532  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1533 
1534  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1535  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1536 
1537  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1538  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1539 
1540  RCP<Matrix> Cij = Teuchos::null;
1541  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1542  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1543  // recursive call
1544 
1545  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1546  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1547  Cij, fos, AHasFixedNnzPerRow);
1548  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1549  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1550  Cij = rcpBopB->getMatrix(i,j);
1551  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1552  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1553  Cij = rcpBopA->getMatrix(i,j);
1554  } else {
1555  Cij = Teuchos::null;
1556  }
1557 
1558  if (!Cij.is_null()) {
1559  if (Cij->isFillComplete())
1560  Cij->resumeFill();
1561  //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
1562  Cij->fillComplete();
1563  bC->setMatrix(i, j, Cij);
1564  } else {
1565  bC->setMatrix(i, j, Teuchos::null);
1566  }
1567  } // loop over columns j
1568  } // loop over rows i
1569 
1570  } // end blocked recursive algorithm
1571  } //MatrixMatrix::TwoMatrixAdd()
1572  }; // end specialization on SC=double, GO=int and NO=EpetraNode
1573 
1574  // specialization MatrixMatrix for SC=double, GO=long long and NO=EptraNode
1575  template<>
1576  class MatrixMatrix<double,int,long long,EpetraNode> {
1577  typedef double Scalar;
1578  typedef int LocalOrdinal;
1579  typedef long long GlobalOrdinal;
1580  typedef EpetraNode Node;
1581 #include "Xpetra_UseShortNames.hpp"
1582 
1583  public:
1584 
1609  static void Multiply(const Matrix& A, bool transposeA,
1610  const Matrix& B, bool transposeB,
1611  Matrix& C,
1612  bool call_FillComplete_on_result = true,
1613  bool doOptimizeStorage = true,
1614  const std::string & label = std::string(),
1615  const RCP<ParameterList>& params = null) {
1616  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
1617  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
1618  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
1619  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
1620  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
1621 
1622  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
1623  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
1624 
1625  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1626 
1627  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
1628 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1629  helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
1630 #else
1631  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
1632 #endif
1633  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
1634 #ifdef HAVE_XPETRA_TPETRA
1635 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1636  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1637  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long, EpetraNode> ETI enabled."));
1638 # else
1639  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
1640  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
1641  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
1642 
1643  //18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
1644  //Previously, Tpetra's matrix matrix multiply did not support fillComplete.
1645  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1646 # endif
1647 #else
1648  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
1649 #endif
1650  }
1651 
1652  if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1653  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
1654  fillParams->set("Optimize Storage",doOptimizeStorage);
1655  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
1656  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
1657  fillParams);
1658  }
1659 
1660  // transfer striding information
1661  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1662  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1663  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
1664  } // end Multiply
1665 
1688  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1689  const Matrix& B, bool transposeB,
1690  RCP<Matrix> C_in,
1691  Teuchos::FancyOStream& fos,
1692  bool doFillComplete = true,
1693  bool doOptimizeStorage = true,
1694  const std::string & label = std::string(),
1695  const RCP<ParameterList>& params = null) {
1696 
1697  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1698  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1699 
1700  // Default case: Xpetra Multiply
1701  RCP<Matrix> C = C_in;
1702 
1703  if (C == Teuchos::null) {
1704  double nnzPerRow = Teuchos::as<double>(0);
1705 
1706 #if 0
1707  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
1708  // For now, follow what ML and Epetra do.
1709  GO numRowsA = A.getGlobalNumRows();
1710  GO numRowsB = B.getGlobalNumRows();
1711  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
1712  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
1713  nnzPerRow *= nnzPerRow;
1714  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1715  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1716  if (totalNnz < minNnz)
1717  totalNnz = minNnz;
1718  nnzPerRow = totalNnz / A.getGlobalNumRows();
1719 
1720  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1721  }
1722 #endif
1723  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1724  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1725 
1726  } else {
1727  C->resumeFill(); // why this is not done inside of Tpetra MxM?
1728  fos << "Reuse C pattern" << std::endl;
1729  }
1730 
1731  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1732 
1733  return C;
1734  }
1735 
1746  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1747  const Matrix& B, bool transposeB,
1748  Teuchos::FancyOStream &fos,
1749  bool callFillCompleteOnResult = true,
1750  bool doOptimizeStorage = true,
1751  const std::string & label = std::string(),
1752  const RCP<ParameterList>& params = null) {
1753  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1754  }
1755 
1756 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1757  // Michael Gee's MLMultiply
1758  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& /* epA */,
1759  const Epetra_CrsMatrix& /* epB */,
1760  Teuchos::FancyOStream& /* fos */) {
1761  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1762  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1763  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1764  }
1765 #endif //ifdef HAVE_XPETRA_EPETRAEXT
1766 
1777  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
1778  const BlockedCrsMatrix& B, bool transposeB,
1779  Teuchos::FancyOStream& fos,
1780  bool doFillComplete = true,
1781  bool doOptimizeStorage = true) {
1782  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1783  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1784 
1785  // Preconditions
1786  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1787  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1788 
1789  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1790  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
1791 
1792  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1793 
1794  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1795  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1796  RCP<Matrix> Cij;
1797 
1798  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1799  RCP<Matrix> crmat1 = A.getMatrix(i,l);
1800  RCP<Matrix> crmat2 = B.getMatrix(l,j);
1801 
1802  if (crmat1.is_null() || crmat2.is_null())
1803  continue;
1804 
1805  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1806  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1807  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
1808  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1809 
1810  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
1811  if (!crop1.is_null())
1812  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
1813  if (!crop2.is_null())
1814  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
1815 
1816  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
1817  "crmat1 does not have global constants");
1818  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
1819  "crmat2 does not have global constants");
1820 
1821  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1822  continue;
1823 
1824  // temporary matrix containing result of local block multiplication
1825  RCP<Matrix> temp = Teuchos::null;
1826 
1827  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1828  temp = Multiply (*crop1, false, *crop2, false, fos);
1829  else {
1830  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1831  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1832  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1833  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1834  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)");
1835  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)");
1836 
1837  // recursive multiplication call
1838  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1839 
1840  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1841  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)");
1842  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)");
1843  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)");
1844  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)");
1845  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)");
1846  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)");
1847  }
1848 
1849  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1850 
1851  RCP<Matrix> addRes = null;
1852  if (Cij.is_null ())
1853  Cij = temp;
1854  else {
1855  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1856  Cij = addRes;
1857  }
1858  }
1859 
1860  if (!Cij.is_null()) {
1861  if (Cij->isFillComplete())
1862  Cij->resumeFill();
1863  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1864  C->setMatrix(i, j, Cij);
1865  } else {
1866  C->setMatrix(i, j, Teuchos::null);
1867  }
1868  }
1869  }
1870 
1871  if (doFillComplete)
1872  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1873 
1874  return C;
1875  } // TwoMatrixMultiplyBlock
1876 
1889  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1890  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
1891  "TwoMatrixAdd: matrix row maps are not the same.");
1892 
1893  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1894 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1895  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1896  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1897 
1898  //FIXME is there a bug if beta=0?
1899  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1900 
1901  if (rv != 0)
1902  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1903  std::ostringstream buf;
1904 #else
1905  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1906 #endif
1907  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1908 #ifdef HAVE_XPETRA_TPETRA
1909 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1910  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1911  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
1912 # else
1913  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1914  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1915 
1916  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1917 # endif
1918 #else
1919  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1920 #endif
1921  }
1922  } //MatrixMatrix::TwoMatrixAdd()
1923 
1924 
1937  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1938  const Matrix& B, bool transposeB, const SC& beta,
1939  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1940  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1941  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1942  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1943  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1944 
1945  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1946  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
1947  "TwoMatrixAdd: matrix row maps are not the same.");
1948  auto lib = A.getRowMap()->lib();
1949  if (lib == Xpetra::UseEpetra) {
1950  #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1951  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1952  const Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(B);
1953  if(C.is_null())
1954  {
1955  size_t maxNzInA = 0;
1956  size_t maxNzInB = 0;
1957  size_t numLocalRows = 0;
1958  if (A.isFillComplete() && B.isFillComplete()) {
1959 
1960  maxNzInA = A.getNodeMaxNumRowEntries();
1961  maxNzInB = B.getNodeMaxNumRowEntries();
1962  numLocalRows = A.getNodeNumRows();
1963  }
1964 
1965  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1966  // first check if either A or B has at most 1 nonzero per row
1967  // the case of both having at most 1 nz per row is handled by the ``else''
1968  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1969 
1970  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1971  for (size_t i = 0; i < numLocalRows; ++i)
1972  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1973 
1974  } else {
1975  for (size_t i = 0; i < numLocalRows; ++i)
1976  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1977  }
1978 
1979  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1980  << ", using static profiling" << std::endl;
1981  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow));
1982 
1983  } else {
1984  // general case
1985  LO maxPossibleNNZ = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
1986  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), maxPossibleNNZ));
1987  }
1988  if (transposeB)
1989  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1990  }
1991  RCP<Epetra_CrsMatrix> epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
1992  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1993 
1994  //FIXME is there a bug if beta=0?
1995  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1996 
1997  if (rv != 0)
1998  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1999  #else
2000  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
2001  #endif
2002  } else if (lib == Xpetra::UseTpetra) {
2003  #ifdef HAVE_XPETRA_TPETRA
2004  # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
2005  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
2006  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
2007  # else
2008  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
2009  using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
2010  const tcrs_matrix_type& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
2011  const tcrs_matrix_type& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
2012  C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
2013  # endif
2014  #else
2015  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
2016  #endif
2017  }
2018 
2020  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
2021  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
2023  }
2024  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
2025  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
2026  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
2027  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2028 
2029  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2030  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2031 
2032  size_t i = 0;
2033  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2034  RCP<Matrix> Cij = Teuchos::null;
2035  if(rcpA != Teuchos::null &&
2036  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2037  // recursive call
2038  TwoMatrixAdd(*rcpA, transposeA, alpha,
2039  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2040  Cij, fos, AHasFixedNnzPerRow);
2041  } else if (rcpA == Teuchos::null &&
2042  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2043  Cij = rcpBopB->getMatrix(i,j);
2044  } else if (rcpA != Teuchos::null &&
2045  rcpBopB->getMatrix(i,j)==Teuchos::null) {
2046  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
2047  } else {
2048  Cij = Teuchos::null;
2049  }
2050 
2051  if (!Cij.is_null()) {
2052  if (Cij->isFillComplete())
2053  Cij->resumeFill();
2054  Cij->fillComplete();
2055  bC->setMatrix(i, j, Cij);
2056  } else {
2057  bC->setMatrix(i, j, Teuchos::null);
2058  }
2059  } // loop over columns j
2060  }
2061  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
2062  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
2063  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2064  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
2065 
2066  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2067  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2068 
2069  size_t j = 0;
2070  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2071  RCP<Matrix> Cij = Teuchos::null;
2072  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2073  rcpB!=Teuchos::null) {
2074  // recursive call
2075  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2076  *rcpB, transposeB, beta,
2077  Cij, fos, AHasFixedNnzPerRow);
2078  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2079  rcpB!=Teuchos::null) {
2080  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
2081  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2082  rcpB==Teuchos::null) {
2083  Cij = rcpBopA->getMatrix(i,j);
2084  } else {
2085  Cij = Teuchos::null;
2086  }
2087 
2088  if (!Cij.is_null()) {
2089  if (Cij->isFillComplete())
2090  Cij->resumeFill();
2091  Cij->fillComplete();
2092  bC->setMatrix(i, j, Cij);
2093  } else {
2094  bC->setMatrix(i, j, Teuchos::null);
2095  }
2096  } // loop over rows i
2097  }
2098  else {
2099  // This is the version for blocked matrices
2100 
2101  // check the compatibility of the blocked operators
2102  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2103  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2104  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)");
2105  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)");
2106  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)");
2107  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)");
2108  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)");
2109  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)");
2110 
2111  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2112  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2113 
2114  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2115  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2116 
2117  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2118  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2119 
2120  RCP<Matrix> Cij = Teuchos::null;
2121  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2122  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2123  // recursive call
2124 
2125  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2126  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2127  Cij, fos, AHasFixedNnzPerRow);
2128  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2129  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2130  Cij = rcpBopB->getMatrix(i,j);
2131  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2132  rcpBopB->getMatrix(i,j)==Teuchos::null) {
2133  Cij = rcpBopA->getMatrix(i,j);
2134  } else {
2135  Cij = Teuchos::null;
2136  }
2137 
2138  if (!Cij.is_null()) {
2139  if (Cij->isFillComplete())
2140  Cij->resumeFill();
2141  //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
2142  Cij->fillComplete();
2143  bC->setMatrix(i, j, Cij);
2144  } else {
2145  bC->setMatrix(i, j, Teuchos::null);
2146  }
2147  } // loop over columns j
2148  } // loop over rows i
2149 
2150  } // end blocked recursive algorithm
2151  } //MatrixMatrix::TwoMatrixAdd()
2152  }; // end specialization on GO=long long and NO=EpetraNode
2153 
2154 #endif // HAVE_XPETRA_EPETRA
2155 
2156 } // end namespace Xpetra
2157 
2158 #define XPETRA_MATRIXMATRIX_SHORT
2159 
2160 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_ */
Xpetra::Helpers::tcrs_matrix_type
Tpetra::CrsMatrix< SC, LO, GO, NO > tcrs_matrix_type
Definition: Xpetra_MatrixMatrix.hpp:218
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition: Xpetra_CrsMatrixWrap_def.hpp:457
Xpetra::Matrix::getNodeMaxNumRowEntries
virtual size_t getNodeMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
Xpetra::Matrix::getNodeNumRows
virtual size_t getNodeNumRows() const =0
Returns the number of matrix rows owned on the calling node.
Xpetra::toString
std::string toString(Xpetra::UnderlyingLib lib)
Convert a Xpetra::UnderlyingLib to a std::string.
Definition: Xpetra_Utils.cpp:51
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:668
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:1353
Xpetra::Matrix::getGlobalNumRows
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
Xpetra::Helpers::tpetraAdd
static Teuchos::RCP< Matrix > tpetraAdd(const tcrs_matrix_type &A, bool transposeA, const typename tcrs_matrix_type::scalar_type alpha, const tcrs_matrix_type &B, bool transposeB, const typename tcrs_matrix_type::scalar_type beta)
Definition: Xpetra_MatrixMatrix.hpp:219
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:1688
Xpetra::BlockedCrsMatrix::isFillComplete
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Definition: Xpetra_BlockedCrsMatrix.hpp:711
GO
GlobalOrdinal GO
Definition: Xpetra_UseShortNamesOrdinal.hpp:139
Xpetra::BlockedCrsMatrix::getRangeMap
RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1073
Xpetra::Helpers::Op2NonConstEpetraCrs
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:145
Xpetra_UseShortNames.hpp
Xpetra::Helpers
Xpetra utility class containing transformation routines between Xpetra::Matrix and Epetra/Tpetra obje...
Definition: Xpetra_MatrixMatrix.hpp:95
Xpetra_EpetraCrsMatrix_fwd.hpp
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory.hpp:234
Xpetra::Matrix::resumeFill
virtual void resumeFill(const RCP< ParameterList > &params=null)=0
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 &quot;in-place&quot;.
Definition: Xpetra_MatrixMatrix.hpp:1777
Xpetra_Matrix.hpp
Xpetra::Helpers::Op2EpetraCrs
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:129
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:1937
Xpetra::MatrixMatrix< double, int, long long, EpetraNode >::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &, const Epetra_CrsMatrix &, Teuchos::FancyOStream &)
Definition: Xpetra_MatrixMatrix.hpp:1758
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:499
Xpetra::Matrix::IsView
bool IsView(const viewLabel_t viewLabel) const
Definition: Xpetra_Matrix.hpp:207
Xpetra::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition: Xpetra_BlockedCrsMatrix.hpp:1332
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1299
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:869
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:949
Xpetra::Helpers::Op2NonConstTpetraCrs
static Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:201
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:443
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:1609
Xpetra::MatrixMatrix::TwoMatrixMultiplyBlock
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply &quot;in-place&quot;.
Definition: Xpetra_MatrixMatrix.hpp:525
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:637
Xpetra::Helpers::epetraExtMult
static void epetraExtMult(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool fillCompleteResult)
Definition: Xpetra_MatrixMatrix.hpp:283
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
Xpetra::Helpers::Op2EpetraCrs
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:101
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1085
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:368
Xpetra::Helpers::Op2NonConstEpetraCrs
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:115
Xpetra_ConfigDefs.hpp
Xpetra::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: Xpetra_Exceptions.hpp:108
Xpetra::BlockedCrsMatrix::getDomainMap
RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1058
SC
Scalar SC
Definition: Xpetra_UseShortNamesScalar.hpp:178
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:1305
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::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::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1296
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:1044
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 &quot;in-place&quot;.
Definition: Xpetra_MatrixMatrix.hpp:1192
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::Helpers::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:175
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:1032
Xpetra::MatrixMatrix::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
Definition: Xpetra_MatrixMatrix.hpp:507
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:1889
Xpetra::Helpers::Op2TpetraCrs
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:186
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:163
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:1746
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1082