Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Xpetra_TripleMatrixMultiply.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_TRIPLEMATRIXMULTIPLY_HPP_
49 #define PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_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 #include "Xpetra_IO.hpp"
62 
63 #ifdef HAVE_XPETRA_TPETRA
64 #include <TpetraExt_TripleMatrixMultiply.hpp>
65 #include <Xpetra_TpetraCrsMatrix.hpp>
66 #include <Tpetra_BlockCrsMatrix.hpp>
67 #include <Tpetra_BlockCrsMatrix_Helpers.hpp>
68 // #include <Xpetra_TpetraMultiVector.hpp>
69 // #include <Xpetra_TpetraVector.hpp>
70 #endif // HAVE_XPETRA_TPETRA
71 
72 namespace Xpetra {
73 
74 template <class Scalar,
75  class LocalOrdinal /*= int*/,
76  class GlobalOrdinal /*= LocalOrdinal*/,
77  class Node /*= Tpetra::KokkosClassic::DefaultNode::DefaultNodeType*/>
78 class TripleMatrixMultiply {
79 #undef XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
80 #include "Xpetra_UseShortNames.hpp"
81 
82  public:
107  static void MultiplyRAP(const Matrix& R, bool transposeR,
108  const Matrix& A, bool transposeA,
109  const Matrix& P, bool transposeP,
110  Matrix& Ac,
111  bool call_FillComplete_on_result = true,
112  bool doOptimizeStorage = true,
113  const std::string& label = std::string(),
114  const RCP<ParameterList>& params = null) {
115  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
116  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
117  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
118  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
119 
120  TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
121  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
122  TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
123 
124  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
125 
126  if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
127  throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
128  } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
129 #ifdef HAVE_XPETRA_TPETRA
130  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
131  if (helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P)) {
132  // All matrices are Crs
133  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(R);
134  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
135  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(P);
136  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpAc = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(Ac);
137 
138  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
139  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
140  Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
141  } else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
142  // All matrices are BlockCrs (except maybe Ac)
143  // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
144  // we'll need to think about refactoring BlockCrs so we can do something smarter here.
145 
146  if (!A.getRowMap()->getComm()->getRank())
147  std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
148 
149  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(R);
150  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
151  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(P);
152  // Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraBlockCrs(Ac);
153 
154  using CRS = Tpetra::CrsMatrix<SC, LO, GO, NO>;
155  RCP<const CRS> Rcrs = Tpetra::convertToCrsMatrix(tpR);
156  RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
157  RCP<const CRS> Pcrs = Tpetra::convertToCrsMatrix(tpP);
158  // RCP<CRS> Accrs = Tpetra::convertToCrsMatrix(tpAc);
159 
160  // FIXME: The lines below only works because we're assuming Ac is Point
161  RCP<CRS> Accrs = Teuchos::rcp(new CRS(Rcrs->getRowMap(), 0));
162  const bool do_fill_complete = true;
163  Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
164 
165  // Temporary output matrix
166  RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*Accrs, A.GetStorageBlockSize());
167  RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
168  RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > Ac_p = Ac_x;
169 
170  // We can now cheat and replace the innards of Ac
171  RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(Teuchos::rcpFromRef(Ac));
172  Ac_w->replaceCrsMatrix(Ac_p);
173  } else {
174  // Mix and match
175  TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
176  }
177 #else
178  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
179 #endif
180  }
181 
182  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
183  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
184  fillParams->set("Optimize Storage", doOptimizeStorage);
185  Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
186  (transposeR) ? R.getDomainMap() : R.getRangeMap(),
187  fillParams);
188  }
189 
190  // transfer striding information
191  RCP<const Map> domainMap = Teuchos::null;
192  RCP<const Map> rangeMap = Teuchos::null;
193 
194  const std::string stridedViewLabel("stridedMaps");
195  const size_t blkSize = 1;
196  std::vector<size_t> stridingInfo(1, blkSize);
197  LocalOrdinal stridedBlockId = -1;
198 
199  if (R.IsView(stridedViewLabel)) {
200  rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
201  } else {
202  rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
203  rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
204  }
205 
206  if (P.IsView(stridedViewLabel)) {
207  domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
208  } else {
209  domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
210  domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
211  }
212  Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
213 
214  } // end Multiply
215 
216 }; // class TripleMatrixMultiply
217 
218 #ifdef HAVE_XPETRA_EPETRA
219 // specialization TripleMatrixMultiply for SC=double, LO=GO=int
220 template <>
221 class TripleMatrixMultiply<double, int, int, EpetraNode> {
222  typedef double Scalar;
223  typedef int LocalOrdinal;
224  typedef int GlobalOrdinal;
225  typedef EpetraNode Node;
226 #include "Xpetra_UseShortNames.hpp"
227 
228  public:
229  static void MultiplyRAP(const Matrix& R, bool transposeR,
230  const Matrix& A, bool transposeA,
231  const Matrix& P, bool transposeP,
232  Matrix& Ac,
233  bool call_FillComplete_on_result = true,
234  bool doOptimizeStorage = true,
235  const std::string& label = std::string(),
236  const RCP<ParameterList>& params = null) {
237  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
238  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
239  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
240  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
241 
242  TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
243  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
244  TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
245 
246  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
247 
248  if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
249  throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
250  } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
251 #ifdef HAVE_XPETRA_TPETRA
252 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
253  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
254  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
255 #else
256  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
257  if (helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P) && helpers::isTpetraCrs(Ac)) {
258  // All matrices are Crs
259  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(R);
260  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
261  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(P);
262  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpAc = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(Ac);
263 
264  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
265  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
266  Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
267  } else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
268  // All matrices are BlockCrs (except maybe Ac)
269  // FIXME: For the moment we're just going to clobber the innards of AC, so no reuse. Once we have a reuse kernel,
270  // we'll need to think about refactoring BlockCrs so we can do something smarter here.
271  if (!A.getRowMap()->getComm()->getRank())
272  std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
273 
274  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(R);
275  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
276  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(P);
277  // Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraBlockCrs(Ac);
278 
279  using CRS = Tpetra::CrsMatrix<SC, LO, GO, NO>;
280  RCP<const CRS> Rcrs = Tpetra::convertToCrsMatrix(tpR);
281  RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
282  RCP<const CRS> Pcrs = Tpetra::convertToCrsMatrix(tpP);
283  // RCP<CRS> Accrs = Tpetra::convertToCrsMatrix(tpAc);
284 
285  // FIXME: The lines below only works because we're assuming Ac is Point
286  RCP<CRS> Accrs = Teuchos::rcp(new CRS(Rcrs->getRowMap(), 0));
287  const bool do_fill_complete = true;
288  Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
289 
290  // Temporary output matrix
291  RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*Accrs, A.GetStorageBlockSize());
292  RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
293  RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > Ac_p = Ac_x;
294 
295  // We can now cheat and replace the innards of Ac
296  RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(Teuchos::rcpFromRef(Ac));
297  Ac_w->replaceCrsMatrix(Ac_p);
298 
299  } else {
300  // Mix and match (not supported)
301  TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
302  }
303 #endif
304 #else
305  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
306 #endif
307  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
308  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
309  fillParams->set("Optimize Storage", doOptimizeStorage);
310  Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
311  (transposeR) ? R.getDomainMap() : R.getRangeMap(),
312  fillParams);
313  }
314 
315  // transfer striding information
316  RCP<const Map> domainMap = Teuchos::null;
317  RCP<const Map> rangeMap = Teuchos::null;
318 
319  const std::string stridedViewLabel("stridedMaps");
320  const size_t blkSize = 1;
321  std::vector<size_t> stridingInfo(1, blkSize);
322  LocalOrdinal stridedBlockId = -1;
323 
324  if (R.IsView(stridedViewLabel)) {
325  rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
326  } else {
327  rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
328  rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
329  }
330 
331  if (P.IsView(stridedViewLabel)) {
332  domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
333  } else {
334  domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
335  domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
336  }
337  Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
338  }
339 
340  } // end Multiply
341 
342 }; // end specialization on SC=double, GO=int and NO=EpetraNode
343 
344 // specialization TripleMatrixMultiply for SC=double, GO=long long and NO=EpetraNode
345 template <>
346 class TripleMatrixMultiply<double, int, long long, EpetraNode> {
347  typedef double Scalar;
348  typedef int LocalOrdinal;
349  typedef long long GlobalOrdinal;
350  typedef EpetraNode Node;
351 #include "Xpetra_UseShortNames.hpp"
352 
353  public:
354  static void MultiplyRAP(const Matrix& R, bool transposeR,
355  const Matrix& A, bool transposeA,
356  const Matrix& P, bool transposeP,
357  Matrix& Ac,
358  bool call_FillComplete_on_result = true,
359  bool doOptimizeStorage = true,
360  const std::string& label = std::string(),
361  const RCP<ParameterList>& params = null) {
362  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
363  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
364  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
365  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
366 
367  TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
368  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
369  TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
370 
371  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
372 
373  if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
374  throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
375  } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
376 #ifdef HAVE_XPETRA_TPETRA
377 #if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
378  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
379  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long,EpetraNode> ETI enabled."));
380 #else
381  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
382  if (helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P)) {
383  // All matrices are Crs
384  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(R);
385  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
386  const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(P);
387  Tpetra::CrsMatrix<SC, LO, GO, NO>& tpAc = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(Ac);
388 
389  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
390  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
391  Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
392  } else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
393  // All matrices are BlockCrs (except maybe Ac)
394  // FIXME: For the moment we're just going to clobber the innards of AC, so no reuse. Once we have a reuse kernel,
395  // we'll need to think about refactoring BlockCrs so we can do something smarter here.
396  if (!A.getRowMap()->getComm()->getRank())
397  std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
398 
399  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(R);
400  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
401  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(P);
402  // Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraBlockCrs(Ac);
403 
404  using CRS = Tpetra::CrsMatrix<SC, LO, GO, NO>;
405  RCP<const CRS> Rcrs = Tpetra::convertToCrsMatrix(tpR);
406  RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
407  RCP<const CRS> Pcrs = Tpetra::convertToCrsMatrix(tpP);
408  // RCP<CRS> Accrs = Tpetra::convertToCrsMatrix(tpAc);
409 
410  // FIXME: The lines below only works because we're assuming Ac is Point
411  RCP<CRS> Accrs = Teuchos::rcp(new CRS(Rcrs->getRowMap(), 0));
412  const bool do_fill_complete = true;
413  Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
414 
415  // Temporary output matrix
416  RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*Accrs, A.GetStorageBlockSize());
417  RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
418  RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > Ac_p = Ac_x;
419 
420  // We can now cheat and replace the innards of Ac
421  RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(Teuchos::rcpFromRef(Ac));
422  Ac_w->replaceCrsMatrix(Ac_p);
423  } else {
424  // Mix and match
425  TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
426  }
427 
428 #endif
429 #else
430  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
431 #endif
432  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
433  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
434  fillParams->set("Optimize Storage", doOptimizeStorage);
435  Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
436  (transposeR) ? R.getDomainMap() : R.getRangeMap(),
437  fillParams);
438  }
439 
440  // transfer striding information
441  RCP<const Map> domainMap = Teuchos::null;
442  RCP<const Map> rangeMap = Teuchos::null;
443 
444  const std::string stridedViewLabel("stridedMaps");
445  const size_t blkSize = 1;
446  std::vector<size_t> stridingInfo(1, blkSize);
447  LocalOrdinal stridedBlockId = -1;
448 
449  if (R.IsView(stridedViewLabel)) {
450  rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
451  } else {
452  rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
453  rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
454  }
455 
456  if (P.IsView(stridedViewLabel)) {
457  domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
458  } else {
459  domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
460  domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
461  }
462  Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
463  }
464 
465  } // end Multiply
466 
467 }; // end specialization on GO=long long and NO=EpetraNode
468 #endif
469 
470 } // end namespace Xpetra
471 
472 #define XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
473 
474 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_HPP_ */
Xpetra::TripleMatrixMultiply< double, int, long long, EpetraNode >::MultiplyRAP
static void MultiplyRAP(const Matrix &R, bool transposeR, const Matrix &A, bool transposeA, const Matrix &P, bool transposeP, Matrix &Ac, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Definition: Xpetra_TripleMatrixMultiply.hpp:354
Xpetra::Helpers::Op2TpetraCrs
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:162
Xpetra::TripleMatrixMultiply::MultiplyRAP
static void MultiplyRAP(const Matrix &R, bool transposeR, const Matrix &A, bool transposeA, const Matrix &P, bool transposeP, Matrix &Ac, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Definition: Xpetra_TripleMatrixMultiply.hpp:107
Xpetra_UseShortNames.hpp
Xpetra::TripleMatrixMultiply< double, int, int, EpetraNode >::MultiplyRAP
static void MultiplyRAP(const Matrix &R, bool transposeR, const Matrix &A, bool transposeA, const Matrix &P, bool transposeP, Matrix &Ac, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Definition: Xpetra_TripleMatrixMultiply.hpp:229
Xpetra::Helpers
Xpetra utility class containing transformation routines between Xpetra::Matrix and Epetra/Tpetra obje...
Definition: Xpetra_MatrixMatrix.hpp:95
Xpetra::Helpers::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:174
Xpetra::Matrix::GetStorageBlockSize
virtual LocalOrdinal GetStorageBlockSize() const =0
Returns the block size of the storage mechanism, which is usually 1, except for Tpetra::BlockCrsMatri...
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:101
Xpetra_Matrix.hpp
Xpetra::Matrix::getColMap
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
Definition: Xpetra_Matrix.hpp:325
Xpetra_IO.hpp
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::IsView
bool IsView(const viewLabel_t viewLabel) const
Definition: Xpetra_Matrix.hpp:202
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
Xpetra::Matrix::CreateView
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
Definition: Xpetra_Matrix.hpp:126
Xpetra::Operator::getRangeMap
virtual const Teuchos::RCP< const map_type > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
Xpetra_MatrixFactory.hpp
Xpetra::EpetraNode
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
Definition: Xpetra_Map_decl.hpp:71
Xpetra::CrsMatrixWrap::replaceCrsMatrix
void replaceCrsMatrix(RCP< CrsMatrix > &M)
Expert only.
Definition: Xpetra_CrsMatrixWrap_def.hpp:495
Xpetra_ConfigDefs.hpp
Xpetra::Helpers::Op2TpetraBlockCrs
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
Definition: Xpetra_MatrixMatrix.hpp:240
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_fwd.hpp:51
Xpetra::Operator::getDomainMap
virtual const Teuchos::RCP< const map_type > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra::StridedMapFactory::Build
static RCP< Xpetra::StridedMap< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
Definition: Xpetra_StridedMapFactory_def.hpp:61
Xpetra::Matrix::getRowMap
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Definition: Xpetra_Matrix.hpp:315
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_fwd.hpp:51