MueLu  Version of the Day
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 #ifdef HAVE_MUELU_EPETRA
50 #include <BelosOperator.hpp>
51 #endif
52 
53 #ifdef HAVE_MUELU_AMGX
54 #include "MueLu_AMGXOperator.hpp"
55 #endif
56 
57 #include <BelosOperatorT.hpp>
58 
59 #include "MueLu_ConfigDefs.hpp"
60 #include "MueLu_Hierarchy.hpp"
61 
62 namespace Belos {
63  using Teuchos::RCP;
64  using Teuchos::rcpFromRef;
65 
66  //
68 
69 
73  class MueLuOpFailure : public BelosError {
74  public:
75  MueLuOpFailure(const std::string& what_arg) : BelosError(what_arg) {}
76  };
77 
78  //TODO: doc
79  template <class Scalar,
80  class LocalOrdinal = int,
81  class GlobalOrdinal = LocalOrdinal,
82  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
83  class MueLuOp :
84  public OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
85 #ifdef HAVE_MUELU_TPETRA
86  , public OperatorT<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
87 #endif
88  {
89  public:
90 
92 
93 
96 #ifdef HAVE_MUELU_AMGX
98 #endif
99  virtual ~MueLuOp() {}
102 
104 
105 
111  void Apply(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {
112 
113  TEUCHOS_TEST_FOR_EXCEPTION(trans!=NOTRANS, MueLuOpFailure,
114  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
115 
116  // This does not matter for Hierarchy, but matters for AMGX
117  y.putScalar(0.0);
118 
119 #ifdef HAVE_MUELU_AMGX
120  if (!AMGX_.is_null()) {
121  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Xpetra::toTpetra(x);
122  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY = Xpetra::toTpetra(y);
123 
124  AMGX_->apply(tX, tY);
125 
126  }
127 #endif
128  if (!Hierarchy_.is_null())
129  Hierarchy_->Iterate(x, y, 1, true);
130  }
132 
133 #ifdef HAVE_MUELU_TPETRA
134  // TO SKIP THE TRAIT IMPLEMENTATION OF XPETRA::MULTIVECTOR
140  void Apply(const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {
141 
142  TEUCHOS_TEST_FOR_EXCEPTION(trans!=NOTRANS, MueLuOpFailure,
143  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
144 
145  //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
146  y.putScalar(0.0);
147 
148 #ifdef HAVE_MUELU_AMGX
149  if (!AMGX_.is_null())
150  AMGX_->apply(x, y);
151 #endif
152 
153  if (!Hierarchy_.is_null()) {
154  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(x);
155 
156  const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
157  Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
158  Hierarchy_->Iterate(tX, tY, 1, true);
159  }
160  }
161 #endif
162 
163  private:
164  RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Hierarchy_;
165 #ifdef HAVE_MUELU_AMGX
166  RCP<MueLu::AMGXOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AMGX_;
167 #endif
168  };
169 
170  template <>
171  class MueLuOp<double, int, int, Xpetra::EpetraNode> :
172  public OperatorT<Xpetra::MultiVector<double, int, int, Xpetra::EpetraNode> >
173 #ifdef HAVE_MUELU_TPETRA
174  // FIXME: guard
175  , public OperatorT<Tpetra::MultiVector<double, int, int, Xpetra::EpetraNode> >
176 #endif
177 #ifdef HAVE_MUELU_EPETRA
178  , public OperatorT<Epetra_MultiVector>
179  , public Belos::Operator<double>
180 #endif
181  {
182  typedef double Scalar;
183  typedef int LocalOrdinal;
184  typedef int GlobalOrdinal;
186 
187  public:
188 
190 #ifdef HAVE_MUELU_AMGX
192 #endif
193  virtual ~MueLuOp() {}
194 
195  void Apply(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {
196 
197  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
198  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
199 
200  //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
201  y.putScalar(0.0);
202 
203 #ifdef HAVE_MUELU_AMGX
204  if (!AMGX_.is_null()) {
205  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Xpetra::toTpetra(x);
206  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY = Xpetra::toTpetra(y);
207 
208  AMGX_->apply(tX, tY);
209  }
210 #endif
211  if (!Hierarchy_.is_null())
212  Hierarchy_->Iterate(x, y, 1, true);
213  }
214 
215 #ifdef HAVE_MUELU_TPETRA
216  void Apply ( const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans=NOTRANS ) const {
217  TEUCHOS_TEST_FOR_EXCEPTION(trans != NOTRANS, MueLuOpFailure,
218  "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");
219 
220  //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
221  y.putScalar(0.0);
222 
223 #ifdef HAVE_MUELU_AMGX
224  if (!AMGX_.is_null())
225  AMGX_->apply(x, y);
226 #endif
227 
228  if (!Hierarchy_.is_null()) {
229  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(x);
230 
231  const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
232  Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
233 
234  tY.putScalar(0.0);
235 
236  Hierarchy_->Iterate(tX, tY, 1, true);
237  }
238  }
239 #endif
240 
241 #ifdef HAVE_MUELU_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 
286 } // namespace Belos
287 
288 #endif // BELOS_MUELU_ADAPTER_HPP
MueLuOp(const RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &H)
void Apply(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &y, ETrans trans=NOTRANS) const
void Apply(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &y, ETrans trans=NOTRANS) const
This routine takes the Tpetra::MultiVector x and applies the operator to it resulting in the Tpetra::...
MueLuOpFailure is thrown when a return value from an MueLu call on an Xpetra::Operator or MueLu::Hier...
RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > AMGX_
void Apply(const Epetra_MultiVector &x, Epetra_MultiVector &y, ETrans trans=NOTRANS) const
This routine takes the Tpetra::MultiVector x and applies the operator to it resulting in the Tpetra::...
void Apply(const Belos::MultiVec< double > &x, Belos::MultiVec< double > &y, ETrans trans=NOTRANS) const
This routine takes the Belos::MultiVec x and applies the operator to it resulting in the Belos::Multi...
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)
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Hierarchy_
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Hierarchy_
MueLuOp(const RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &H)
Default constructor.
MueLuOp(const RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
MueLuOp(const RCP< MueLu::AMGXOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
void Apply(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &x, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &y, ETrans trans=NOTRANS) const
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Adapter for AmgX library from Nvidia.