Xpetra Version of the Day
Loading...
Searching...
No Matches
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
72namespace 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:
83
95 FillComplete() has already been called on C or not. If it has,
96 then C's graph must already contain all nonzero locations that
97 will be produced when forming the product A*B. On exit,
98 C.FillComplete() will have been called, unless the last argument
99 to this function is specified to be false.
100 @param call_FillComplete_on_result Optional argument, defaults to true.
101 Power users may specify this argument to be false if they *DON'T*
102 want this function to call C.FillComplete. (It is often useful
103 to allow this function to call C.FillComplete, in cases where
104 one or both of the input matrices are rectangular and it is not
105 trivial to know which maps to use for the domain- and range-maps.)
106
107*/
108 static void MultiplyRAP(const Matrix& R, bool transposeR,
109 const Matrix& A, bool transposeA,
110 const Matrix& P, bool transposeP,
111 Matrix& Ac,
112 bool call_FillComplete_on_result = true,
113 bool doOptimizeStorage = true,
114 const std::string & label = std::string(),
115 const RCP<ParameterList>& params = null) {
116
117 TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
118 Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
119 TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
120 Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
121
122 TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
123 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
124 TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
125
126 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
127
128 if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
129 throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
130 } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
131#ifdef HAVE_XPETRA_TPETRA
132 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
133 if(helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P)) {
134 // All matrices are Crs
135 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(R);
136 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
137 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(P);
138 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(Ac);
139
140 // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
141 // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
142 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
143 }
144 else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
145 // All matrices are BlockCrs (except maybe Ac)
146 // FIXME: For the moment we're just going to clobber the innards of Ac, so no reuse. Once we have a reuse kernel,
147 // we'll need to think about refactoring BlockCrs so we can do something smarter here.
148
149 if(!A.getRowMap()->getComm()->getRank())
150 std::cout<<"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
151
152 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(R);
153 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(A);
154 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(P);
155 // Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraBlockCrs(Ac);
156
157 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
158 RCP<const CRS> Rcrs = Tpetra::convertToCrsMatrix(tpR);
159 RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
160 RCP<const CRS> Pcrs = Tpetra::convertToCrsMatrix(tpP);
161 // RCP<CRS> Accrs = Tpetra::convertToCrsMatrix(tpAc);
162
163 // FIXME: The lines below only works because we're assuming Ac is Point
164 RCP<CRS> Accrs = Teuchos::rcp(new CRS(Rcrs->getRowMap(),0));
165 const bool do_fill_complete=true;
166 Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
167
168 // Temporary output matrix
169 RCP<Tpetra::BlockCrsMatrix<SC,LO,GO,NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*Accrs,A.GetStorageBlockSize());
170 RCP<Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO>(Ac_t));
171 RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > Ac_p = Ac_x;
172
173 // We can now cheat and replace the innards of Ac
174 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(Teuchos::rcpFromRef(Ac));
175 Ac_w->replaceCrsMatrix(Ac_p);
176 }
177 else {
178 // Mix and match
179 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
180 }
181#else
182 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
183#endif
184 }
185
186 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
187 RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
188 fillParams->set("Optimize Storage", doOptimizeStorage);
189 Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
190 (transposeR) ? R.getDomainMap() : R.getRangeMap(),
191 fillParams);
192 }
193
194 // transfer striding information
195 RCP<const Map> domainMap = Teuchos::null;
196 RCP<const Map> rangeMap = Teuchos::null;
197
198 const std::string stridedViewLabel("stridedMaps");
199 const size_t blkSize = 1;
200 std::vector<size_t> stridingInfo(1, blkSize);
201 LocalOrdinal stridedBlockId = -1;
202
203 if (R.IsView(stridedViewLabel)) {
204 rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
205 } else {
206 rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
207 rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
208 }
209
210 if (P.IsView(stridedViewLabel)) {
211 domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
212 } else {
213 domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
214 domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
215 }
216 Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
217
218 } // end Multiply
219
220 }; // class TripleMatrixMultiply
221
222#ifdef HAVE_XPETRA_EPETRA
223 // specialization TripleMatrixMultiply for SC=double, LO=GO=int
224 template <>
225 class TripleMatrixMultiply<double,int,int,EpetraNode> {
226 typedef double Scalar;
227 typedef int LocalOrdinal;
228 typedef int GlobalOrdinal;
229 typedef EpetraNode Node;
230#include "Xpetra_UseShortNames.hpp"
231
232 public:
233
234 static void MultiplyRAP(const Matrix& R, bool transposeR,
235 const Matrix& A, bool transposeA,
236 const Matrix& P, bool transposeP,
237 Matrix& Ac,
238 bool call_FillComplete_on_result = true,
239 bool doOptimizeStorage = true,
240 const std::string & label = std::string(),
241 const RCP<ParameterList>& params = null) {
242
243 TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
244 Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
245 TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
246 Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
247
248 TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
249 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
250 TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
251
252 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
253
254 if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
255 throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
256 } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
257#ifdef HAVE_XPETRA_TPETRA
258# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
259 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
260 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
261# else
262 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
263 if(helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P) && helpers::isTpetraCrs(Ac)) {
264 // All matrices are Crs
265 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(R);
266 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
267 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(P);
268 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(Ac);
269
270 // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
271 // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
272 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
273 }
274 else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
275 // All matrices are BlockCrs (except maybe Ac)
276 // FIXME: For the moment we're just going to clobber the innards of AC, so no reuse. Once we have a reuse kernel,
277 // we'll need to think about refactoring BlockCrs so we can do something smarter here.
278 if(!A.getRowMap()->getComm()->getRank())
279 std::cout<<"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
280
281 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(R);
282 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(A);
283 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(P);
284 // Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraBlockCrs(Ac);
285
286 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
287 RCP<const CRS> Rcrs = Tpetra::convertToCrsMatrix(tpR);
288 RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
289 RCP<const CRS> Pcrs = Tpetra::convertToCrsMatrix(tpP);
290 // RCP<CRS> Accrs = Tpetra::convertToCrsMatrix(tpAc);
291
292 // FIXME: The lines below only works because we're assuming Ac is Point
293 RCP<CRS> Accrs = Teuchos::rcp(new CRS(Rcrs->getRowMap(),0));
294 const bool do_fill_complete=true;
295 Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
296
297 // Temporary output matrix
298 RCP<Tpetra::BlockCrsMatrix<SC,LO,GO,NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*Accrs,A.GetStorageBlockSize());
299 RCP<Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO>(Ac_t));
300 RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > Ac_p = Ac_x;
301
302 // We can now cheat and replace the innards of Ac
303 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(Teuchos::rcpFromRef(Ac));
304 Ac_w->replaceCrsMatrix(Ac_p);
305
306 }
307 else {
308 // Mix and match (not supported)
309 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
310 }
311# endif
312#else
313 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
314#endif
315 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
316 RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
317 fillParams->set("Optimize Storage", doOptimizeStorage);
318 Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
319 (transposeR) ? R.getDomainMap() : R.getRangeMap(),
320 fillParams);
321 }
322
323 // transfer striding information
324 RCP<const Map> domainMap = Teuchos::null;
325 RCP<const Map> rangeMap = Teuchos::null;
326
327 const std::string stridedViewLabel("stridedMaps");
328 const size_t blkSize = 1;
329 std::vector<size_t> stridingInfo(1, blkSize);
330 LocalOrdinal stridedBlockId = -1;
331
332 if (R.IsView(stridedViewLabel)) {
333 rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
334 } else {
335 rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
336 rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
337 }
338
339 if (P.IsView(stridedViewLabel)) {
340 domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
341 } else {
342 domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
343 domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
344 }
345 Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
346 }
347
348 } // end Multiply
349
350 }; // end specialization on SC=double, GO=int and NO=EpetraNode
351
352 // specialization TripleMatrixMultiply for SC=double, GO=long long and NO=EpetraNode
353 template <>
354 class TripleMatrixMultiply<double,int,long long,EpetraNode> {
355 typedef double Scalar;
356 typedef int LocalOrdinal;
357 typedef long long GlobalOrdinal;
358 typedef EpetraNode Node;
359#include "Xpetra_UseShortNames.hpp"
360
361 public:
362
363 static void MultiplyRAP(const Matrix& R, bool transposeR,
364 const Matrix& A, bool transposeA,
365 const Matrix& P, bool transposeP,
366 Matrix& Ac,
367 bool call_FillComplete_on_result = true,
368 bool doOptimizeStorage = true,
369 const std::string & label = std::string(),
370 const RCP<ParameterList>& params = null) {
371
372 TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
373 Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
374 TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
375 Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
376
377 TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
378 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
379 TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
380
381 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
382
383 if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
384 throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
385 } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
386#ifdef HAVE_XPETRA_TPETRA
387# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
388 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
389 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long,EpetraNode> ETI enabled."));
390# else
391 using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
392 if(helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P)) {
393 // All matrices are Crs
394 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(R);
395 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
396 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(P);
397 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(Ac);
398
399 // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
400 // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
401 Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
402 }
403 else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
404 // All matrices are BlockCrs (except maybe Ac)
405 // FIXME: For the moment we're just going to clobber the innards of AC, so no reuse. Once we have a reuse kernel,
406 // we'll need to think about refactoring BlockCrs so we can do something smarter here.
407 if(!A.getRowMap()->getComm()->getRank())
408 std::cout<<"WARNING: Using inefficient BlockCrs Multiply Placeholder"<<std::endl;
409
410 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(R);
411 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(A);
412 const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraBlockCrs(P);
413 // Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraBlockCrs(Ac);
414
415 using CRS=Tpetra::CrsMatrix<SC,LO,GO,NO>;
416 RCP<const CRS> Rcrs = Tpetra::convertToCrsMatrix(tpR);
417 RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
418 RCP<const CRS> Pcrs = Tpetra::convertToCrsMatrix(tpP);
419 // RCP<CRS> Accrs = Tpetra::convertToCrsMatrix(tpAc);
420
421 // FIXME: The lines below only works because we're assuming Ac is Point
422 RCP<CRS> Accrs = Teuchos::rcp(new CRS(Rcrs->getRowMap(),0));
423 const bool do_fill_complete=true;
424 Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
425
426 // Temporary output matrix
427 RCP<Tpetra::BlockCrsMatrix<SC,LO,GO,NO> > Ac_t = Tpetra::convertToBlockCrsMatrix(*Accrs,A.GetStorageBlockSize());
428 RCP<Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO>(Ac_t));
429 RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > Ac_p = Ac_x;
430
431 // We can now cheat and replace the innards of Ac
432 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(Teuchos::rcpFromRef(Ac));
433 Ac_w->replaceCrsMatrix(Ac_p);
434 }
435 else {
436 // Mix and match
437 TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
438 }
439
440# endif
441#else
442 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
443#endif
444 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
445 RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
446 fillParams->set("Optimize Storage", doOptimizeStorage);
447 Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
448 (transposeR) ? R.getDomainMap() : R.getRangeMap(),
449 fillParams);
450 }
451
452 // transfer striding information
453 RCP<const Map> domainMap = Teuchos::null;
454 RCP<const Map> rangeMap = Teuchos::null;
455
456 const std::string stridedViewLabel("stridedMaps");
457 const size_t blkSize = 1;
458 std::vector<size_t> stridingInfo(1, blkSize);
459 LocalOrdinal stridedBlockId = -1;
460
461 if (R.IsView(stridedViewLabel)) {
462 rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
463 } else {
464 rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
465 rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
466 }
467
468 if (P.IsView(stridedViewLabel)) {
469 domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
470 } else {
471 domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
472 domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
473 }
474 Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
475 }
476
477 } // end Multiply
478
479 }; // end specialization on GO=long long and NO=EpetraNode
480#endif
481
482} // end namespace Xpetra
483
484#define XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
485
486#endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_HPP_ */
Xpetra_ConfigDefs.hpp
Xpetra_IO.hpp
Xpetra_MatrixFactory.hpp
Xpetra_Matrix.hpp
Xpetra_UseShortNames.hpp
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition Xpetra_Exceptions.hpp:102
Xpetra::Matrix
Xpetra-specific matrix class.
Definition Xpetra_Matrix.hpp:97
Xpetra::Matrix::CreateView
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
Definition Xpetra_Matrix.hpp:126
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode.
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:328
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:205
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:318
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::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....
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::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:234
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:363
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:108
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra
Xpetra namespace
Definition Xpetra_BlockedCrsMatrix.hpp:89
Xpetra::EpetraNode
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
Definition Xpetra_Map_decl.hpp:71
Xpetra::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)