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 // ***********************************************************************
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
61 #include "MueLu_AMGXOperator.hpp"
62 #endif
63 
64 #include <BelosOperatorT.hpp>
65 
66 #include "MueLu_ConfigDefs.hpp"
67 #include "MueLu_Hierarchy.hpp"
68 
69 namespace 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
116  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 
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 
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:
182 #ifdef HAVE_MUELU_AMGX
184 #endif
185  };
186 
187 #ifdef HAVE_XPETRA_EPETRA
188 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
189 
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
327 
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
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLuOpFailure is thrown when a return value from an MueLu call on an Xpetra::Operator or MueLu::Hier...
MueLuOp derives from Belos::OperatorT and administrates a MueLu::Hierarchy. It implements the apply c...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::DefaultNode Node
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::...
RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > AMGX_
MueLuOpFailure(const std::string &what_arg)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Hierarchy_
MueLuOp(const RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &H)
Default constructor.
virtual ~MueLuOp()
Destructor.
MueLuOp(const RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Adapter for AmgX library from Nvidia.
bool is_null() const