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_Matrix.hpp"
54 #include "Xpetra_CrsMatrixWrap.hpp"
55 
56 #include "Xpetra_Map.hpp"
57 #include "Xpetra_StridedMap.hpp"
59 #include "Xpetra_MapExtractor.hpp"
60 #include "Xpetra_MatrixFactory.hpp"
62 
63 #ifdef HAVE_XPETRA_EPETRA
65 #endif
66 
67 #ifdef HAVE_XPETRA_EPETRAEXT
68 #include <EpetraExt_MatrixMatrix.h>
69 #include <EpetraExt_RowMatrixOut.h>
70 #include <Epetra_RowMatrixTransposer.h>
71 #endif // HAVE_XPETRA_EPETRAEXT
72 
73 #ifdef HAVE_XPETRA_TPETRA
74 #include <TpetraExt_MatrixMatrix.hpp>
75 #include <MatrixMarket_Tpetra.hpp>
77 #include <Xpetra_TpetraVector.hpp>
79 #endif // HAVE_XPETRA_TPETRA
80 
81 namespace Xpetra {
82 
89  template <class Scalar,
90  class LocalOrdinal = int,
91  class GlobalOrdinal = LocalOrdinal,
93  class Helpers {
94 #include "Xpetra_UseShortNames.hpp"
95 
96  public:
97 
98 #ifdef HAVE_XPETRA_EPETRA
101  // Get the underlying Epetra Mtx
103  if (crsOp == Teuchos::null)
104  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
106  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
107  if (tmp_ECrsMtx == Teuchos::null)
108  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed"));
109  A = tmp_ECrsMtx->getEpetra_CrsMatrix();
110  return A;
111  } //Op2EpetraCrs
112 
115  // Get the underlying Epetra Mtx
117  if (crsOp == Teuchos::null)
118  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
120  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
121  if (tmp_ECrsMtx == Teuchos::null)
122  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed"));
123  A = tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
124  return A;
125  } //Op2NonConstEpetraCrs
126 
127  static const Epetra_CrsMatrix & Op2EpetraCrs(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
129  // Get the underlying Epetra Mtx
130  try {
133  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
134  if (tmp_ECrsMtx == Teuchos::null)
135  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed"));
136  A = tmp_ECrsMtx->getEpetra_CrsMatrix();
137  return *A;
138  } catch(...) {
139  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
140  }
141  } //Op2EpetraCrs
142 
145  // Get the underlying Epetra Mtx
146  try {
149  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > &tmp_ECrsMtx = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(tmp_CrsMtx);
150  if (tmp_ECrsMtx == Teuchos::null)
151  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed"));
152  A = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(tmp_ECrsMtx->getEpetra_CrsMatrix());
153  return *A;
154  } catch(...) {
155  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
156  }
157  } //Op2NonConstEpetraCrs
158 #endif // HAVE_XPETRA_EPETRA
159 
160 #ifdef HAVE_XPETRA_TPETRA
163  // Get the underlying Tpetra Mtx
165  if (crsOp == Teuchos::null)
166  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
169  if (tmp_ECrsMtx == Teuchos::null)
170  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed"));
171  A = tmp_ECrsMtx->getTpetra_CrsMatrix();
172  return A;
173  } //Op2TpetraCrs
174 
177  // Get the underlying Tpetra Mtx
179  if (crsOp == Teuchos::null)
180  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
183  if (tmp_ECrsMtx == Teuchos::null)
184  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed"));
185  A = tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
186  return A;
187  } //Op2NonConstTpetraCrs
188 
189  static const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op2TpetraCrs(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
191  // Get the underlying Tpetra Mtx
192  try{
196  if (tmp_TCrsMtx == Teuchos::null)
197  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed"));
198  A = tmp_TCrsMtx->getTpetra_CrsMatrix();
199  return *A;
200  } catch (...) {
201  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
202  }
203  } //Op2TpetraCrs
204 
205  static Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op2NonConstTpetraCrs(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Op) {
207  // Get the underlying Tpetra Mtx
208  try{
212  if (tmp_TCrsMtx == Teuchos::null)
213  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed"));
214  A = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tmp_TCrsMtx->getTpetra_CrsMatrix());
215  return *A;
216  } catch (...) {
217  throw(Xpetra::Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed"));
218  }
219  } //Op2NonConstTpetraCrs
220 
221 #endif // HAVE_XPETRA_TPETRA
222 
223  };
224 
225  template <class Scalar,
226  class LocalOrdinal /*= int*/,
227  class GlobalOrdinal /*= LocalOrdinal*/,
228  class Node /*= KokkosClassic::DefaultNode::DefaultNodeType*/>
229  class MatrixMatrix {
230 #undef XPETRA_MATRIXMATRIX_SHORT
231 #include "Xpetra_UseShortNames.hpp"
232 
233  public:
234 
259  static void Multiply(
261  bool transposeA,
263  bool transposeB,
265  bool call_FillComplete_on_result = true,
266  bool doOptimizeStorage = true,
267  const std::string & label = std::string()) {
268 
269  if(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false) {
270  std::string msg = "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A";
272  }
273  else if(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false) {
274  std::string msg = "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A";
276  }
277 
278 
279  if (!A.isFillComplete())
280  throw(Xpetra::Exceptions::RuntimeError("A is not fill-completed"));
281  if (!B.isFillComplete())
282  throw(Xpetra::Exceptions::RuntimeError("B is not fill-completed"));
283 
284  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
285 
286  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
287 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
291 
292  int i = EpetraExt::MatrixMatrix::Multiply(epA,transposeA,epB,transposeB,epC,haveMultiplyDoFillComplete);
293  if (haveMultiplyDoFillComplete) {
294  // Due to Epetra wrapper intricacies, we need to explicitly call
295  // fillComplete on Xpetra matrix here. Specifically, EpetraCrsMatrix
296  // only keeps an internal variable to check whether we are in resumed
297  // state or not, but never touches the underlying Epetra object. As
298  // such, we need to explicitly update the state of Xpetra matrix to
299  // that of Epetra one afterwords
300  C.fillComplete();
301  }
302 
303  if (i != 0) {
304  std::ostringstream buf;
305  buf << i;
306  std::string msg = "EpetraExt::MatrixMatrix::Multiply return value of " + buf.str();
307  throw(Exceptions::RuntimeError(msg));
308  }
309 #else
310  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
311 
312 #endif
313  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
314 #ifdef HAVE_XPETRA_TPETRA
315  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & tpA = Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Op2TpetraCrs(A);
316  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & tpB = Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Op2TpetraCrs(B);
317  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> & tpC = Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Op2NonConstTpetraCrs(C);
318 
319  //18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
320  //Previously, Tpetra's matrix matrix multiply did not support fillComplete.
321  Tpetra::MatrixMatrix::Multiply(tpA,transposeA,tpB,transposeB,tpC,haveMultiplyDoFillComplete,label);
322 #else
323  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
324 #endif
325  }
326 
327  if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
329  params->set("Optimize Storage",doOptimizeStorage);
330  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
331  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
332  params);
333  }
334 
335  // transfer striding information
338  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
339  } // end Multiply
340 
364  const Matrix& B, bool transposeB,
365  RCP<Matrix> C_in,
367  bool doFillComplete = true,
368  bool doOptimizeStorage = true,
369  const std::string & label = std::string()) {
370 
371  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
372  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
373 
374  // Default case: Xpetra Multiply
375  RCP<Matrix> C = C_in;
376 
377  if (C == Teuchos::null) {
378  double nnzPerRow = Teuchos::as<double>(0);
379 
380  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
381  // For now, follow what ML and Epetra do.
382  GO numRowsA = A.getGlobalNumRows();
383  GO numRowsB = B.getGlobalNumRows();
384  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
385  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
386  nnzPerRow *= nnzPerRow;
387  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
388  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
389  if (totalNnz < minNnz)
390  totalNnz = minNnz;
391  nnzPerRow = totalNnz / A.getGlobalNumRows();
392 
393  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
394  }
395 
396  if (transposeA) C = MatrixFactory::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
397  else C = MatrixFactory::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
398 
399  } else {
400  C->resumeFill(); // why this is not done inside of Tpetra MxM?
401  fos << "Reuse C pattern" << std::endl;
402  }
403 
404  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage,label); // call Multiply routine from above
405 
406  return C;
407  }
408 
420  bool transposeA,
421  const Matrix & B,
422  bool transposeB,
424  bool callFillCompleteOnResult = true,
425  bool doOptimizeStorage = true,
426  const std::string & label = std::string()){
427  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage,label);
428  }
429 
430 #ifdef HAVE_XPETRA_EPETRAEXT
431  // Michael Gee's MLMultiply
432  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
433  const Epetra_CrsMatrix& epB,
434  Teuchos::FancyOStream& fos) {
435  throw(Xpetra::Exceptions::RuntimeError("MLTwoMatrixMultiply only available for SC=double and GO=LO=int."));
436  return Teuchos::null;
437  }
438 #endif //ifdef HAVE_XPETRA_EPETRAEXT
439 
451  BlockedCrsMatrix& B, bool transposeB,
453  bool doFillComplete = true,
454  bool doOptimizeStorage = true) {
455  if (transposeA || transposeB)
456  throw Exceptions::RuntimeError("TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
457 
458  // Preconditions
459  if (!A.isFillComplete())
460  throw Exceptions::RuntimeError("A is not fill-completed");
461  if (!B.isFillComplete())
462  throw Exceptions::RuntimeError("B is not fill-completed");
463 
466 
467  RCP<BlockedCrsMatrix> C = rcp(new BlockedCrsMatrix(rgmapextractor, domapextractor, 33 /* TODO fix me */));
468 
469  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
470  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
471  RCP<Matrix> Cij;
472 
473  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
474  RCP<CrsMatrix> crmat1 = A.getMatrix(i,l);
475  RCP<CrsMatrix> crmat2 = B.getMatrix(l,j);
476 
477  if (crmat1.is_null() || crmat2.is_null()) {
478  continue;
479  }
480 
481  RCP<CrsMatrixWrap> crop1 = rcp(new CrsMatrixWrap(crmat1));
482  RCP<CrsMatrixWrap> crop2 = rcp(new CrsMatrixWrap(crmat2));
483 
484  RCP<Matrix> temp = Multiply (*crop1, false, *crop2, false, fos);
485 
486  RCP<Matrix> addRes = null;
487  if (Cij.is_null ())
488  Cij = temp;
489  else {
490  Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
491  Cij = addRes;
492  }
493  }
494 
495  if (!Cij.is_null()) {
496  if (Cij->isFillComplete())
497  Cij->resumeFill();
498  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
499 
500  RCP<CrsMatrixWrap> crsCij = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Cij);
502  "MatrixFactory failed in generating a CrsMatrixWrap." );
503 
504  RCP<CrsMatrix> crsMatCij = crsCij->getCrsMatrix();
505 
506  C->setMatrix(i, j, crsMatCij);
507 
508  } else {
509  C->setMatrix(i, j, Teuchos::null);
510  }
511  }
512  }
513 
514  if (doFillComplete)
515  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
516 
517  return C;
518  } // TwoMatrixMultiplyBlock
519 
532  static void TwoMatrixAdd(const Matrix& A, bool transposeA, SC alpha, Matrix& B, SC beta) {
533  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
534  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
535 
536  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
537  throw Exceptions::RuntimeError("TwoMatrixAdd for Epetra matrices needs <double,int,int> for Scalar, LocalOrdinal and GlobalOrdinal.");
538  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
539 #ifdef HAVE_XPETRA_TPETRA
540  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpA = Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Op2TpetraCrs(A);
541  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpB = Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Op2NonConstTpetraCrs(B);
542 
543  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
544 #else
545  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
546 #endif
547  }
548  } //MatrixMatrix::TwoMatrixAdd()
549 
550 
563  static void TwoMatrixAdd(const Matrix& A, bool transposeA, const SC& alpha,
564  const Matrix& B, bool transposeB, const SC& beta,
565  RCP<Matrix>& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
566  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
567  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
568 
569  if (C == Teuchos::null) {
570  if (!A.isFillComplete() || !B.isFillComplete())
571  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Global statistics are not available for estimates.");
572 
573  size_t maxNzInA = A.getGlobalMaxNumRowEntries();
574  size_t maxNzInB = B.getGlobalMaxNumRowEntries();
575  size_t numLocalRows = A.getNodeNumRows();
576 
577  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
578  // first check if either A or B has at most 1 nonzero per row
579  // the case of both having at most 1 nz per row is handled by the ``else''
580  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
581 
582  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
583  for (size_t i = 0; i < numLocalRows; ++i)
584  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
585 
586  } else {
587  for (size_t i = 0; i < numLocalRows; ++i)
588  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
589  }
590 
591  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
592  << ", using static profiling" << std::endl;
593  C = rcp(new CrsMatrixWrap(A.getRowMap(), exactNnzPerRow, Xpetra::StaticProfile));
594 
595  } else {
596  // general case
597  double nnzPerRowInA = Teuchos::as<double>(A.getGlobalNumEntries()) / A.getGlobalNumRows();
598  double nnzPerRowInB = Teuchos::as<double>(B.getGlobalNumEntries()) / B.getGlobalNumRows();
599  LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
600 
602  //Use static profiling (more efficient) if the estimate is at least as big as the max
603  //possible nnz's in any single row of the result.
604  Xpetra::ProfileType pft = (maxPossible) > nnzToAllocate ? Xpetra::DynamicProfile : Xpetra::StaticProfile;
605 
606  fos << "nnzPerRowInA = " << nnzPerRowInA << ", nnzPerRowInB = " << nnzPerRowInB << std::endl;
607  fos << "MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
608  << ", max possible nnz per row in sum = " << maxPossible
609  << ", using " << (pft == Xpetra::DynamicProfile ? "dynamic" : "static" ) << " profiling"
610  << std::endl;
611  C = rcp(new CrsMatrixWrap(A.getRowMap(), nnzToAllocate, pft));
612  }
613  if (transposeB)
614  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
615  }
616 
617  if (C->getRowMap()->lib() == Xpetra::UseEpetra) {
618  throw Exceptions::RuntimeError("MatrixMatrix::Add for Epetra only available with Scalar = double, LO = GO = int.");
619  } else if (C->getRowMap()->lib() == Xpetra::UseTpetra) {
620 #ifdef HAVE_XPETRA_TPETRA
621  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpA =
623  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpB =
627 
628  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
629 #else
630  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
631 #endif
632  }
633 
635  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
636  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
638 
639  } //MatrixMatrix::TwoMatrixAdd()
640 
641 
642  }; // class MatrixMatrix
643 
644 
645 #ifdef HAVE_XPETRA_EPETRA
646  // specialization MatrixMatrix for SC=double, LO=GO=int
647  template<>
648  class MatrixMatrix<double,int,int,EpetraNode> {
649  typedef double SC;
650  typedef int LO;
651  typedef int GO;
652  typedef EpetraNode NO;
653 
657 
658  public:
659 
684  static void Multiply(
686  bool transposeA,
688  bool transposeB,
690  bool call_FillComplete_on_result = true,
691  bool doOptimizeStorage = true,
692  const std::string & label = std::string()) {
693  if(transposeA == false && C.getRowMap()->isSameAs(*A.getRowMap()) == false) {
694  std::string msg = "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as row map of A";
696  }
697  else if(transposeA == true && C.getRowMap()->isSameAs(*A.getDomainMap()) == false) {
698  std::string msg = "XpetraExt::MatrixMatrix::Multiply: row map of C is not same as domain map of A";
700  }
701 
702 
703  if (!A.isFillComplete())
704  throw(Xpetra::Exceptions::RuntimeError("A is not fill-completed"));
705  if (!B.isFillComplete())
706  throw(Xpetra::Exceptions::RuntimeError("B is not fill-completed"));
707 
708  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
709 
710  if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
711 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
712  Epetra_CrsMatrix & epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(A);
713  Epetra_CrsMatrix & epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
714  Epetra_CrsMatrix & epC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(C);
715 
716  int i = EpetraExt::MatrixMatrix::Multiply(epA,transposeA,epB,transposeB,epC,haveMultiplyDoFillComplete);
717  if (haveMultiplyDoFillComplete) {
718  // Due to Epetra wrapper intricacies, we need to explicitly call
719  // fillComplete on Xpetra matrix here. Specifically, EpetraCrsMatrix
720  // only keeps an internal variable to check whether we are in resumed
721  // state or not, but never touches the underlying Epetra object. As
722  // such, we need to explicitly update the state of Xpetra matrix to
723  // that of Epetra one afterwords
724  C.fillComplete();
725  }
726 
727  if (i != 0) {
728  std::ostringstream buf;
729  buf << i;
730  std::string msg = "EpetraExt::MatrixMatrix::Multiply return value of " + buf.str();
731  throw(Exceptions::RuntimeError(msg));
732  }
733 
734 #else
735  throw(Xpetra::Exceptions::RuntimeError("Xpetra::MatrixMatrix::Multiply requires EpetraExt to be compiled."));
736 #endif
737  } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
738 #ifdef HAVE_XPETRA_TPETRA
739 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
740  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
741  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
742 # else
743  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
744  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
745  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
746 
747  //18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
748  //Previously, Tpetra's matrix matrix multiply did not support fillComplete.
749  Tpetra::MatrixMatrix::Multiply(tpA,transposeA,tpB,transposeB,tpC,haveMultiplyDoFillComplete,label);
750 # endif
751 #else
752  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
753 #endif
754  }
755 
756  if(call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
758  params->set("Optimize Storage",doOptimizeStorage);
759  C.fillComplete((transposeB) ? B.getRangeMap() : B.getDomainMap(),
760  (transposeA) ? A.getDomainMap() : A.getRangeMap(),
761  params);
762  }
763 
764  // transfer striding information
765  RCP<Xpetra::Matrix<SC,LO,GO,NO> > rcpA = Teuchos::rcp_const_cast<Xpetra::Matrix<SC,LO,GO,NO> >(Teuchos::rcpFromRef(A));
766  RCP<Xpetra::Matrix<SC,LO,GO,NO> > rcpB = Teuchos::rcp_const_cast<Xpetra::Matrix<SC,LO,GO,NO> >(Teuchos::rcpFromRef(B));
767  C.CreateView("stridedMaps", rcpA, transposeA, rcpB, transposeB); // TODO use references instead of RCPs
768  } // end Multiply
769 
792  static RCP<Xpetra::Matrix<SC,LO,GO,NO> > Multiply(const Matrix& A, bool transposeA,
793  const Matrix& B, bool transposeB,
794  RCP<Matrix> C_in,
796  bool doFillComplete = true,
797  bool doOptimizeStorage = true,
798  const std::string & label = std::string()) {
799 
800  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
801  TEUCHOS_TEST_FOR_EXCEPTION(!B.isFillComplete(), Exceptions::RuntimeError, "B is not fill-completed");
802 
803  // Optimization using ML Multiply when available and requested
804  // This feature is currently not supported. We would have to introduce the HAVE_XPETRA_ML_MMM flag
805 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT) && defined(HAVE_XPETRA_ML_MMM)
806  if (B.getDomainMap()->lib() == Xpetra::UseEpetra && !transposeA && !transposeB) {
809  RCP<Epetra_CrsMatrix> epC = MLTwoMatrixMultiply(*epA, *epB, fos);
810 
811  RCP<Matrix> C = Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO> (epC);
812  if (doFillComplete) {
814  params->set("Optimize Storage", doOptimizeStorage);
815  C->fillComplete(B.getDomainMap(), A.getRangeMap(), params);
816  }
817 
818  // Fill strided maps information
819  // This is necessary since the ML matrix matrix multiplication routine has no handling for this
820  // TODO: move this call to MLMultiply...
821  C->CreateView("stridedMaps", rcpFromRef(A), transposeA, rcpFromRef(B), transposeB);
822 
823  return C;
824  }
825 #endif // EPETRA + EPETRAEXT + ML
826 
827  // Default case: Xpetra Multiply
828  RCP<Matrix> C = C_in;
829 
830  if (C == Teuchos::null) {
831  double nnzPerRow = Teuchos::as<double>(0);
832 
833  if (A.getDomainMap()->lib() == Xpetra::UseTpetra) {
834  // For now, follow what ML and Epetra do.
835  GO numRowsA = A.getGlobalNumRows();
836  GO numRowsB = B.getGlobalNumRows();
837  nnzPerRow = sqrt(Teuchos::as<double>(A.getGlobalNumEntries())/numRowsA) +
838  sqrt(Teuchos::as<double>(B.getGlobalNumEntries())/numRowsB) - 1;
839  nnzPerRow *= nnzPerRow;
840  double totalNnz = nnzPerRow * A.getGlobalNumRows() * 0.75 + 100;
841  double minNnz = Teuchos::as<double>(1.2 * A.getGlobalNumEntries());
842  if (totalNnz < minNnz)
843  totalNnz = minNnz;
844  nnzPerRow = totalNnz / A.getGlobalNumRows();
845 
846  fos << "Matrix product nnz per row estimate = " << Teuchos::as<LO>(nnzPerRow) << std::endl;
847  }
848 
849  if (transposeA) C = Xpetra::MatrixFactory<SC,LO,GO,NO>::Build(A.getDomainMap(), Teuchos::as<LO>(nnzPerRow));
850  else C = Xpetra::MatrixFactory<SC,LO,GO,NO>::Build(A.getRowMap(), Teuchos::as<LO>(nnzPerRow));
851 
852  } else {
853  C->resumeFill(); // why this is not done inside of Tpetra MxM?
854  fos << "Reuse C pattern" << std::endl;
855  }
856 
857  Multiply(A, transposeA, B, transposeB, *C, doFillComplete, doOptimizeStorage,label); // call Multiply routine from above
858 
859  return C;
860  }
861 
872  static RCP<Xpetra::Matrix<SC,LO,GO,NO> > Multiply(const Matrix & A,
873  bool transposeA,
874  const Matrix & B,
875  bool transposeB,
877  bool callFillCompleteOnResult = true,
878  bool doOptimizeStorage = true,
879  const std::string & label = std::string()){
880  return Multiply(A, transposeA, B, transposeB, Teuchos::null, fos, callFillCompleteOnResult, doOptimizeStorage,label);
881  }
882 
883 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
884  // Michael Gee's MLMultiply
885  static RCP<Epetra_CrsMatrix> MLTwoMatrixMultiply(const Epetra_CrsMatrix& epA,
886  const Epetra_CrsMatrix& epB,
887  Teuchos::FancyOStream& fos) {
888 #if defined(HAVE_XPETRA_ML_MMM) // Note: this is currently not supported
889  ML_Comm* comm;
890  ML_Comm_Create(&comm);
891  fos << "****** USING ML's MATRIX MATRIX MULTIPLY ******" << std::endl;
892 #ifdef HAVE_MPI
893  // ML_Comm uses MPI_COMM_WORLD, so try to use the same communicator as epA.
894  const Epetra_MpiComm * Mcomm = dynamic_cast<const Epetra_MpiComm*>(&(epA.Comm()));
895  if (Mcomm)
896  ML_Comm_Set_UsrComm(comm,Mcomm->GetMpiComm());
897 # endif
898  //in order to use ML, there must be no indices missing from the matrix column maps.
899  EpetraExt::CrsMatrix_SolverMap Atransform;
900  EpetraExt::CrsMatrix_SolverMap Btransform;
901  const Epetra_CrsMatrix& A = Atransform(const_cast<Epetra_CrsMatrix&>(epA));
902  const Epetra_CrsMatrix& B = Btransform(const_cast<Epetra_CrsMatrix&>(epB));
903 
904  if (!A.Filled()) throw Exceptions::RuntimeError("A has to be FillCompleted");
905  if (!B.Filled()) throw Exceptions::RuntimeError("B has to be FillCompleted");
906 
907  // create ML operators from EpetraCrsMatrix
908  ML_Operator* ml_As = ML_Operator_Create(comm);
909  ML_Operator* ml_Bs = ML_Operator_Create(comm);
910  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?
911  ML_Operator_WrapEpetraCrsMatrix(const_cast<Epetra_CrsMatrix*>(&B),ml_Bs);
912  ML_Operator* ml_AtimesB = ML_Operator_Create(comm);
913  {
914  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("ML_2matmult kernel"));
915  ML_2matmult(ml_As,ml_Bs,ml_AtimesB,ML_CSR_MATRIX); // do NOT use ML_EpetraCRS_MATRIX!!!
916  }
917  ML_Operator_Destroy(&ml_As);
918  ML_Operator_Destroy(&ml_Bs);
919 
920  // For ml_AtimesB we have to reconstruct the column map in global indexing,
921  // The following is going down to the salt-mines of ML ...
922  // note: we use integers, since ML only works for Epetra...
923  int N_local = ml_AtimesB->invec_leng;
924  ML_CommInfoOP* getrow_comm = ml_AtimesB->getrow->pre_comm;
925  if (!getrow_comm) throw(Exceptions::RuntimeError("ML_Operator does not have a CommInfo"));
926  ML_Comm* comm_AB = ml_AtimesB->comm; // get comm object
927  if (N_local != B.DomainMap().NumMyElements())
928  throw(Exceptions::RuntimeError("Mismatch in local row dimension between ML and Epetra"));
929  int N_rcvd = 0;
930  int N_send = 0;
931  int flag = 0;
932  for (int i = 0; i < getrow_comm->N_neighbors; i++) {
933  N_rcvd += (getrow_comm->neighbors)[i].N_rcv;
934  N_send += (getrow_comm->neighbors)[i].N_send;
935  if ( ((getrow_comm->neighbors)[i].N_rcv != 0) &&
936  ((getrow_comm->neighbors)[i].rcv_list != NULL) ) flag = 1;
937  }
938  // For some unknown reason, ML likes to have stuff one larger than
939  // neccessary...
940  std::vector<double> dtemp(N_local + N_rcvd + 1); // "double" vector for comm function
941  std::vector<int> cmap (N_local + N_rcvd + 1); // vector for gids
942  for (int i = 0; i < N_local; ++i) {
943  cmap[i] = B.DomainMap().GID(i);
944  dtemp[i] = (double) cmap[i];
945  }
946  ML_cheap_exchange_bdry(&dtemp[0],getrow_comm,N_local,N_send,comm_AB); // do communication
947  if (flag) { // process received data
948  int count = N_local;
949  const int neighbors = getrow_comm->N_neighbors;
950  for (int i = 0; i < neighbors; i++) {
951  const int nrcv = getrow_comm->neighbors[i].N_rcv;
952  for (int j = 0; j < nrcv; j++)
953  cmap[getrow_comm->neighbors[i].rcv_list[j]] = (int) dtemp[count++];
954  }
955  } else {
956  for (int i = 0; i < N_local+N_rcvd; ++i)
957  cmap[i] = (int)dtemp[i];
958  }
959  dtemp.clear(); // free double array
960 
961  // we can now determine a matching column map for the result
962  Epetra_Map gcmap(-1, N_local+N_rcvd, &cmap[0], B.ColMap().IndexBase(), A.Comm());
963 
964  int allocated = 0;
965  int rowlength;
966  double* val = NULL;
967  int* bindx = NULL;
968 
969  const int myrowlength = A.RowMap().NumMyElements();
970  const Epetra_Map& rowmap = A.RowMap();
971 
972  // Determine the maximum bandwith for the result matrix.
973  // replaces the old, very(!) memory-consuming guess:
974  // int guessnpr = A.MaxNumEntries()*B.MaxNumEntries();
975  int educatedguess = 0;
976  for (int i = 0; i < myrowlength; ++i) {
977  // get local row
978  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
979  if (rowlength>educatedguess)
980  educatedguess = rowlength;
981  }
982 
983  // allocate our result matrix and fill it
984  RCP<Epetra_CrsMatrix> result = rcp(new Epetra_CrsMatrix(::Copy, A.RangeMap(), gcmap, educatedguess, false));
985 
986  std::vector<int> gcid(educatedguess);
987  for (int i = 0; i < myrowlength; ++i) {
988  const int grid = rowmap.GID(i);
989  // get local row
990  ML_get_matrix_row(ml_AtimesB, 1, &i, &allocated, &bindx, &val, &rowlength, 0);
991  if (!rowlength) continue;
992  if ((int)gcid.size() < rowlength) gcid.resize(rowlength);
993  for (int j = 0; j < rowlength; ++j) {
994  gcid[j] = gcmap.GID(bindx[j]);
995  if (gcid[j] < 0)
996  throw Exceptions::RuntimeError("Error: cannot find gcid!");
997  }
998  int err = result->InsertGlobalValues(grid, rowlength, val, &gcid[0]);
999  if (err != 0 && err != 1) {
1000  std::ostringstream errStr;
1001  errStr << "Epetra_CrsMatrix::InsertGlobalValues returned err=" << err;
1002  throw Exceptions::RuntimeError(errStr.str());
1003  }
1004  }
1005  // free memory
1006  if (bindx) ML_free(bindx);
1007  if (val) ML_free(val);
1008  ML_Operator_Destroy(&ml_AtimesB);
1009  ML_Comm_Destroy(&comm);
1010 
1011  return result;
1012 #else // no MUELU_ML
1014  "No ML multiplication available. This feature is currently not supported by Xpetra.");
1015  return Teuchos::null;
1016 #endif
1017  }
1018 #endif //ifdef HAVE_XPETRA_EPETRAEXT
1019 
1031  Xpetra::BlockedCrsMatrix<SC,LO,GO,NO>& B, bool transposeB,
1032  Teuchos::FancyOStream& fos,
1033  bool doFillComplete = true,
1034  bool doOptimizeStorage = true) {
1035  if (transposeA || transposeB)
1036  throw Exceptions::RuntimeError("TwoMatrixMultiply for BlockedCrsMatrix not implemented for transposeA==true or transposeB==true");
1037 
1038  // Preconditions
1039  if (!A.isFillComplete())
1040  throw Exceptions::RuntimeError("A is not fill-completed");
1041  if (!B.isFillComplete())
1042  throw Exceptions::RuntimeError("B is not fill-completed");
1043 
1046 
1048  rcp(new Xpetra::BlockedCrsMatrix<SC,LO,GO,NO>(rgmapextractor, domapextractor, 33 /* TODO fix me */));
1049 
1050  for (size_t i = 0; i < A.Rows(); ++i) { // loop over all block rows of A
1051  for (size_t j = 0; j < B.Cols(); ++j) { // loop over all block columns of B
1052  RCP<Matrix> Cij = Teuchos::null;
1053 
1054  for (size_t l = 0; l < B.Rows(); ++l) { // loop for calculating entry C_{ij}
1057 
1058  if (crmat1.is_null() || crmat2.is_null()) {
1059  continue;
1060  }
1061 
1064 
1065  RCP<Xpetra::Matrix<SC,LO,GO,NO> > temp = Multiply (*crop1, false, *crop2, false, fos);
1066 
1067  RCP<Matrix> addRes = Teuchos::null;
1068  if (Cij.is_null ())
1069  Cij = temp;
1070  else {
1071  Xpetra::MatrixMatrix<SC,LO,GO,NO>::TwoMatrixAdd (*temp, false, 1.0, *Cij, false, 1.0, addRes, fos);
1072  Cij = addRes;
1073  }
1074  }
1075 
1076  if (!Cij.is_null()) {
1077  if (Cij->isFillComplete())
1078  Cij->resumeFill();
1079  Cij->fillComplete(B.getDomainMap(j), A.getRangeMap(i));
1080 
1081  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > crsCij = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(Cij);
1083  "MatrixFactory failed in generating a CrsMatrixWrap." );
1084 
1085  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > crsMatCij = crsCij->getCrsMatrix();
1086 
1087  C->setMatrix(i, j, crsMatCij);
1088 
1089  } else {
1090  C->setMatrix(i, j, Teuchos::null);
1091  }
1092  }
1093  }
1094 
1095  if (doFillComplete)
1096  C->fillComplete(); // call default fillComplete for BlockCrsMatrixWrap objects
1097 
1098  return C;
1099  } // TwoMatrixMultiplyBlock
1100 
1113  static void TwoMatrixAdd(const Xpetra::Matrix<SC,LO,GO,NO>& A, bool transposeA, SC alpha, Xpetra::Matrix<SC,LO,GO,NO>& B, SC beta) {
1114  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
1115  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
1116 
1117  if (A.getRowMap()->lib() == Xpetra::UseEpetra) {
1118 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1119  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1120  Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstEpetraCrs(B);
1121 
1122  //FIXME is there a bug if beta=0?
1123  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, beta);
1124 
1125  if (rv != 0)
1126  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value " + Teuchos::toString(rv));
1127  std::ostringstream buf;
1128 #else
1129  throw Exceptions::RuntimeError("Xpetra must be compiled with EpetraExt.");
1130 #endif
1131  } else if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
1132 #ifdef HAVE_XPETRA_TPETRA
1133 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1134  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1135  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1136 # else
1137  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
1138  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
1139 
1140  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, beta);
1141 # endif
1142 #else
1143  throw Exceptions::RuntimeError("Xpetra must be compiled with Tpetra.");
1144 #endif
1145  }
1146  } //MatrixMatrix::TwoMatrixAdd()
1147 
1148 
1161  static void TwoMatrixAdd(const Xpetra::Matrix<SC,LO,GO,NO>& A, bool transposeA, const SC& alpha,
1162  const Xpetra::Matrix<SC,LO,GO,NO>& B, bool transposeB, const SC& beta,
1163  RCP<Xpetra::Matrix<SC,LO,GO,NO> >& C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow = false) {
1164  if (!(A.getRowMap()->isSameAs(*(B.getRowMap()))))
1165  throw Exceptions::Incompatible("TwoMatrixAdd: matrix row maps are not the same.");
1166 
1167  if (C == Teuchos::null) {
1168  if (!A.isFillComplete() || !B.isFillComplete())
1169  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Global statistics are not available for estimates.");
1170 
1171  size_t maxNzInA = A.getGlobalMaxNumRowEntries();
1172  size_t maxNzInB = B.getGlobalMaxNumRowEntries();
1173  size_t numLocalRows = A.getNodeNumRows();
1174 
1175  if (maxNzInA == 1 || maxNzInB == 1 || AHasFixedNnzPerRow) {
1176  // first check if either A or B has at most 1 nonzero per row
1177  // the case of both having at most 1 nz per row is handled by the ``else''
1178  Teuchos::ArrayRCP<size_t> exactNnzPerRow(numLocalRows);
1179 
1180  if ((maxNzInA == 1 && maxNzInB > 1) || AHasFixedNnzPerRow) {
1181  for (size_t i = 0; i < numLocalRows; ++i)
1182  exactNnzPerRow[i] = B.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInA;
1183 
1184  } else {
1185  for (size_t i = 0; i < numLocalRows; ++i)
1186  exactNnzPerRow[i] = A.getNumEntriesInLocalRow(Teuchos::as<LO>(i)) + maxNzInB;
1187  }
1188 
1189  fos << "MatrixMatrix::TwoMatrixAdd : special case detected (one matrix has a fixed nnz per row)"
1190  << ", using static profiling" << std::endl;
1192 
1193  } else {
1194  // general case
1195  double nnzPerRowInA = Teuchos::as<double>(A.getGlobalNumEntries()) / A.getGlobalNumRows();
1196  double nnzPerRowInB = Teuchos::as<double>(B.getGlobalNumEntries()) / B.getGlobalNumRows();
1197  LO nnzToAllocate = Teuchos::as<LO>( (nnzPerRowInA + nnzPerRowInB) * 1.5) + Teuchos::as<LO>(1);
1198 
1199  LO maxPossible = A.getGlobalMaxNumRowEntries() + B.getGlobalMaxNumRowEntries();
1200  //Use static profiling (more efficient) if the estimate is at least as big as the max
1201  //possible nnz's in any single row of the result.
1202  Xpetra::ProfileType pft = (maxPossible) > nnzToAllocate ? Xpetra::DynamicProfile : Xpetra::StaticProfile;
1203 
1204  fos << "nnzPerRowInA = " << nnzPerRowInA << ", nnzPerRowInB = " << nnzPerRowInB << std::endl;
1205  fos << "MatrixMatrix::TwoMatrixAdd : space allocated per row = " << nnzToAllocate
1206  << ", max possible nnz per row in sum = " << maxPossible
1207  << ", using " << (pft == Xpetra::DynamicProfile ? "dynamic" : "static" ) << " profiling"
1208  << std::endl;
1209  C = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(A.getRowMap(), nnzToAllocate, pft));
1210  }
1211  if (transposeB)
1212  fos << "MatrixMatrix::TwoMatrixAdd : ** WARNING ** estimate could be badly wrong because second summand is transposed" << std::endl;
1213  }
1214 
1215  if (C->getRowMap()->lib() == Xpetra::UseEpetra) {
1216 #if defined(HAVE_XPETRA_EPETRA) && defined(HAVE_XPETRA_EPETRAEXT)
1217  const Epetra_CrsMatrix& epA = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(A);
1218  const Epetra_CrsMatrix& epB = Xpetra::Helpers<SC,LO,GO,NO>::Op2EpetraCrs(B);
1220  Epetra_CrsMatrix* ref2epC = &*epC; //to avoid a compiler error...
1221 
1222  //FIXME is there a bug if beta=0?
1223  int rv = EpetraExt::MatrixMatrix::Add(epA, transposeA, alpha, epB, transposeB, beta, ref2epC);
1224 
1225  if (rv != 0)
1226  throw Exceptions::RuntimeError("EpetraExt::MatrixMatrix::Add return value of " + Teuchos::toString(rv));
1227 #else
1228  throw Exceptions::RuntimeError("MueLu must be compile with EpetraExt.");
1229 #endif
1230  } else if (C->getRowMap()->lib() == Xpetra::UseTpetra) {
1231 #ifdef HAVE_XPETRA_TPETRA
1232 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
1233  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
1234  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=int enabled."));
1235 # else
1236  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA =
1238  const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB =
1242 
1243  Tpetra::MatrixMatrix::Add(tpA, transposeA, alpha, tpB, transposeB, beta, tpC);
1244 # endif
1245 #else
1246  throw Exceptions::RuntimeError("Xpetra must be compile with Tpetra.");
1247 #endif
1248  }
1249 
1251  if (A.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(A));
1252  if (B.IsView("stridedMaps")) C->CreateView("stridedMaps", rcpFromRef(B));
1254 
1255  } //MatrixMatrix::TwoMatrixAdd()
1256  };
1257 #endif // HAVE_XPETRA_EPETRA
1258 
1259 } // end namespace Xpetra
1260 
1261 #define XPETRA_MATRIXMATRIX_SHORT
1262 
1263 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_MATRIXMATRIX_HPP_ */
virtual size_t getNodeNumRows() const =0
Returns the number of matrix rows owned on the calling node.
Xpetra::MultiVector< double, int, int, NO > MultiVector
RCP< const MapExtractor > getRangeMapExtractor()
Returns map extractor class for range map.
static void TwoMatrixAdd(const Matrix &A, bool transposeA, const SC &alpha, const Matrix &B, bool transposeB, const SC &beta, RCP< Matrix > &C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow=false)
Helper function to calculate C = alpha*A + beta*B.
static void Multiply(const Xpetra::Matrix< SC, LO, GO, NO > &A, bool transposeA, const Xpetra::Matrix< SC, LO, GO, NO > &B, bool transposeB, Xpetra::Matrix< SC, LO, GO, NO > &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string())
RCP< CrsMatrix > getCrsMatrix() const
virtual global_size_t getGlobalNumRows() const =0
Returns the number of global rows in this matrix.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > 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())
Helper function to do matrix-matrix multiply.
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
GlobalOrdinal GO
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static const Epetra_CrsMatrix & Op2EpetraCrs(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
RCP< const Map > getRangeMap() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Xpetra utility class containing transformation routines between Xpetra::Matrix and Epetra/Tpetra obje...
Xpetra namespace
Teuchos::RCP< CrsMatrix > getMatrix(size_t r, size_t c) const
return block (r,c)
static Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > & Op2NonConstTpetraCrs(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Exception throws to report errors in the internal logical of the program.
LocalOrdinal LO
virtual size_t getGlobalMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on all nodes.
static void Multiply(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool transposeA, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B, bool transposeB, Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string())
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
static RCP< Xpetra::Matrix< SC, LO, GO, NO > > 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())
Helper function to do matrix-matrix multiply.
static RCP< Time > getNewTimer(const std::string &name)
Exception indicating invalid cast attempted.
static RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(BlockedCrsMatrix &A, bool transposeA, BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool IsView(const viewLabel_t viewLabel) const
RCP< const MapExtractor > getDomainMapExtractor()
Returns map extractor for domain map.
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
virtual size_t Cols() const
number of column blocks
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
Returns the current number of entries on this node in the specified local row.
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
static void TwoMatrixAdd(const Xpetra::Matrix< SC, LO, GO, NO > &A, bool transposeA, SC alpha, Xpetra::Matrix< SC, LO, GO, NO > &B, SC beta)
Helper function to calculate B = alpha*A + beta*B.
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Xpetra::BlockedCrsMatrix< SC, LO, GO, NO > > TwoMatrixMultiplyBlock(Xpetra::BlockedCrsMatrix< SC, LO, GO, NO > &A, bool transposeA, Xpetra::BlockedCrsMatrix< SC, LO, GO, NO > &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
Helper function to do matrix-matrix multiply "in-place".
static void TwoMatrixAdd(const Xpetra::Matrix< SC, LO, GO, NO > &A, bool transposeA, const SC &alpha, const Xpetra::Matrix< SC, LO, GO, NO > &B, bool transposeB, const SC &beta, RCP< Xpetra::Matrix< SC, LO, GO, NO > > &C, Teuchos::FancyOStream &fos, bool AHasFixedNnzPerRow=false)
Helper function to calculate C = alpha*A + beta*B.
Exception throws to report incompatible objects (like maps).
RCP< const Map > getDomainMap() const
Returns the Map associated with the full domain of this operator. This will be null until fillComplet...
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > 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())
Helper function to do matrix-matrix multiply.
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Copy
virtual global_size_t getGlobalNumEntries() const =0
Returns the global number of entries in this matrix.
Concrete implementation of Xpetra::Matrix.
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
virtual size_t Rows() const
number of row blocks
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor specifying the number of non-zeros for all rows.
static RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
static RCP< Epetra_CrsMatrix > MLTwoMatrixMultiply(const Epetra_CrsMatrix &epA, const Epetra_CrsMatrix &epB, Teuchos::FancyOStream &fos)
static RCP< Xpetra::Matrix< SC, LO, GO, NO > > 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())
Helper function to do matrix-matrix multiply.
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
static const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > & Op2TpetraCrs(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Xpetra-specific matrix class.
std::string toString(const T &t)
bool is_null() const