Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Xpetra_EpetraOperator.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
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 XPETRA_EPETRAOPERATOR_HPP
47 #define XPETRA_EPETRAOPERATOR_HPP
48 
50 
51 #include <Epetra_Operator.h>
52 #include <Epetra_Map.h>
53 
54 #include "Xpetra_Map.hpp"
55 #include "Xpetra_EpetraMap.hpp"
56 #include "Xpetra_MultiVector.hpp"
58 #include "Xpetra_Operator.hpp"
59 
60 #include "Xpetra_Utils.hpp"
61 
62 namespace Xpetra {
63 
64 template <class EpetraGlobalOrdinal, class Node>
65 class EpetraOperator : public Operator<double, int, EpetraGlobalOrdinal, Node> {
66  typedef double Scalar;
67  typedef int LocalOrdinal;
68  typedef EpetraGlobalOrdinal GlobalOrdinal;
69 
70  public:
72 
74  virtual const Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > getDomainMap() const {
75  XPETRA_MONITOR("EpetraOperator::getDomainMap()");
76  return toXpetra<GlobalOrdinal, Node>(op_->OperatorDomainMap());
77  }
78 
80  virtual const Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > getRangeMap() const {
81  XPETRA_MONITOR("EpetraOperator::getRangeMap()");
82  return toXpetra<GlobalOrdinal, Node>(op_->OperatorRangeMap());
83  }
84 
86 
91  virtual void
94  Teuchos::ETransp mode = Teuchos::NO_TRANS,
95  Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
96  Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const {
97  XPETRA_MONITOR("EpetraOperator::apply");
98 
99  XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, X, eX, "Xpetra::EpetraOperator->apply(): cannot cast input to Xpetra::EpetraMultiVectorT");
100  XPETRA_DYNAMIC_CAST(EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, Y, eY, "Xpetra::EpetraOperator->apply(): cannot cast input to Xpetra::EpetraMultiVectorT");
101 
102  TEUCHOS_TEST_FOR_EXCEPTION((mode != Teuchos::NO_TRANS) && (mode != Teuchos::TRANS), Exceptions::NotImplemented,
103  "Xpetra::EpetraOperator->apply(): can only accept mode == NO_TRANS or mode == TRANS");
104  TEUCHOS_TEST_FOR_EXCEPTION(mode == Teuchos::TRANS && !hasTransposeApply(), Exceptions::RuntimeError,
105  "Xpetra::EpetraOperator->apply(): cannot apply transpose as underlying Epetra operator does not support it");
106 
107  // Helper vector for string A*X
108  RCP<Epetra_MultiVector> epY = eY.getEpetra_MultiVector();
109  RCP<Epetra_MultiVector> tmp = Teuchos::rcp(new Epetra_MultiVector(*epY));
110  tmp->PutScalar(0.0);
111 
112  op_->SetUseTranspose(mode == Teuchos::TRANS);
113  XPETRA_ERR_CHECK(op_->Apply(*eX.getEpetra_MultiVector(), *tmp));
114 
115  // calculate alpha * A * x + beta * y
116  XPETRA_ERR_CHECK(eY.getEpetra_MultiVector()->Update(alpha, *tmp, beta));
117  }
118 
120  virtual bool hasTransposeApply() const {
121  // We do not currently use transpose, try setting it
122  int err = op_->SetUseTranspose(true);
123  op_->SetUseTranspose(false);
124  return (err == 0);
125  }
126 
128 
130 
131 
133  std::string description() const {
134  XPETRA_MONITOR("EpetraOperator::description");
135  return "Epetra_Operator";
136  }
137 
139  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
140  XPETRA_MONITOR("EpetraOperator::describe");
141  out << "Epetra_Operator" << std::endl;
142  }
143 
145 
147 
148 
150  EpetraOperator(const Teuchos::RCP<Epetra_Operator> &op)
151  : op_(op) {} // TODO removed const
152 
157  using STS = Teuchos::ScalarTraits<Scalar>;
158  R.update(STS::one(), B, STS::zero());
159  this->apply(X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
160  }
161 
163 
164  private:
166  RCP<Epetra_Operator> op_;
167 
168 }; // EpetraOperator class
169 
170 template <class EpetraGlobalOrdinal, class Node>
171 class EpetraInverseOperator : public Operator<double, int, EpetraGlobalOrdinal, Node> {
172  typedef double Scalar;
173  typedef int LocalOrdinal;
174  typedef EpetraGlobalOrdinal GlobalOrdinal;
175 
176  public:
178 
180  virtual Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getDomainMap() const {
181  XPETRA_MONITOR("EpetraOperator::getDomainMap()");
182  return toXpetra<GlobalOrdinal, Node>(op_->OperatorDomainMap());
183  }
184 
186  virtual Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getRangeMap() const {
187  XPETRA_MONITOR("EpetraOperator::getRangeMap()");
188  return toXpetra<GlobalOrdinal, Node>(op_->OperatorRangeMap());
189  }
190 
192 
197  virtual void
200  Teuchos::ETransp mode = Teuchos::NO_TRANS,
201  Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
202  Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const {
203  XPETRA_MONITOR("EpetraOperator::apply");
204 
205  XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, X, eX, "Xpetra::EpetraOperator->apply(): cannot cast input to Xpetra::EpetraMultiVectorT");
206  XPETRA_DYNAMIC_CAST(EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, Y, eY, "Xpetra::EpetraOperator->apply(): cannot cast input to Xpetra::EpetraMultiVectorT");
207 
208  TEUCHOS_TEST_FOR_EXCEPTION((mode != Teuchos::NO_TRANS) && (mode != Teuchos::TRANS), Exceptions::NotImplemented,
209  "Xpetra::EpetraOperator->apply(): can only accept mode == NO_TRANS or mode == TRANS");
210  TEUCHOS_TEST_FOR_EXCEPTION(mode == Teuchos::TRANS && !hasTransposeApply(), Exceptions::RuntimeError,
211  "Xpetra::EpetraOperator->apply(): cannot apply transpose as underlying Epetra operator does not support it");
212 
213  // Helper vector for string A*X
214  RCP<Epetra_MultiVector> epY = eY.getEpetra_MultiVector();
215  RCP<Epetra_MultiVector> tmp = Teuchos::rcp(new Epetra_MultiVector(*epY));
216  tmp->PutScalar(0.0);
217 
218  op_->SetUseTranspose(mode == Teuchos::TRANS);
219  XPETRA_ERR_CHECK(op_->ApplyInverse(*eX.getEpetra_MultiVector(), *tmp));
220 
221  // calculate alpha * A * x + beta * y
222  XPETRA_ERR_CHECK(eY.getEpetra_MultiVector()->Update(alpha, *tmp, beta));
223  }
224 
226  virtual bool hasTransposeApply() const {
227  // We do not currently use transpose, try setting it
228  int err = op_->SetUseTranspose(true);
229  op_->SetUseTranspose(false);
230  return (err == 0);
231  }
232 
234 
236 
237 
239  std::string description() const {
240  XPETRA_MONITOR("EpetraOperator::description");
241  return "Epetra_Operator";
242  }
243 
245  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
246  XPETRA_MONITOR("EpetraOperator::describe");
247  out << "Epetra_Operator" << std::endl;
248  }
249 
251 
253 
254 
256  EpetraInverseOperator(const Teuchos::RCP<Epetra_Operator> &op)
257  : op_(op) {} // TODO removed const
258 
263  using STS = Teuchos::ScalarTraits<Scalar>;
264  R.update(STS::one(), B, STS::zero());
265  this->apply(X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
266  }
267 
269 
270  private:
272  RCP<Epetra_Operator> op_;
273 };
274 
275 } // namespace Xpetra
276 
277 #endif // XPETRA_EPETRAOPERATOR_HPP
std::string description() const
A simple one-line description of this object.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with the given verbosity level to a FancyOStream.
RCP< Epetra_Operator > op_
The Tpetra::Operator which this class wraps.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
The Map associated with the range of this operator, which must be compatible with Y...
Exception throws to report errors in the internal logical of the program.
EpetraInverseOperator(const Teuchos::RCP< Epetra_Operator > &op)
EpetraOperator constructor to wrap a Epetra_Operator object.
virtual void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta)=0
Update multi-vector values with scaled values of A, this = beta*this + alpha*A.
RCP< Epetra_Operator > op_
The Tpetra::Operator which this class wraps.
#define XPETRA_ERR_CHECK(arg)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with the given verbosity level to a FancyOStream.
void residual(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &R) const
Compute a residual R = B - (*this) * X.
EpetraOperator(const Teuchos::RCP< Epetra_Operator > &op)
EpetraOperator constructor to wrap a Epetra_Operator object.
virtual bool hasTransposeApply() const
Whether this operator supports applying the transpose or conjugate transpose.
Exception throws when you call an unimplemented method of Xpetra.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
The Map associated with the domain of this operator, which must be compatible with X...
#define XPETRA_DYNAMIC_CAST(type, obj, newObj, exceptionMsg)
virtual void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
Computes the operator-multivector application.
virtual bool hasTransposeApply() const
Whether this operator supports applying the transpose or conjugate transpose.
EpetraGlobalOrdinal GlobalOrdinal
virtual const Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
The Map associated with the domain of this operator, which must be compatible with X...
virtual const Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
The Map associated with the range of this operator, which must be compatible with Y...
#define XPETRA_MONITOR(funcName)
std::string description() const
A simple one-line description of this object.
virtual void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
Computes the operator-multivector application.
void residual(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &R) const
Compute a residual R = B - (*this) * X.