Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DiagonalQuadraticResponseOnlyModelEvaluator_def.hpp
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Thyra: Interfaces and Support for Abstract Numerical Algorithms
6// Copyright (2004) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
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 Roscoe A. Bartlett (bartlettra@ornl.gov)
39//
40// ***********************************************************************
41// @HEADER
42*/
43
44#ifndef OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
45#define OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
46
47
48#include "Thyra_DiagonalQuadraticResponseOnlyModelEvaluator_decl.hpp"
49#include "Thyra_DiagonalScalarProd.hpp"
50#include "Thyra_VectorStdOps.hpp"
51#include "Thyra_DefaultSpmdVectorSpace.hpp"
52#include "Thyra_DetachedSpmdVectorView.hpp"
53#include "Teuchos_DefaultComm.hpp"
54#include "Teuchos_CommHelpers.hpp"
55#include "Teuchos_Assert.hpp"
56
57
58namespace Thyra {
59
60
61//
62// DiagonalQuadraticResponseOnlyModelEvaluator
63//
64
65
66// Constructors, Initialization, Misc.
67
68
69template<class Scalar>
71 const int localDim,
72 const RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm
73 )
74 :Np_(1), Ng_(1), comm_(comm), localDim_(localDim),
75 nonlinearTermFactor_(0.0), g_offset_(0.0)
76{
77
78 typedef ScalarTraits<Scalar> ST;
79 using Thyra::createMember;
80
81 TEUCHOS_ASSERT( localDim > 0 );
82
83 // Get the comm
84 if (is_null(comm_)) {
86 }
87
88 // Locally replicated space for g
89 g_space_ = Thyra::locallyReplicatedDefaultSpmdVectorSpace<Scalar>(comm_, 1);
90
91 // Distributed space for p
92 p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm_, localDim, -1);
93
94 // Default solution
95 const RCP<Thyra::VectorBase<Scalar> > ps = createMember<Scalar>(p_space_);
96 V_S(ps.ptr(), ST::zero());
97 ps_ = ps;
98
99 // Default diagonal
100 const RCP<Thyra::VectorBase<Scalar> > diag = createMember<Scalar>(p_space_);
101 V_S(diag.ptr(), ST::one());
102 diag_ = diag;
103 diag_bar_ = diag;
104
105 // Default inner product
106 const RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_);
107 V_S(s_bar.ptr(), ST::one());
108 s_bar_ = s_bar;
109
110 // Default response offset
111 g_offset_ = ST::zero();
112
113}
114
115
116template<class Scalar>
118 const RCP<const Thyra::VectorBase<Scalar> > &ps)
119{
120 ps_ = ps.assert_not_null();
121}
122
123
124template<class Scalar>
130
131
132template<class Scalar>
134 const RCP<const Thyra::VectorBase<Scalar> > &diag)
135{
136 diag_ = diag;
137 diag_bar_ = diag;
138}
139
140
141template<class Scalar>
143 const RCP<const Thyra::VectorBase<Scalar> > &diag_bar)
144{
145
147 using Teuchos::rcp_dynamic_cast;
148 using Thyra::createMember;
149 using Thyra::ele_wise_divide;
150 using Thyra::V_S;
151
152 diag_bar_ = diag_bar.assert_not_null();
153
154 // Reset the scalar product for p_space!
155
156 RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_->clone());
157 // NOTE: We have to clone the vector space in order to avoid creating a
158 // circular reference between the space and the vector that defines the
159 // scalar product for the vector space.
160
161 // s_bar[i] = diag[i] / diag_bar[i]
162 V_S( s_bar.ptr(), ST::zero() );
163 ele_wise_divide( ST::one(), *diag_, *diag_bar_, s_bar.ptr() );
164 s_bar_ = s_bar;
165
167 rcp_dynamic_cast<Thyra::ScalarProdVectorSpaceBase<Scalar> >(p_space_, true);
168 sp_p_space->setScalarProd(diagonalScalarProd<Scalar>(s_bar_));
169
170}
171
172
173template<class Scalar>
175 const Scalar &nonlinearTermFactor)
176{
177 nonlinearTermFactor_ = nonlinearTermFactor;
178}
179
180
181template<class Scalar>
183 const Scalar &g_offset)
184{
185 g_offset_ = g_offset;
186}
187
188
189// Public functions overridden from ModelEvaulator
190
191
192template<class Scalar>
197
198
199template<class Scalar>
204
205
206template<class Scalar>
209{
210#ifdef TEUCHOS_DEBUG
212#else
213 (void)l;
214#endif
215 return p_space_;
216}
217
218
219template<class Scalar>
222{
223#ifdef TEUCHOS_DEBUG
225#else
226 (void)j;
227#endif
228 return g_space_;
229}
230
231
232template<class Scalar>
235{
236 typedef Thyra::ModelEvaluatorBase MEB;
237 MEB::InArgsSetup<Scalar> inArgs;
238 inArgs.setModelEvalDescription(this->description());
239 inArgs.set_Np(Np_);
240 return inArgs;
241}
242
243
244// Private functions overridden from ModelEvaulatorDefaultBase
245
246
247template<class Scalar>
250{
251 typedef Thyra::ModelEvaluatorBase MEB;
252 MEB::OutArgsSetup<Scalar> outArgs;
253 outArgs.setModelEvalDescription(this->description());
254 outArgs.set_Np_Ng(Np_,Ng_);
255 outArgs.setSupports(MEB::OUT_ARG_DgDp, 0 ,0, MEB::DERIV_TRANS_MV_BY_ROW);
256 return outArgs;
257}
258
259
260template<class Scalar>
261void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::evalModelImpl(
264 ) const
265{
266
267 using Teuchos::as;
268 using Teuchos::outArg;
270 using Thyra::get_mv;
273 typedef Thyra::ModelEvaluatorBase MEB;
274
275 const ConstDetachedSpmdVectorView<Scalar> p(inArgs.get_p(0));
276 const ConstDetachedSpmdVectorView<Scalar> ps(ps_);
277 const ConstDetachedSpmdVectorView<Scalar> diag(diag_);
278 const ConstDetachedSpmdVectorView<Scalar> s_bar(s_bar_);
279
280 // g(p)
281 if (!is_null(outArgs.get_g(0))) {
282 Scalar g_val = ST::zero();
283 for (Ordinal i = 0; i < p.subDim(); ++i) {
284 const Scalar p_ps = p[i] - ps[i];
285 g_val += diag[i] * p_ps*p_ps;
286 if (nonlinearTermFactor_ != ST::zero()) {
287 g_val += nonlinearTermFactor_ * p_ps * p_ps * p_ps;
288 }
289 }
290 Scalar global_g_val;
292 g_val, outArg(global_g_val) );
293 DetachedSpmdVectorView<Scalar>(outArgs.get_g(0))[0] =
294 as<Scalar>(0.5) * global_g_val + g_offset_;
295 }
296
297 // DgDp[i]
298 if (!outArgs.get_DgDp(0,0).isEmpty()) {
299 const RCP<Thyra::MultiVectorBase<Scalar> > DgDp_trans_mv =
300 get_mv<Scalar>(outArgs.get_DgDp(0,0), "DgDp^T", MEB::DERIV_TRANS_MV_BY_ROW);
301 const DetachedSpmdVectorView<Scalar> DgDp_grad(DgDp_trans_mv->col(0));
302 for (Thyra::Ordinal i = 0; i < p.subDim(); ++i) {
303 const Scalar p_ps = p[i] - ps[i];
304 Scalar DgDp_grad_i = diag[i] * p_ps;
305 if (nonlinearTermFactor_ != ST::zero()) {
306 DgDp_grad_i += as<Scalar>(1.5) * nonlinearTermFactor_ * p_ps * p_ps;
307 }
308 DgDp_grad[i] = DgDp_grad_i / s_bar[i];
309
310 }
311 }
312
313}
314
315
316} // namespace Thyra
317
318
319#endif // OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Ptr< T > ptr() const
Create an explicit detached non-mutable (const) view of all of the local elements on this process of ...
Create an explicit detached mutable (non-const) view of all of the local elements on this process of ...
void setDiagonalVector(const RCP< const VectorBase< Scalar > > &diag)
Set the diagonal vector diag.
void setSolutionVector(const RCP< const VectorBase< Scalar > > &ps)
Set the solution vector ps .
void setNonlinearTermFactor(const Scalar &nonlinearTermFactor)
Set nonlinear term factory.
DiagonalQuadraticResponseOnlyModelEvaluator(const int localDim, const RCP< const Teuchos::Comm< Ordinal > > &comm=Teuchos::null)
void setDiagonalBarVector(const RCP< const VectorBase< Scalar > > &diag_bar)
Set the diagonal vector diag_bar.
const RCP< const VectorBase< Scalar > > getSolutionVector() const
Get the solution vector ps .
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
RCP< const VectorBase< Scalar > > get_p(int l) const
Get p(l) where 0 <= l && l < this->Np().
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object.
Derivative< Scalar > get_DgDp(int j, int l) const
Precondition: supports(OUT_ARG_DgDp,j,l)==true.
Evaluation< VectorBase< Scalar > > get_g(int j) const
Precondition: supports(OUT_ARG_g)==true..
Base subclass for ModelEvaluator that defines some basic types.
Abstract interface for finite-dimensional dense vectors.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
TypeTo as(const TypeFrom &t)
T_To & dyn_cast(T_From &from)