MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosMueLuAdapter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_MUELU_ADAPTER_HPP
11 #define BELOS_MUELU_ADAPTER_HPP
12 
13 // TAW: 3/4/2016: we use the Xpetra macros
14 // These are available and Xpetra is prerequisite for MueLu
15 #ifdef HAVE_XPETRA_EPETRA
16 #include <Epetra_config.h>
17 #include <BelosOperator.hpp>
18 #endif
19 
20 //#ifdef HAVE_XPETRA_TPETRA
21 //#include "TpetraCore_config.h"
22 //#endif
23 
24 #ifdef HAVE_MUELU_AMGX
25 #include "MueLu_AMGXOperator.hpp"
26 #endif
27 
28 #include <BelosOperatorT.hpp>
29 
30 #include "MueLu_ConfigDefs.hpp"
31 #include "MueLu_Hierarchy.hpp"
32 
33 namespace Belos {
34 using Teuchos::RCP;
35 using Teuchos::rcpFromRef;
36 
37 //
39 
40 
44 class MueLuOpFailure : public BelosError {
45  public:
46  MueLuOpFailure(const std::string& what_arg)
47  : BelosError(what_arg) {}
48 };
49 
61 template <class Scalar,
62  class LocalOrdinal = int,
63  class GlobalOrdinal = LocalOrdinal,
64  class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
65 class MueLuOp : public OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
66 #ifdef HAVE_XPETRA_TPETRA
67  ,
68  public OperatorT<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
69 #endif
70 {
71  public:
73 
74 
76  MueLuOp(const RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& H)
77  : Hierarchy_(H) {}
78 #ifdef HAVE_MUELU_AMGX
79  MueLuOp(const RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& A)
80  : AMGX_(A) {}
81 #endif
82  virtual ~MueLuOp() {}
85 
87 
88 
94  void Apply(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS) const {
95  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
96  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
97 
98  // This does not matter for Hierarchy, but matters for AMGX
99  y.putScalar(0.0);
100 
101 #ifdef HAVE_MUELU_AMGX
102  if (!AMGX_.is_null()) {
103  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Xpetra::toTpetra(x);
104  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY = Xpetra::toTpetra(y);
105 
106  AMGX_->apply(tX, tY);
107  }
108 #endif
109  if (!Hierarchy_.is_null())
110  Hierarchy_->Iterate(x, y, 1, true);
111  }
113 
114 #ifdef HAVE_XPETRA_TPETRA
115  // TO SKIP THE TRAIT IMPLEMENTATION OF XPETRA::MULTIVECTOR
121  void Apply(const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS) const {
122  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
123  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
124 
125  // FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
126  y.putScalar(0.0);
127 
128 #ifdef HAVE_MUELU_AMGX
129  if (!AMGX_.is_null())
130  AMGX_->apply(x, y);
131 #endif
132 
133  if (!Hierarchy_.is_null()) {
134  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(x);
135 
136  const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
137  Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
138  Hierarchy_->Iterate(tX, tY, 1, true);
139  }
140  }
141 #endif
142 
143  private:
144  RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Hierarchy_;
145 #ifdef HAVE_MUELU_AMGX
146  RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AMGX_;
147 #endif
148 };
149 
150 #ifdef HAVE_XPETRA_EPETRA
151 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
152 
163 template <>
164 class MueLuOp<double, int, int, Xpetra::EpetraNode> : public OperatorT<Xpetra::MultiVector<double, int, int, Xpetra::EpetraNode> >
165 #ifdef HAVE_XPETRA_TPETRA
166 // check whether Tpetra is instantiated on double,int,int,EpetraNode
167 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
168  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
169  ,
170  public OperatorT<Tpetra::MultiVector<double, int, int, Xpetra::EpetraNode> >
171 #endif
172 #endif
173 #ifdef HAVE_XPETRA_EPETRA
174  ,
175  public OperatorT<Epetra_MultiVector>,
176  public Belos::Operator<double>
177 #endif
178 {
179  typedef double Scalar;
180  typedef int LocalOrdinal;
181  typedef int GlobalOrdinal;
182  typedef Xpetra::EpetraNode Node;
183 
184  public:
185  MueLuOp(const RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& H)
186  : Hierarchy_(H) {}
187 #ifdef HAVE_MUELU_AMGX
188  MueLuOp(const RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& A)
189  : AMGX_(A) {}
190 #endif
191  virtual ~MueLuOp() {}
192 
193  void Apply(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS) const {
194  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
195  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
196 
197  // FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
198  y.putScalar(0.0);
199 
200 #ifdef HAVE_MUELU_AMGX
201  if (!AMGX_.is_null()) {
202  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Xpetra::toTpetra(x);
203  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY = Xpetra::toTpetra(y);
204 
205  AMGX_->apply(tX, tY);
206  }
207 #endif
208  if (!Hierarchy_.is_null())
209  Hierarchy_->Iterate(x, y, 1, true);
210  }
211 
212 #ifdef HAVE_XPETRA_TPETRA
213 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
214  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
215  void Apply(const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS) const {
216  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
217  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
218 
219  // FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
220  y.putScalar(0.0);
221 
222 #ifdef HAVE_MUELU_AMGX
223  if (!AMGX_.is_null())
224  AMGX_->apply(x, y);
225 #endif
226 
227  if (!Hierarchy_.is_null()) {
228  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(x);
229 
230  const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
231  Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
232 
233  tY.putScalar(0.0);
234 
235  Hierarchy_->Iterate(tX, tY, 1, true);
236  }
237  }
238 #endif
239 #endif
240 
241 #ifdef HAVE_XPETRA_EPETRA
242  // TO SKIP THE TRAIT IMPLEMENTATION OF XPETRA::MULTIVECTOR
248  void Apply(const Epetra_MultiVector& x, Epetra_MultiVector& y, ETrans trans = NOTRANS) const {
249  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
250  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
251 
252  Epetra_MultiVector& temp_x = const_cast<Epetra_MultiVector&>(x);
253 
254  const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
255  Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node> tY(rcpFromRef(y));
256 
257  // FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate().
258  tY.putScalar(0.0);
259 
260  Hierarchy_->Iterate(tX, tY, 1, true);
261  }
262 
268  void Apply(const Belos::MultiVec<double>& x, Belos::MultiVec<double>& y, ETrans trans = NOTRANS) const {
269  const Epetra_MultiVector* vec_x = dynamic_cast<const Epetra_MultiVector*>(&x);
270  Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector*>(&y);
271 
272  TEUCHOS_TEST_FOR_EXCEPTION(vec_x == NULL || vec_y == NULL, MueLuOpFailure,
273  "Belos::MueLuOp::Apply, x and/or y cannot be dynamic cast to an Epetra_MultiVector.");
274 
275  Apply(*vec_x, *vec_y, trans);
276  }
277 #endif
278 
279  private:
280  RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Hierarchy_;
281 #ifdef HAVE_MUELU_AMGX
282  RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AMGX_;
283 #endif
284 };
285 #endif // !EPETRA_NO_32BIT_GLOBAL_INDICES
286 #endif // HAVE_XPETRA_EPETRA
287 
288 #ifdef HAVE_XPETRA_EPETRA
289 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
290 
301 template <>
302 class MueLuOp<double, int, long long, Xpetra::EpetraNode> : public OperatorT<Xpetra::MultiVector<double, int, long long, Xpetra::EpetraNode> >
303 #ifdef HAVE_XPETRA_TPETRA
304 // check whether Tpetra is instantiated on double,int,int,EpetraNode
305 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
306  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
307  ,
308  public OperatorT<Tpetra::MultiVector<double, int, long long, Xpetra::EpetraNode> >
309 #endif
310 #endif
311 #ifdef HAVE_XPETRA_EPETRA
312  ,
313  public OperatorT<Epetra_MultiVector>,
314  public Belos::Operator<double>
315 #endif
316 {
317  typedef double Scalar;
318  typedef int LocalOrdinal;
319  typedef long long GlobalOrdinal;
320  typedef Xpetra::EpetraNode Node;
321 
322  public:
323  MueLuOp(const RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& H)
324  : Hierarchy_(H) {}
325 #ifdef HAVE_MUELU_AMGX
326  MueLuOp(const RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& A)
327  : AMGX_(A) {}
328 #endif
329  virtual ~MueLuOp() {}
330 
331  void Apply(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS) const {
332  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
333  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
334 
335  // FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
336  y.putScalar(0.0);
337 
338 #ifdef HAVE_MUELU_AMGX
339  if (!AMGX_.is_null()) {
340  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Xpetra::toTpetra(x);
341  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY = Xpetra::toTpetra(y);
342 
343  AMGX_->apply(tX, tY);
344  }
345 #endif
346  if (!Hierarchy_.is_null())
347  Hierarchy_->Iterate(x, y, 1, true);
348  }
349 
350 #ifdef HAVE_XPETRA_TPETRA
351 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
352  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
353  void Apply(const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS) const {
354  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
355  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
356 
357  // FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
358  y.putScalar(0.0);
359 
360 #ifdef HAVE_MUELU_AMGX
361  if (!AMGX_.is_null())
362  AMGX_->apply(x, y);
363 #endif
364 
365  if (!Hierarchy_.is_null()) {
366  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(x);
367 
368  const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
369  Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
370 
371  tY.putScalar(0.0);
372 
373  Hierarchy_->Iterate(tX, tY, 1, true);
374  }
375  }
376 #endif
377 #endif
378 
379 #ifdef HAVE_XPETRA_EPETRA
380  // TO SKIP THE TRAIT IMPLEMENTATION OF XPETRA::MULTIVECTOR
386  void Apply(const Epetra_MultiVector& x, Epetra_MultiVector& y, ETrans trans = NOTRANS) const {
387  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
388  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
389 
390  Epetra_MultiVector& temp_x = const_cast<Epetra_MultiVector&>(x);
391 
392  const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
393  Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node> tY(rcpFromRef(y));
394 
395  // FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate().
396  tY.putScalar(0.0);
397 
398  Hierarchy_->Iterate(tX, tY, 1, true);
399  }
400 
406  void Apply(const Belos::MultiVec<double>& x, Belos::MultiVec<double>& y, ETrans trans = NOTRANS) const {
407  const Epetra_MultiVector* vec_x = dynamic_cast<const Epetra_MultiVector*>(&x);
408  Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector*>(&y);
409 
410  TEUCHOS_TEST_FOR_EXCEPTION(vec_x == NULL || vec_y == NULL, MueLuOpFailure,
411  "Belos::MueLuOp::Apply, x and/or y cannot be dynamic cast to an Epetra_MultiVector.");
412 
413  Apply(*vec_x, *vec_y, trans);
414  }
415 #endif
416 
417  private:
418  RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Hierarchy_;
419 #ifdef HAVE_MUELU_AMGX
420  RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AMGX_;
421 #endif
422 };
423 #endif // !EPETRA_NO_64BIT_GLOBAL_INDICES
424 #endif // HAVE_XPETRA_EPETRA
425 } // namespace Belos
426 
427 #endif // BELOS_MUELU_ADAPTER_HPP
Epetra_MultiVector
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:13
Belos::MueLuOpFailure
MueLuOpFailure is thrown when a return value from an MueLu call on an Xpetra::Operator or MueLu::Hier...
Definition: BelosMueLuAdapter.hpp:44
Tpetra::MultiVector::putScalar
void putScalar(const Scalar &value)
Belos::MueLuOp
MueLuOp derives from Belos::OperatorT and administrates a MueLu::Hierarchy. It implements the apply c...
Definition: BelosMueLuAdapter.hpp:65
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Xpetra::MultiVector::putScalar
virtual void putScalar(const Scalar &value)=0
MueLu_ConfigDefs.hpp
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:15
Belos::MueLuOp::Apply
void Apply(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &y, ETrans trans=NOTRANS) const
This routine takes the Xpetra::MultiVector x and applies the operator to it resulting in the Xpetra::...
Definition: BelosMueLuAdapter.hpp:94
Belos::MueLuOp::AMGX_
RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > AMGX_
Definition: BelosMueLuAdapter.hpp:146
Belos::MueLuOpFailure::MueLuOpFailure
MueLuOpFailure(const std::string &what_arg)
Definition: BelosMueLuAdapter.hpp:46
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:12
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:14
Belos::MueLuOp::Hierarchy_
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Hierarchy_
Definition: BelosMueLuAdapter.hpp:144
MueLu_AMGXOperator.hpp
Belos::MueLuOp::MueLuOp
MueLuOp(const RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &H)
Default constructor.
Definition: BelosMueLuAdapter.hpp:76
Belos::MueLuOp::~MueLuOp
virtual ~MueLuOp()
Destructor.
Definition: BelosMueLuAdapter.hpp:83
Xpetra::toTpetra
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
Belos::MueLuOp::MueLuOp
MueLuOp(const RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
Definition: BelosMueLuAdapter.hpp:79
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition: MueLu_Hierarchy_decl.hpp:63
MueLu::AMGXOperator
Adapter for AmgX library from Nvidia.
Definition: MueLu_AMGXOperator_fwd.hpp:18
Teuchos::RCP::is_null
bool is_null() const