MueLu Version of the Day
Loading...
Searching...
No Matches
BelosMueLuAdapter.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
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//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef BELOS_MUELU_ADAPTER_HPP
47#define BELOS_MUELU_ADAPTER_HPP
48
49// TAW: 3/4/2016: we use the Xpetra macros
50// These are available and Xpetra is prerequisite for MueLu
51#ifdef HAVE_XPETRA_EPETRA
52#include <Epetra_config.h>
53#include <BelosOperator.hpp>
54#endif
55
56//#ifdef HAVE_XPETRA_TPETRA
57//#include "TpetraCore_config.h"
58//#endif
59
60#ifdef HAVE_MUELU_AMGX
62#endif
63
64#include <BelosOperatorT.hpp>
65
66#include "MueLu_ConfigDefs.hpp"
67#include "MueLu_Hierarchy.hpp"
68
69namespace Belos {
70 using Teuchos::RCP;
71 using Teuchos::rcpFromRef;
72
73 //
75
76
80 class MueLuOpFailure : public BelosError {
81 public:
82 MueLuOpFailure(const std::string& what_arg) : BelosError(what_arg) {}
83 };
84
96 template <class Scalar,
97 class LocalOrdinal = int,
99 class Node = KokkosClassic::DefaultNode::DefaultNodeType>
100 class MueLuOp :
101 public OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
102#ifdef HAVE_XPETRA_TPETRA
103 , public OperatorT<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
104#endif
105 {
106 public:
107
109
110
113#ifdef HAVE_MUELU_AMGX
115#endif
117 virtual ~MueLuOp() {}
119
121
122
128 void Apply(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {
129
130 TEUCHOS_TEST_FOR_EXCEPTION(trans!=NOTRANS, MueLuOpFailure,
131 "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
132
133 // This does not matter for Hierarchy, but matters for AMGX
134 y.putScalar(0.0);
135
136#ifdef HAVE_MUELU_AMGX
137 if (!AMGX_.is_null()) {
138 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Xpetra::toTpetra(x);
139 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY = Xpetra::toTpetra(y);
140
141 AMGX_->apply(tX, tY);
142
143 }
144#endif
145 if (!Hierarchy_.is_null())
146 Hierarchy_->Iterate(x, y, 1, true);
147 }
149
150#ifdef HAVE_XPETRA_TPETRA
151 // TO SKIP THE TRAIT IMPLEMENTATION OF XPETRA::MULTIVECTOR
157 void Apply(const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {
158
159 TEUCHOS_TEST_FOR_EXCEPTION(trans!=NOTRANS, MueLuOpFailure,
160 "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
161
162 //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
163 y.putScalar(0.0);
164
165#ifdef HAVE_MUELU_AMGX
166 if (!AMGX_.is_null())
167 AMGX_->apply(x, y);
168#endif
169
170 if (!Hierarchy_.is_null()) {
171 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(x);
172
173 const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
174 Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
175 Hierarchy_->Iterate(tX, tY, 1, true);
176 }
177 }
178#endif
179
180 private:
181 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Hierarchy_;
182#ifdef HAVE_MUELU_AMGX
183 RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AMGX_;
184#endif
185 };
186
187#ifdef HAVE_XPETRA_EPETRA
188#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
200 template <>
201 class MueLuOp<double, int, int, Xpetra::EpetraNode> :
202 public OperatorT<Xpetra::MultiVector<double, int, int, Xpetra::EpetraNode> >
203#ifdef HAVE_XPETRA_TPETRA
204 // check whether Tpetra is instantiated on double,int,int,EpetraNode
205#if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
206 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
207 , public OperatorT<Tpetra::MultiVector<double, int, int, Xpetra::EpetraNode> >
208#endif
209#endif
210#ifdef HAVE_XPETRA_EPETRA
211 , public OperatorT<Epetra_MultiVector>
212 , public Belos::Operator<double>
213#endif
214 {
215 typedef double Scalar;
216 typedef int LocalOrdinal;
217 typedef int GlobalOrdinal;
218 typedef Xpetra::EpetraNode Node;
219
220 public:
221
223#ifdef HAVE_MUELU_AMGX
225#endif
226 virtual ~MueLuOp() {}
227
228 void Apply(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {
229
230 TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
231 "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
232
233 //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
234 y.putScalar(0.0);
235
236#ifdef HAVE_MUELU_AMGX
237 if (!AMGX_.is_null()) {
238 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Xpetra::toTpetra(x);
239 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY = Xpetra::toTpetra(y);
240
241 AMGX_->apply(tX, tY);
242 }
243#endif
244 if (!Hierarchy_.is_null())
245 Hierarchy_->Iterate(x, y, 1, true);
246 }
247
248#ifdef HAVE_XPETRA_TPETRA
249#if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
250 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
251 void Apply ( const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans=NOTRANS ) const {
252 TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
253 "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
254
255 //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
256 y.putScalar(0.0);
257
258#ifdef HAVE_MUELU_AMGX
259 if (!AMGX_.is_null())
260 AMGX_->apply(x, y);
261#endif
262
263 if (!Hierarchy_.is_null()) {
264 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(x);
265
266 const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
267 Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
268
269 tY.putScalar(0.0);
270
271 Hierarchy_->Iterate(tX, tY, 1, true);
272 }
273 }
274#endif
275#endif
276
277#ifdef HAVE_XPETRA_EPETRA
278 // TO SKIP THE TRAIT IMPLEMENTATION OF XPETRA::MULTIVECTOR
284 void Apply(const Epetra_MultiVector& x, Epetra_MultiVector& y, ETrans trans = NOTRANS) const {
285 TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
286 "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
287
288 Epetra_MultiVector& temp_x = const_cast<Epetra_MultiVector&>(x);
289
290 const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> tX(rcpFromRef(temp_x));
291 Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> tY(rcpFromRef(y));
292
293 //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate().
294 tY.putScalar(0.0);
295
296 Hierarchy_->Iterate(tX, tY, 1, true);
297 }
298
304 void Apply(const Belos::MultiVec<double>& x, Belos::MultiVec<double>& y, ETrans trans = NOTRANS ) const {
305 const Epetra_MultiVector* vec_x = dynamic_cast<const Epetra_MultiVector*>(&x);
306 Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector*>(&y);
307
308 TEUCHOS_TEST_FOR_EXCEPTION(vec_x==NULL || vec_y==NULL, MueLuOpFailure,
309 "Belos::MueLuOp::Apply, x and/or y cannot be dynamic cast to an Epetra_MultiVector.");
310
311 Apply(*vec_x, *vec_y, trans);
312 }
313#endif
314
315 private:
316 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Hierarchy_;
317#ifdef HAVE_MUELU_AMGX
318 RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AMGX_;
319#endif
320 };
321#endif // !EPETRA_NO_32BIT_GLOBAL_INDICES
322#endif // HAVE_XPETRA_EPETRA
323
324
325#ifdef HAVE_XPETRA_EPETRA
326#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
338 template <>
339 class MueLuOp<double, int, long long, Xpetra::EpetraNode> :
340 public OperatorT<Xpetra::MultiVector<double, int, long long, Xpetra::EpetraNode> >
341#ifdef HAVE_XPETRA_TPETRA
342 // check whether Tpetra is instantiated on double,int,int,EpetraNode
343#if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
344 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
345 , public OperatorT<Tpetra::MultiVector<double, int, long long, Xpetra::EpetraNode> >
346#endif
347#endif
348#ifdef HAVE_XPETRA_EPETRA
349 , public OperatorT<Epetra_MultiVector>
350 , public Belos::Operator<double>
351#endif
352 {
353 typedef double Scalar;
354 typedef int LocalOrdinal;
355 typedef long long GlobalOrdinal;
356 typedef Xpetra::EpetraNode Node;
357
358 public:
359
361#ifdef HAVE_MUELU_AMGX
363#endif
364 virtual ~MueLuOp() {}
365
366 void Apply(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {
367
368 TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
369 "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
370
371 //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
372 y.putScalar(0.0);
373
374#ifdef HAVE_MUELU_AMGX
375 if (!AMGX_.is_null()) {
376 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Xpetra::toTpetra(x);
377 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY = Xpetra::toTpetra(y);
378
379 AMGX_->apply(tX, tY);
380 }
381#endif
382 if (!Hierarchy_.is_null())
383 Hierarchy_->Iterate(x, y, 1, true);
384 }
385
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 void Apply ( const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans=NOTRANS ) const {
390 TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
391 "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
392
393 //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
394 y.putScalar(0.0);
395
396#ifdef HAVE_MUELU_AMGX
397 if (!AMGX_.is_null())
398 AMGX_->apply(x, y);
399#endif
400
401 if (!Hierarchy_.is_null()) {
402 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(x);
403
404 const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
405 Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
406
407 tY.putScalar(0.0);
408
409 Hierarchy_->Iterate(tX, tY, 1, true);
410 }
411 }
412#endif
413#endif
414
415#ifdef HAVE_XPETRA_EPETRA
416 // TO SKIP THE TRAIT IMPLEMENTATION OF XPETRA::MULTIVECTOR
422 void Apply(const Epetra_MultiVector& x, Epetra_MultiVector& y, ETrans trans = NOTRANS) const {
423 TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
424 "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
425
426 Epetra_MultiVector& temp_x = const_cast<Epetra_MultiVector&>(x);
427
428 const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> tX(rcpFromRef(temp_x));
429 Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> tY(rcpFromRef(y));
430
431 //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate().
432 tY.putScalar(0.0);
433
434 Hierarchy_->Iterate(tX, tY, 1, true);
435 }
436
442 void Apply(const Belos::MultiVec<double>& x, Belos::MultiVec<double>& y, ETrans trans = NOTRANS ) const {
443 const Epetra_MultiVector* vec_x = dynamic_cast<const Epetra_MultiVector*>(&x);
444 Epetra_MultiVector* vec_y = dynamic_cast<Epetra_MultiVector*>(&y);
445
446 TEUCHOS_TEST_FOR_EXCEPTION(vec_x==NULL || vec_y==NULL, MueLuOpFailure,
447 "Belos::MueLuOp::Apply, x and/or y cannot be dynamic cast to an Epetra_MultiVector.");
448
449 Apply(*vec_x, *vec_y, trans);
450 }
451#endif
452
453 private:
454 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Hierarchy_;
455#ifdef HAVE_MUELU_AMGX
456 RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AMGX_;
457#endif
458 };
459#endif // !EPETRA_NO_64BIT_GLOBAL_INDICES
460#endif // HAVE_XPETRA_EPETRA
461} // namespace Belos
462
463#endif // BELOS_MUELU_ADAPTER_HPP
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
MueLuOpFailure is thrown when a return value from an MueLu call on an Xpetra::Operator or MueLu::Hier...
MueLuOpFailure(const std::string &what_arg)
MueLuOp derives from Belos::OperatorT and administrates a MueLu::Hierarchy. It implements the apply c...
RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > AMGX_
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Hierarchy_
virtual ~MueLuOp()
Destructor.
MueLuOp(const RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &H)
Default constructor.
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::...
MueLuOp(const RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
Adapter for AmgX library from Nvidia.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.