Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Xpetra_MatrixMatrix.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 
48 #ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_
49 #define PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_
50 
51 #include "Xpetra_ConfigDefs.hpp"
52 
53 #include "Xpetra_BlockedCrsMatrix.hpp"
54 #include "Xpetra_CrsMatrixWrap.hpp"
55 #include "Xpetra_MapExtractor.hpp"
56 #include "Xpetra_Map.hpp"
57 #include "Xpetra_MatrixFactory.hpp"
58 #include "Xpetra_Matrix.hpp"
59 #include "Xpetra_StridedMapFactory.hpp"
60 #include "Xpetra_StridedMap.hpp"
61 
62 #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_TpetraBlockCrsMatrix.hpp>
78 #include <Tpetra_BlockCrsMatrix_Helpers.hpp>
79 #include <Xpetra_TpetraMultiVector.hpp>
80 #include <Xpetra_TpetraVector.hpp>
81 #endif // HAVE_XPETRA_TPETRA
82 
83 namespace Xpetra {
84 
91 template <class Scalar,
92  class LocalOrdinal = int,
93  class GlobalOrdinal = LocalOrdinal,
94  class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
95 class Helpers {
96 #include "Xpetra_UseShortNames.hpp"
97 
98  public:
99 #ifdef HAVE_XPETRA_EPETRA
100  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<Matrix> Op) {
101  // Get the underlying Epetra Mtx
102  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
103  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast,
104  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
105 
106  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
107  const RCP<const EpetraCrsMatrixT<GO, NO> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO, NO> >(tmp_CrsMtx);
108  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast,
109  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
110 
111  return tmp_ECrsMtx->getEpetra_CrsMatrix();
112  }
113 
114  static RCP<Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) {
115  RCP<Epetra_CrsMatrix> A;
116  // Get the underlying Epetra Mtx
117  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
118  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast,
119  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
120 
121  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
122  const RCP<const EpetraCrsMatrixT<GO, NO> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO, NO> >(tmp_CrsMtx);
123  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
124 
125  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
126  }
127 
128  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) {
129  // Get the underlying Epetra Mtx
130  try {
131  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
132  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
133  const RCP<const EpetraCrsMatrixT<GO, NO> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO, NO> >(tmp_CrsMtx);
134  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast,
135  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
136 
137  return *tmp_ECrsMtx->getEpetra_CrsMatrix();
138 
139  } catch (...) {
140  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
141  }
142  }
143 
144  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(const Matrix& Op) {
145  RCP<Epetra_CrsMatrix> A;
146  // Get the underlying Epetra Mtx
147  try {
148  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
149  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
150  const RCP<const EpetraCrsMatrixT<GO, NO> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO, NO> >(tmp_CrsMtx);
151  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
152 
153  return *Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
154 
155  } catch (...) {
156  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
157  }
158  }
159 #endif // HAVE_XPETRA_EPETRA
160 
161 #ifdef HAVE_XPETRA_TPETRA
162  static RCP<const Tpetra::CrsMatrix<SC, LO, GO, NO> > Op2TpetraCrs(RCP<Matrix> Op) {
163  // Get the underlying Tpetra Mtx
164  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
165  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
166 
167  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
168  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >(tmp_CrsMtx);
169  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
170 
171  return tmp_ECrsMtx->getTpetra_CrsMatrix();
172  }
173 
174  static RCP<Tpetra::CrsMatrix<SC, LO, GO, NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) {
175  // Get the underlying Tpetra Mtx
176  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
177  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
178  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
179  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >(tmp_CrsMtx);
180  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
181 
182  return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
183  }
184 
185  static const Tpetra::CrsMatrix<SC, LO, GO, NO>& Op2TpetraCrs(const Matrix& Op) {
186  // Get the underlying Tpetra Mtx
187  try {
188  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
189  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
190  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >& tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >(tmp_CrsMtx);
191  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
192 
193  return *tmp_TCrsMtx->getTpetra_CrsMatrix();
194 
195  } catch (...) {
196  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
197  }
198  }
199 
200  static Tpetra::CrsMatrix<SC, LO, GO, NO>& Op2NonConstTpetraCrs(const Matrix& Op) {
201  // Get the underlying Tpetra Mtx
202  try {
203  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
204  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
205  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >& tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >(tmp_CrsMtx);
206  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
207 
208  return *Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_TCrsMtx->getTpetra_CrsMatrix());
209 
210  } catch (...) {
211  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
212  }
213  }
214 
215  static bool isTpetraCrs(RCP<Matrix> Op) {
216  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
217  if (crsOp == Teuchos::null) return false;
218  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
219  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >(tmp_CrsMtx);
220  if (tmp_ECrsMtx == Teuchos::null)
221  return false;
222  else
223  return true;
224  }
225 
226  static bool isTpetraCrs(const Matrix& Op) {
227  try {
228  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
229  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
230  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO> >(tmp_CrsMtx);
231  if (tmp_ECrsMtx == Teuchos::null)
232  return false;
233  else
234  return true;
235  } catch (...) {
236  return false;
237  }
238  }
239 
240  static RCP<const Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Op2TpetraBlockCrs(RCP<Matrix> Op) {
241  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
242  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
243 
244  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
245  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
246  TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
247  return tmp_BlockCrs->getTpetra_BlockCrsMatrix();
248  }
249 
250  static RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Op2NonConstTpetraBlockCrs(RCP<Matrix> Op) {
251  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
252  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
253 
254  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
255  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
256  TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
257  return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
258  }
259 
260  static const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& Op2TpetraBlockCrs(const Matrix& Op) {
261  try {
262  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
263  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
264  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
265  TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
266  return *tmp_BlockCrs->getTpetra_BlockCrsMatrix();
267  } catch (...) {
268  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
269  }
270  }
271 
272  static Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& Op2NonConstTpetraBlockCrs(const Matrix& Op) {
273  try {
274  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
275  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
276  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
277  TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
278  return *tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
279  } catch (...) {
280  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
281  }
282  }
283 
284  static bool isTpetraBlockCrs(RCP<Matrix> Op) {
285  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
286  if (crsOp == Teuchos::null) return false;
287  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
288  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
289  if (tmp_BlockCrs == Teuchos::null)
290  return false;
291  else
292  return true;
293  }
294 
295  static bool isTpetraBlockCrs(const Matrix& Op) {
296  try {
297  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
298  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
299  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
300  if (tmp_BlockCrs == Teuchos::null)
301  return false;
302  else
303  return true;
304  } catch (...) {
305  return false;
306  }
307  }
308 #else // HAVE_XPETRA_TPETRA
309  static bool isTpetraCrs(const Matrix& Op) {
310  return false;
311  }
312 
313  static bool isTpetraBlockCrs(const Matrix& Op) {
314  return false;
315  }
316 
317 #endif // HAVE_XPETRA_TPETRA
318 
319 #ifdef HAVE_XPETRA_TPETRA
320  using tcrs_matrix_type = Tpetra::CrsMatrix<SC, LO, GO, NO>;
321  static Teuchos::RCP<Matrix> tpetraAdd(
322  const tcrs_matrix_type& A, bool transposeA, const typename tcrs_matrix_type::scalar_type alpha,
323  const tcrs_matrix_type& B, bool transposeB, const typename tcrs_matrix_type::scalar_type beta) {
324  using Teuchos::Array;
325  using Teuchos::ParameterList;
326  using Teuchos::RCP;
327  using Teuchos::rcp;
328  using Teuchos::rcp_implicit_cast;
329  using Teuchos::rcpFromRef;
330  using XTCrsType = Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>;
331  using CrsType = Xpetra::CrsMatrix<SC, LO, GO, NO>;
332  using CrsWrap = Xpetra::CrsMatrixWrap<SC, LO, GO, NO>;
333  using transposer_type = Tpetra::RowMatrixTransposer<SC, LO, GO, NO>;
334  using import_type = Tpetra::Import<LO, GO, NO>;
335  RCP<const tcrs_matrix_type> Aprime = rcpFromRef(A);
336  if (transposeA)
337  Aprime = transposer_type(Aprime).createTranspose();
338  // Decide whether the fast code path can be taken.
339  if (A.isFillComplete() && B.isFillComplete()) {
340  RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), 0));
341  RCP<ParameterList> addParams = rcp(new ParameterList);
342  addParams->set("Call fillComplete", false);
343  // passing A after B means C will have the same domain/range map as A (or A^T if transposeA)
344  Tpetra::MatrixMatrix::add<SC, LO, GO, NO>(beta, transposeB, B, alpha, false, *Aprime, *C, Teuchos::null, Teuchos::null, addParams);
345  return rcp_implicit_cast<Matrix>(rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C))))));
346  } else {
347  // Slow case - one or both operands are non-fill complete.
348  // TODO: deprecate this.
349  // Need to compute the explicit transpose before add if transposeA and/or transposeB.
350  // This is to count the number of entries to allocate per row in the final sum.
351  RCP<const tcrs_matrix_type> Bprime = rcpFromRef(B);
352  if (transposeB)
353  Bprime = transposer_type(Bprime).createTranspose();
354  // Use Aprime's row map as C's row map.
355  if (!(Aprime->getRowMap()->isSameAs(*(Bprime->getRowMap())))) {
356  auto import = rcp(new import_type(Bprime->getRowMap(), Aprime->getRowMap()));
357  Bprime = Tpetra::importAndFillCompleteCrsMatrix<tcrs_matrix_type>(Bprime, *import, Aprime->getDomainMap(), Aprime->getRangeMap());
358  }
359  // Count the entries to allocate for C in each local row.
360  LO numLocalRows = Aprime->getLocalNumRows();
361  Array<size_t> allocPerRow(numLocalRows);
362  // 0 is always the minimum LID
363  for (LO i = 0; i < numLocalRows; i++) {
364  allocPerRow[i] = Aprime->getNumEntriesInLocalRow(i) + Bprime->getNumEntriesInLocalRow(i);
365  }
366  // Construct C
367  RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), allocPerRow()));
368  // Compute the sum in C (already took care of transposes)
369  Tpetra::MatrixMatrix::Add<SC, LO, GO, NO>(
370  *Aprime, false, alpha,
371  *Bprime, false, beta,
372  C);
373  return rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C)))));
374  }
375  }
376 #endif
377 
378 #ifdef HAVE_XPETRA_EPETRAEXT
379  static void epetraExtMult(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Matrix& C, bool fillCompleteResult) {
380  Epetra_CrsMatrix& epA = Op2NonConstEpetraCrs(A);
381  Epetra_CrsMatrix& epB = Op2NonConstEpetraCrs(B);
382  Epetra_CrsMatrix& epC = Op2NonConstEpetraCrs(C);
383  // Want a static profile (possibly fill complete) matrix as the result.
384  // But, EpetraExt Multiply needs C to be dynamic profile, so compute the product in a temporary DynamicProfile matrix.
385  const Epetra_Map& Crowmap = epC.RowMap();
386  int errCode = 0;
387  Epetra_CrsMatrix Ctemp(::Copy, Crowmap, 0, false);
388  if (fillCompleteResult) {
389  errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, true);
390  if (!errCode) {
391  epC = Ctemp;
392  C.fillComplete();
393  }
394  } else {
395  errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, false);
396  if (!errCode) {
397  int numLocalRows = Crowmap.NumMyElements();
398  long long* globalElementList = nullptr;
399  Crowmap.MyGlobalElementsPtr(globalElementList);
400  Teuchos::Array<int> entriesPerRow(numLocalRows, 0);
401  for (int i = 0; i < numLocalRows; i++) {
402  entriesPerRow[i] = Ctemp.NumGlobalEntries(globalElementList[i]);
403  }
404  // know how many entries to allocate in epC (which must be static profile)
405  epC = Epetra_CrsMatrix(::Copy, Crowmap, entriesPerRow.data(), true);
406  for (int i = 0; i < numLocalRows; i++) {
407  int gid = globalElementList[i];
408  int numEntries;
409  double* values;
410  int* inds;
411  Ctemp.ExtractGlobalRowView(gid, numEntries, values, inds);
412  epC.InsertGlobalValues(gid, numEntries, values, inds);
413  }
414  }
415  }
416  if (errCode) {
417  std::ostringstream buf;
418  buf << errCode;
419  std::string msg = "EpetraExt::MatrixMatrix::Multiply returned nonzero error code " + buf.str();
420  throw(Exceptions::RuntimeError(msg));
421  }
422  }
423 #endif
424 };
425 
426 template <class Scalar,
427  class LocalOrdinal /*= int*/,
428  class GlobalOrdinal /*= LocalOrdinal*/,
429  class Node /*= Tpetra::KokkosClassic::DefaultNode::DefaultNodeType*/>
430 class MatrixMatrix {
431 #undef XPETRA_MATRIXMATRIX_SHORT
432 #include "Xpetra_UseShortNames.hpp"
433 
434  public:
459  static void Multiply(const Matrix& A, bool transposeA,
460  const Matrix& B, bool transposeB,
461  Matrix& C,
462  bool call_FillComplete_on_result = true,
463  bool doOptimizeStorage = true,
464  const std::string& label = std::string(),
465  const RCP<ParameterList>& params = null) {
466  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
467  Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
468  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
469  Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
470 
471  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
472  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
473 
474  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
475 
476  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
477 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
478  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
479 #else
480  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
481 
482 #endif
483  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
484 #ifdef HAVE_XPETRA_TPETRA
485  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
486  if (helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
487  // All matrices are Crs
488  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = helpers::Op2TpetraCrs(A);
489  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = helpers::Op2TpetraCrs(B);
490  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpC = helpers::Op2NonConstTpetraCrs(C);
491 
492  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
493  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
494  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
495  } else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
496  // All matrices are BlockCrs (except maybe Ac)
497  // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
498  // we'll need to think about refactoring BlockCrs so we can do something smarter here.
499  if (!A.getRowMap()->getComm()->getRank())
500  std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
501 
502  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
503  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(B);
504  using CRS = Tpetra::CrsMatrix<SC, LO, GO, NO>;
505  RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
506  RCP<const CRS> Bcrs = Tpetra::convertToCrsMatrix(tpB);
507 
508  // We need the global constants to do the copy back to BlockCrs
509  RCP<ParameterList> new_params;
510  if (!params.is_null()) {
511  new_params = rcp(new Teuchos::ParameterList(*params));
512  new_params->set("compute global constants", true);
513  }
514 
515  // FIXME: The lines below only works because we're assuming Ac is Point
516  RCP<CRS> tempAc = Teuchos::rcp(new CRS(Acrs->getRowMap(), 0));
517  Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
518 
519  // Temporary output matrix
520  RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*tempAc, A.GetStorageBlockSize());
521  RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
522  RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > Ac_p = Ac_x;
523 
524  // We can now cheat and replace the innards of Ac
525  RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(Teuchos::rcpFromRef(C));
526  Ac_w->replaceCrsMatrix(Ac_p);
527  } else {
528  // Mix and match
529  TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
530  }
531 #else
532  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
533 #endif
534  }
535 
536  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
537  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
538  fillParams->set("Optimize Storage", doOptimizeStorage);
539  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
540  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
541  fillParams);
542  }
543 
544  // transfer striding information
545  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
546  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
547  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
548  } // end Multiply
549 
572  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, RCP<Matrix> C_in,
573  Teuchos::FancyOStream& fos,
574  bool doFillComplete = true,
575  bool doOptimizeStorage = true,
576  const std::string& label = std::string(),
577  const RCP<ParameterList>& params = null) {
578  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
579  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
580 
581  // Default case: Xpetra Multiply
582  RCP<Matrix> C = C_in;
583 
584  if (C == Teuchos::null) {
585  double nnzPerRow = Teuchos::as<double>(0);
586 
587 #if 0
588  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
589  // For now, follow what ML and Epetra do.
590  GO numRowsA = A.getGlobalNumRows();
591  GO numRowsB = B.getGlobalNumRows();
592  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
593  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
594  nnzPerRow *= nnzPerRow;
595  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
596  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
597  if (totalNnz < minNnz)
598  totalNnz = minNnz;
599  nnzPerRow = totalNnz / A.getGlobalNumRows();
600 
601  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
602  }
603 #endif
604  if (transposeA)
605  C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
606  else
607  C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
608 
609  } else {
610  C->resumeFill(); // why this is not done inside of Tpetra MxM?
611  fos << "Reuse C pattern" << std::endl;
612  }
613 
614  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
615 
616  return C;
617  }
618 
629  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Teuchos::FancyOStream& fos,
630  bool callFillCompleteOnResult = true, bool doOptimizeStorage = true, const std::string& label = std::string(),
631  const RCP<ParameterList>& params = null) {
632  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
633  }
634 
635 #ifdef HAVE_XPETRA_EPETRAEXT
636  // Michael Gee's MLMultiply
637  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
638  const Epetra_CrsMatrix& epB,
639  Teuchos::FancyOStream& fos) {
640  throw(Xpetra::Exceptions::RuntimeError("MLTwoMatrixMultiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
641  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
642  }
643 #endif // ifdef HAVE_XPETRA_EPETRAEXT
644 
655  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
656  const BlockedCrsMatrix& B, bool transposeB,
657  Teuchos::FancyOStream& fos,
658  bool doFillComplete = true,
659  bool doOptimizeStorage = true) {
660  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
661  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
662 
663  // Preconditions
664  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
665  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
666 
667  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
668  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
669 
670  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
671 
672  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
673  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
674  RCP<Matrix> Cij;
675 
676  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
677  RCP<Matrix> crmat1 = A.getMatrix(i, l);
678  RCP<Matrix> crmat2 = B.getMatrix(l, j);
679 
680  if (crmat1.is_null() || crmat2.is_null())
681  continue;
682 
683  // try unwrapping 1x1 blocked matrix
684  {
685  auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
686  auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
687 
688  if (unwrap1.is_null() != unwrap2.is_null()) {
689  if (unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
690  crmat1 = unwrap1->getCrsMatrix();
691  if (unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
692  crmat2 = unwrap2->getCrsMatrix();
693  }
694  }
695 
696  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
697  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
698  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
699  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
700 
701  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
702  if (!crop1.is_null())
703  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
704  if (!crop2.is_null())
705  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
706 
707  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
708  "crmat1 does not have global constants");
709  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
710  "crmat2 does not have global constants");
711 
712  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
713  continue;
714 
715  // temporary matrix containing result of local block multiplication
716  RCP<Matrix> temp = Teuchos::null;
717 
718  if (crop1 != Teuchos::null && crop2 != Teuchos::null)
719  temp = Multiply(*crop1, false, *crop2, false, fos);
720  else {
721  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
722  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
723  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
724  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
725  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)");
726  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)");
727 
728  // recursive multiplication call
729  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
730 
731  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
732  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)");
733  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)");
734  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)");
735  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)");
736  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)");
737  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)");
738  }
739 
740  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
741 
742  RCP<Matrix> addRes = null;
743  if (Cij.is_null())
744  Cij = temp;
745  else {
746  MatrixMatrix::TwoMatrixAdd(*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
747  Cij = addRes;
748  }
749  }
750 
751  if (!Cij.is_null()) {
752  if (Cij->isFillComplete())
753  Cij->resumeFill();
754  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
755  C->setMatrix(i, j, Cij);
756  } else {
757  C->setMatrix(i, j, Teuchos::null);
758  }
759  }
760  }
761 
762  if (doFillComplete)
763  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
764 
765  return C;
766  } // TwoMatrixMultiplyBlock
767 
780  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
781  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
782  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
783 
784  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
785  throw Exceptions::RuntimeError("TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
786  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
787 #ifdef HAVE_XPETRA_TPETRA
788  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
789  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(B);
790 
791  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
792 #else
793  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
794 #endif
795  }
796  } // MatrixMatrix::TwoMatrixAdd()
797 
812  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
813  const Matrix& B, bool transposeB, const SC& beta,
814  RCP<Matrix>& C, Teuchos::FancyOStream& fos, bool AHasFixedNnzPerRow = false) {
815  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
816  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
817  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
818  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
819  // We have to distinguish 4 cases:
820  // both matrices are CrsMatrixWrap based, only one of them is CrsMatrixWrap based, or none.
821 
822  // C can be null, so just use A to get the lib
823  auto lib = A.getRowMap()->lib();
824 
825  // Both matrices are CrsMatrixWrap
826  if (rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
827  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
828  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
829  if (lib == Xpetra::UseEpetra) {
830  throw Exceptions::RuntimeError("MatrixMatrix::Add for Epetra only available with Scalar = double, LO = GO = int.");
831  } else if (lib == Xpetra::UseTpetra) {
832 #ifdef HAVE_XPETRA_TPETRA
833  using tcrs_matrix_type = Tpetra::CrsMatrix<SC, LO, GO, NO>;
834  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
835  const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
836  const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
837  C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
838 #else
839  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
840 #endif
841  }
843  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
844  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
846  }
847  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
848  else if (rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
849  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
850  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
851 
852  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
853  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
854 
855  size_t i = 0;
856  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
857  RCP<Matrix> Cij = Teuchos::null;
858  if (rcpA != Teuchos::null &&
859  rcpBopB->getMatrix(i, j) != Teuchos::null) {
860  // recursive call
861  TwoMatrixAdd(*rcpA, transposeA, alpha,
862  *(rcpBopB->getMatrix(i, j)), transposeB, beta,
863  Cij, fos, AHasFixedNnzPerRow);
864  } else if (rcpA == Teuchos::null &&
865  rcpBopB->getMatrix(i, j) != Teuchos::null) {
866  Cij = rcpBopB->getMatrix(i, j);
867  } else if (rcpA != Teuchos::null &&
868  rcpBopB->getMatrix(i, j) == Teuchos::null) {
869  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
870  } else {
871  Cij = Teuchos::null;
872  }
873 
874  if (!Cij.is_null()) {
875  if (Cij->isFillComplete())
876  Cij->resumeFill();
877  Cij->fillComplete();
878  bC->setMatrix(i, j, Cij);
879  } else {
880  bC->setMatrix(i, j, Teuchos::null);
881  }
882  } // loop over columns j
883  }
884  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
885  else if (rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
886  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
887  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
888 
889  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
890  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
891  size_t j = 0;
892  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
893  RCP<Matrix> Cij = Teuchos::null;
894  if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
895  rcpB != Teuchos::null) {
896  // recursive call
897  TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
898  *rcpB, transposeB, beta,
899  Cij, fos, AHasFixedNnzPerRow);
900  } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
901  rcpB != Teuchos::null) {
902  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
903  } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
904  rcpB == Teuchos::null) {
905  Cij = rcpBopA->getMatrix(i, j);
906  } else {
907  Cij = Teuchos::null;
908  }
909 
910  if (!Cij.is_null()) {
911  if (Cij->isFillComplete())
912  Cij->resumeFill();
913  Cij->fillComplete();
914  bC->setMatrix(i, j, Cij);
915  } else {
916  bC->setMatrix(i, j, Teuchos::null);
917  }
918  } // loop over rows i
919  } else {
920  // This is the version for blocked matrices
921 
922  // check the compatibility of the blocked operators
923  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
924  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
925  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)");
926  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)");
927  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)");
928  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)");
929  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)");
930  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)");
931 
932  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
933  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
934 
935  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
936  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
937  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
938  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
939 
940  RCP<Matrix> Cij = Teuchos::null;
941  if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
942  rcpBopB->getMatrix(i, j) != Teuchos::null) {
943  // recursive call
944  TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
945  *(rcpBopB->getMatrix(i, j)), transposeB, beta,
946  Cij, fos, AHasFixedNnzPerRow);
947  } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
948  rcpBopB->getMatrix(i, j) != Teuchos::null) {
949  Cij = rcpBopB->getMatrix(i, j);
950  } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
951  rcpBopB->getMatrix(i, j) == Teuchos::null) {
952  Cij = rcpBopA->getMatrix(i, j);
953  } else {
954  Cij = Teuchos::null;
955  }
956 
957  if (!Cij.is_null()) {
958  if (Cij->isFillComplete())
959  Cij->resumeFill();
960  Cij->fillComplete();
961  bC->setMatrix(i, j, Cij);
962  } else {
963  bC->setMatrix(i, j, Teuchos::null);
964  }
965  } // loop over columns j
966  } // loop over rows i
967 
968  } // end blocked recursive algorithm
969  } // MatrixMatrix::TwoMatrixAdd()
970 
971 }; // class MatrixMatrix
972 
973 #ifdef HAVE_XPETRA_EPETRA
974 // specialization MatrixMatrix for SC=double, LO=GO=int
975 template <>
976 class MatrixMatrix<double, int, int, EpetraNode> {
977  typedef double Scalar;
978  typedef int LocalOrdinal;
979  typedef int GlobalOrdinal;
980  typedef EpetraNode Node;
981 #include "Xpetra_UseShortNames.hpp"
982 
983  public:
1008  static void Multiply(const Matrix& A, bool transposeA,
1009  const Matrix& B, bool transposeB,
1010  Matrix& C,
1011  bool call_FillComplete_on_result = true,
1012  bool doOptimizeStorage = true,
1013  const std::string& label = std::string(),
1014  const RCP<ParameterList>& params = null) {
1015  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
1016  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
1017  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
1018  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
1019 
1020  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
1021  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
1022 
1023  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1024 
1025  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
1026 
1027  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
1028 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1029  helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
1030 #else
1031  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
1032 #endif
1033  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
1034 #ifdef HAVE_XPETRA_TPETRA
1035 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1036  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1037  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
1038 #else
1039  if (helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
1040  // All matrices are Crs
1041  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = helpers::Op2TpetraCrs(A);
1042  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = helpers::Op2TpetraCrs(B);
1043  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpC = helpers::Op2NonConstTpetraCrs(C);
1044 
1045  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
1046  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
1047  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1048  } else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
1049  // All matrices are BlockCrs (except maybe Ac)
1050  // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
1051  // we'll need to think about refactoring BlockCrs so we can do something smarter here.
1052 
1053  if (!A.getRowMap()->getComm()->getRank())
1054  std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
1055 
1056  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
1057  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(B);
1058  using CRS = Tpetra::CrsMatrix<SC, LO, GO, NO>;
1059  RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
1060  RCP<const CRS> Bcrs = Tpetra::convertToCrsMatrix(tpB);
1061 
1062  // We need the global constants to do the copy back to BlockCrs
1063  RCP<ParameterList> new_params;
1064  if (!params.is_null()) {
1065  new_params = rcp(new Teuchos::ParameterList(*params));
1066  new_params->set("compute global constants", true);
1067  }
1068 
1069  // FIXME: The lines below only works because we're assuming Ac is Point
1070  RCP<CRS> tempAc = Teuchos::rcp(new CRS(Acrs->getRowMap(), 0));
1071  Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
1072 
1073  // Temporary output matrix
1074  RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*tempAc, A.GetStorageBlockSize());
1075  RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
1076  RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > Ac_p = Ac_x;
1077 
1078  // We can now cheat and replace the innards of Ac
1079  RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(Teuchos::rcpFromRef(C));
1080  Ac_w->replaceCrsMatrix(Ac_p);
1081  } else {
1082  // Mix and match
1083  TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
1084  }
1085 #endif
1086 #else
1087  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
1088 #endif
1089  }
1090 
1091  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1092  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
1093  fillParams->set("Optimize Storage", doOptimizeStorage);
1094  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
1095  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
1096  fillParams);
1097  }
1098 
1099  // transfer striding information
1100  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1101  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1102  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
1103  } // end Multiply
1104 
1127  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1128  const Matrix& B, bool transposeB,
1129  RCP<Matrix> C_in,
1130  Teuchos::FancyOStream& fos,
1131  bool doFillComplete = true,
1132  bool doOptimizeStorage = true,
1133  const std::string& label = std::string(),
1134  const RCP<ParameterList>& params = null) {
1135  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1136  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1137 
1138  // Optimization using ML Multiply when available and requested
1139  // This feature is currently not supported. We would have to introduce the HAVE_XPETRA_ML_MMM flag
1140 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) && defined(HAVE_XPETRA_ML_MMM)
1141  if (B.getDomainMap()->lib() == Xpetra::UseEpetra && !transposeA && !transposeB) {
1142  RCP<const Epetra_CrsMatrix> epA = Xpetra::Helpers<SC, LO, GO, NO>::Op2EpetraCrs(rcpFromRef(A));
1143  RCP<const Epetra_CrsMatrix> epB = Xpetra::Helpers<SC, LO, GO, NO>::Op2EpetraCrs(rcpFromRef(B));
1144  RCP<Epetra_CrsMatrix> epC = MLTwoMatrixMultiply(*epA, *epB, fos);
1145 
1146  RCP<Matrix> C = Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC, LO, GO, NO>(epC);
1147  if (doFillComplete) {
1148  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
1149  fillParams->set("Optimize Storage", doOptimizeStorage);
1150  C->fillComplete(B.getDomainMap(), A.getRangeMap(), fillParams);
1151  }
1152 
1153  // Fill strided maps information
1154  // This is necessary since the ML matrix matrix multiplication routine has no handling for this
1155  // TODO: move this call to MLMultiply...
1156  C->CreateView("stridedMaps", rcpFromRef(A), transposeA, rcpFromRef(B), transposeB);
1157 
1158  return C;
1159  }
1160 #endif // EPETRA + EPETRAEXT + ML
1161 
1162  // Default case: Xpetra Multiply
1163  RCP<Matrix> C = C_in;
1164 
1165  if (C == Teuchos::null) {
1166  double nnzPerRow = Teuchos::as<double>(0);
1167 
1168 #if 0
1169  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
1170  // For now, follow what ML and Epetra do.
1171  GO numRowsA = A.getGlobalNumRows();
1172  GO numRowsB = B.getGlobalNumRows();
1173  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
1174  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
1175  nnzPerRow *= nnzPerRow;
1176  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1177  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1178  if (totalNnz < minNnz)
1179  totalNnz = minNnz;
1180  nnzPerRow = totalNnz / A.getGlobalNumRows();
1181 
1182  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1183  }
1184 #endif
1185 
1186  if (transposeA)
1187  C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1188  else
1189  C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1190 
1191  } else {
1192  C->resumeFill(); // why this is not done inside of Tpetra MxM?
1193  fos << "Reuse C pattern" << std::endl;
1194  }
1195 
1196  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1197 
1198  return C;
1199  }
1200 
1211  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1212  const Matrix& B, bool transposeB,
1213  Teuchos::FancyOStream& fos,
1214  bool callFillCompleteOnResult = true,
1215  bool doOptimizeStorage = true,
1216  const std::string& label = std::string(),
1217  const RCP<ParameterList>& params = null) {
1218  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1219  }
1220 
1221 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1222  // Michael Gee's MLMultiply
1223  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
1224  const Epetra_CrsMatrix& epB,
1225  Teuchos::FancyOStream& fos) {
1226 #if defined(HAVE_XPETRA_ML_MMM) // Note: this is currently not supported
1227  ML_Comm* comm;
1228  ML_Comm_Create(&comm);
1229  fos << "****** USING ML's MATRIX MATRIX MULTIPLY ******" << std::endl;
1230 #ifdef HAVE_MPI
1231  // ML_Comm uses MPI_COMM_WORLD, so try to use the same communicator as epA.
1232  const Epetra_MpiComm* Mcomm = dynamic_cast<const Epetra_MpiComm*>(&(epA.Comm()));
1233  if (Mcomm)
1234  ML_Comm_Set_UsrComm(comm, Mcomm->GetMpiComm());
1235 #endif
1236  // in order to use ML, there must be no indices missing from the matrix column maps.
1237  EpetraExt::CrsMatrix_SolverMap Atransform;
1238  EpetraExt::CrsMatrix_SolverMap Btransform;
1239  const Epetra_CrsMatrix& A = Atransform(const_cast<Epetra_CrsMatrix&>(epA));
1240  const Epetra_CrsMatrix& B = Btransform(const_cast<Epetra_CrsMatrix&>(epB));
1241 
1242  if (!A.Filled()) throw Exceptions::RuntimeError("A has to be FillCompleted");
1243  if (!B.Filled()) throw Exceptions::RuntimeError("B has to be FillCompleted");
1244 
1245  // create ML operators from EpetraCrsMatrix
1246  ML_Operator* ml_As = ML_Operator_Create(comm);
1247  ML_Operator* ml_Bs = ML_Operator_Create(comm);
1248  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?
1249  ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&B), ml_Bs);
1250  ML_Operator* ml_AtimesB = ML_Operator_Create(comm);
1251  {
1252  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("ML_2matmult kernel"));
1253  ML_2matmult(ml_As, ml_Bs, ml_AtimesB, ML_CSR_MATRIX); // do NOT use ML_EpetraCRS_MATRIX!!!
1254  }
1255  ML_Operator_Destroy(&ml_As);
1256  ML_Operator_Destroy(&ml_Bs);
1257 
1258  // For ml_AtimesB we have to reconstruct the column map in global indexing,
1259  // The following is going down to the salt-mines of ML ...
1260  // note: we use integers, since ML only works for Epetra...
1261  int N_local = ml_AtimesB->invec_leng;
1262  ML_CommInfoOP* getrow_comm = ml_AtimesB->getrow->pre_comm;
1263  if (!getrow_comm) throw(Exceptions::RuntimeError("ML_Operator does not have a CommInfo"));
1264  ML_Comm* comm_AB = ml_AtimesB->comm; // get comm object
1265  if (N_local != B.DomainMap().NumMyElements())
1266  throw(Exceptions::RuntimeError("Mismatch in local row dimension between ML and Epetra"));
1267  int N_rcvd = 0;
1268  int N_send = 0;
1269  int flag = 0;
1270  for (int i = 0; i < getrow_comm->N_neighbors; i++) {
1271  N_rcvd += (getrow_comm->neighbors)[i].N_rcv;
1272  N_send += (getrow_comm->neighbors)[i].N_send;
1273  if (((getrow_comm->neighbors)[i].N_rcv != 0) &&
1274  ((getrow_comm->neighbors)[i].rcv_list != NULL)) flag = 1;
1275  }
1276  // For some unknown reason, ML likes to have stuff one larger than
1277  // neccessary...
1278  std::vector<double> dtemp(N_local + N_rcvd + 1); // "double" vector for comm function
1279  std::vector<int> cmap(N_local + N_rcvd + 1); // vector for gids
1280  for (int i = 0; i < N_local; ++i) {
1281  cmap[i] = B.DomainMap().GID(i);
1282  dtemp[i] = (double)cmap[i];
1283  }
1284  ML_cheap_exchange_bdry(&dtemp[0], getrow_comm, N_local, N_send, comm_AB); // do communication
1285  if (flag) { // process received data
1286  int count = N_local;
1287  const int neighbors = getrow_comm->N_neighbors;
1288  for (int i = 0; i < neighbors; i++) {
1289  const int nrcv = getrow_comm->neighbors[i].N_rcv;
1290  for (int j = 0; j < nrcv; j++)
1291  cmap[getrow_comm->neighbors[i].rcv_list[j]] = (int)dtemp[count++];
1292  }
1293  } else {
1294  for (int i = 0; i < N_local + N_rcvd; ++i)
1295  cmap[i] = (int)dtemp[i];
1296  }
1297  dtemp.clear(); // free double array
1298 
1299  // we can now determine a matching column map for the result
1300  Epetra_Map gcmap(-1, N_local + N_rcvd, &cmap[0], B.ColMap().IndexBase(), A.Comm());
1301 
1302  int allocated = 0;
1303  int rowlength;
1304  double* val = NULL;
1305  int* bindx = NULL;
1306 
1307  const int myrowlength = A.RowMap().NumMyElements();
1308  const Epetra_Map& rowmap = A.RowMap();
1309 
1310  // Determine the maximum bandwith for the result matrix.
1311  // replaces the old, very(!) memory-consuming guess:
1312  // int guessnpr = A.MaxNumEntries()*B.MaxNumEntries();
1313  int educatedguess = 0;
1314  for (int i = 0; i < myrowlength; ++i) {
1315  // get local row
1316  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1317  if (rowlength > educatedguess)
1318  educatedguess = rowlength;
1319  }
1320 
1321  // allocate our result matrix and fill it
1322  RCP<Epetra_CrsMatrix> result = rcp(new Epetra_CrsMatrix(::Copy, A.RangeMap(), gcmap, educatedguess, false));
1323 
1324  std::vector<int> gcid(educatedguess);
1325  for (int i = 0; i < myrowlength; ++i) {
1326  const int grid = rowmap.GID(i);
1327  // get local row
1328  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1329  if (!rowlength) continue;
1330  if ((int)gcid.size() < rowlength) gcid.resize(rowlength);
1331  for (int j = 0; j < rowlength; ++j) {
1332  gcid[j] = gcmap.GID(bindx[j]);
1333  if (gcid[j] < 0)
1334  throw Exceptions::RuntimeError("Error: cannot find gcid!");
1335  }
1336  int err = result->InsertGlobalValues(grid, rowlength, val, &gcid[0]);
1337  if (err != 0 && err != 1) {
1338  std::ostringstream errStr;
1339  errStr << "Epetra_CrsMatrix::InsertGlobalValues returned err=" << err;
1340  throw Exceptions::RuntimeError(errStr.str());
1341  }
1342  }
1343  // free memory
1344  if (bindx) ML_free(bindx);
1345  if (val) ML_free(val);
1346  ML_Operator_Destroy(&ml_AtimesB);
1347  ML_Comm_Destroy(&comm);
1348 
1349  return result;
1350 #else // no MUELU_ML
1351  (void)epA;
1352  (void)epB;
1353  (void)fos;
1354  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1355  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1356  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1357 #endif
1358  }
1359 #endif // ifdef HAVE_XPETRA_EPETRAEXT
1360 
1371  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
1372  const BlockedCrsMatrix& B, bool transposeB,
1373  Teuchos::FancyOStream& fos,
1374  bool doFillComplete = true,
1375  bool doOptimizeStorage = true) {
1376  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1377  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1378 
1379  // Preconditions
1380  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1381  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1382 
1383  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1384  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
1385 
1386  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1387 
1388  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1389  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1390  RCP<Matrix> Cij;
1391 
1392  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1393  RCP<Matrix> crmat1 = A.getMatrix(i, l);
1394  RCP<Matrix> crmat2 = B.getMatrix(l, j);
1395 
1396  if (crmat1.is_null() || crmat2.is_null())
1397  continue;
1398 
1399  // try unwrapping 1x1 blocked matrix
1400  {
1401  auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1402  auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1403 
1404  if (unwrap1.is_null() != unwrap2.is_null()) {
1405  if (unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
1406  crmat1 = unwrap1->getCrsMatrix();
1407  if (unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
1408  crmat2 = unwrap2->getCrsMatrix();
1409  }
1410  }
1411 
1412  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1413  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1414  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
1415  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1416 
1417  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
1418  if (!crop1.is_null())
1419  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
1420  if (!crop2.is_null())
1421  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
1422 
1423  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
1424  "crmat1 does not have global constants");
1425  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
1426  "crmat2 does not have global constants");
1427 
1428  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1429  continue;
1430 
1431  // temporary matrix containing result of local block multiplication
1432  RCP<Matrix> temp = Teuchos::null;
1433 
1434  if (crop1 != Teuchos::null && crop2 != Teuchos::null)
1435  temp = Multiply(*crop1, false, *crop2, false, fos);
1436  else {
1437  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1438  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1439  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1440  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1441  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)");
1442  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)");
1443 
1444  // recursive multiplication call
1445  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1446 
1447  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1448  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)");
1449  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)");
1450  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)");
1451  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)");
1452  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)");
1453  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)");
1454  }
1455 
1456  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1457 
1458  RCP<Matrix> addRes = null;
1459  if (Cij.is_null())
1460  Cij = temp;
1461  else {
1462  MatrixMatrix::TwoMatrixAdd(*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1463  Cij = addRes;
1464  }
1465  }
1466 
1467  if (!Cij.is_null()) {
1468  if (Cij->isFillComplete())
1469  Cij->resumeFill();
1470  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1471  C->setMatrix(i, j, Cij);
1472  } else {
1473  C->setMatrix(i, j, Teuchos::null);
1474  }
1475  }
1476  }
1477 
1478  if (doFillComplete)
1479  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1480 
1481  return C;
1482  } // TwoMatrixMultiplyBlock
1483 
1496  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1497  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*B.getRowMap())), Exceptions::Incompatible,
1498  "TwoMatrixAdd: matrix row maps are not the same.");
1499 
1500  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1501 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1502  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC, LO, GO, NO>::Op2EpetraCrs(A);
1503  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstEpetraCrs(B);
1504 
1505  // FIXME is there a bug if beta=0?
1506  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1507 
1508  if (rv != 0)
1509  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1510  std::ostringstream buf;
1511 #else
1512  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1513 #endif
1514  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1515 #ifdef HAVE_XPETRA_TPETRA
1516 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1517  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1518  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1519 #else
1520  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
1521  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(B);
1522 
1523  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1524 #endif
1525 #else
1526  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1527 #endif
1528  }
1529  } // MatrixMatrix::TwoMatrixAdd()
1530 
1545  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1546  const Matrix& B, bool transposeB, const SC& beta,
1547  RCP<Matrix>& C, Teuchos::FancyOStream& fos, bool AHasFixedNnzPerRow = false) {
1548  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
1549  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1550  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1551  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1552  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1553 
1554  if (rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1555  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
1556  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
1557 
1558  auto lib = A.getRowMap()->lib();
1559  if (lib == Xpetra::UseEpetra) {
1560 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1561  const Epetra_CrsMatrix& epA = helpers::Op2EpetraCrs(A);
1562  const Epetra_CrsMatrix& epB = helpers::Op2EpetraCrs(B);
1563  if (C.is_null()) {
1564  size_t maxNzInA = 0;
1565  size_t maxNzInB = 0;
1566  size_t numLocalRows = 0;
1567  if (A.isFillComplete() && B.isFillComplete()) {
1568  maxNzInA = A.getLocalMaxNumRowEntries();
1569  maxNzInB = B.getLocalMaxNumRowEntries();
1570  numLocalRows = A.getLocalNumRows();
1571  }
1572 
1573  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1574  // first check if either A or B has at most 1 nonzero per row
1575  // the case of both having at most 1 nz per row is handled by the ``else''
1576  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1577 
1578  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1579  for (size_t i = 0; i < numLocalRows; ++i)
1580  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1581 
1582  } else {
1583  for (size_t i = 0; i < numLocalRows; ++i)
1584  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1585  }
1586 
1587  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1588  << ", using static profiling" << std::endl;
1589  C = rcp(new Xpetra::CrsMatrixWrap<SC, LO, GO, NO>(A.getRowMap(), exactNnzPerRow));
1590 
1591  } else {
1592  // general case
1593  LO maxPossibleNNZ = A.getLocalMaxNumRowEntries() + B.getLocalMaxNumRowEntries();
1594  C = rcp(new Xpetra::CrsMatrixWrap<SC, LO, GO, NO>(A.getRowMap(), maxPossibleNNZ));
1595  }
1596  if (transposeB)
1597  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1598  }
1599  RCP<Epetra_CrsMatrix> epC = helpers::Op2NonConstEpetraCrs(C);
1600  Epetra_CrsMatrix* ref2epC = &*epC; // to avoid a compiler error...
1601 
1602  // FIXME is there a bug if beta=0?
1603  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1604 
1605  if (rv != 0)
1606  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1607 #else
1608  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1609 #endif
1610  } else if (lib == Xpetra::UseTpetra) {
1611 #ifdef HAVE_XPETRA_TPETRA
1612 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1613  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1614  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1615 #else
1616  using tcrs_matrix_type = Tpetra::CrsMatrix<SC, LO, GO, NO>;
1617  const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
1618  const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
1619  C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
1620 #endif
1621 #else
1622  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
1623 #endif
1624  }
1625 
1627  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
1628  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
1630  }
1631  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
1632  else if (rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1633  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1634  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1635 
1636  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1637  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1638 
1639  size_t i = 0;
1640  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1641  RCP<Matrix> Cij = Teuchos::null;
1642  if (rcpA != Teuchos::null &&
1643  rcpBopB->getMatrix(i, j) != Teuchos::null) {
1644  // recursive call
1645  TwoMatrixAdd(*rcpA, transposeA, alpha,
1646  *(rcpBopB->getMatrix(i, j)), transposeB, beta,
1647  Cij, fos, AHasFixedNnzPerRow);
1648  } else if (rcpA == Teuchos::null &&
1649  rcpBopB->getMatrix(i, j) != Teuchos::null) {
1650  Cij = rcpBopB->getMatrix(i, j);
1651  } else if (rcpA != Teuchos::null &&
1652  rcpBopB->getMatrix(i, j) == Teuchos::null) {
1653  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
1654  } else {
1655  Cij = Teuchos::null;
1656  }
1657 
1658  if (!Cij.is_null()) {
1659  if (Cij->isFillComplete())
1660  Cij->resumeFill();
1661  Cij->fillComplete();
1662  bC->setMatrix(i, j, Cij);
1663  } else {
1664  bC->setMatrix(i, j, Teuchos::null);
1665  }
1666  } // loop over columns j
1667  }
1668  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
1669  else if (rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
1670  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1671  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
1672 
1673  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1674  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1675 
1676  size_t j = 0;
1677  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1678  RCP<Matrix> Cij = Teuchos::null;
1679  if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
1680  rcpB != Teuchos::null) {
1681  // recursive call
1682  TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
1683  *rcpB, transposeB, beta,
1684  Cij, fos, AHasFixedNnzPerRow);
1685  } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
1686  rcpB != Teuchos::null) {
1687  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
1688  } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
1689  rcpB == Teuchos::null) {
1690  Cij = rcpBopA->getMatrix(i, j);
1691  } else {
1692  Cij = Teuchos::null;
1693  }
1694 
1695  if (!Cij.is_null()) {
1696  if (Cij->isFillComplete())
1697  Cij->resumeFill();
1698  Cij->fillComplete();
1699  bC->setMatrix(i, j, Cij);
1700  } else {
1701  bC->setMatrix(i, j, Teuchos::null);
1702  }
1703  } // loop over rows i
1704  } else {
1705  // This is the version for blocked matrices
1706 
1707  // check the compatibility of the blocked operators
1708  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1709  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1710  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)");
1711  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)");
1712  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)");
1713  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)");
1714  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)");
1715  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)");
1716 
1717  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1718  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1719 
1720  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1721  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1722 
1723  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1724  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1725 
1726  RCP<Matrix> Cij = Teuchos::null;
1727  if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
1728  rcpBopB->getMatrix(i, j) != Teuchos::null) {
1729  // recursive call
1730 
1731  TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
1732  *(rcpBopB->getMatrix(i, j)), transposeB, beta,
1733  Cij, fos, AHasFixedNnzPerRow);
1734  } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
1735  rcpBopB->getMatrix(i, j) != Teuchos::null) {
1736  Cij = rcpBopB->getMatrix(i, j);
1737  } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
1738  rcpBopB->getMatrix(i, j) == Teuchos::null) {
1739  Cij = rcpBopA->getMatrix(i, j);
1740  } else {
1741  Cij = Teuchos::null;
1742  }
1743 
1744  if (!Cij.is_null()) {
1745  if (Cij->isFillComplete())
1746  Cij->resumeFill();
1747  // Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
1748  Cij->fillComplete();
1749  bC->setMatrix(i, j, Cij);
1750  } else {
1751  bC->setMatrix(i, j, Teuchos::null);
1752  }
1753  } // loop over columns j
1754  } // loop over rows i
1755 
1756  } // end blocked recursive algorithm
1757  } // MatrixMatrix::TwoMatrixAdd()
1758 }; // end specialization on SC=double, GO=int and NO=EpetraNode
1759 
1760 // specialization MatrixMatrix for SC=double, GO=long long and NO=EptraNode
1761 template <>
1762 class MatrixMatrix<double, int, long long, EpetraNode> {
1763  typedef double Scalar;
1764  typedef int LocalOrdinal;
1765  typedef long long GlobalOrdinal;
1766  typedef EpetraNode Node;
1767 #include "Xpetra_UseShortNames.hpp"
1768 
1769  public:
1794  static void Multiply(const Matrix& A, bool transposeA,
1795  const Matrix& B, bool transposeB,
1796  Matrix& C,
1797  bool call_FillComplete_on_result = true,
1798  bool doOptimizeStorage = true,
1799  const std::string& label = std::string(),
1800  const RCP<ParameterList>& params = null) {
1801  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
1802  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
1803  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
1804  TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
1805  Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
1806 
1807  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError, "A is not fill-completed");
1808  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Xpetra::Exceptions::RuntimeError, "B is not fill-completed");
1809 
1810  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1811 
1812  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
1813 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1814  helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
1815 #else
1816  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
1817 #endif
1818  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
1819 #ifdef HAVE_XPETRA_TPETRA
1820 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1821  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1822  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long, EpetraNode> ETI enabled."));
1823 #else
1824  if (helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
1825  // All matrices are Crs
1826  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = helpers::Op2TpetraCrs(A);
1827  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = helpers::Op2TpetraCrs(B);
1828  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpC = helpers::Op2NonConstTpetraCrs(C);
1829 
1830  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
1831  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
1832  Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1833  } else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
1834  // All matrices are BlockCrs (except maybe Ac)
1835  // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
1836  // we'll need to think about refactoring BlockCrs so we can do something smarter here.
1837 
1838  if (!A.getRowMap()->getComm()->getRank())
1839  std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
1840 
1841  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
1842  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(B);
1843  using CRS = Tpetra::CrsMatrix<SC, LO, GO, NO>;
1844  RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
1845  RCP<const CRS> Bcrs = Tpetra::convertToCrsMatrix(tpB);
1846 
1847  // We need the global constants to do the copy back to BlockCrs
1848  RCP<ParameterList> new_params;
1849  if (!params.is_null()) {
1850  new_params = rcp(new Teuchos::ParameterList(*params));
1851  new_params->set("compute global constants", true);
1852  }
1853 
1854  // FIXME: The lines below only works because we're assuming Ac is Point
1855  RCP<CRS> tempAc = Teuchos::rcp(new CRS(Acrs->getRowMap(), 0));
1856  Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
1857 
1858  // Temporary output matrix
1859  RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*tempAc, A.GetStorageBlockSize());
1860  RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
1861  RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > Ac_p = Ac_x;
1862 
1863  // We can now cheat and replace the innards of Ac
1864  RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(Teuchos::rcpFromRef(C));
1865  Ac_w->replaceCrsMatrix(Ac_p);
1866  } else {
1867  // Mix and match
1868  TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
1869  }
1870 
1871 #endif
1872 #else
1873  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
1874 #endif
1875  }
1876 
1877  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1878  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
1879  fillParams->set("Optimize Storage", doOptimizeStorage);
1880  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
1881  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
1882  fillParams);
1883  }
1884 
1885  // transfer striding information
1886  RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1887  RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1888  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
1889  } // end Multiply
1890 
1913  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1914  const Matrix& B, bool transposeB,
1915  RCP<Matrix> C_in,
1916  Teuchos::FancyOStream& fos,
1917  bool doFillComplete = true,
1918  bool doOptimizeStorage = true,
1919  const std::string& label = std::string(),
1920  const RCP<ParameterList>& params = null) {
1921  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1922  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1923 
1924  // Default case: Xpetra Multiply
1925  RCP<Matrix> C = C_in;
1926 
1927  if (C == Teuchos::null) {
1928  double nnzPerRow = Teuchos::as<double>(0);
1929 
1930 #if 0
1931  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
1932  // For now, follow what ML and Epetra do.
1933  GO numRowsA = A.getGlobalNumRows();
1934  GO numRowsB = B.getGlobalNumRows();
1935  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
1936  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
1937  nnzPerRow *= nnzPerRow;
1938  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1939  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1940  if (totalNnz < minNnz)
1941  totalNnz = minNnz;
1942  nnzPerRow = totalNnz / A.getGlobalNumRows();
1943 
1944  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1945  }
1946 #endif
1947  if (transposeA)
1948  C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1949  else
1950  C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1951 
1952  } else {
1953  C->resumeFill(); // why this is not done inside of Tpetra MxM?
1954  fos << "Reuse C pattern" << std::endl;
1955  }
1956 
1957  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1958 
1959  return C;
1960  }
1961 
1972  static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1973  const Matrix& B, bool transposeB,
1974  Teuchos::FancyOStream& fos,
1975  bool callFillCompleteOnResult = true,
1976  bool doOptimizeStorage = true,
1977  const std::string& label = std::string(),
1978  const RCP<ParameterList>& params = null) {
1979  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1980  }
1981 
1982 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1983  // Michael Gee's MLMultiply
1984  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& /* epA */,
1985  const Epetra_CrsMatrix& /* epB */,
1986  Teuchos::FancyOStream& /* fos */) {
1987  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
1988  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1989  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1990  }
1991 #endif // ifdef HAVE_XPETRA_EPETRAEXT
1992 
2003  static RCP<BlockedCrsMatrix> TwoMatrixMultiplyBlock(const BlockedCrsMatrix& A, bool transposeA,
2004  const BlockedCrsMatrix& B, bool transposeB,
2005  Teuchos::FancyOStream& fos,
2006  bool doFillComplete = true,
2007  bool doOptimizeStorage = true) {
2008  TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
2009  "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
2010 
2011  // Preconditions
2012  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
2013  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
2014 
2015  RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
2016  RCP<const MapExtractor> domapextractor = B.getDomainMapExtractor();
2017 
2018  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2019 
2020  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
2021  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
2022  RCP<Matrix> Cij;
2023 
2024  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
2025  RCP<Matrix> crmat1 = A.getMatrix(i, l);
2026  RCP<Matrix> crmat2 = B.getMatrix(l, j);
2027 
2028  if (crmat1.is_null() || crmat2.is_null())
2029  continue;
2030 
2031  // try unwrapping 1x1 blocked matrix
2032  {
2033  auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
2034  auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
2035 
2036  if (unwrap1.is_null() != unwrap2.is_null()) {
2037  if (unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
2038  crmat1 = unwrap1->getCrsMatrix();
2039  if (unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
2040  crmat2 = unwrap2->getCrsMatrix();
2041  }
2042  }
2043 
2044  RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
2045  RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
2046  TEUCHOS_TEST_FOR_EXCEPTION(crop1.is_null() != crop2.is_null(), Xpetra::Exceptions::RuntimeError,
2047  "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
2048 
2049  // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
2050  if (!crop1.is_null())
2051  Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
2052  if (!crop2.is_null())
2053  Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
2054 
2055  TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
2056  "crmat1 does not have global constants");
2057  TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
2058  "crmat2 does not have global constants");
2059 
2060  if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
2061  continue;
2062 
2063  // temporary matrix containing result of local block multiplication
2064  RCP<Matrix> temp = Teuchos::null;
2065 
2066  if (crop1 != Teuchos::null && crop2 != Teuchos::null)
2067  temp = Multiply(*crop1, false, *crop2, false, fos);
2068  else {
2069  RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
2070  RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
2071  TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
2072  TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
2073  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)");
2074  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)");
2075 
2076  // recursive multiplication call
2077  temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
2078 
2079  RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
2080  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)");
2081  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)");
2082  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)");
2083  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)");
2084  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)");
2085  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)");
2086  }
2087 
2088  TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
2089 
2090  RCP<Matrix> addRes = null;
2091  if (Cij.is_null())
2092  Cij = temp;
2093  else {
2094  MatrixMatrix::TwoMatrixAdd(*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
2095  Cij = addRes;
2096  }
2097  }
2098 
2099  if (!Cij.is_null()) {
2100  if (Cij->isFillComplete())
2101  Cij->resumeFill();
2102  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
2103  C->setMatrix(i, j, Cij);
2104  } else {
2105  C->setMatrix(i, j, Teuchos::null);
2106  }
2107  }
2108  }
2109 
2110  if (doFillComplete)
2111  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
2112 
2113  return C;
2114  } // TwoMatrixMultiplyBlock
2115 
2128  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
2129  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
2130  "TwoMatrixAdd: matrix row maps are not the same.");
2131 
2132  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
2133 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
2134  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC, LO, GO, NO>::Op2EpetraCrs(A);
2135  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstEpetraCrs(B);
2136 
2137  // FIXME is there a bug if beta=0?
2138  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
2139 
2140  if (rv != 0)
2141  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
2142  std::ostringstream buf;
2143 #else
2144  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
2145 #endif
2146  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
2147 #ifdef HAVE_XPETRA_TPETRA
2148 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
2149  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
2150  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
2151 #else
2152  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
2153  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(B);
2154 
2155  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
2156 #endif
2157 #else
2158  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
2159 #endif
2160  }
2161  } // MatrixMatrix::TwoMatrixAdd()
2162 
2177  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
2178  const Matrix& B, bool transposeB, const SC& beta,
2179  RCP<Matrix>& C, Teuchos::FancyOStream& fos, bool AHasFixedNnzPerRow = false) {
2180  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
2181  RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
2182  RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
2183  RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
2184 
2185  if (rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
2186  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRowMap()->isSameAs(*(B.getRowMap()))), Exceptions::Incompatible,
2187  "TwoMatrixAdd: matrix row maps are not the same.");
2188  auto lib = A.getRowMap()->lib();
2189  if (lib == Xpetra::UseEpetra) {
2190 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
2191  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC, LO, GO, NO>::Op2EpetraCrs(A);
2192  const Epetra_CrsMatrix& epB = Xpetra::Helpers<SC, LO, GO, NO>::Op2EpetraCrs(B);
2193  if (C.is_null()) {
2194  size_t maxNzInA = 0;
2195  size_t maxNzInB = 0;
2196  size_t numLocalRows = 0;
2197  if (A.isFillComplete() && B.isFillComplete()) {
2198  maxNzInA = A.getLocalMaxNumRowEntries();
2199  maxNzInB = B.getLocalMaxNumRowEntries();
2200  numLocalRows = A.getLocalNumRows();
2201  }
2202 
2203  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
2204  // first check if either A or B has at most 1 nonzero per row
2205  // the case of both having at most 1 nz per row is handled by the ``else''
2206  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
2207 
2208  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
2209  for (size_t i = 0; i < numLocalRows; ++i)
2210  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
2211 
2212  } else {
2213  for (size_t i = 0; i < numLocalRows; ++i)
2214  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
2215  }
2216 
2217  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
2218  << ", using static profiling" << std::endl;
2219  C = rcp(new Xpetra::CrsMatrixWrap<SC, LO, GO, NO>(A.getRowMap(), exactNnzPerRow));
2220 
2221  } else {
2222  // general case
2223  LO maxPossibleNNZ = A.getLocalMaxNumRowEntries() + B.getLocalMaxNumRowEntries();
2224  C = rcp(new Xpetra::CrsMatrixWrap<SC, LO, GO, NO>(A.getRowMap(), maxPossibleNNZ));
2225  }
2226  if (transposeB)
2227  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
2228  }
2229  RCP<Epetra_CrsMatrix> epC = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstEpetraCrs(C);
2230  Epetra_CrsMatrix* ref2epC = &*epC; // to avoid a compiler error...
2231 
2232  // FIXME is there a bug if beta=0?
2233  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
2234 
2235  if (rv != 0)
2236  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
2237 #else
2238  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
2239 #endif
2240  } else if (lib == Xpetra::UseTpetra) {
2241 #ifdef HAVE_XPETRA_TPETRA
2242 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
2243  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
2244  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
2245 #else
2246  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
2247  using tcrs_matrix_type = Tpetra::CrsMatrix<SC, LO, GO, NO>;
2248  const tcrs_matrix_type& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
2249  const tcrs_matrix_type& tpB = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(B);
2250  C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
2251 #endif
2252 #else
2253  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
2254 #endif
2255  }
2256 
2258  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
2259  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
2261  }
2262  // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
2263  else if (rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
2264  RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
2265  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2266 
2267  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2268  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2269 
2270  size_t i = 0;
2271  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2272  RCP<Matrix> Cij = Teuchos::null;
2273  if (rcpA != Teuchos::null &&
2274  rcpBopB->getMatrix(i, j) != Teuchos::null) {
2275  // recursive call
2276  TwoMatrixAdd(*rcpA, transposeA, alpha,
2277  *(rcpBopB->getMatrix(i, j)), transposeB, beta,
2278  Cij, fos, AHasFixedNnzPerRow);
2279  } else if (rcpA == Teuchos::null &&
2280  rcpBopB->getMatrix(i, j) != Teuchos::null) {
2281  Cij = rcpBopB->getMatrix(i, j);
2282  } else if (rcpA != Teuchos::null &&
2283  rcpBopB->getMatrix(i, j) == Teuchos::null) {
2284  Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
2285  } else {
2286  Cij = Teuchos::null;
2287  }
2288 
2289  if (!Cij.is_null()) {
2290  if (Cij->isFillComplete())
2291  Cij->resumeFill();
2292  Cij->fillComplete();
2293  bC->setMatrix(i, j, Cij);
2294  } else {
2295  bC->setMatrix(i, j, Teuchos::null);
2296  }
2297  } // loop over columns j
2298  }
2299  // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
2300  else if (rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
2301  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2302  RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
2303 
2304  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2305  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2306 
2307  size_t j = 0;
2308  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2309  RCP<Matrix> Cij = Teuchos::null;
2310  if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
2311  rcpB != Teuchos::null) {
2312  // recursive call
2313  TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
2314  *rcpB, transposeB, beta,
2315  Cij, fos, AHasFixedNnzPerRow);
2316  } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
2317  rcpB != Teuchos::null) {
2318  Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
2319  } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
2320  rcpB == Teuchos::null) {
2321  Cij = rcpBopA->getMatrix(i, j);
2322  } else {
2323  Cij = Teuchos::null;
2324  }
2325 
2326  if (!Cij.is_null()) {
2327  if (Cij->isFillComplete())
2328  Cij->resumeFill();
2329  Cij->fillComplete();
2330  bC->setMatrix(i, j, Cij);
2331  } else {
2332  bC->setMatrix(i, j, Teuchos::null);
2333  }
2334  } // loop over rows i
2335  } else {
2336  // This is the version for blocked matrices
2337 
2338  // check the compatibility of the blocked operators
2339  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null() == true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2340  TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null() == true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2341  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)");
2342  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)");
2343  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)");
2344  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)");
2345  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)");
2346  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)");
2347 
2348  RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2349  RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2350 
2351  C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2352  RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2353 
2354  for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2355  for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2356 
2357  RCP<Matrix> Cij = Teuchos::null;
2358  if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
2359  rcpBopB->getMatrix(i, j) != Teuchos::null) {
2360  // recursive call
2361 
2362  TwoMatrixAdd(*(rcpBopA->getMatrix(i, j)), transposeA, alpha,
2363  *(rcpBopB->getMatrix(i, j)), transposeB, beta,
2364  Cij, fos, AHasFixedNnzPerRow);
2365  } else if (rcpBopA->getMatrix(i, j) == Teuchos::null &&
2366  rcpBopB->getMatrix(i, j) != Teuchos::null) {
2367  Cij = rcpBopB->getMatrix(i, j);
2368  } else if (rcpBopA->getMatrix(i, j) != Teuchos::null &&
2369  rcpBopB->getMatrix(i, j) == Teuchos::null) {
2370  Cij = rcpBopA->getMatrix(i, j);
2371  } else {
2372  Cij = Teuchos::null;
2373  }
2374 
2375  if (!Cij.is_null()) {
2376  if (Cij->isFillComplete())
2377  Cij->resumeFill();
2378  // Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
2379  Cij->fillComplete();
2380  bC->setMatrix(i, j, Cij);
2381  } else {
2382  bC->setMatrix(i, j, Teuchos::null);
2383  }
2384  } // loop over columns j
2385  } // loop over rows i
2386 
2387  } // end blocked recursive algorithm
2388  } // MatrixMatrix::TwoMatrixAdd()
2389 }; // end specialization on GO=long long and NO=EpetraNode
2390 
2391 #endif // HAVE_XPETRA_EPETRA
2392 
2393 } // end namespace Xpetra
2394 
2395 #define XPETRA_MATRIXMATRIX_SHORT
2396 
2397 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_ */
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition: Xpetra_CrsMatrixWrap_def.hpp:470
Xpetra::Helpers::isTpetraBlockCrs
static bool isTpetraBlockCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:295
Xpetra::Helpers::Op2NonConstTpetraCrs
static Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:200
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:812
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:1545
Xpetra::Helpers::Op2TpetraCrs
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:185
Xpetra::Helpers::Op2TpetraCrs
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:162
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:321
Xpetra::Helpers::isTpetraCrs
static bool isTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:215
Xpetra::BlockedCrsMatrix::getRangeMap
const RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1099
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:1913
Xpetra::BlockedCrsMatrix::isFillComplete
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Definition: Xpetra_BlockedCrsMatrix.hpp:714
Xpetra::Helpers::Op2NonConstEpetraCrs
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:144
Xpetra_UseShortNames.hpp
Xpetra::Helpers
Xpetra utility class containing transformation routines between Xpetra::Matrix and Epetra/Tpetra obje...
Definition: Xpetra_MatrixMatrix.hpp:95
Xpetra_EpetraCrsMatrix_fwd.hpp
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory.hpp:235
Xpetra::Matrix::resumeFill
virtual void resumeFill(const RCP< ParameterList > &params=null)=0
Xpetra::Helpers::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:174
Xpetra::Matrix::GetStorageBlockSize
virtual LocalOrdinal GetStorageBlockSize() const =0
Returns the block size of the storage mechanism, which is usually 1, except for Tpetra::BlockCrsMatri...
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:101
Xpetra::MatrixMatrix< double, int, long long, EpetraNode >::TwoMatrixMultiplyBlock
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply &quot;in-place&quot;.
Definition: Xpetra_MatrixMatrix.hpp:2003
Xpetra::BlockedCrsMatrix::getCrsMatrix
Teuchos::RCP< Matrix > getCrsMatrix() const
return unwrap 1x1 blocked operators
Definition: Xpetra_BlockedCrsMatrix.hpp:1350
Xpetra_Matrix.hpp
Xpetra::Helpers::isTpetraBlockCrs
static bool isTpetraBlockCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:284
Xpetra::Helpers::Op2EpetraCrs
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:128
Xpetra::Helpers::Op2TpetraBlockCrs
static const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > & Op2TpetraBlockCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:260
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:2177
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:1984
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:629
Xpetra::Matrix::IsView
bool IsView(const viewLabel_t viewLabel) const
Definition: Xpetra_Matrix.hpp:202
Xpetra::BlockedCrsMatrix::getMatrix
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
Definition: Xpetra_BlockedCrsMatrix.hpp:1380
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Xpetra::BlockedCrsMatrix::Cols
virtual size_t Cols() const
number of column blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1344
Xpetra::Matrix::CreateView
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
Definition: Xpetra_Matrix.hpp:126
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:1008
Xpetra::Helpers::Op2NonConstTpetraBlockCrs
static RCP< Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraBlockCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:250
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:1127
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:572
Xpetra::Operator::getRangeMap
virtual const Teuchos::RCP< const map_type > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
Xpetra::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:1794
Xpetra::MatrixMatrix::TwoMatrixMultiplyBlock
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply &quot;in-place&quot;.
Definition: Xpetra_MatrixMatrix.hpp:655
Xpetra::Matrix::getLocalNumRows
virtual size_t getLocalNumRows() const =0
Returns the number of matrix rows owned on the calling node.
Xpetra::Helpers::isTpetraCrs
static bool isTpetraCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:226
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:780
Xpetra::Helpers::tcrs_matrix_type
Tpetra::CrsMatrix< SC, LO, GO, NO > tcrs_matrix_type
Definition: Xpetra_MatrixMatrix.hpp:320
Xpetra::Helpers::epetraExtMult
static void epetraExtMult(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool fillCompleteResult)
Definition: Xpetra_MatrixMatrix.hpp:379
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
Xpetra::Helpers::Op2EpetraCrs
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:100
Xpetra::BlockedCrsMatrix::getDomainMapExtractor
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1123
Xpetra::EpetraNode
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
Definition: Xpetra_Map_decl.hpp:71
Xpetra::BlockedCrsMatrix::getDomainMap
const RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
Definition: Xpetra_BlockedCrsMatrix.hpp:1069
Xpetra::CrsMatrixWrap::replaceCrsMatrix
void replaceCrsMatrix(RCP< CrsMatrix > &M)
Expert only.
Definition: Xpetra_CrsMatrixWrap_def.hpp:495
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:459
Xpetra::Helpers::Op2NonConstEpetraCrs
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:114
Xpetra_ConfigDefs.hpp
Xpetra::Helpers::Op2TpetraBlockCrs
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:240
Xpetra::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: Xpetra_Exceptions.hpp:108
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:1496
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::Rows
virtual size_t Rows() const
number of row blocks
Definition: Xpetra_BlockedCrsMatrix.hpp:1338
Xpetra::Operator::getDomainMap
virtual const Teuchos::RCP< const map_type > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra::Matrix::getLocalMaxNumRowEntries
virtual size_t getLocalMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
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:1223
Xpetra::Helpers::Op2NonConstTpetraBlockCrs
static Tpetra::BlockCrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraBlockCrs(const Matrix &Op)
Definition: Xpetra_MatrixMatrix.hpp:272
Xpetra::MatrixMatrix< double, int, int, EpetraNode >::TwoMatrixMultiplyBlock
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply &quot;in-place&quot;.
Definition: Xpetra_MatrixMatrix.hpp:1371
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:315
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:1211
Xpetra::MatrixMatrix::MLTwoMatrixMultiply
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
Definition: Xpetra_MatrixMatrix.hpp:637
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:2128
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_fwd.hpp:51
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:1972
Xpetra::BlockedCrsMatrix::getRangeMapExtractor
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Definition: Xpetra_BlockedCrsMatrix.hpp:1117