Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Xpetra_Helpers_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Xpetra: A linear algebra interface package
4 //
5 // Copyright 2012 NTESS and the Xpetra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef XPETRA_HELPERS_DEF_HPP
11 #define XPETRA_HELPERS_DEF_HPP
12 
13 #include "Xpetra_Helpers_decl.hpp"
14 
15 namespace Xpetra {
16 
17 #ifdef HAVE_XPETRA_EPETRA
18 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
19 RCP<const Epetra_CrsMatrix> Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2EpetraCrs(RCP<Matrix> Op) {
20  // Get the underlying Epetra Mtx
21  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
22  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast,
23  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
24 
25  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
26  const RCP<const EpetraCrsMatrixT<GO, NO>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO, NO>>(tmp_CrsMtx);
27  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast,
28  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
29 
30  return tmp_ECrsMtx->getEpetra_CrsMatrix();
31 }
32 
33 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
34 RCP<Epetra_CrsMatrix> Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(RCP<Matrix> Op) {
35  RCP<Epetra_CrsMatrix> A;
36  // Get the underlying Epetra Mtx
37  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
38  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast,
39  "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
40 
41  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
42  const RCP<const EpetraCrsMatrixT<GO, NO>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO, NO>>(tmp_CrsMtx);
43  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
44 
45  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
46 }
47 
48 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
49 const Epetra_CrsMatrix& Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2EpetraCrs(const Matrix& Op) {
50  // Get the underlying Epetra Mtx
51  try {
52  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
53  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
54  const RCP<const EpetraCrsMatrixT<GO, NO>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO, NO>>(tmp_CrsMtx);
55  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast,
56  "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
57 
58  return *tmp_ECrsMtx->getEpetra_CrsMatrix();
59 
60  } catch (...) {
61  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
62  }
63 }
64 
65 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
66 Epetra_CrsMatrix& Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(const Matrix& Op) {
67  RCP<Epetra_CrsMatrix> A;
68  // Get the underlying Epetra Mtx
69  try {
70  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
71  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
72  const RCP<const EpetraCrsMatrixT<GO, NO>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO, NO>>(tmp_CrsMtx);
73  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
74 
75  return *Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
76 
77  } catch (...) {
78  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
79  }
80 }
81 #endif // HAVE_XPETRA_EPETRA
82 
83 #ifdef HAVE_XPETRA_TPETRA
84 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
85 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(RCP<Matrix> Op) {
86  // Get the underlying Tpetra Mtx
87  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
88  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
89 
90  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
91  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>(tmp_CrsMtx);
92  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
93 
94  return tmp_ECrsMtx->getTpetra_CrsMatrix();
95 }
96 
97 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
98 RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraCrs(RCP<Matrix> Op) {
99  // Get the underlying Tpetra Mtx
100  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
101  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
102  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
103  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>(tmp_CrsMtx);
104  TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
105 
106  return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
107 }
108 
109 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
110 const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(const Matrix& Op) {
111  // Get the underlying Tpetra Mtx
112  try {
113  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
114  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
115  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>& tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>(tmp_CrsMtx);
116  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
117 
118  return *tmp_TCrsMtx->getTpetra_CrsMatrix();
119 
120  } catch (...) {
121  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
122  }
123 }
124 
125 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
126 Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraCrs(const Matrix& Op) {
127  // Get the underlying Tpetra Mtx
128  try {
129  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
130  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
131  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>& tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>(tmp_CrsMtx);
132  TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
133 
134  return *Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tmp_TCrsMtx->getTpetra_CrsMatrix());
135 
136  } catch (...) {
137  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
138  }
139 }
140 
141 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142 bool Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetraCrs(RCP<Matrix> Op) {
143  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
144  if (crsOp == Teuchos::null) return false;
145  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
146  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>(tmp_CrsMtx);
147  if (tmp_ECrsMtx == Teuchos::null)
148  return false;
149  else
150  return true;
151 }
152 
153 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
154 bool Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetraCrs(const Matrix& Op) {
155  try {
156  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
157  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
158  const RCP<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>& tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>>(tmp_CrsMtx);
159  if (tmp_ECrsMtx == Teuchos::null)
160  return false;
161  else
162  return true;
163  } catch (...) {
164  return false;
165  }
166 }
167 
168 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
169 RCP<const Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraBlockCrs(RCP<Matrix> Op) {
170  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
171  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
172 
173  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
174  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
175  TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
176  return tmp_BlockCrs->getTpetra_BlockCrsMatrix();
177 }
178 
179 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
180 RCP<Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraBlockCrs(RCP<Matrix> Op) {
181  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
182  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
183 
184  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
185  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
186  TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
187  return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
188 }
189 
190 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
191 const Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraBlockCrs(const Matrix& Op) {
192  try {
193  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
194  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
195  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
196  TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
197  return *tmp_BlockCrs->getTpetra_BlockCrsMatrix();
198  } catch (...) {
199  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
200  }
201 }
202 
203 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
204 Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraBlockCrs(const Matrix& Op) {
205  try {
206  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
207  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
208  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
209  TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
210  return *tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
211  } catch (...) {
212  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
213  }
214 }
215 
216 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
217 bool Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetraBlockCrs(RCP<Matrix> Op) {
218  RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
219  if (crsOp == Teuchos::null) return false;
220  RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
221  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
222  if (tmp_BlockCrs == Teuchos::null)
223  return false;
224  else
225  return true;
226 }
227 
228 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
229 bool Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetraBlockCrs(const Matrix& Op) {
230  try {
231  const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
232  RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
233  RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
234  if (tmp_BlockCrs == Teuchos::null)
235  return false;
236  else
237  return true;
238  } catch (...) {
239  return false;
240  }
241 }
242 #else // HAVE_XPETRA_TPETRA
243 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
244 bool Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetraCrs(const Matrix& Op) {
245  return false;
246 }
247 
248 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
249 bool Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetraBlockCrs(const Matrix& Op) {
250  return false;
251 }
252 
253 #endif // HAVE_XPETRA_TPETRA
254 
255 #ifdef HAVE_XPETRA_TPETRA
256 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
257 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::tpetraAdd(
258  const tcrs_matrix_type& A, bool transposeA, const typename tcrs_matrix_type::scalar_type alpha,
259  const tcrs_matrix_type& B, bool transposeB, const typename tcrs_matrix_type::scalar_type beta) {
260  using Teuchos::Array;
261  using Teuchos::ParameterList;
262  using Teuchos::RCP;
263  using Teuchos::rcp;
264  using Teuchos::rcp_implicit_cast;
265  using Teuchos::rcpFromRef;
266  using XTCrsType = Xpetra::TpetraCrsMatrix<SC, LO, GO, NO>;
267  using CrsType = Xpetra::CrsMatrix<SC, LO, GO, NO>;
268  using CrsWrap = Xpetra::CrsMatrixWrap<SC, LO, GO, NO>;
269  using transposer_type = Tpetra::RowMatrixTransposer<SC, LO, GO, NO>;
270  using import_type = Tpetra::Import<LO, GO, NO>;
271  RCP<const tcrs_matrix_type> Aprime = rcpFromRef(A);
272  if (transposeA)
273  Aprime = transposer_type(Aprime).createTranspose();
274  // Decide whether the fast code path can be taken.
275  if (A.isFillComplete() && B.isFillComplete()) {
276  RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), 0));
277  RCP<ParameterList> addParams = rcp(new ParameterList);
278  addParams->set("Call fillComplete", false);
279  // passing A after B means C will have the same domain/range map as A (or A^T if transposeA)
280  Tpetra::MatrixMatrix::add<SC, LO, GO, NO>(beta, transposeB, B, alpha, false, *Aprime, *C, Teuchos::null, Teuchos::null, addParams);
281  return rcp_implicit_cast<Matrix>(rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C))))));
282  } else {
283  // Slow case - one or both operands are non-fill complete.
284  // TODO: deprecate this.
285  // Need to compute the explicit transpose before add if transposeA and/or transposeB.
286  // This is to count the number of entries to allocate per row in the final sum.
287  RCP<const tcrs_matrix_type> Bprime = rcpFromRef(B);
288  if (transposeB)
289  Bprime = transposer_type(Bprime).createTranspose();
290  // Use Aprime's row map as C's row map.
291  if (!(Aprime->getRowMap()->isSameAs(*(Bprime->getRowMap())))) {
292  auto import = rcp(new import_type(Bprime->getRowMap(), Aprime->getRowMap()));
293  Bprime = Tpetra::importAndFillCompleteCrsMatrix<tcrs_matrix_type>(Bprime, *import, Aprime->getDomainMap(), Aprime->getRangeMap());
294  }
295  // Count the entries to allocate for C in each local row.
296  LO numLocalRows = Aprime->getLocalNumRows();
297  Array<size_t> allocPerRow(numLocalRows);
298  // 0 is always the minimum LID
299  for (LO i = 0; i < numLocalRows; i++) {
300  allocPerRow[i] = Aprime->getNumEntriesInLocalRow(i) + Bprime->getNumEntriesInLocalRow(i);
301  }
302  // Construct C
303  RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), allocPerRow()));
304  // Compute the sum in C (already took care of transposes)
305  Tpetra::MatrixMatrix::Add<SC, LO, GO, NO>(
306  *Aprime, false, alpha,
307  *Bprime, false, beta,
308  C);
309  return rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C)))));
310  }
311 }
312 #endif
313 
314 #ifdef HAVE_XPETRA_EPETRAEXT
315 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
316 void Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>::epetraExtMult(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Matrix& C, bool fillCompleteResult) {
317  Epetra_CrsMatrix& epA = Op2NonConstEpetraCrs(A);
318  Epetra_CrsMatrix& epB = Op2NonConstEpetraCrs(B);
319  Epetra_CrsMatrix& epC = Op2NonConstEpetraCrs(C);
320  // Want a static profile (possibly fill complete) matrix as the result.
321  // But, EpetraExt Multiply needs C to be dynamic profile, so compute the product in a temporary DynamicProfile matrix.
322  const Epetra_Map& Crowmap = epC.RowMap();
323  int errCode = 0;
324  Epetra_CrsMatrix Ctemp(::Copy, Crowmap, 0, false);
325  if (fillCompleteResult) {
326  errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, true);
327  if (!errCode) {
328  epC = Ctemp;
329  C.fillComplete();
330  }
331  } else {
332  errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, false);
333  if (!errCode) {
334  int numLocalRows = Crowmap.NumMyElements();
335  long long* globalElementList = nullptr;
336  Crowmap.MyGlobalElementsPtr(globalElementList);
337  Teuchos::Array<int> entriesPerRow(numLocalRows, 0);
338  for (int i = 0; i < numLocalRows; i++) {
339  entriesPerRow[i] = Ctemp.NumGlobalEntries(globalElementList[i]);
340  }
341  // know how many entries to allocate in epC (which must be static profile)
342  epC = Epetra_CrsMatrix(::Copy, Crowmap, entriesPerRow.data(), true);
343  for (int i = 0; i < numLocalRows; i++) {
344  int gid = globalElementList[i];
345  int numEntries;
346  double* values;
347  int* inds;
348  Ctemp.ExtractGlobalRowView(gid, numEntries, values, inds);
349  epC.InsertGlobalValues(gid, numEntries, values, inds);
350  }
351  }
352  }
353  if (errCode) {
354  std::ostringstream buf;
355  buf << errCode;
356  std::string msg = "EpetraExt::MatrixMatrix::Multiply returned nonzero error code " + buf.str();
357  throw(Exceptions::RuntimeError(msg));
358  }
359 }
360 #endif
361 
362 } // namespace Xpetra
363 
364 #endif
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Definition: Xpetra_CrsMatrixWrap_def.hpp:433
Xpetra_Helpers_decl.hpp
Xpetra::Helpers::isTpetraBlockCrs
static bool isTpetraBlockCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:217
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:65
Xpetra::Helpers::Op2TpetraCrs
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:85
Xpetra::Helpers::isTpetraCrs
static bool isTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:142
Xpetra::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: Xpetra_Exceptions.hpp:50
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::Helpers::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:98
Xpetra::Helpers::tcrs_matrix_type
Tpetra::CrsMatrix< SC, LO, GO, NO > tcrs_matrix_type
Definition: Xpetra_Helpers_decl.hpp:102
Xpetra::Helpers::Op2NonConstTpetraBlockCrs
static RCP< Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraBlockCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:180
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_Helpers_def.hpp:257
Xpetra::Helpers::Op2EpetraCrs
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:19
Xpetra::Helpers::epetraExtMult
static void epetraExtMult(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool fillCompleteResult)
Definition: Xpetra_Helpers_def.hpp:316
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_decl.hpp:46
Xpetra::Helpers::Op2TpetraBlockCrs
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:169
Xpetra::Helpers::Op2NonConstEpetraCrs
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_Helpers_def.hpp:34
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_decl.hpp:60