Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_TpetraExplicitAdjointModelEvaluator.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef Thyra_TpetraExplicitAdjointModelEvaluator_hpp
43#define Thyra_TpetraExplicitAdjointModelEvaluator_hpp
44
45#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
46#include "Thyra_TpetraLinearOp.hpp"
47#include "Tpetra_RowMatrixTransposer.hpp"
48
49namespace Thyra {
50
61 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal=LocalOrdinal, typename Node=KokkosClassic::DefaultNode::DefaultNodeType>
63 public ModelEvaluatorDelegatorBase<Scalar>{
64public:
65
70
75
78
80 {
81 // This ME should use the same InArgs as the underlying model. However
82 // we can't just use it's InArgs directly because the description won't
83 // match (which is checked in debug builds). Instead create a new
84 // InArgsSetup initialized by the underlying model's createInArgs() and
85 // set the description appropriately.
87 this->getUnderlyingModel()->createInArgs();
89 return inArgs;
90 }
91
94 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> TO;
95 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TCM;
96 typedef Tpetra::RowMatrixTransposer<Scalar,LocalOrdinal,GlobalOrdinal,Node> TRMT;
97 if (thyra_fwd_op == Teuchos::null)
98 thyra_fwd_op = this->getUnderlyingModel()->create_W_op();
99 RCP<TLO> thyra_tpetra_fwd_op =
100 Teuchos::rcp_dynamic_cast<TLO>(thyra_fwd_op,true);
101 RCP<TO> tpetra_fwd_op = thyra_tpetra_fwd_op->getTpetraOperator();
102 RCP<TCM> tpetra_fwd_mat =
103 Teuchos::rcp_dynamic_cast<TCM>(tpetra_fwd_op,true);
104 TRMT transposer(tpetra_fwd_mat);
105 RCP<TCM> tpetra_trans_mat = transposer.createTranspose();
106 return tpetraLinearOp(thyra_op->range(), thyra_op->domain(),
107 tpetra_trans_mat);
108 }
109
110private:
111
112 mutable RCP< Thyra::LinearOpBase<Scalar> > thyra_fwd_op;
113
114 ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const
115 {
116 typedef ModelEvaluatorBase MEB;
117 MEB::OutArgs<Scalar> model_outArgs =
118 this->getUnderlyingModel()->createOutArgs();
119 MEB::OutArgsSetup<Scalar> outArgs;
120 outArgs.setModelEvalDescription(this->description());
121 outArgs.setSupports(MEB::OUT_ARG_f); //All models must support f, apparently
122 outArgs.setSupports(MEB::OUT_ARG_W_op);
123 TEUCHOS_ASSERT(model_outArgs.supports(MEB::OUT_ARG_W_op));
124 return outArgs;
125 }
126
127 void evalModelImpl(
128 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
129 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
130 {
131 typedef ModelEvaluatorBase MEB;
132 typedef TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> TLO;
133 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> TO;
134 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TCM;
135 typedef Tpetra::RowMatrixTransposer<Scalar,LocalOrdinal,GlobalOrdinal,Node> TRMT;
136
137 MEB::OutArgs<Scalar> model_outArgs =
138 this->getUnderlyingModel()->createOutArgs();
139
140 if (model_outArgs.supports(MEB::OUT_ARG_W_op) &&
141 outArgs.get_W_op() != Teuchos::null) {
142 // Compute underlying W_op.
143 if (thyra_fwd_op == Teuchos::null)
144 thyra_fwd_op = this->getUnderlyingModel()->create_W_op();
145 model_outArgs->set_W_op(thyra_fwd_op);
146 this->getUnderlyingModel()->evalModel(inArgs, model_outArgs);
147
148 // Transpose W_op
149 // Unfortunately, because of the horrendous design of the Tpetra
150 // RowMatrixTransposer, this creates a new transposed matrix each time
151 RCP<TLO> thyra_tpetra_fwd_op =
152 Teuchos::rcp_dynamic_cast<TLO>(thyra_fwd_op,true);
153 RCP<TO> tpetra_fwd_op = thyra_tpetra_fwd_op->getTpetraOperator();
154 RCP<TCM> tpetra_fwd_mat =
155 Teuchos::rcp_dynamic_cast<TCM>(tpetra_fwd_op,true);
156 TRMT transposer(tpetra_fwd_mat);
157 RCP<TCM> tpetra_trans_mat = transposer.createTranspose();
158
159 // Copy transposed matrix into our outArg
160 RCP<LOB> thyra_adj_op = outArgs.get_W_op();
161 RCP<TLO> thyra_tpetra_adj_op =
162 Teuchos::rcp_dynamic_cast<TLO>(thyra_adj_op,true);
163 RCP<TO> tpetra_adj_op = thyra_tpetra_adj_op->getTpetraOperator();
164 RCP<TCM> tpetra_adj_mat =
165 Teuchos::rcp_dynamic_cast<TCM>(tpetra_adj_op,true);
166 *tpetra_adj_mat = *tpetra_trans_mat;
167 }
168 }
169
170};
171
172template <typename Scalar>
173RCP<TpetraExplicitAdjointModelEvaluator<Scalar> >
174tpetraExplicitAdjointModelEvaluator(
175 const RCP<const ModelEvaluator<Scalar> >& model)
176{
177 return Teuchos::rcp(new TpetraExplicitAdjointModelEvaluator<Scalar>(model));
178}
179
180template <typename Scalar>
181RCP<TpetraExplicitAdjointModelEvaluator<Scalar> >
182tpetraExplicitAdjointModelEvaluator(
183 const RCP<ModelEvaluator<Scalar> >& model)
184{
185 return Teuchos::rcp(new TpetraExplicitAdjointModelEvaluator<Scalar>(model));
186}
187
188} // namespace Thyra
189
190#endif
virtual std::string description() const
Protected subclass of InArgs that only ModelEvaluator subclasses can access to set up the selection o...
void setModelEvalDescription(const std::string &modelEvalDescription)
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object.
Base subclass for ModelEvaluator that defines some basic types.
This is a base class that delegetes almost all function to a wrapped model evaluator object.
virtual RCP< const ModelEvaluator< Scalar > > getUnderlyingModel() const
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
A model evaluator decorator for computing an explicit adjoint.
TpetraExplicitAdjointModelEvaluator(const RCP< const ModelEvaluator< Scalar > > &model)
Constructor.
TpetraExplicitAdjointModelEvaluator(const RCP< ModelEvaluator< Scalar > > &model)
Constructor.
virtual ~TpetraExplicitAdjointModelEvaluator()=default
Destructor.
Concrete Thyra::LinearOpBase subclass for Tpetra::Operator.
#define TEUCHOS_ASSERT(assertion_test)
T_To & dyn_cast(T_From &from)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)