MueLu Version of the Day
Loading...
Searching...
No Matches
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
33namespace Belos {
34using Teuchos::RCP;
35using Teuchos::rcpFromRef;
36
37//
39
40
44class MueLuOpFailure : public BelosError {
45 public:
46 MueLuOpFailure(const std::string& what_arg)
47 : BelosError(what_arg) {}
48};
49
61template <class Scalar,
62 class LocalOrdinal = int,
63 class GlobalOrdinal = LocalOrdinal,
64 class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
65class 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
78#ifdef HAVE_MUELU_AMGX
81#endif
83 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
163template <>
164class 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
301template <>
302class 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
MueLu_AMGXOperator.hpp
MueLu_ConfigDefs.hpp
EpetraNode
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
Definition MueLu_ETI_3arg.hpp:29
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition MueLu_UseDefaultTypes.hpp:13
Scalar
MueLu::DefaultScalar Scalar
Definition MueLu_UseDefaultTypes.hpp:12
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition MueLu_UseDefaultTypes.hpp:14
Node
MueLu::DefaultNode Node
Definition MueLu_UseDefaultTypes.hpp:15
Belos::MueLuOpFailure
MueLuOpFailure is thrown when a return value from an MueLu call on an Xpetra::Operator or MueLu::Hier...
Definition BelosMueLuAdapter.hpp:44
Belos::MueLuOpFailure::MueLuOpFailure
MueLuOpFailure(const std::string &what_arg)
Definition BelosMueLuAdapter.hpp:46
Belos::MueLuOp
MueLuOp derives from Belos::OperatorT and administrates a MueLu::Hierarchy. It implements the apply c...
Definition BelosMueLuAdapter.hpp:70
Belos::MueLuOp::AMGX_
RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > AMGX_
Definition BelosMueLuAdapter.hpp:146
Belos::MueLuOp::Hierarchy_
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Hierarchy_
Definition BelosMueLuAdapter.hpp:144
Belos::MueLuOp::~MueLuOp
virtual ~MueLuOp()
Destructor.
Definition BelosMueLuAdapter.hpp:83
Belos::MueLuOp::MueLuOp
MueLuOp(const RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &H)
Default constructor.
Definition BelosMueLuAdapter.hpp:76
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::MueLuOp
MueLuOp(const RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
Definition BelosMueLuAdapter.hpp:79
Epetra_MultiVector
MueLu::AMGXOperator
Adapter for AmgX library from Nvidia.
Definition MueLu_AMGXOperator_decl.hpp:44
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition MueLu_Hierarchy_decl.hpp:63