MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_XpetraOperator_decl.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 MUELU_XPETRAOPERATOR_DECL_HPP
47 #define MUELU_XPETRAOPERATOR_DECL_HPP
48 
49 #include "MueLu_ConfigDefs.hpp"
50 
51 #include <Xpetra_Operator.hpp>
52 #include <Xpetra_MultiVector.hpp>
53 #include "MueLu_Level.hpp"
54 #include "MueLu_Hierarchy_decl.hpp"
55 
56 namespace MueLu {
57 
60  template <class Scalar = DefaultScalar,
63  class Node = DefaultNode>
64  class XpetraOperator : public Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
65  protected:
67  public:
68 
70 
71 
74 
76  virtual ~XpetraOperator() { }
77 
79 
82  typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> Matrix;
83 
84  RCP<Matrix> A = Hierarchy_->GetLevel(0)->template Get<RCP<Matrix> >("A");
85  return A->getDomainMap();
86  }
87 
90  typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> Matrix;
91 
92  RCP<Matrix> A = Hierarchy_->GetLevel(0)->template Get<RCP<Matrix> >("A");
93  return A->getRangeMap();
94  }
95 
97 
101  void apply(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
102  Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
103  Teuchos::ETransp /* mode */ = Teuchos::NO_TRANS,
105  Scalar /* beta */ = Teuchos::ScalarTraits<Scalar>::one()) const{
106  try {
107 #ifdef HAVE_MUELU_DEBUG
108  typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> Matrix;
109  RCP<Matrix> A = Hierarchy_->GetLevel(0)->template Get<RCP<Matrix> >("A");
110 
111  // X is supposed to live in the range map of the operator (const rhs = B)
113  Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getRangeMap(),X.getNumVectors());
115  Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getDomainMap(),Y.getNumVectors());
116  TEUCHOS_TEST_FOR_EXCEPTION(A->getRangeMap()->isSameAs(*(Xop->getMap())) == false, std::logic_error,
117  "MueLu::XpetraOperator::apply: map of X is incompatible with range map of A");
118  TEUCHOS_TEST_FOR_EXCEPTION(A->getDomainMap()->isSameAs(*(Yop->getMap())) == false, std::logic_error,
119  "MueLu::XpetraOperator::apply: map of Y is incompatible with domain map of A");
120 #endif
121 
123  Hierarchy_->Iterate(X, Y, 1, true);
124  } catch (std::exception& e) {
125  //FIXME add message and rethrow
126  std::cerr << "Caught an exception in MueLu::XpetraOperator::apply():" << std::endl
127  << e.what() << std::endl;
128  }
129  }
130 
132  bool hasTransposeApply() const { return false; }
133 
135  void residual(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & X,
136  const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & B,
137  Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & R) const {
138  using STS = Teuchos::ScalarTraits<Scalar>;
139  R.update(STS::one(),B,STS::zero());
140  this->apply (X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
141  }
142 
144 
145 
148 
150 
151  private:
153  };
154 
155 } // namespace
156 
157 #endif // MUELU_XPETRAOPERATOR_DECL_HPP
void apply(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp=Teuchos::NO_TRANS, Scalar=Teuchos::ScalarTraits< Scalar >::one(), Scalar=Teuchos::ScalarTraits< Scalar >::one()) const
Returns in Y the result of a Xpetra::Operator applied to a Xpetra::MultiVector X. ...
MueLu::DefaultLocalOrdinal LocalOrdinal
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
MueLu::DefaultNode Node
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetHierarchy() const
Direct access to the underlying MueLu::Hierarchy.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Hierarchy_
void residual(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &R) const
Compute a residual R = B - (*this) * X.
XpetraOperator(const RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &H)
Constructor.
Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
virtual ~XpetraOperator()
Destructor.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.