Xpetra Version of the Day
Loading...
Searching...
No Matches
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
54#include "Xpetra_CrsMatrixWrap.hpp"
55#include "Xpetra_MapExtractor.hpp"
56#include "Xpetra_Map.hpp"
58#include "Xpetra_Matrix.hpp"
59#include "Xpetra_StridedMapFactory.hpp"
60#include "Xpetra_StridedMap.hpp"
61
62#ifdef HAVE_XPETRA_EPETRA
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
83namespace Xpetra {
84
91 template <class Scalar,
92 class LocalOrdinal = int,
93 class GlobalOrdinal = LocalOrdinal,
95 class Helpers {
97
98 public:
99
100#ifdef HAVE_XPETRA_EPETRA
102 // Get the underlying Epetra Mtx
103 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
105 "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
106
107 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
108 const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
109 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast,
110 "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
111
112 return tmp_ECrsMtx->getEpetra_CrsMatrix();
113 }
114
117 // Get the underlying Epetra Mtx
118 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
119 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast,
120 "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
121
122 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
123 const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
124 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
125
126 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
127 }
128
129 static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) {
130 // Get the underlying Epetra Mtx
131 try {
132 const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
133 RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
134 const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
135 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast,
136 "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
137
138 return *tmp_ECrsMtx->getEpetra_CrsMatrix();
139
140 } catch(...) {
141 throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
142 }
143 }
144
147 // Get the underlying Epetra Mtx
148 try {
149 const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
150 RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
151 const RCP<const EpetraCrsMatrixT<GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const EpetraCrsMatrixT<GO,NO> >(tmp_CrsMtx);
152 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
153
154 return *Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
155
156 } catch(...) {
157 throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
158 }
159 }
160#endif // HAVE_XPETRA_EPETRA
161
162#ifdef HAVE_XPETRA_TPETRA
164 // Get the underlying Tpetra Mtx
165 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
166 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
167
168 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
169 const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
170 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
171
172 return tmp_ECrsMtx->getTpetra_CrsMatrix();
173 }
174
176 // Get the underlying Tpetra Mtx
177 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
178 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
179 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
180 const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
181 TEUCHOS_TEST_FOR_EXCEPTION(tmp_ECrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
182
183 return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
184 }
185
186 static const Tpetra::CrsMatrix<SC,LO,GO,NO> & Op2TpetraCrs(const Matrix& Op) {
187 // Get the underlying Tpetra Mtx
188 try{
189 const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
190 RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
191 const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
192 TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
193
194 return *tmp_TCrsMtx->getTpetra_CrsMatrix();
195
196 } catch (...) {
197 throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
198 }
199 }
200
201 static Tpetra::CrsMatrix<SC,LO,GO,NO> & Op2NonConstTpetraCrs(const Matrix& Op) {
202 // Get the underlying Tpetra Mtx
203 try{
204 const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap& >(Op);
205 RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
206 const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_TCrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
207 TEUCHOS_TEST_FOR_EXCEPTION(tmp_TCrsMtx == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
208
209 return *Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_TCrsMtx->getTpetra_CrsMatrix());
210
211 } catch (...) {
212 throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
213 }
214 }
215
216 static bool isTpetraCrs(RCP<Matrix> Op) {
217 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
218 if(crsOp == Teuchos::null) return false;
219 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
220 const RCP<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> >(tmp_CrsMtx);
221 if(tmp_ECrsMtx == Teuchos::null) return false;
222 else return true;
223 }
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) return false;
232 else return true;
233 } catch(...) {
234 return false;
235 }
236 }
237
239 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
240 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
241
242 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
243 RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs= Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
244 TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
245 return tmp_BlockCrs->getTpetra_BlockCrsMatrix();
246 }
247
249 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
250 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
251
252 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
253 RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs= Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
254 TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
255 return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
256 }
257
258 static const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & Op2TpetraBlockCrs(const Matrix& Op) {
259 try {
260 const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
261 RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
262 RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs= Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
263 TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
264 return *tmp_BlockCrs->getTpetra_BlockCrsMatrix();
265 } catch(...) {
266 throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
267 }
268 }
269
270 static Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & Op2NonConstTpetraBlockCrs(const Matrix& Op) {
271 try {
272 const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
273 RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
274 RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs= Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
275 TEUCHOS_TEST_FOR_EXCEPTION(tmp_BlockCrs == Teuchos::null, Xpetra::Exceptions::BadCast, "Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
276 return *tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
277 } catch(...) {
278 throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
279 }
280 }
281
283 RCP<const CrsMatrixWrap> crsOp = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(Op);
284 if(crsOp == Teuchos::null) return false;
285 RCP<const CrsMatrix> tmp_CrsMtx = crsOp->getCrsMatrix();
286 RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs= Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
287 if(tmp_BlockCrs == Teuchos::null) return false;
288 else return true;
289 }
290
291 static bool isTpetraBlockCrs(const Matrix& Op) {
292 try {
293 const CrsMatrixWrap& crsOp = dynamic_cast<const CrsMatrixWrap&>(Op);
294 RCP<const CrsMatrix> tmp_CrsMtx = crsOp.getCrsMatrix();
295 RCP<const TpetraBlockCrsMatrix> tmp_BlockCrs= Teuchos::rcp_dynamic_cast<const TpetraBlockCrsMatrix>(tmp_CrsMtx);
296 if(tmp_BlockCrs == Teuchos::null) return false;
297 else return true;
298 } catch(...) {
299 return false;
300 }
301 }
302#else // HAVE_XPETRA_TPETRA
303 static bool isTpetraCrs(const Matrix& Op) {
304 return false;
305 }
306
307 static bool isTpetraBlockCrs(const Matrix& Op) {
308 return false;
309 }
310
311#endif // HAVE_XPETRA_TPETRA
312
313#ifdef HAVE_XPETRA_TPETRA
314 using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
316 const tcrs_matrix_type& A, bool transposeA, const typename tcrs_matrix_type::scalar_type alpha,
317 const tcrs_matrix_type& B, bool transposeB, const typename tcrs_matrix_type::scalar_type beta)
318 {
319 using Teuchos::Array;
320 using Teuchos::RCP;
321 using Teuchos::rcp;
322 using Teuchos::rcp_implicit_cast;
323 using Teuchos::rcpFromRef;
325 using XTCrsType = Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>;
326 using CrsType = Xpetra::CrsMatrix<SC,LO,GO,NO>;
328 using transposer_type = Tpetra::RowMatrixTransposer<SC,LO,GO,NO>;
329 using import_type = Tpetra::Import<LO,GO,NO>;
330 RCP<const tcrs_matrix_type> Aprime = rcpFromRef(A);
331 if(transposeA)
332 Aprime = transposer_type(Aprime).createTranspose();
333 //Decide whether the fast code path can be taken.
334 if(A.isFillComplete() && B.isFillComplete())
335 {
336 RCP<tcrs_matrix_type> C = rcp(new tcrs_matrix_type(Aprime->getRowMap(), 0));
337 RCP<ParameterList> addParams = rcp(new ParameterList);
338 addParams->set("Call fillComplete", false);
339 //passing A after B means C will have the same domain/range map as A (or A^T if transposeA)
340 Tpetra::MatrixMatrix::add<SC,LO,GO,NO>(beta, transposeB, B, alpha, false, *Aprime, *C, Teuchos::null, Teuchos::null, addParams);
341 return rcp_implicit_cast<Matrix>(rcp(new CrsWrap(rcp_implicit_cast<CrsType>(rcp(new XTCrsType(C))))));
342 }
343 else
344 {
345 //Slow case - one or both operands are non-fill complete.
346 //TODO: deprecate this.
347 //Need to compute the explicit transpose before add if transposeA and/or transposeB.
348 //This is to count the number of entries to allocate per row in the final sum.
349 RCP<const tcrs_matrix_type> Bprime = rcpFromRef(B);
350 if(transposeB)
351 Bprime = transposer_type(Bprime).createTranspose();
352 //Use Aprime's row map as C's row map.
353 if(!(Aprime->getRowMap()->isSameAs(*(Bprime->getRowMap()))))
354 {
355 auto import = rcp(new import_type(Bprime->getRowMap(), Aprime->getRowMap()));
356 Bprime = Tpetra::importAndFillCompleteCrsMatrix<tcrs_matrix_type>(Bprime, *import, Aprime->getDomainMap(), Aprime->getRangeMap());
357 }
358 //Count the entries to allocate for C in each local row.
359 LO numLocalRows = Aprime->getLocalNumRows();
360 Array<size_t> allocPerRow(numLocalRows);
361 //0 is always the minimum LID
362 for(LO i = 0; i < numLocalRows; i++)
363 {
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 {
384 //Want a static profile (possibly fill complete) matrix as the result.
385 //But, EpetraExt Multiply needs C to be dynamic profile, so compute the product in a temporary DynamicProfile matrix.
386 const Epetra_Map& Crowmap = epC.RowMap();
387 int errCode = 0;
388 Epetra_CrsMatrix Ctemp(::Copy, Crowmap, 0, false);
389 if(fillCompleteResult) {
390 errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, true);
391 if(!errCode) {
392 epC = Ctemp;
393 C.fillComplete();
394 }
395 }
396 else {
397 errCode = EpetraExt::MatrixMatrix::Multiply(epA, transposeA, epB, transposeB, Ctemp, false);
398 if(!errCode) {
399 int numLocalRows = Crowmap.NumMyElements();
400 long long* globalElementList = nullptr;
401 Crowmap.MyGlobalElementsPtr(globalElementList);
402 Teuchos::Array<int> entriesPerRow(numLocalRows, 0);
403 for(int i = 0; i < numLocalRows; i++)
404 {
405 entriesPerRow[i] = Ctemp.NumGlobalEntries(globalElementList[i]);
406 }
407 //know how many entries to allocate in epC (which must be static profile)
408 epC = Epetra_CrsMatrix(::Copy, Crowmap, entriesPerRow.data(), true);
409 for(int i = 0; i < numLocalRows; i++)
410 {
411 int gid = globalElementList[i];
412 int numEntries;
413 double* values;
414 int* inds;
415 Ctemp.ExtractGlobalRowView(gid, numEntries, values, inds);
416 epC.InsertGlobalValues(gid, numEntries, values, inds);
417 }
418 }
419 }
420 if(errCode) {
421 std::ostringstream buf;
422 buf << errCode;
423 std::string msg = "EpetraExt::MatrixMatrix::Multiply returned nonzero error code " + buf.str();
424 throw(Exceptions::RuntimeError(msg));
425 }
426 }
427#endif
428 };
429
430 template <class Scalar,
431 class LocalOrdinal /*= int*/,
432 class GlobalOrdinal /*= LocalOrdinal*/,
433 class Node /*= KokkosClassic::DefaultNode::DefaultNodeType*/>
435#undef XPETRA_MATRIXMATRIX_SHORT
437
438 public:
439
464 static void Multiply(const Matrix& A, bool transposeA,
465 const Matrix& B, bool transposeB,
466 Matrix& C,
467 bool call_FillComplete_on_result = true,
468 bool doOptimizeStorage = true,
469 const std::string & label = std::string(),
470 const RCP<ParameterList>& params = null) {
471
472 TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
473 Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
474 TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
475 Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
476
479
480 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
481
482 if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
483#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
484 throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
485#else
486 throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
487
488#endif
489 } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
490#ifdef HAVE_XPETRA_TPETRA
491 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
492 if(helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
493 // All matrices are Crs
494 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
495 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
496 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
497
498 // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
499 // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
500 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
501 }
502 else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
503 // All matrices are BlockCrs (except maybe Ac)
504 // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
505 // we'll need to think about refactoring BlockCrs so we can do something smarter here.
506 if(!A.getRowMap()->getComm()->getRank())
507 std::cout<<"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
508
509 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(A);
510 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(B);
511 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
512 RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
513 RCP<const CRS> Bcrs = Tpetra::convertToCrsMatrix(tpB);
514
515 // We need the global constants to do the copy back to BlockCrs
516 RCP<ParameterList> new_params;
517 if(!params.is_null()) {
518 new_params = rcp(new Teuchos::ParameterList(*params));
519 new_params->set("compute global constants",true);
520 }
521
522 // FIXME: The lines below only works because we're assuming Ac is Point
523 RCP<CRS> tempAc = Teuchos::rcp(new CRS(Acrs->getRowMap(),0));
524 Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
525
526 // Temporary output matrix
527 RCP<Tpetra::BlockCrsMatrix<SC,LO,GO,NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*tempAc,A.GetStorageBlockSize());
530
531 // We can now cheat and replace the innards of Ac
532 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(Teuchos::rcpFromRef(C));
533 Ac_w->replaceCrsMatrix(Ac_p);
534 }
535 else {
536 // Mix and match
537 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
538 }
539#else
540 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
541#endif
542 }
543
544 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
546 fillParams->set("Optimize Storage", doOptimizeStorage);
547 C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
548 (transposeA) ? A.getDomainMap() : A.getRangeMap(),
549 fillParams);
550 }
551
552 // transfer striding information
553 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
554 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
555 C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
556 } // end Multiply
557
580 static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, RCP<Matrix> C_in,
582 bool doFillComplete = true,
583 bool doOptimizeStorage = true,
584 const std::string & label = std::string(),
585 const RCP<ParameterList>& params = null) {
586
589
590 // Default case: Xpetra Multiply
591 RCP<Matrix> C = C_in;
592
593 if (C == Teuchos::null) {
594 double nnzPerRow = Teuchos::as<double>(0);
595
596#if 0
597 if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
598 // For now, follow what ML and Epetra do.
599 GO numRowsA = A.getGlobalNumRows();
600 GO numRowsB = B.getGlobalNumRows();
601 nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
602 sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
603 nnzPerRow *= nnzPerRow;
604 double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
605 double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
606 if (totalNnz < minNnz)
607 totalNnz = minNnz;
608 nnzPerRow = totalNnz / A.getGlobalNumRows();
609
610 fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
611 }
612#endif
613 if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
614 else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
615
616 } else {
617 C->resumeFill(); // why this is not done inside of Tpetra MxM?
618 fos << "Reuse C pattern" << std::endl;
619 }
620
621 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
622
623 return C;
624 }
625
636 static RCP<Matrix> Multiply(const Matrix& A, bool transposeA, const Matrix& B, bool transposeB, Teuchos::FancyOStream &fos,
637 bool callFillCompleteOnResult = true, bool doOptimizeStorage = true, const std::string& label = std::string(),
638 const RCP<ParameterList>& params = null) {
639 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
640 }
641
642#ifdef HAVE_XPETRA_EPETRAEXT
643 // Michael Gee's MLMultiply
645 const Epetra_CrsMatrix& epB,
647 throw(Xpetra::Exceptions::RuntimeError("MLTwoMatrixMultiply only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)"));
648 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
649 }
650#endif //ifdef HAVE_XPETRA_EPETRAEXT
651
663 const BlockedCrsMatrix& B, bool transposeB,
665 bool doFillComplete = true,
666 bool doOptimizeStorage = true) {
668 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
669
670 // Preconditions
673
676
677 RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
678
679 for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
680 for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
681 RCP<Matrix> Cij;
682
683 for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
684 RCP<Matrix> crmat1 = A.getMatrix(i,l);
685 RCP<Matrix> crmat2 = B.getMatrix(l,j);
686
687 if (crmat1.is_null() || crmat2.is_null())
688 continue;
689
690 // try unwrapping 1x1 blocked matrix
691 {
692 auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
693 auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
694
695 if(unwrap1.is_null() != unwrap2.is_null()) {
696 if(unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
697 crmat1 = unwrap1->getCrsMatrix();
698 if(unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
699 crmat2 = unwrap2->getCrsMatrix();
700 }
701 }
702
703 RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
704 RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
706 "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
707
708 // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
709 if (!crop1.is_null())
710 Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
711 if (!crop2.is_null())
712 Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
713
714 TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
715 "crmat1 does not have global constants");
716 TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
717 "crmat2 does not have global constants");
718
719 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
720 continue;
721
722 // temporary matrix containing result of local block multiplication
723 RCP<Matrix> temp = Teuchos::null;
724
725 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
726 temp = Multiply (*crop1, false, *crop2, false, fos);
727 else {
728 RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
729 RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
730 TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
731 TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
732 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)");
733 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)");
734
735 // recursive multiplication call
736 temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
737
738 RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
739 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)");
740 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)");
741 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)");
742 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)");
743 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)");
744 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)");
745 }
746
747 TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
748
749 RCP<Matrix> addRes = null;
750 if (Cij.is_null ())
751 Cij = temp;
752 else {
753 MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
754 Cij = addRes;
755 }
756 }
757
758 if (!Cij.is_null()) {
759 if (Cij->isFillComplete())
760 Cij->resumeFill();
761 Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
762 C->setMatrix(i, j, Cij);
763 } else {
764 C->setMatrix(i, j, Teuchos::null);
765 }
766 }
767 }
768
769 if (doFillComplete)
770 C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
771
772 return C;
773 } // TwoMatrixMultiplyBlock
774
787 static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
788 if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
789 throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
790
791 if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
792 throw Exceptions::RuntimeError("TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
793 } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
794#ifdef HAVE_XPETRA_TPETRA
795 const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
796 Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
797
798 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
799#else
800 throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
801#endif
802 }
803 } //MatrixMatrix::TwoMatrixAdd()
804
805
820 static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
821 const Matrix& B, bool transposeB, const SC& beta,
822 RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
823
824 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
825 RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
826 RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
827 RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
828 // We have to distinguish 4 cases:
829 // both matrices are CrsMatrixWrap based, only one of them is CrsMatrixWrap based, or none.
830
831 // C can be null, so just use A to get the lib
832 auto lib = A.getRowMap()->lib();
833
834 // Both matrices are CrsMatrixWrap
835 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
836 if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
837 throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
838 if (lib == Xpetra::UseEpetra) {
839 throw Exceptions::RuntimeError("MatrixMatrix::Add for Epetra only available with Scalar = double, LO = GO = int.");
840 } else if (lib == Xpetra::UseTpetra) {
841 #ifdef HAVE_XPETRA_TPETRA
842 using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
843 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
844 const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
845 const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
846 C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
847 #else
848 throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
849 #endif
850 }
852 if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
853 if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
855 }
856 // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
857 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
858 RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
859 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
860
861 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
862 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
863
864 size_t i = 0;
865 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
866 RCP<Matrix> Cij = Teuchos::null;
867 if(rcpA != Teuchos::null &&
868 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
869 // recursive call
870 TwoMatrixAdd(*rcpA, transposeA, alpha,
871 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
872 Cij, fos, AHasFixedNnzPerRow);
873 } else if (rcpA == Teuchos::null &&
874 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
875 Cij = rcpBopB->getMatrix(i,j);
876 } else if (rcpA != Teuchos::null &&
877 rcpBopB->getMatrix(i,j)==Teuchos::null) {
878 Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
879 } else {
880 Cij = Teuchos::null;
881 }
882
883 if (!Cij.is_null()) {
884 if (Cij->isFillComplete())
885 Cij->resumeFill();
886 Cij->fillComplete();
887 bC->setMatrix(i, j, Cij);
888 } else {
889 bC->setMatrix(i, j, Teuchos::null);
890 }
891 } // loop over columns j
892 }
893 // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
894 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
895 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
896 RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
897
898 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
899 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
900 size_t j = 0;
901 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
902 RCP<Matrix> Cij = Teuchos::null;
903 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
904 rcpB!=Teuchos::null) {
905 // recursive call
906 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
907 *rcpB, transposeB, beta,
908 Cij, fos, AHasFixedNnzPerRow);
909 } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
910 rcpB!=Teuchos::null) {
911 Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
912 } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
913 rcpB==Teuchos::null) {
914 Cij = rcpBopA->getMatrix(i,j);
915 } else {
916 Cij = Teuchos::null;
917 }
918
919 if (!Cij.is_null()) {
920 if (Cij->isFillComplete())
921 Cij->resumeFill();
922 Cij->fillComplete();
923 bC->setMatrix(i, j, Cij);
924 } else {
925 bC->setMatrix(i, j, Teuchos::null);
926 }
927 } // loop over rows i
928 }
929 else {
930 // This is the version for blocked matrices
931
932 // check the compatibility of the blocked operators
933 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
934 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
935 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)");
936 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)");
937 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)");
938 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)");
939 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)");
940 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)");
941
942 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
943 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
944
945 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
946 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
947 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
948 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
949
950 RCP<Matrix> Cij = Teuchos::null;
951 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
952 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
953 // recursive call
954 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
955 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
956 Cij, fos, AHasFixedNnzPerRow);
957 } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
958 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
959 Cij = rcpBopB->getMatrix(i,j);
960 } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
961 rcpBopB->getMatrix(i,j)==Teuchos::null) {
962 Cij = rcpBopA->getMatrix(i,j);
963 } else {
964 Cij = Teuchos::null;
965 }
966
967 if (!Cij.is_null()) {
968 if (Cij->isFillComplete())
969 Cij->resumeFill();
970 Cij->fillComplete();
971 bC->setMatrix(i, j, Cij);
972 } else {
973 bC->setMatrix(i, j, Teuchos::null);
974 }
975 } // loop over columns j
976 } // loop over rows i
977
978 } // end blocked recursive algorithm
979 } //MatrixMatrix::TwoMatrixAdd()
980
981
982 }; // class MatrixMatrix
983
984
985#ifdef HAVE_XPETRA_EPETRA
986 // specialization MatrixMatrix for SC=double, LO=GO=int
987 template<>
988 class MatrixMatrix<double,int,int,EpetraNode> {
989 typedef double Scalar;
990 typedef int LocalOrdinal;
991 typedef int GlobalOrdinal;
994
995 public:
996
1021 static void Multiply(const Matrix& A, bool transposeA,
1022 const Matrix& B, bool transposeB,
1023 Matrix& C,
1024 bool call_FillComplete_on_result = true,
1025 bool doOptimizeStorage = true,
1026 const std::string & label = std::string(),
1027 const RCP<ParameterList>& params = null) {
1028 TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
1029 Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
1030 TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
1031 Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
1032
1035
1036 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1037
1038 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
1039
1040 if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
1041#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1042 helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
1043#else
1044 throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
1045#endif
1046 } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
1047#ifdef HAVE_XPETRA_TPETRA
1048# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1049 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1050 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
1051# else
1052 if(helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
1053 // All matrices are Crs
1054 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
1055 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
1056 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
1057
1058 // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
1059 // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
1060 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1061 }
1062 else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
1063 // All matrices are BlockCrs (except maybe Ac)
1064 // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
1065 // we'll need to think about refactoring BlockCrs so we can do something smarter here.
1066
1067 if(!A.getRowMap()->getComm()->getRank())
1068 std::cout<<"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
1069
1070 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(A);
1071 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(B);
1072 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
1073 RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
1074 RCP<const CRS> Bcrs = Tpetra::convertToCrsMatrix(tpB);
1075
1076 // We need the global constants to do the copy back to BlockCrs
1077 RCP<ParameterList> new_params;
1078 if(!params.is_null()) {
1079 new_params = rcp(new Teuchos::ParameterList(*params));
1080 new_params->set("compute global constants",true);
1081 }
1082
1083 // FIXME: The lines below only works because we're assuming Ac is Point
1084 RCP<CRS> tempAc = Teuchos::rcp(new CRS(Acrs->getRowMap(),0));
1085 Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
1086
1087 // Temporary output matrix
1088 RCP<Tpetra::BlockCrsMatrix<SC,LO,GO,NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*tempAc,A.GetStorageBlockSize());
1091
1092 // We can now cheat and replace the innards of Ac
1093 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(Teuchos::rcpFromRef(C));
1094 Ac_w->replaceCrsMatrix(Ac_p);
1095 }
1096 else {
1097 // Mix and match
1098 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
1099 }
1100# endif
1101#else
1102 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
1103#endif
1104 }
1105
1106 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1108 fillParams->set("Optimize Storage",doOptimizeStorage);
1109 C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
1110 (transposeA) ? A.getDomainMap() : A.getRangeMap(),
1111 fillParams);
1112 }
1113
1114 // transfer striding information
1115 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1116 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1117 C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
1118 } // end Multiply
1119
1142 static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1143 const Matrix& B, bool transposeB,
1144 RCP<Matrix> C_in,
1146 bool doFillComplete = true,
1147 bool doOptimizeStorage = true,
1148 const std::string & label = std::string(),
1149 const RCP<ParameterList>& params = null) {
1150
1151 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1152 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1153
1154 // Optimization using ML Multiply when available and requested
1155 // This feature is currently not supported. We would have to introduce the HAVE_XPETRA_ML_MMM flag
1156#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) && defined(HAVE_XPETRA_ML_MMM)
1157 if (B.getDomainMap()->lib() == Xpetra::UseEpetra && !transposeA && !transposeB) {
1160 RCP<Epetra_CrsMatrix> epC = MLTwoMatrixMultiply(*epA, *epB, fos);
1161
1162 RCP<Matrix> C = Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO> (epC);
1163 if (doFillComplete) {
1165 fillParams->set("Optimize Storage", doOptimizeStorage);
1166 C->fillComplete(B.getDomainMap(), A.getRangeMap(), fillParams);
1167 }
1168
1169 // Fill strided maps information
1170 // This is necessary since the ML matrix matrix multiplication routine has no handling for this
1171 // TODO: move this call to MLMultiply...
1172 C->CreateView("stridedMaps", rcpFromRef(A), transposeA, rcpFromRef(B), transposeB);
1173
1174 return C;
1175 }
1176#endif // EPETRA + EPETRAEXT + ML
1177
1178 // Default case: Xpetra Multiply
1179 RCP<Matrix> C = C_in;
1180
1181 if (C == Teuchos::null) {
1182 double nnzPerRow = Teuchos::as<double>(0);
1183
1184#if 0
1185 if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
1186 // For now, follow what ML and Epetra do.
1187 GO numRowsA = A.getGlobalNumRows();
1188 GO numRowsB = B.getGlobalNumRows();
1189 nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
1190 sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
1191 nnzPerRow *= nnzPerRow;
1192 double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1193 double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1194 if (totalNnz < minNnz)
1195 totalNnz = minNnz;
1196 nnzPerRow = totalNnz / A.getGlobalNumRows();
1197
1198 fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1199 }
1200#endif
1201
1202 if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1203 else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1204
1205 } else {
1206 C->resumeFill(); // why this is not done inside of Tpetra MxM?
1207 fos << "Reuse C pattern" << std::endl;
1208 }
1209
1210 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1211
1212 return C;
1213 }
1214
1225 static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1226 const Matrix& B, bool transposeB,
1228 bool callFillCompleteOnResult = true,
1229 bool doOptimizeStorage = true,
1230 const std::string & label = std::string(),
1231 const RCP<ParameterList>& params = null) {
1232 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
1233 }
1234
1235#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1236 // Michael Gee's MLMultiply
1238 const Epetra_CrsMatrix& epB,
1239 Teuchos::FancyOStream& fos) {
1240#if defined(HAVE_XPETRA_ML_MMM) // Note: this is currently not supported
1241 ML_Comm* comm;
1242 ML_Comm_Create(&comm);
1243 fos << "****** USING ML's MATRIX MATRIX MULTIPLY ******" << std::endl;
1244#ifdef HAVE_MPI
1245 // ML_Comm uses MPI_COMM_WORLD, so try to use the same communicator as epA.
1246 const Epetra_MpiComm * Mcomm = dynamic_cast<const Epetra_MpiComm*>(&(epA.Comm()));
1247 if (Mcomm)
1248 ML_Comm_Set_UsrComm(comm,Mcomm->GetMpiComm());
1249# endif
1250 //in order to use ML, there must be no indices missing from the matrix column maps.
1251 EpetraExt::CrsMatrix_SolverMap Atransform;
1252 EpetraExt::CrsMatrix_SolverMap Btransform;
1253 const Epetra_CrsMatrix& A = Atransform(const_cast<Epetra_CrsMatrix&>(epA));
1254 const Epetra_CrsMatrix& B = Btransform(const_cast<Epetra_CrsMatrix&>(epB));
1255
1256 if (!A.Filled()) throw Exceptions::RuntimeError("A has to be FillCompleted");
1257 if (!B.Filled()) throw Exceptions::RuntimeError("B has to be FillCompleted");
1258
1259 // create ML operators from EpetraCrsMatrix
1260 ML_Operator* ml_As = ML_Operator_Create(comm);
1261 ML_Operator* ml_Bs = ML_Operator_Create(comm);
1262 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?
1263 ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&B),ml_Bs);
1264 ML_Operator* ml_AtimesB = ML_Operator_Create(comm);
1265 {
1267 ML_2matmult(ml_As,ml_Bs,ml_AtimesB,ML_CSR_MATRIX); // do NOT use ML_EpetraCRS_MATRIX!!!
1268 }
1269 ML_Operator_Destroy(&ml_As);
1270 ML_Operator_Destroy(&ml_Bs);
1271
1272 // For ml_AtimesB we have to reconstruct the column map in global indexing,
1273 // The following is going down to the salt-mines of ML ...
1274 // note: we use integers, since ML only works for Epetra...
1275 int N_local = ml_AtimesB->invec_leng;
1276 ML_CommInfoOP* getrow_comm = ml_AtimesB->getrow->pre_comm;
1277 if (!getrow_comm) throw(Exceptions::RuntimeError("ML_Operator does not have a CommInfo"));
1278 ML_Comm* comm_AB = ml_AtimesB->comm; // get comm object
1279 if (N_local != B.DomainMap().NumMyElements())
1280 throw(Exceptions::RuntimeError("Mismatch in local row dimension between ML and Epetra"));
1281 int N_rcvd = 0;
1282 int N_send = 0;
1283 int flag = 0;
1284 for (int i = 0; i < getrow_comm->N_neighbors; i++) {
1285 N_rcvd += (getrow_comm->neighbors)[i].N_rcv;
1286 N_send += (getrow_comm->neighbors)[i].N_send;
1287 if ( ((getrow_comm->neighbors)[i].N_rcv != 0) &&
1288 ((getrow_comm->neighbors)[i].rcv_list != NULL) ) flag = 1;
1289 }
1290 // For some unknown reason, ML likes to have stuff one larger than
1291 // neccessary...
1292 std::vector<double> dtemp(N_local + N_rcvd + 1); // "double" vector for comm function
1293 std::vector<int> cmap (N_local + N_rcvd + 1); // vector for gids
1294 for (int i = 0; i < N_local; ++i) {
1295 cmap[i] = B.DomainMap().GID(i);
1296 dtemp[i] = (double) cmap[i];
1297 }
1298 ML_cheap_exchange_bdry(&dtemp[0],getrow_comm,N_local,N_send,comm_AB); // do communication
1299 if (flag) { // process received data
1300 int count = N_local;
1301 const int neighbors = getrow_comm->N_neighbors;
1302 for (int i = 0; i < neighbors; i++) {
1303 const int nrcv = getrow_comm->neighbors[i].N_rcv;
1304 for (int j = 0; j < nrcv; j++)
1305 cmap[getrow_comm->neighbors[i].rcv_list[j]] = (int) dtemp[count++];
1306 }
1307 } else {
1308 for (int i = 0; i < N_local+N_rcvd; ++i)
1309 cmap[i] = (int)dtemp[i];
1310 }
1311 dtemp.clear(); // free double array
1312
1313 // we can now determine a matching column map for the result
1314 Epetra_Map gcmap(-1, N_local+N_rcvd, &cmap[0], B.ColMap().IndexBase(), A.Comm());
1315
1316 int allocated = 0;
1317 int rowlength;
1318 double* val = NULL;
1319 int* bindx = NULL;
1320
1321 const int myrowlength = A.RowMap().NumMyElements();
1322 const Epetra_Map& rowmap = A.RowMap();
1323
1324 // Determine the maximum bandwith for the result matrix.
1325 // replaces the old, very(!) memory-consuming guess:
1326 // int guessnpr = A.MaxNumEntries()*B.MaxNumEntries();
1327 int educatedguess = 0;
1328 for (int i = 0; i < myrowlength; ++i) {
1329 // get local row
1330 ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1331 if (rowlength>educatedguess)
1332 educatedguess = rowlength;
1333 }
1334
1335 // allocate our result matrix and fill it
1336 RCP<Epetra_CrsMatrix> result = rcp(new Epetra_CrsMatrix(::Copy, A.RangeMap(), gcmap, educatedguess, false));
1337
1338 std::vector<int> gcid(educatedguess);
1339 for (int i = 0; i < myrowlength; ++i) {
1340 const int grid = rowmap.GID(i);
1341 // get local row
1342 ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
1343 if (!rowlength) continue;
1344 if ((int)gcid.size() < rowlength) gcid.resize(rowlength);
1345 for (int j = 0; j < rowlength; ++j) {
1346 gcid[j] = gcmap.GID(bindx[j]);
1347 if (gcid[j] < 0)
1348 throw Exceptions::RuntimeError("Error: cannot find gcid!");
1349 }
1350 int err = result->InsertGlobalValues(grid, rowlength, val, &gcid[0]);
1351 if (err != 0 && err != 1) {
1352 std::ostringstream errStr;
1353 errStr << "Epetra_CrsMatrix::InsertGlobalValues returned err=" << err;
1354 throw Exceptions::RuntimeError(errStr.str());
1355 }
1356 }
1357 // free memory
1358 if (bindx) ML_free(bindx);
1359 if (val) ML_free(val);
1360 ML_Operator_Destroy(&ml_AtimesB);
1361 ML_Comm_Destroy(&comm);
1362
1363 return result;
1364#else // no MUELU_ML
1365 (void)epA;
1366 (void)epB;
1367 (void)fos;
1369 "No ML multiplication available. This feature is currently not supported by Xpetra.");
1370 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
1371#endif
1372 }
1373#endif //ifdef HAVE_XPETRA_EPETRAEXT
1374
1386 const BlockedCrsMatrix& B, bool transposeB,
1388 bool doFillComplete = true,
1389 bool doOptimizeStorage = true) {
1390 TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
1391 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1392
1393 // Preconditions
1394 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1395 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1396
1397 RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
1399
1400 RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1401
1402 for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1403 for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1404 RCP<Matrix> Cij;
1405
1406 for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1407 RCP<Matrix> crmat1 = A.getMatrix(i,l);
1408 RCP<Matrix> crmat2 = B.getMatrix(l,j);
1409
1410 if (crmat1.is_null() || crmat2.is_null())
1411 continue;
1412
1413 // try unwrapping 1x1 blocked matrix
1414 {
1415 auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1416 auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1417
1418 if(unwrap1.is_null() != unwrap2.is_null()) {
1419 if(unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
1420 crmat1 = unwrap1->getCrsMatrix();
1421 if(unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
1422 crmat2 = unwrap2->getCrsMatrix();
1423 }
1424 }
1425
1426 RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
1427 RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
1429 "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
1430
1431 // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
1432 if (!crop1.is_null())
1433 Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
1434 if (!crop2.is_null())
1435 Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
1436
1437 TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
1438 "crmat1 does not have global constants");
1439 TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
1440 "crmat2 does not have global constants");
1441
1442 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
1443 continue;
1444
1445
1446 // temporary matrix containing result of local block multiplication
1447 RCP<Matrix> temp = Teuchos::null;
1448
1449 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
1450 temp = Multiply (*crop1, false, *crop2, false, fos);
1451 else {
1452 RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
1453 RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
1454 TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1455 TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
1456 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)");
1457 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)");
1458
1459 // recursive multiplication call
1460 temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
1461
1462 RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
1463 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)");
1464 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)");
1465 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)");
1466 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)");
1467 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)");
1468 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)");
1469 }
1470
1471 TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
1472
1473 RCP<Matrix> addRes = null;
1474 if (Cij.is_null ())
1475 Cij = temp;
1476 else {
1477 MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1478 Cij = addRes;
1479 }
1480 }
1481
1482 if (!Cij.is_null()) {
1483 if (Cij->isFillComplete())
1484 Cij->resumeFill();
1485 Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1486 C->setMatrix(i, j, Cij);
1487 } else {
1488 C->setMatrix(i, j, Teuchos::null);
1489 }
1490 }
1491 }
1492
1493 if (doFillComplete)
1494 C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1495
1496 return C;
1497 } // TwoMatrixMultiplyBlock
1498
1511 static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
1513 "TwoMatrixAdd: matrix row maps are not the same.");
1514
1515 if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1516#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1519
1520 //FIXME is there a bug if beta=0?
1521 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1522
1523 if (rv != 0)
1524 throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1525 std::ostringstream buf;
1526#else
1527 throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1528#endif
1529 } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1530#ifdef HAVE_XPETRA_TPETRA
1531# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1532 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1533 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1534# else
1535 const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1536 Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1537
1538 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1539# endif
1540#else
1541 throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1542#endif
1543 }
1544 } //MatrixMatrix::TwoMatrixAdd()
1545
1546
1561 static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
1562 const Matrix& B, bool transposeB, const SC& beta,
1563 RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1564 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
1565 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1566 RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
1567 RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
1568 RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
1569
1570
1571 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
1572
1573
1574 if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
1575 throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
1576
1577 auto lib = A.getRowMap()->lib();
1578 if (lib == Xpetra::UseEpetra) {
1579 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1580 const Epetra_CrsMatrix& epA = helpers::Op2EpetraCrs(A);
1581 const Epetra_CrsMatrix& epB = helpers::Op2EpetraCrs(B);
1582 if(C.is_null())
1583 {
1584 size_t maxNzInA = 0;
1585 size_t maxNzInB = 0;
1586 size_t numLocalRows = 0;
1587 if (A.isFillComplete() && B.isFillComplete()) {
1588
1589 maxNzInA = A.getLocalMaxNumRowEntries();
1590 maxNzInB = B.getLocalMaxNumRowEntries();
1591 numLocalRows = A.getLocalNumRows();
1592 }
1593
1594 if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1595 // first check if either A or B has at most 1 nonzero per row
1596 // the case of both having at most 1 nz per row is handled by the ``else''
1597 Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1598
1599 if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1600 for (size_t i = 0; i < numLocalRows; ++i)
1601 exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1602
1603 } else {
1604 for (size_t i = 0; i < numLocalRows; ++i)
1605 exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1606 }
1607
1608 fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1609 << ", using static profiling" << std::endl;
1610 C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow));
1611
1612 } else {
1613 // general case
1614 LO maxPossibleNNZ = A.getLocalMaxNumRowEntries() + B.getLocalMaxNumRowEntries();
1615 C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), maxPossibleNNZ));
1616 }
1617 if (transposeB)
1618 fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1619 }
1620 RCP<Epetra_CrsMatrix> epC = helpers::Op2NonConstEpetraCrs(C);
1621 Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1622
1623 //FIXME is there a bug if beta=0?
1624 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1625
1626 if (rv != 0)
1627 throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1628 #else
1629 throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1630 #endif
1631 } else if (lib == Xpetra::UseTpetra) {
1632 #ifdef HAVE_XPETRA_TPETRA
1633 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1634 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1635 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1636 # else
1637 using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
1638 const tcrs_matrix_type& tpA = helpers::Op2TpetraCrs(A);
1639 const tcrs_matrix_type& tpB = helpers::Op2TpetraCrs(B);
1640 C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
1641 # endif
1642 #else
1643 throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
1644 #endif
1645 }
1646
1648 if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
1649 if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
1651 }
1652 // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
1653 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
1654 RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
1655 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1656
1657 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1658 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1659
1660 size_t i = 0;
1661 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1662 RCP<Matrix> Cij = Teuchos::null;
1663 if(rcpA != Teuchos::null &&
1664 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1665 // recursive call
1666 TwoMatrixAdd(*rcpA, transposeA, alpha,
1667 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1668 Cij, fos, AHasFixedNnzPerRow);
1669 } else if (rcpA == Teuchos::null &&
1670 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1671 Cij = rcpBopB->getMatrix(i,j);
1672 } else if (rcpA != Teuchos::null &&
1673 rcpBopB->getMatrix(i,j)==Teuchos::null) {
1674 Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
1675 } else {
1676 Cij = Teuchos::null;
1677 }
1678
1679 if (!Cij.is_null()) {
1680 if (Cij->isFillComplete())
1681 Cij->resumeFill();
1682 Cij->fillComplete();
1683 bC->setMatrix(i, j, Cij);
1684 } else {
1685 bC->setMatrix(i, j, Teuchos::null);
1686 }
1687 } // loop over columns j
1688 }
1689 // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
1690 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
1691 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1692 RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
1693
1694 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1695 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1696
1697 size_t j = 0;
1698 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1699 RCP<Matrix> Cij = Teuchos::null;
1700 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1701 rcpB!=Teuchos::null) {
1702 // recursive call
1703 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1704 *rcpB, transposeB, beta,
1705 Cij, fos, AHasFixedNnzPerRow);
1706 } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1707 rcpB!=Teuchos::null) {
1708 Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
1709 } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1710 rcpB==Teuchos::null) {
1711 Cij = rcpBopA->getMatrix(i,j);
1712 } else {
1713 Cij = Teuchos::null;
1714 }
1715
1716 if (!Cij.is_null()) {
1717 if (Cij->isFillComplete())
1718 Cij->resumeFill();
1719 Cij->fillComplete();
1720 bC->setMatrix(i, j, Cij);
1721 } else {
1722 bC->setMatrix(i, j, Teuchos::null);
1723 }
1724 } // loop over rows i
1725 }
1726 else {
1727 // This is the version for blocked matrices
1728
1729 // check the compatibility of the blocked operators
1730 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1731 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
1732 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)");
1733 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)");
1734 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)");
1735 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)");
1736 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)");
1737 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)");
1738
1739 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
1740 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
1741
1742 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1743 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
1744
1745 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
1746 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
1747
1748 RCP<Matrix> Cij = Teuchos::null;
1749 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
1750 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1751 // recursive call
1752
1753 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
1754 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
1755 Cij, fos, AHasFixedNnzPerRow);
1756 } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
1757 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
1758 Cij = rcpBopB->getMatrix(i,j);
1759 } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
1760 rcpBopB->getMatrix(i,j)==Teuchos::null) {
1761 Cij = rcpBopA->getMatrix(i,j);
1762 } else {
1763 Cij = Teuchos::null;
1764 }
1765
1766 if (!Cij.is_null()) {
1767 if (Cij->isFillComplete())
1768 Cij->resumeFill();
1769 //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
1770 Cij->fillComplete();
1771 bC->setMatrix(i, j, Cij);
1772 } else {
1773 bC->setMatrix(i, j, Teuchos::null);
1774 }
1775 } // loop over columns j
1776 } // loop over rows i
1777
1778 } // end blocked recursive algorithm
1779 } //MatrixMatrix::TwoMatrixAdd()
1780 }; // end specialization on SC=double, GO=int and NO=EpetraNode
1781
1782 // specialization MatrixMatrix for SC=double, GO=long long and NO=EptraNode
1783 template<>
1784 class MatrixMatrix<double,int,long long,EpetraNode> {
1785 typedef double Scalar;
1786 typedef int LocalOrdinal;
1787 typedef long long GlobalOrdinal;
1789#include "Xpetra_UseShortNames.hpp"
1790
1791 public:
1792
1817 static void Multiply(const Matrix& A, bool transposeA,
1818 const Matrix& B, bool transposeB,
1819 Matrix& C,
1820 bool call_FillComplete_on_result = true,
1821 bool doOptimizeStorage = true,
1822 const std::string & label = std::string(),
1823 const RCP<ParameterList>& params = null) {
1824 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
1825 TEUCHOS_TEST_FOR_EXCEPTION(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false,
1826 Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A");
1827 TEUCHOS_TEST_FOR_EXCEPTION(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false,
1828 Xpetra::Exceptions::RuntimeError, "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A");
1829
1832
1833 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
1834
1835 if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
1836#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1837 helpers::epetraExtMult(A, transposeA, B, transposeB, C, haveMultiplyDoFillComplete);
1838#else
1839 throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
1840#endif
1841 } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
1842#ifdef HAVE_XPETRA_TPETRA
1843# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
1844 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
1845 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long, EpetraNode> ETI enabled."));
1846# else
1847 if(helpers::isTpetraCrs(A) && helpers::isTpetraCrs(B) && helpers::isTpetraCrs(C)) {
1848 // All matrices are Crs
1849 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = helpers::Op2TpetraCrs(A);
1850 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = helpers::Op2TpetraCrs(B);
1851 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = helpers::Op2NonConstTpetraCrs(C);
1852
1853 // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
1854 // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
1855 Tpetra::MatrixMatrix::Multiply(tpA, transposeA, tpB, transposeB, tpC, haveMultiplyDoFillComplete, label, params);
1856 }
1857 else if (helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(B)) {
1858 // All matrices are BlockCrs (except maybe Ac)
1859 // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
1860 // we'll need to think about refactoring BlockCrs so we can do something smarter here.
1861
1862 if(!A.getRowMap()->getComm()->getRank())
1863 std::cout<<"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
1864
1865 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(A);
1866 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(B);
1867 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
1868 RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
1869 RCP<const CRS> Bcrs = Tpetra::convertToCrsMatrix(tpB);
1870
1871 // We need the global constants to do the copy back to BlockCrs
1872 RCP<ParameterList> new_params;
1873 if(!params.is_null()) {
1874 new_params = rcp(new Teuchos::ParameterList(*params));
1875 new_params->set("compute global constants",true);
1876 }
1877
1878 // FIXME: The lines below only works because we're assuming Ac is Point
1879 RCP<CRS> tempAc = Teuchos::rcp(new CRS(Acrs->getRowMap(),0));
1880 Tpetra::MatrixMatrix::Multiply(*Acrs, transposeA, *Bcrs, transposeB, *tempAc, haveMultiplyDoFillComplete, label, new_params);
1881
1882 // Temporary output matrix
1883 RCP<Tpetra::BlockCrsMatrix<SC,LO,GO,NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*tempAc,A.GetStorageBlockSize());
1886
1887 // We can now cheat and replace the innards of Ac
1888 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(Teuchos::rcpFromRef(C));
1889 Ac_w->replaceCrsMatrix(Ac_p);
1890 }
1891 else {
1892 // Mix and match
1893 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
1894 }
1895
1896# endif
1897#else
1898 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
1899#endif
1900 }
1901
1902 if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
1904 fillParams->set("Optimize Storage",doOptimizeStorage);
1905 C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
1906 (transposeA) ? A.getDomainMap() : A.getRangeMap(),
1907 fillParams);
1908 }
1909
1910 // transfer striding information
1911 RCP<Matrix> rcpA = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(A));
1912 RCP<Matrix> rcpB = Teuchos::rcp_const_cast<Matrix>(Teuchos::rcpFromRef(B));
1913 C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
1914 } // end Multiply
1915
1938 static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1939 const Matrix& B, bool transposeB,
1940 RCP<Matrix> C_in,
1942 bool doFillComplete = true,
1943 bool doOptimizeStorage = true,
1944 const std::string & label = std::string(),
1945 const RCP<ParameterList>& params = null) {
1946
1947 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
1948 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
1949
1950 // Default case: Xpetra Multiply
1951 RCP<Matrix> C = C_in;
1952
1953 if (C == Teuchos::null) {
1954 double nnzPerRow = Teuchos::as<double>(0);
1955
1956#if 0
1957 if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
1958 // For now, follow what ML and Epetra do.
1959 GO numRowsA = A.getGlobalNumRows();
1960 GO numRowsB = B.getGlobalNumRows();
1961 nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
1962 sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
1963 nnzPerRow *= nnzPerRow;
1964 double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
1965 double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
1966 if (totalNnz < minNnz)
1967 totalNnz = minNnz;
1968 nnzPerRow = totalNnz / A.getGlobalNumRows();
1969
1970 fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
1971 }
1972#endif
1973 if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
1974 else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
1975
1976 } else {
1977 C->resumeFill(); // why this is not done inside of Tpetra MxM?
1978 fos << "Reuse C pattern" << std::endl;
1979 }
1980
1981 Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage, label, params); // call Multiply routine from above
1982
1983 return C;
1984 }
1985
1996 static RCP<Matrix> Multiply(const Matrix& A, bool transposeA,
1997 const Matrix& B, bool transposeB,
1999 bool callFillCompleteOnResult = true,
2000 bool doOptimizeStorage = true,
2001 const std::string & label = std::string(),
2002 const RCP<ParameterList>& params = null) {
2003 return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage, label, params);
2004 }
2005
2006#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
2007 // Michael Gee's MLMultiply
2009 const Epetra_CrsMatrix& /* epB */,
2010 Teuchos::FancyOStream& /* fos */) {
2012 "No ML multiplication available. This feature is currently not supported by Xpetra.");
2013 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
2014 }
2015#endif //ifdef HAVE_XPETRA_EPETRAEXT
2016
2028 const BlockedCrsMatrix& B, bool transposeB,
2030 bool doFillComplete = true,
2031 bool doOptimizeStorage = true) {
2032 TEUCHOS_TEST_FOR_EXCEPTION(transposeA || transposeB, Exceptions::RuntimeError,
2033 "TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
2034
2035 // Preconditions
2036 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
2037 TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
2038
2039 RCP<const MapExtractor> rgmapextractor = A.getRangeMapExtractor();
2041
2042 RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2043
2044 for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
2045 for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
2046 RCP<Matrix> Cij;
2047
2048 for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
2049 RCP<Matrix> crmat1 = A.getMatrix(i,l);
2050 RCP<Matrix> crmat2 = B.getMatrix(l,j);
2051
2052 if (crmat1.is_null() || crmat2.is_null())
2053 continue;
2054
2055 // try unwrapping 1x1 blocked matrix
2056 {
2057 auto unwrap1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
2058 auto unwrap2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
2059
2060 if(unwrap1.is_null() != unwrap2.is_null()) {
2061 if(unwrap1 != Teuchos::null && unwrap1->Rows() == 1 && unwrap1->Cols() == 1)
2062 crmat1 = unwrap1->getCrsMatrix();
2063 if(unwrap2 != Teuchos::null && unwrap2->Rows() == 1 && unwrap2->Cols() == 1)
2064 crmat2 = unwrap2->getCrsMatrix();
2065 }
2066 }
2067
2068 RCP<CrsMatrixWrap> crop1 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat1);
2069 RCP<CrsMatrixWrap> crop2 = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(crmat2);
2071 "A and B must be either both (compatible) BlockedCrsMatrix objects or both CrsMatrixWrap objects.");
2072
2073 // Forcibly compute the global constants if we don't have them (only works for real CrsMatrices, not nested blocks)
2074 if (!crop1.is_null())
2075 Teuchos::rcp_const_cast<CrsGraph>(crmat1->getCrsGraph())->computeGlobalConstants();
2076 if (!crop2.is_null())
2077 Teuchos::rcp_const_cast<CrsGraph>(crmat2->getCrsGraph())->computeGlobalConstants();
2078
2079 TEUCHOS_TEST_FOR_EXCEPTION(!crmat1->haveGlobalConstants(), Exceptions::RuntimeError,
2080 "crmat1 does not have global constants");
2081 TEUCHOS_TEST_FOR_EXCEPTION(!crmat2->haveGlobalConstants(), Exceptions::RuntimeError,
2082 "crmat2 does not have global constants");
2083
2084 if (crmat1->getGlobalNumEntries() == 0 || crmat2->getGlobalNumEntries() == 0)
2085 continue;
2086
2087 // temporary matrix containing result of local block multiplication
2088 RCP<Matrix> temp = Teuchos::null;
2089
2090 if(crop1 != Teuchos::null && crop2 != Teuchos::null)
2091 temp = Multiply (*crop1, false, *crop2, false, fos);
2092 else {
2093 RCP<BlockedCrsMatrix> bop1 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat1);
2094 RCP<BlockedCrsMatrix> bop2 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(crmat2);
2095 TEUCHOS_TEST_FOR_EXCEPTION(bop1.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
2096 TEUCHOS_TEST_FOR_EXCEPTION(bop2.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixMultiplyBlock)");
2097 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)");
2098 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)");
2099
2100 // recursive multiplication call
2101 temp = TwoMatrixMultiplyBlock(*bop1, transposeA, *bop2, transposeB, fos, doFillComplete, doOptimizeStorage);
2102
2103 RCP<BlockedCrsMatrix> btemp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(temp);
2104 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)");
2105 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)");
2106 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)");
2107 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)");
2108 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)");
2109 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)");
2110 }
2111
2112 TEUCHOS_TEST_FOR_EXCEPTION(temp->isFillComplete() == false, Xpetra::Exceptions::RuntimeError, "Local block is not filled. (TwoMatrixMultiplyBlock)");
2113
2114 RCP<Matrix> addRes = null;
2115 if (Cij.is_null ())
2116 Cij = temp;
2117 else {
2118 MatrixMatrix::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
2119 Cij = addRes;
2120 }
2121 }
2122
2123 if (!Cij.is_null()) {
2124 if (Cij->isFillComplete())
2125 Cij->resumeFill();
2126 Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
2127 C->setMatrix(i, j, Cij);
2128 } else {
2129 C->setMatrix(i, j, Teuchos::null);
2130 }
2131 }
2132 }
2133
2134 if (doFillComplete)
2135 C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
2136
2137 return C;
2138 } // TwoMatrixMultiplyBlock
2139
2152 static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
2154 "TwoMatrixAdd: matrix row maps are not the same.");
2155
2156 if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
2157#if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
2160
2161 //FIXME is there a bug if beta=0?
2162 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
2163
2164 if (rv != 0)
2165 throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
2166 std::ostringstream buf;
2167#else
2168 throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
2169#endif
2170 } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
2171#ifdef HAVE_XPETRA_TPETRA
2172# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
2173 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
2174 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
2175# else
2176 const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
2177 Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
2178
2179 Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
2180# endif
2181#else
2182 throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
2183#endif
2184 }
2185 } //MatrixMatrix::TwoMatrixAdd()
2186
2187
2202 static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
2203 const Matrix& B, bool transposeB, const SC& beta,
2204 RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
2205 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
2206 RCP<const Matrix> rcpB = Teuchos::rcpFromRef(B);
2207 RCP<const BlockedCrsMatrix> rcpBopA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpA);
2208 RCP<const BlockedCrsMatrix> rcpBopB = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpB);
2209
2210 if(rcpBopA == Teuchos::null && rcpBopB == Teuchos::null) {
2212 "TwoMatrixAdd: matrix row maps are not the same.");
2213 auto lib = A.getRowMap()->lib();
2214 if (lib == Xpetra::UseEpetra) {
2215 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
2218 if(C.is_null())
2219 {
2220 size_t maxNzInA = 0;
2221 size_t maxNzInB = 0;
2222 size_t numLocalRows = 0;
2223 if (A.isFillComplete() && B.isFillComplete()) {
2224
2225 maxNzInA = A.getLocalMaxNumRowEntries();
2226 maxNzInB = B.getLocalMaxNumRowEntries();
2227 numLocalRows = A.getLocalNumRows();
2228 }
2229
2230 if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
2231 // first check if either A or B has at most 1 nonzero per row
2232 // the case of both having at most 1 nz per row is handled by the ``else''
2233 Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
2234
2235 if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
2236 for (size_t i = 0; i < numLocalRows; ++i)
2237 exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
2238
2239 } else {
2240 for (size_t i = 0; i < numLocalRows; ++i)
2241 exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
2242 }
2243
2244 fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
2245 << ", using static profiling" << std::endl;
2246 C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), exactNnzPerRow));
2247
2248 } else {
2249 // general case
2250 LO maxPossibleNNZ = A.getLocalMaxNumRowEntries() + B.getLocalMaxNumRowEntries();
2251 C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), maxPossibleNNZ));
2252 }
2253 if (transposeB)
2254 fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
2255 }
2257 Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
2258
2259 //FIXME is there a bug if beta=0?
2260 int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
2261
2262 if (rv != 0)
2263 throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
2264 #else
2265 throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
2266 #endif
2267 } else if (lib == Xpetra::UseTpetra) {
2268 #ifdef HAVE_XPETRA_TPETRA
2269 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
2270 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
2271 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=long long enabled."));
2272 # else
2273 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
2274 using tcrs_matrix_type = Tpetra::CrsMatrix<SC,LO,GO,NO>;
2275 const tcrs_matrix_type& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
2276 const tcrs_matrix_type& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
2277 C = helpers::tpetraAdd(tpA, transposeA, alpha, tpB, transposeB, beta);
2278 # endif
2279 #else
2280 throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
2281 #endif
2282 }
2283
2285 if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
2286 if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
2288 }
2289 // the first matrix is of type CrsMatrixWrap, the second is a blocked operator
2290 else if(rcpBopA == Teuchos::null && rcpBopB != Teuchos::null) {
2291 RCP<const MapExtractor> rgmapextractor = rcpBopB->getRangeMapExtractor();
2292 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2293
2294 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2295 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2296
2297 size_t i = 0;
2298 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2299 RCP<Matrix> Cij = Teuchos::null;
2300 if(rcpA != Teuchos::null &&
2301 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2302 // recursive call
2303 TwoMatrixAdd(*rcpA, transposeA, alpha,
2304 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2305 Cij, fos, AHasFixedNnzPerRow);
2306 } else if (rcpA == Teuchos::null &&
2307 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2308 Cij = rcpBopB->getMatrix(i,j);
2309 } else if (rcpA != Teuchos::null &&
2310 rcpBopB->getMatrix(i,j)==Teuchos::null) {
2311 Cij = Teuchos::rcp_const_cast<Matrix>(rcpA);
2312 } else {
2313 Cij = Teuchos::null;
2314 }
2315
2316 if (!Cij.is_null()) {
2317 if (Cij->isFillComplete())
2318 Cij->resumeFill();
2319 Cij->fillComplete();
2320 bC->setMatrix(i, j, Cij);
2321 } else {
2322 bC->setMatrix(i, j, Teuchos::null);
2323 }
2324 } // loop over columns j
2325 }
2326 // the second matrix is of type CrsMatrixWrap, the first is a blocked operator
2327 else if(rcpBopA != Teuchos::null && rcpBopB == Teuchos::null) {
2328 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2329 RCP<const MapExtractor> domapextractor = rcpBopA->getDomainMapExtractor();
2330
2331 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2332 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2333
2334 size_t j = 0;
2335 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2336 RCP<Matrix> Cij = Teuchos::null;
2337 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2338 rcpB!=Teuchos::null) {
2339 // recursive call
2340 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2341 *rcpB, transposeB, beta,
2342 Cij, fos, AHasFixedNnzPerRow);
2343 } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2344 rcpB!=Teuchos::null) {
2345 Cij = Teuchos::rcp_const_cast<Matrix>(rcpB);
2346 } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2347 rcpB==Teuchos::null) {
2348 Cij = rcpBopA->getMatrix(i,j);
2349 } else {
2350 Cij = Teuchos::null;
2351 }
2352
2353 if (!Cij.is_null()) {
2354 if (Cij->isFillComplete())
2355 Cij->resumeFill();
2356 Cij->fillComplete();
2357 bC->setMatrix(i, j, Cij);
2358 } else {
2359 bC->setMatrix(i, j, Teuchos::null);
2360 }
2361 } // loop over rows i
2362 }
2363 else {
2364 // This is the version for blocked matrices
2365
2366 // check the compatibility of the blocked operators
2367 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopA.is_null()==true, Xpetra::Exceptions::BadCast, "A is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2368 TEUCHOS_TEST_FOR_EXCEPTION(rcpBopB.is_null()==true, Xpetra::Exceptions::BadCast, "B is not a BlockedCrsMatrix. (TwoMatrixAdd)");
2369 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)");
2370 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)");
2371 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)");
2372 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)");
2373 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)");
2374 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)");
2375
2376 RCP<const MapExtractor> rgmapextractor = rcpBopA->getRangeMapExtractor();
2377 RCP<const MapExtractor> domapextractor = rcpBopB->getDomainMapExtractor();
2378
2379 C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
2380 RCP<BlockedCrsMatrix> bC = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(C);
2381
2382 for (size_t i = 0; i < rcpBopA->Rows(); ++i) { // loop over all block rows of A
2383 for (size_t j = 0; j < rcpBopB->Cols(); ++j) { // loop over all block columns of B
2384
2385 RCP<Matrix> Cij = Teuchos::null;
2386 if(rcpBopA->getMatrix(i,j)!=Teuchos::null &&
2387 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2388 // recursive call
2389
2390 TwoMatrixAdd(*(rcpBopA->getMatrix(i,j)), transposeA, alpha,
2391 *(rcpBopB->getMatrix(i,j)), transposeB, beta,
2392 Cij, fos, AHasFixedNnzPerRow);
2393 } else if (rcpBopA->getMatrix(i,j) == Teuchos::null &&
2394 rcpBopB->getMatrix(i,j)!=Teuchos::null) {
2395 Cij = rcpBopB->getMatrix(i,j);
2396 } else if (rcpBopA->getMatrix(i,j) != Teuchos::null &&
2397 rcpBopB->getMatrix(i,j)==Teuchos::null) {
2398 Cij = rcpBopA->getMatrix(i,j);
2399 } else {
2400 Cij = Teuchos::null;
2401 }
2402
2403 if (!Cij.is_null()) {
2404 if (Cij->isFillComplete())
2405 Cij->resumeFill();
2406 //Cij->fillComplete(rcpBopA->getDomainMap(j), rcpBopA->getRangeMap(i));
2407 Cij->fillComplete();
2408 bC->setMatrix(i, j, Cij);
2409 } else {
2410 bC->setMatrix(i, j, Teuchos::null);
2411 }
2412 } // loop over columns j
2413 } // loop over rows i
2414
2415 } // end blocked recursive algorithm
2416 } //MatrixMatrix::TwoMatrixAdd()
2417 }; // end specialization on GO=long long and NO=EpetraNode
2418
2419#endif // HAVE_XPETRA_EPETRA
2420
2421} // end namespace Xpetra
2422
2423#define XPETRA_MATRIXMATRIX_SHORT
2424
2425#endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_ */
LocalOrdinal LO
GlobalOrdinal GO
int GID(int LID) const
int IndexBase() const
int NumMyElements() const
bool Filled() const
int ExtractGlobalRowView(int GlobalRow, int &NumEntries, double *&Values, int *&Indices) const
int NumGlobalEntries(long long Row) const
const Epetra_Map & RowMap() const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
const Epetra_Comm & Comm() const
const Epetra_Map & RangeMap() const
const Epetra_Map & ColMap() const
const Epetra_Map & DomainMap() const
MPI_Comm GetMpiComm() const
bool is_null() const
static RCP< Time > getNewTimer(const std::string &name)
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual size_t Rows() const
number of row blocks
RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode.
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
virtual size_t Cols() const
number of column blocks
Concrete implementation of Xpetra::Matrix.
RCP< CrsMatrix > getCrsMatrix() const
Exception indicating invalid cast attempted.
Exception throws to report incompatible objects (like maps).
Exception throws to report errors in the internal logical of the program.
Tpetra::CrsMatrix< SC, LO, GO, NO > tcrs_matrix_type
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
static Tpetra::BlockCrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraBlockCrs(const Matrix &Op)
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
static bool isTpetraCrs(const Matrix &Op)
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Matrix &Op)
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)
static bool isTpetraBlockCrs(RCP< Matrix > Op)
static bool isTpetraBlockCrs(const Matrix &Op)
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Matrix > Op)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
static bool isTpetraCrs(RCP< Matrix > Op)
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
static Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraCrs(const Matrix &Op)
static RCP< Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraBlockCrs(RCP< Matrix > Op)
static const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > & Op2TpetraBlockCrs(const Matrix &Op)
static void epetraExtMult(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool fillCompleteResult)
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
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.
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
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.
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
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.
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)
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
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.
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
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.
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.
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &, const Epetra_CrsMatrix &, Teuchos::FancyOStream &)
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)
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
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.
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
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.
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
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.
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)
Xpetra-specific matrix class.
virtual global_size_t getGlobalNumEntries() const =0
Returns the global number of entries in this matrix.
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode.
virtual size_t getLocalNumRows() const =0
Returns the number of matrix rows owned on the calling node.
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.
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
Returns the current number of entries on this node in the specified local row.
bool IsView(const viewLabel_t viewLabel) const
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
virtual size_t getLocalMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
virtual LocalOrdinal GetStorageBlockSize() const =0
Returns the block size of the storage mechanism, which is usually 1, except for Tpetra::BlockCrsMatri...
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y....
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X....
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
std::string toString(const T &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace