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"
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 /*= KokkosClassic::DefaultNode::DefaultNodeType*/>
79#undef XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
81
82 public:
83
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
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());
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) {
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;
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
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());
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) {
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;
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
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());
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) {
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_ */
Exception throws to report errors in the internal logical of the program.
static RCP< const Tpetra::BlockCrsMatrix< SC, LO, GO, NO > > Op2TpetraBlockCrs(RCP< Matrix > Op)
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Xpetra-specific matrix class.
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode.
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
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.
bool IsView(const viewLabel_t viewLabel) const
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
virtual LocalOrdinal GetStorageBlockSize() const =0
Returns the block size of the storage mechanism, which is usually 1, except for Tpetra::BlockCrsMatri...
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y....
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X....
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.
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)
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)
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace