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 
62 #ifdef HAVE_XPETRA_TPETRA
63 #include <TpetraExt_TripleMatrixMultiply.hpp>
64 #include <Xpetra_TpetraCrsMatrix.hpp>
65 // #include <Xpetra_TpetraMultiVector.hpp>
66 // #include <Xpetra_TpetraVector.hpp>
67 #endif // HAVE_XPETRA_TPETRA
68 
69 namespace Xpetra {
70 
71  template <class Scalar,
72  class LocalOrdinal /*= int*/,
73  class GlobalOrdinal /*= LocalOrdinal*/,
74  class Node /*= KokkosClassic::DefaultNode::DefaultNodeType*/>
76 #undef XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
77 #include "Xpetra_UseShortNames.hpp"
78 
79  public:
80 
105  static void MultiplyRAP(const Matrix& R, bool transposeR,
106  const Matrix& A, bool transposeA,
107  const Matrix& P, bool transposeP,
108  Matrix& Ac,
109  bool call_FillComplete_on_result = true,
110  bool doOptimizeStorage = true,
111  const std::string & label = std::string(),
112  const RCP<ParameterList>& params = null) {
113 
114  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
115  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
116  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
117  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
118 
119  TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
120  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
121  TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
122 
123  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
124 
125  if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
126  throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
127  } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
128 #ifdef HAVE_XPETRA_TPETRA
129  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(R);
130  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
131  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(P);
132  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(Ac);
133 
134  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
135  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
136  Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
137 #else
138  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
139 #endif
140  }
141 
142  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
143  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
144  fillParams->set("Optimize Storage", doOptimizeStorage);
145  Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
146  (transposeR) ? R.getDomainMap() : R.getRangeMap(),
147  fillParams);
148  }
149 
150  // transfer striding information
151  RCP<const Map> domainMap = Teuchos::null;
152  RCP<const Map> rangeMap = Teuchos::null;
153 
154  const std::string stridedViewLabel("stridedMaps");
155  const size_t blkSize = 1;
156  std::vector<size_t> stridingInfo(1, blkSize);
157  LocalOrdinal stridedBlockId = -1;
158 
159  if (R.IsView(stridedViewLabel)) {
160  rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
161  } else {
162  rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
163  rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
164  }
165 
166  if (P.IsView(stridedViewLabel)) {
167  domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
168  } else {
169  domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
170  domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
171  }
172  Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
173 
174  } // end Multiply
175 
176  }; // class TripleMatrixMultiply
177 
178 #ifdef HAVE_XPETRA_EPETRA
179  // specialization TripleMatrixMultiply for SC=double, LO=GO=int
180  template <>
181  class TripleMatrixMultiply<double,int,int,EpetraNode> {
182  typedef double Scalar;
183  typedef int LocalOrdinal;
184  typedef int GlobalOrdinal;
185  typedef EpetraNode Node;
186 #include "Xpetra_UseShortNames.hpp"
187 
188  public:
189 
190  static void MultiplyRAP(const Matrix& R, bool transposeR,
191  const Matrix& A, bool transposeA,
192  const Matrix& P, bool transposeP,
193  Matrix& Ac,
194  bool call_FillComplete_on_result = true,
195  bool doOptimizeStorage = true,
196  const std::string & label = std::string(),
197  const RCP<ParameterList>& params = null) {
198 
199  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
200  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
201  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
202  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
203 
204  TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
205  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
206  TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
207 
208  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
209 
210  if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
211  throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
212  } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
213 #ifdef HAVE_XPETRA_TPETRA
214 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
215  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
216  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
217 # else
218  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(R);
219  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
220  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(P);
221  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(Ac);
222 
223  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
224  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
225  Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
226 # endif
227 #else
228  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
229 #endif
230  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
231  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
232  fillParams->set("Optimize Storage", doOptimizeStorage);
233  Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
234  (transposeR) ? R.getDomainMap() : R.getRangeMap(),
235  fillParams);
236  }
237 
238  // transfer striding information
239  RCP<const Map> domainMap = Teuchos::null;
240  RCP<const Map> rangeMap = Teuchos::null;
241 
242  const std::string stridedViewLabel("stridedMaps");
243  const size_t blkSize = 1;
244  std::vector<size_t> stridingInfo(1, blkSize);
245  LocalOrdinal stridedBlockId = -1;
246 
247  if (R.IsView(stridedViewLabel)) {
248  rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
249  } else {
250  rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
251  rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
252  }
253 
254  if (P.IsView(stridedViewLabel)) {
255  domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
256  } else {
257  domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
258  domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
259  }
260  Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
261  }
262 
263  } // end Multiply
264 
265  }; // end specialization on SC=double, GO=int and NO=EpetraNode
266 
267  // specialization TripleMatrixMultiply for SC=double, GO=long long and NO=EpetraNode
268  template <>
269  class TripleMatrixMultiply<double,int,long long,EpetraNode> {
270  typedef double Scalar;
271  typedef int LocalOrdinal;
272  typedef long long GlobalOrdinal;
273  typedef EpetraNode Node;
274 #include "Xpetra_UseShortNames.hpp"
275 
276  public:
277 
278  static void MultiplyRAP(const Matrix& R, bool transposeR,
279  const Matrix& A, bool transposeA,
280  const Matrix& P, bool transposeP,
281  Matrix& Ac,
282  bool call_FillComplete_on_result = true,
283  bool doOptimizeStorage = true,
284  const std::string & label = std::string(),
285  const RCP<ParameterList>& params = null) {
286 
287  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
288  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
289  TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
290  Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
291 
292  TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
293  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
294  TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
295 
296  bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
297 
298  if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
299  throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
300  } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
301 #ifdef HAVE_XPETRA_TPETRA
302 # if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
303  (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
304  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long,EpetraNode> ETI enabled."));
305 # else
306  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpR = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(R);
307  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
308  const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpP = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(P);
309  Tpetra::CrsMatrix<SC,LO,GO,NO> & tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(Ac);
310 
311  // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
312  // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
313  Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
314 # endif
315 #else
316  throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
317 #endif
318  if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
319  RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
320  fillParams->set("Optimize Storage", doOptimizeStorage);
321  Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
322  (transposeR) ? R.getDomainMap() : R.getRangeMap(),
323  fillParams);
324  }
325 
326  // transfer striding information
327  RCP<const Map> domainMap = Teuchos::null;
328  RCP<const Map> rangeMap = Teuchos::null;
329 
330  const std::string stridedViewLabel("stridedMaps");
331  const size_t blkSize = 1;
332  std::vector<size_t> stridingInfo(1, blkSize);
333  LocalOrdinal stridedBlockId = -1;
334 
335  if (R.IsView(stridedViewLabel)) {
336  rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
337  } else {
338  rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
339  rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
340  }
341 
342  if (P.IsView(stridedViewLabel)) {
343  domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
344  } else {
345  domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
346  domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
347  }
348  Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
349  }
350 
351  } // end Multiply
352 
353  }; // end specialization on GO=long long and NO=EpetraNode
354 #endif
355 
356 } // end namespace Xpetra
357 
358 #define XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
359 
360 #endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_HPP_ */
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)
Exception throws to report errors in the internal logical of the program.
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 bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
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.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra-specific matrix class.
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)