Xpetra  Version of the Day
Xpetra_MatrixMatrix.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 
48 #ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_
49 #define PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_
50 
51 #include "Xpetra_ConfigDefs.hpp"
52 
53 #include "Xpetra_BlockedCrsMatrix.hpp"
54 #include "Xpetra_CrsMatrixWrap.hpp"
55 #include "Xpetra_MapExtractor.hpp"
56 #include "Xpetra_Map.hpp"
57 #include "Xpetra_MatrixFactory.hpp"
58 #include "Xpetra_Matrix.hpp"
59 #include "Xpetra_StridedMapFactory.hpp"
60 #include "Xpetra_StridedMap.hpp"
61 
62 #ifdef HAVE_XPETRA_EPETRA
63 #include <Xpetra_EpetraCrsMatrix_fwd.hpp>
64 #endif
65 
66 #ifdef HAVE_XPETRA_EPETRAEXT
67 #include <EpetraExt_MatrixMatrix.h>
68 #include <EpetraExt_RowMatrixOut.h>
69 #include <Epetra_RowMatrixTransposer.h>
70 #endif // HAVE_XPETRA_EPETRAEXT
71 
72 #ifdef HAVE_XPETRA_TPETRA
73 #include <TpetraExt_MatrixMatrix.hpp>
74 #include <Tpetra_RowMatrixTransposer.hpp>
75 #include <MatrixMarket_Tpetra.hpp>
76 #include <Xpetra_TpetraCrsMatrix.hpp>
77 #include <Xpetra_TpetraMultiVector.hpp>
78 #include <Xpetra_TpetraVector.hpp>
79 #endif // HAVE_XPETRA_TPETRA
80 
81 namespace Xpetra {
82 
89  template <class Scalar,
90  class LocalOrdinal = int,
91  class GlobalOrdinal = LocalOrdinal,
92  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
93  class Helpers {
94 #include "Xpetra_UseShortNames.hpp"
95 
96  public:
97 
98 #ifdef HAVE_XPETRA_EPETRA
99  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<Matrix> Op) {
100  // Get the underlying Epetra Mtx
101  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
102  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast,
103  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
104 
105  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
106  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
107  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast,
108  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
109 
110  return tmp_ECrsMtx->getEpetra_CrsMatrix();
111  }
112 
113  static RCP<Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) {
114  RCP<Epetra_CrsMatrix> A;
115  // Get the underlying Epetra Mtx
116  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
117  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast,
118  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
119 
120  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
121  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
122  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
123 
124  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
125  }
126 
127  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) {
128  // Get the underlying Epetra Mtx
129  try {
130  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
131  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
132  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
133  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast,
134  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
135 
136  return *tmp_ECrsMtx->getEpetra_CrsMatrix();
137 
138  } catch(...) {
139  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
140  }
141  }
142 
143  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(const Matrix& Op) {
144  RCP<Epetra_CrsMatrix> A;
145  // Get the underlying Epetra Mtx
146  try {
147  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
148  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
149  const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
150  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
151 
152  return *Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
153 
154  } catch(...) {
155  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
156  }
157  }
158 #endif // HAVE_XPETRA_EPETRA
159 
160 #ifdef HAVE_XPETRA_TPETRA
161  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<Matrix> Op) {
162  // Get the underlying Tpetra Mtx
163  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
164  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
165 
166  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
167  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
168  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
169 
170  return tmp_ECrsMtx->getTpetra_CrsMatrix();
171  }
172 
173  static RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) {
174  // Get the underlying Tpetra Mtx
175  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
176  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
177  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
178  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
179  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
180 
181  return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
182  }
183 
184  static const Tpetra::CrsMatrix<SC,LO,GO,NO> & Op2TpetraCrs(const Matrix& Op) {
185  // Get the underlying Tpetra Mtx
186  try{
187  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
188  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
189  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
190  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
191 
192  return *tmp_TCrsMtx->getTpetra_CrsMatrix();
193 
194  } catch (...) {
195  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
196  }
197  }
198 
199  static Tpetra::CrsMatrix<SC,LO,GO,NO> & Op2NonConstTpetraCrs(const Matrix& Op) {
200  // Get the underlying Tpetra Mtx
201  try{
202  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap& >(Op);
203  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
204  const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
205  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
206 
207  return *Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_TCrsMtx->getTpetra_CrsMatrix());
208 
209  } catch (...) {
210  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
211  }
212  }
213 #endif // HAVE_XPETRA_TPETRA
214 
215 #ifdef HAVE_XPETRA_TPETRA
216  using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
217  static Teuchos::RCP<Matrix> tpetraAdd(
218  const tcrs_matrix_type& A, bool transposeA, const typename tcrs_matrix_type::scalar_type alpha,
219  const tcrs_matrix_type& B, bool transposeB, const typename tcrs_matrix_type::scalar_type beta)
220  {
221  using Teuchos::Array;
222  using Teuchos::RCP;
223  using Teuchos::rcp;
224  using Teuchos::rcp_implicit_cast;
225  using Teuchos::rcpFromRef;
226  using Teuchos::ParameterList;
227  using XTCrsType = Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>;
228  using CrsType = Xpetra::CrsMatrix<SC,LO,GO,NO>;
229  using CrsWrap = Xpetra::CrsMatrixWrap<SC,LO,GO,NO>;
230  using transposer_type = Tpetra::RowMatrixTransposer<SC,LO,GO,NO>;
231  using import_type = Tpetra::Import<LO,GO,NO>;
232  RCP<const tcrs_matrix_type> Aprime = rcpFromRef(A);
233  if(transposeA)
234  Aprime = transposer_type(Aprime).createTranspose();
235  //Decide whether the fast code path can be taken.
236  if(A.isFillComplete() && B.isFillComplete())
237  {
238  RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), 0));
239  RCP<ParameterList> addParams = rcp(new ParameterList);
240  addParams->set("Call fillComplete", false);
241  //passing A after B means C will have the same domain/range map as A (or A^T if transposeA)
242  Tpetra::MatrixMatrix::add<SC,LO,GO,NO>(beta, transposeB, B, alpha, false, *Aprime, *C, Teuchos::null, Teuchos::null, addParams);
243  return rcp_implicit_cast<Matrix>(rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C))))));
244  }
245  else
246  {
247  //Slow case - one or both operands are non-fill complete.
248  //TODO: deprecate this.
249  //Need to compute the explicit transpose before add if transposeA and/or transposeB.
250  //This is to count the number of entries to allocate per row in the final sum.
251  RCP<const tcrs_matrix_type> Bprime = rcpFromRef(B);
252  if(transposeB)
253  Bprime = transposer_type(Bprime).createTranspose();
254  //Use Aprime's row map as C's row map.
255  if(!(Aprime->getRowMap()->isSameAs(*(Bprime->getRowMap()))))
256  {
257  auto import = rcp(new import_type(Bprime->getRowMap(), Aprime->getRowMap()));
258  Bprime = Tpetra::importAndFillCompleteCrsMatrix<tcrs_matrix_type>(Bprime, *import, Aprime->getDomainMap(), Aprime->getRangeMap());
259  }
260  //Count the entries to allocate for C in each local row.
261  LO numLocalRows = Aprime->getNodeNumRows();
262  Array<size_t> allocPerRow(numLocalRows);
263  //0 is always the minimum LID
264  for(LO i = 0; i < numLocalRows; i++)
265  {
266  allocPerRow[i] = Aprime->getNumEntriesInLocalRow(i) + Bprime->getNumEntriesInLocalRow(i);
267  }
268  //Construct C
269  RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), allocPerRow(), Tpetra::StaticProfile));
270  //Compute the sum in C (already took care of transposes)
271  Tpetra::MatrixMatrix::Add<SC,LO,GO,NO>(
272  *Aprime, false, alpha,
273  *Bprime, false, beta,
274  C);
275  return rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C)))));
276  }
277  }
278 #endif
279 
280 #ifdef HAVE_XPETRA_EPETRAEXT
281  static void epetraExtMult(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Matrix& C, bool fillCompleteResult)
282  {
283  Epetra_CrsMatrix& epA = Op2NonConstEpetraCrs(A);
284  Epetra_CrsMatrix& epB = Op2NonConstEpetraCrs(B);
285  Epetra_CrsMatrix& epC = Op2NonConstEpetraCrs(C);
286  //Want a static profile (possibly fill complete) matrix as the result.
287  //But, EpetraExt Multiply needs C to be dynamic profile, so compute the product in a temporary DynamicProfile matrix.
288  const Epetra_Map& Crowmap = epC.RowMap();
289  int errCode = 0;
290  Epetra_CrsMatrix Ctemp(::Copy, Crowmap, 0, false);
291  if(fillCompleteResult) {
292  errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, true);
293  if(!errCode) {
294  epC = Ctemp;
295  C.fillComplete();
296  }
297  }
298  else {
299  errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, false);
300  if(!errCode) {
301  int numLocalRows = Crowmap.NumMyElements();
302  long long* globalElementList = nullptr;
303  Crowmap.MyGlobalElementsPtr(globalElementList);
304  Teuchos::Array<int> entriesPerRow(numLocalRows, 0);
305  for(int i = 0; i < numLocalRows; i++)
306  {
307  entriesPerRow[i] = Ctemp.NumGlobalEntries(globalElementList[i]);
308  }
309  //know how many entries to allocate in epC (which must be static profile)
310  epC = Epetra_CrsMatrix(::Copy, Crowmap, entriesPerRow.data(), true);
311  for(int i = 0; i < numLocalRows; i++)
312  {
313  int gid = globalElementList[i];
314  int numEntries;
315  double* values;
316  int* inds;
317  Ctemp.ExtractGlobalRowView(gid, numEntries, values, inds);
318  epC.InsertGlobalValues(gid, numEntries, values, inds);
319  }
320  }
321  }
322  if(errCode) {
323  std::ostringstream buf;
324  buf << errCode;
325  std::string msg = "EpetraExt::MatrixMatrix::Multiply returned nonzero error code " + buf.str();
326  throw(Exceptions::RuntimeError(msg));
327  }
328  }
329 #endif
330  };
331 
332  template <class Scalar,
333  class LocalOrdinal /*= int*/,
334  class GlobalOrdinal /*= LocalOrdinal*/,
335  class Node /*= KokkosClassic::DefaultNode::DefaultNodeType*/>
336  class MatrixMatrix {
337 #undef XPETRA_MATRIXMATRIX_SHORT
338 #include "Xpetra_UseShortNames.hpp"
339 
340  public:
341 
366  static void Multiply(const Matrix& A, bool transposeA,
367  const Matrix& B, bool transposeB,
368  Matrix& C,
369  bool call_FillComplete_on_result = true,
370  bool doOptimizeStorage = true,
371  const std::string & label = std::string(),
372  const RCP<ParameterList>& params = null) {
373 
374  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
375  Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
376  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
377  Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
378 
379  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
380  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
381 
382  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
383 
384  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
385 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
386  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
387 #else
388  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
389 
390 #endif
391  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
392 #ifdef HAVE_XPETRA_TPETRA
393  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
394  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
395  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
396 
397  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
398  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
399  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
400 #else
401  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
402 #endif
403  }
404 
405  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
406  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
407  fillParams->set("Optimize Storage", doOptimizeStorage);
408  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
409  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
410  fillParams);
411  }
412 
413  // transfer striding information
414  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
415  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
416  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
417  } // end Multiply
418 
441  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, RCP<Matrix> C_in,
442  Teuchos::FancyOStream& fos,
443  bool doFillComplete = true,
444  bool doOptimizeStorage = true,
445  const std::string & label = std::string(),
446  const RCP<ParameterList>& params = null) {
447 
448  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
449  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
450 
451  // Default case: Xpetra Multiply
452  RCP<Matrix> C = C_in;
453 
454  if (C == Teuchos::null) {
455  double nnzPerRow = Teuchos::as<double>(0);
456 
457 #if 0
458  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
459  // For now, follow what ML and Epetra do.
460  GO numRowsA = A.getGlobalNumRows();
461  GO numRowsB = B.getGlobalNumRows();
462  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
463  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
464  nnzPerRow *= nnzPerRow;
465  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
466  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
467  if (totalNnz < minNnz)
468  totalNnz = minNnz;
469  nnzPerRow = totalNnz / A.getGlobalNumRows();
470 
471  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
472  }
473 #endif
474  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
475  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
476 
477  } else {
478  C->resumeFill(); // why this is not done inside of Tpetra MxM?
479  fos << "Reuse C pattern" << std::endl;
480  }
481 
482  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
483 
484  return C;
485  }
486 
497  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Teuchos::FancyOStream &fos,
498  bool callFillCompleteOnResult = true, bool doOptimizeStorage = true, const std::string& label = std::string(),
499  const RCP<ParameterList>& params = null) {
500  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
501  }
502 
503 #ifdef HAVE_XPETRA_EPETRAEXT
504  // Michael Gee's MLMultiply
505  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
506  const Epetra_CrsMatrix& epB,
507  Teuchos::FancyOStream& fos) {
508  throw(Xpetra::Exceptions::RuntimeError("MLTwoMatrixMultiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
509  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
510  }
511 #endif //ifdef HAVE_XPETRA_EPETRAEXT
512 
523  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
524  const BlockedCrsMatrix& B, bool transposeB,
525  Teuchos::FancyOStream& fos,
526  bool doFillComplete = true,
527  bool doOptimizeStorage = true) {
528  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
529  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
530 
531  // Preconditions
532  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
533  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
534 
535  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
536  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
537 
538  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
539 
540  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
541  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
542  RCP<Matrix> Cij;
543 
544  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
545  RCP<Matrix> crmat1 = A.getMatrix(i,l);
546  RCP<Matrix> crmat2 = B.getMatrix(l,j);
547 
548  if (crmat1.is_null() || crmat2.is_null())
549  continue;
550 
551  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
552  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
553  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
554  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
555 
556  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
557  if (!crop1.is_null())
558  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
559  if (!crop2.is_null())
560  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
561 
562  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
563  "crmat1 does not have global constants");
564  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
565  "crmat2 does not have global constants");
566 
567  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
568  continue;
569 
570  // temporary matrix containing result of local block multiplication
571  RCP<Matrix> temp = Teuchos::null;
572 
573  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
574  temp = Multiply (*crop1, false, *crop2, false, fos);
575  else {
576  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
577  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
578  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
579  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
580  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)");
581  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)");
582 
583  // recursive multiplication call
584  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
585 
586  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
587  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)");
588  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)");
589  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)");
590  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)");
591  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)");
592  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)");
593  }
594 
595  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
596 
597  RCP<Matrix> addRes = null;
598  if (Cij.is_null ())
599  Cij = temp;
600  else {
601  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
602  Cij = addRes;
603  }
604  }
605 
606  if (!Cij.is_null()) {
607  if (Cij->isFillComplete())
608  Cij->resumeFill();
609  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
610  C->setMatrix(i, j, Cij);
611  } else {
612  C->setMatrix(i, j, Teuchos::null);
613  }
614  }
615  }
616 
617  if (doFillComplete)
618  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
619 
620  return C;
621  } // TwoMatrixMultiplyBlock
622 
635  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
636  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
637  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
638 
639  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
640  throw Exceptions::RuntimeError("TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
641  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
642 #ifdef HAVE_XPETRA_TPETRA
643  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
644  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
645 
646  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
647 #else
648  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
649 #endif
650  }
651  } //MatrixMatrix::TwoMatrixAdd()
652 
653 
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 
1355  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1356  const Matrix& B, bool transposeB, const SC& beta,
1357  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1358  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
1359  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1360  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1361  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1362  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1363 
1364 
1365  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1366 
1367 
1368  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
1369  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
1370 
1371  auto lib = A.getRowMap()->lib();
1372  if (lib == Xpetra::UseEpetra) {
1373  #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1374  const Epetra_CrsMatrix& epA = helpers::Op2EpetraCrs(A);
1375  const Epetra_CrsMatrix& epB = helpers::Op2EpetraCrs(B);
1376  if(C.is_null())
1377  {
1378  size_t maxNzInA = 0;
1379  size_t maxNzInB = 0;
1380  size_t numLocalRows = 0;
1381  if (A.isFillComplete() && B.isFillComplete()) {
1382 
1383  maxNzInA = A.getNodeMaxNumRowEntries();
1384  maxNzInB = B.getNodeMaxNumRowEntries();
1385  numLocalRows = A.getNodeNumRows();
1386  }
1387 
1388  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1389  // first check if either A or B has at most 1 nonzero per row
1390  // the case of both having at most 1 nz per row is handled by the ``else''
1391  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1392 
1393  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1394  for (size_t i = 0; i < numLocalRows; ++i)
1395  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1396 
1397  } else {
1398  for (size_t i = 0; i < numLocalRows; ++i)
1399  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1400  }
1401 
1402  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1403  << ", using static profiling" << std::endl;
1404  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow));
1405 
1406  } else {
1407  // general case
1408  LO maxPossibleNNZ = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
1409  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), maxPossibleNNZ));
1410  }
1411  if (transposeB)
1412  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1413  }
1414  RCP<Epetra_CrsMatrix> epC = helpers::Op2NonConstEpetraCrs(C);
1415  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1416 
1417  //FIXME is there a bug if beta=0?
1418  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1419 
1420  if (rv != 0)
1421  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1422  #else
1423  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1424  #endif
1425  } else if (lib == Xpetra::UseTpetra) {
1426  #ifdef HAVE_XPETRA_TPETRA
1427  # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1428  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1429  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1430  # else
1431  using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
1432  const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
1433  const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
1434  C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
1435  # endif
1436  #else
1437  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
1438  #endif
1439  }
1440 
1442  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
1443  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
1445  }
1446  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
1447  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1448  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1449  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1450 
1451  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1452  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1453 
1454  size_t i = 0;
1455  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1456  RCP<Matrix> Cij = Teuchos::null;
1457  if(rcpA != Teuchos::null &&
1458  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1459  // recursive call
1460  TwoMatrixAdd(*rcpA, transposeA, alpha,
1461  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1462  Cij, fos, AHasFixedNnzPerRow);
1463  } else if (rcpA == Teuchos::null &&
1464  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1465  Cij = rcpBopB->getMatrix(i,j);
1466  } else if (rcpA != Teuchos::null &&
1467  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1468  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
1469  } else {
1470  Cij = Teuchos::null;
1471  }
1472 
1473  if (!Cij.is_null()) {
1474  if (Cij->isFillComplete())
1475  Cij->resumeFill();
1476  Cij->fillComplete();
1477  bC->setMatrix(i, j, Cij);
1478  } else {
1479  bC->setMatrix(i, j, Teuchos::null);
1480  }
1481  } // loop over columns j
1482  }
1483  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
1484  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
1485  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1486  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
1487 
1488  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1489  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1490 
1491  size_t j = 0;
1492  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1493  RCP<Matrix> Cij = Teuchos::null;
1494  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1495  rcpB!=Teuchos::null) {
1496  // recursive call
1497  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1498  *rcpB, transposeB, beta,
1499  Cij, fos, AHasFixedNnzPerRow);
1500  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1501  rcpB!=Teuchos::null) {
1502  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
1503  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1504  rcpB==Teuchos::null) {
1505  Cij = rcpBopA->getMatrix(i,j);
1506  } else {
1507  Cij = Teuchos::null;
1508  }
1509 
1510  if (!Cij.is_null()) {
1511  if (Cij->isFillComplete())
1512  Cij->resumeFill();
1513  Cij->fillComplete();
1514  bC->setMatrix(i, j, Cij);
1515  } else {
1516  bC->setMatrix(i, j, Teuchos::null);
1517  }
1518  } // loop over rows i
1519  }
1520  else {
1521  // This is the version for blocked matrices
1522 
1523  // check the compatibility of the blocked operators
1524  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1525  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1526  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)");
1527  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)");
1528  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)");
1529  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)");
1530  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)");
1531  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)");
1532 
1533  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1534  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1535 
1536  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1537  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1538 
1539  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1540  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1541 
1542  RCP<Matrix> Cij = Teuchos::null;
1543  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1544  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1545  // recursive call
1546 
1547  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1548  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1549  Cij, fos, AHasFixedNnzPerRow);
1550  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1551  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1552  Cij = rcpBopB->getMatrix(i,j);
1553  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1554  rcpBopB->getMatrix(i,j)==Teuchos::null) {
1555  Cij = rcpBopA->getMatrix(i,j);
1556  } else {
1557  Cij = Teuchos::null;
1558  }
1559 
1560  if (!Cij.is_null()) {
1561  if (Cij->isFillComplete())
1562  Cij->resumeFill();
1563  //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
1564  Cij->fillComplete();
1565  bC->setMatrix(i, j, Cij);
1566  } else {
1567  bC->setMatrix(i, j, Teuchos::null);
1568  }
1569  } // loop over columns j
1570  } // loop over rows i
1571 
1572  } // end blocked recursive algorithm
1573  } //MatrixMatrix::TwoMatrixAdd()
1574  }; // end specialization on SC=double, GO=int and NO=EpetraNode
1575 
1576  // specialization MatrixMatrix for SC=double, GO=long long and NO=EptraNode
1577  template<>
1578  class MatrixMatrix<double,int,long long,EpetraNode> {
1579  typedef double Scalar;
1580  typedef int LocalOrdinal;
1581  typedef long long GlobalOrdinal;
1582  typedef EpetraNode Node;
1583 #include "Xpetra_UseShortNames.hpp"
1584 
1585  public:
1586 
1611  static void Multiply(const Matrix& A, bool transposeA,
1612  const Matrix& B, bool transposeB,
1613  Matrix& C,
1614  bool call_FillComplete_on_result = true,
1615  bool doOptimizeStorage = true,
1616  const std::string & label = std::string(),
1617  const RCP<ParameterList>& params = null) {
1618  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
1619  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
1620  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
1621  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
1622  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
1623 
1624  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
1625  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
1626 
1627  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1628 
1629  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
1630 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1631  helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
1632 #else
1633  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
1634 #endif
1635  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
1636 #ifdef HAVE_XPETRA_TPETRA
1637 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1638  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1639  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long, EpetraNode> ETI enabled."));
1640 # else
1641  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
1642  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
1643  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
1644 
1645  //18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
1646  //Previously, Tpetra's matrix matrix multiply did not support fillComplete.
1647  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1648 # endif
1649 #else
1650  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
1651 #endif
1652  }
1653 
1654  if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1655  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
1656  fillParams->set("Optimize Storage",doOptimizeStorage);
1657  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
1658  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
1659  fillParams);
1660  }
1661 
1662  // transfer striding information
1663  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1664  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1665  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
1666  } // end Multiply
1667 
1690  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1691  const Matrix& B, bool transposeB,
1692  RCP<Matrix> C_in,
1693  Teuchos::FancyOStream& fos,
1694  bool doFillComplete = true,
1695  bool doOptimizeStorage = true,
1696  const std::string & label = std::string(),
1697  const RCP<ParameterList>& params = null) {
1698 
1699  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1700  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1701 
1702  // Default case: Xpetra Multiply
1703  RCP<Matrix> C = C_in;
1704 
1705  if (C == Teuchos::null) {
1706  double nnzPerRow = Teuchos::as<double>(0);
1707 
1708 #if 0
1709  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
1710  // For now, follow what ML and Epetra do.
1711  GO numRowsA = A.getGlobalNumRows();
1712  GO numRowsB = B.getGlobalNumRows();
1713  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
1714  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
1715  nnzPerRow *= nnzPerRow;
1716  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1717  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1718  if (totalNnz < minNnz)
1719  totalNnz = minNnz;
1720  nnzPerRow = totalNnz / A.getGlobalNumRows();
1721 
1722  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1723  }
1724 #endif
1725  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1726  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1727 
1728  } else {
1729  C->resumeFill(); // why this is not done inside of Tpetra MxM?
1730  fos << "Reuse C pattern" << std::endl;
1731  }
1732 
1733  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1734 
1735  return C;
1736  }
1737 
1748  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1749  const Matrix& B, bool transposeB,
1750  Teuchos::FancyOStream &fos,
1751  bool callFillCompleteOnResult = true,
1752  bool doOptimizeStorage = true,
1753  const std::string & label = std::string(),
1754  const RCP<ParameterList>& params = null) {
1755  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1756  }
1757 
1758 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1759  // Michael Gee's MLMultiply
1760  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& /* epA */,
1761  const Epetra_CrsMatrix& /* epB */,
1762  Teuchos::FancyOStream& /* fos */) {
1763  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1764  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1765  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1766  }
1767 #endif //ifdef HAVE_XPETRA_EPETRAEXT
1768 
1779  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
1780  const BlockedCrsMatrix& B, bool transposeB,
1781  Teuchos::FancyOStream& fos,
1782  bool doFillComplete = true,
1783  bool doOptimizeStorage = true) {
1784  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1785  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1786 
1787  // Preconditions
1788  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1789  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1790 
1791  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1792  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
1793 
1794  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1795 
1796  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1797  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1798  RCP<Matrix> Cij;
1799 
1800  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1801  RCP<Matrix> crmat1 = A.getMatrix(i,l);
1802  RCP<Matrix> crmat2 = B.getMatrix(l,j);
1803 
1804  if (crmat1.is_null() || crmat2.is_null())
1805  continue;
1806 
1807  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1808  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1809  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
1810  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1811 
1812  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
1813  if (!crop1.is_null())
1814  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
1815  if (!crop2.is_null())
1816  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
1817 
1818  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
1819  "crmat1 does not have global constants");
1820  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
1821  "crmat2 does not have global constants");
1822 
1823  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1824  continue;
1825 
1826  // temporary matrix containing result of local block multiplication
1827  RCP<Matrix> temp = Teuchos::null;
1828 
1829  if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1830  temp = Multiply (*crop1, false, *crop2, false, fos);
1831  else {
1832  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1833  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1834  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1835  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1836  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)");
1837  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)");
1838 
1839  // recursive multiplication call
1840  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1841 
1842  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1843  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)");
1844  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)");
1845  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)");
1846  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)");
1847  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)");
1848  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)");
1849  }
1850 
1851  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1852 
1853  RCP<Matrix> addRes = null;
1854  if (Cij.is_null ())
1855  Cij = temp;
1856  else {
1857  MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1858  Cij = addRes;
1859  }
1860  }
1861 
1862  if (!Cij.is_null()) {
1863  if (Cij->isFillComplete())
1864  Cij->resumeFill();
1865  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1866  C->setMatrix(i, j, Cij);
1867  } else {
1868  C->setMatrix(i, j, Teuchos::null);
1869  }
1870  }
1871  }
1872 
1873  if (doFillComplete)
1874  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1875 
1876  return C;
1877  } // TwoMatrixMultiplyBlock
1878 
1891  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1892  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
1893  "TwoMatrixAdd: matrix row maps are not the same.");
1894 
1895  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1896 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1897  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1898  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1899 
1900  //FIXME is there a bug if beta=0?
1901  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1902 
1903  if (rv != 0)
1904  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1905  std::ostringstream buf;
1906 #else
1907  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1908 #endif
1909  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1910 #ifdef HAVE_XPETRA_TPETRA
1911 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1912  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1913  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
1914 # else
1915  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1916  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1917 
1918  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1919 # endif
1920 #else
1921  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1922 #endif
1923  }
1924  } //MatrixMatrix::TwoMatrixAdd()
1925 
1926 
1941  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1942  const Matrix& B, bool transposeB, const SC& beta,
1943  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1944  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1945  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1946  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1947  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1948 
1949  if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1950  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
1951  "TwoMatrixAdd: matrix row maps are not the same.");
1952  auto lib = A.getRowMap()->lib();
1953  if (lib == Xpetra::UseEpetra) {
1954  #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1955  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1956  const Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(B);
1957  if(C.is_null())
1958  {
1959  size_t maxNzInA = 0;
1960  size_t maxNzInB = 0;
1961  size_t numLocalRows = 0;
1962  if (A.isFillComplete() && B.isFillComplete()) {
1963 
1964  maxNzInA = A.getNodeMaxNumRowEntries();
1965  maxNzInB = B.getNodeMaxNumRowEntries();
1966  numLocalRows = A.getNodeNumRows();
1967  }
1968 
1969  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1970  // first check if either A or B has at most 1 nonzero per row
1971  // the case of both having at most 1 nz per row is handled by the ``else''
1972  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1973 
1974  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1975  for (size_t i = 0; i < numLocalRows; ++i)
1976  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1977 
1978  } else {
1979  for (size_t i = 0; i < numLocalRows; ++i)
1980  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1981  }
1982 
1983  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1984  << ", using static profiling" << std::endl;
1985  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow));
1986 
1987  } else {
1988  // general case
1989  LO maxPossibleNNZ = A.getNodeMaxNumRowEntries() + B.getNodeMaxNumRowEntries();
1990  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), maxPossibleNNZ));
1991  }
1992  if (transposeB)
1993  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1994  }
1995  RCP<Epetra_CrsMatrix> epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
1996  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1997 
1998  //FIXME is there a bug if beta=0?
1999  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
2000 
2001  if (rv != 0)
2002  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
2003  #else
2004  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
2005  #endif
2006  } else if (lib == Xpetra::UseTpetra) {
2007  #ifdef HAVE_XPETRA_TPETRA
2008  # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
2009  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
2010  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
2011  # else
2012  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
2013  using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
2014  const tcrs_matrix_type& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
2015  const tcrs_matrix_type& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
2016  C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
2017  # endif
2018  #else
2019  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
2020  #endif
2021  }
2022 
2024  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
2025  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
2027  }
2028  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
2029  else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
2030  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
2031  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2032 
2033  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2034  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2035 
2036  size_t i = 0;
2037  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2038  RCP<Matrix> Cij = Teuchos::null;
2039  if(rcpA != Teuchos::null &&
2040  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2041  // recursive call
2042  TwoMatrixAdd(*rcpA, transposeA, alpha,
2043  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2044  Cij, fos, AHasFixedNnzPerRow);
2045  } else if (rcpA == Teuchos::null &&
2046  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2047  Cij = rcpBopB->getMatrix(i,j);
2048  } else if (rcpA != Teuchos::null &&
2049  rcpBopB->getMatrix(i,j)==Teuchos::null) {
2050  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
2051  } else {
2052  Cij = Teuchos::null;
2053  }
2054 
2055  if (!Cij.is_null()) {
2056  if (Cij->isFillComplete())
2057  Cij->resumeFill();
2058  Cij->fillComplete();
2059  bC->setMatrix(i, j, Cij);
2060  } else {
2061  bC->setMatrix(i, j, Teuchos::null);
2062  }
2063  } // loop over columns j
2064  }
2065  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
2066  else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
2067  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2068  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
2069 
2070  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2071  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2072 
2073  size_t j = 0;
2074  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2075  RCP<Matrix> Cij = Teuchos::null;
2076  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2077  rcpB!=Teuchos::null) {
2078  // recursive call
2079  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2080  *rcpB, transposeB, beta,
2081  Cij, fos, AHasFixedNnzPerRow);
2082  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2083  rcpB!=Teuchos::null) {
2084  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
2085  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2086  rcpB==Teuchos::null) {
2087  Cij = rcpBopA->getMatrix(i,j);
2088  } else {
2089  Cij = Teuchos::null;
2090  }
2091 
2092  if (!Cij.is_null()) {
2093  if (Cij->isFillComplete())
2094  Cij->resumeFill();
2095  Cij->fillComplete();
2096  bC->setMatrix(i, j, Cij);
2097  } else {
2098  bC->setMatrix(i, j, Teuchos::null);
2099  }
2100  } // loop over rows i
2101  }
2102  else {
2103  // This is the version for blocked matrices
2104 
2105  // check the compatibility of the blocked operators
2106  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2107  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2108  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)");
2109  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)");
2110  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)");
2111  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)");
2112  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)");
2113  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)");
2114 
2115  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2116  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2117 
2118  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2119  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2120 
2121  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2122  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2123 
2124  RCP<Matrix> Cij = Teuchos::null;
2125  if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2126  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2127  // recursive call
2128 
2129  TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2130  *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2131  Cij, fos, AHasFixedNnzPerRow);
2132  } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2133  rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2134  Cij = rcpBopB->getMatrix(i,j);
2135  } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2136  rcpBopB->getMatrix(i,j)==Teuchos::null) {
2137  Cij = rcpBopA->getMatrix(i,j);
2138  } else {
2139  Cij = Teuchos::null;
2140  }
2141 
2142  if (!Cij.is_null()) {
2143  if (Cij->isFillComplete())
2144  Cij->resumeFill();
2145  //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
2146  Cij->fillComplete();
2147  bC->setMatrix(i, j, Cij);
2148  } else {
2149  bC->setMatrix(i, j, Teuchos::null);
2150  }
2151  } // loop over columns j
2152  } // loop over rows i
2153 
2154  } // end blocked recursive algorithm
2155  } //MatrixMatrix::TwoMatrixAdd()
2156  }; // end specialization on GO=long long and NO=EpetraNode
2157 
2158 #endif // HAVE_XPETRA_EPETRA
2159 
2160 } // end namespace Xpetra
2161 
2162 #define XPETRA_MATRIXMATRIX_SHORT
2163 
2164 #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:216
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:1355
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:217
Xpetra::Matrix::getRowMap
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Xpetra_Matrix.hpp:320
Xpetra::BlockedCrsMatrix::isFillComplete
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Definition: Xpetra_BlockedCrsMatrix.hpp:713
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:1690
GO
GlobalOrdinal GO
Definition: Xpetra_UseShortNamesOrdinal.hpp:139
Xpetra::Helpers::Op2NonConstEpetraCrs
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:143
Xpetra_UseShortNames.hpp
Xpetra::Helpers
Xpetra utility class containing transformation routines between Xpetra::Matrix and Epetra/Tpetra obje...
Definition: Xpetra_MatrixMatrix.hpp:93
Xpetra
Xpetra namespace
Definition: Xpetra_BlockedCrsMatrix.hpp:89
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 "in-place".
Definition: Xpetra_MatrixMatrix.hpp:1779
Xpetra_Matrix.hpp
Xpetra::Helpers::Op2EpetraCrs
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:127
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1090
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:1941
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:1760
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:497
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Xpetra::Matrix::CreateView
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
Definition: Xpetra_Matrix.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:199
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:441
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:1611
Xpetra::MatrixMatrix::TwoMatrixMultiplyBlock
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
Definition: Xpetra_MatrixMatrix.hpp:523
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:635
Xpetra::Helpers::epetraExtMult
static void epetraExtMult(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool fillCompleteResult)
Definition: Xpetra_MatrixMatrix.hpp:281
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
Xpetra::Helpers::Op2EpetraCrs
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:99
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:366
Xpetra::Helpers::Op2NonConstEpetraCrs
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:113
Xpetra_ConfigDefs.hpp
Xpetra::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: Xpetra_Exceptions.hpp:108
SC
Scalar SC
Definition: Xpetra_UseShortNamesScalar.hpp:182
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::BlockedCrsMatrix::getRangeMap
RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1081
Xpetra::Matrix::getGlobalNumEntries
virtual global_size_t getGlobalNumEntries() const =0
Returns the global number of entries in this matrix.
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_fwd.hpp:51
Xpetra::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition: Xpetra_BlockedCrsMatrix.hpp:1340
Xpetra::Operator::getRangeMap
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
Xpetra::Matrix::IsView
bool IsView(const viewLabel_t viewLabel) const
Definition: Xpetra_Matrix.hpp:207
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
Definition: Xpetra_MatrixMatrix.hpp:1044
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1093
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::TwoMatrixMultiplyBlock
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
Definition: Xpetra_MatrixMatrix.hpp:1192
Xpetra::BlockedCrsMatrix::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1304
Xpetra::Helpers::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:173
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1307
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition: Xpetra_CrsMatrixWrap_def.hpp:482
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::BlockedCrsMatrix::getDomainMap
RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1066
Xpetra::MatrixMatrix::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
Definition: Xpetra_MatrixMatrix.hpp:505
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:1891
Xpetra::Helpers::Op2TpetraCrs
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:184
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:161
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:1748