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