Thyra  Version of the Day
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 
58 namespace Thyra {
59 
60 
61 //
62 // DiagonalQuadraticResponseOnlyModelEvaluator
63 //
64 
65 
66 // Constructors, Initialization, Misc.
67 
68 
69 template<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_)) {
85  comm_ = Teuchos::DefaultComm<Thyra::Ordinal>::getComm();
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 
116 template<class Scalar>
118  const RCP<const Thyra::VectorBase<Scalar> > &ps)
119 {
120  ps_ = ps.assert_not_null();
121 }
122 
123 
124 template<class Scalar>
125 const RCP<const Thyra::VectorBase<Scalar> >
127 {
128  return ps_;
129 }
130 
131 
132 template<class Scalar>
134  const RCP<const Thyra::VectorBase<Scalar> > &diag)
135 {
136  diag_ = diag;
137  diag_bar_ = diag;
138 }
139 
140 
141 template<class Scalar>
143  const RCP<const Thyra::VectorBase<Scalar> > &diag_bar)
144 {
145 
146  typedef Teuchos::ScalarTraits<Scalar> ST;
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 
166  const RCP<Thyra::ScalarProdVectorSpaceBase<Scalar> > sp_p_space =
167  rcp_dynamic_cast<Thyra::ScalarProdVectorSpaceBase<Scalar> >(p_space_, true);
168  sp_p_space->setScalarProd(diagonalScalarProd<Scalar>(s_bar_));
169 
170 }
171 
172 
173 template<class Scalar>
175  const Scalar &nonlinearTermFactor)
176 {
177  nonlinearTermFactor_ = nonlinearTermFactor;
178 }
179 
180 
181 template<class Scalar>
183  const Scalar &g_offset)
184 {
185  g_offset_ = g_offset;
186 }
187 
188 
189 // Public functions overridden from ModelEvaulator
190 
191 
192 template<class Scalar>
194 {
195  return Np_;
196 }
197 
198 
199 template<class Scalar>
201 {
202  return Ng_;
203 }
204 
205 
206 template<class Scalar>
207 RCP<const Thyra::VectorSpaceBase<Scalar> >
209 {
210 #ifdef TEUCHOS_DEBUG
211  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
212 #endif
213  return p_space_;
214 }
215 
216 
217 template<class Scalar>
218 RCP<const Thyra::VectorSpaceBase<Scalar> >
220 {
221 #ifdef TEUCHOS_DEBUG
222  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
223 #endif
224  return g_space_;
225 }
226 
227 
228 template<class Scalar>
231 {
232  typedef Thyra::ModelEvaluatorBase MEB;
233  MEB::InArgsSetup<Scalar> inArgs;
234  inArgs.setModelEvalDescription(this->description());
235  inArgs.set_Np(Np_);
236  return inArgs;
237 }
238 
239 
240 // Private functions overridden from ModelEvaulatorDefaultBase
241 
242 
243 template<class Scalar>
246 {
247  typedef Thyra::ModelEvaluatorBase MEB;
248  MEB::OutArgsSetup<Scalar> outArgs;
249  outArgs.setModelEvalDescription(this->description());
250  outArgs.set_Np_Ng(Np_,Ng_);
251  outArgs.setSupports(MEB::OUT_ARG_DgDp, 0 ,0, MEB::DERIV_TRANS_MV_BY_ROW);
252  return outArgs;
253 }
254 
255 
256 template<class Scalar>
257 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::evalModelImpl(
260  ) const
261 {
262 
263  using Teuchos::as;
264  using Teuchos::outArg;
265  typedef Teuchos::ScalarTraits<Scalar> ST;
266  using Thyra::get_mv;
269  typedef Thyra::ModelEvaluatorBase MEB;
270 
271  const ConstDetachedSpmdVectorView<Scalar> p(inArgs.get_p(0));
272  const ConstDetachedSpmdVectorView<Scalar> ps(ps_);
273  const ConstDetachedSpmdVectorView<Scalar> diag(diag_);
274  const ConstDetachedSpmdVectorView<Scalar> s_bar(s_bar_);
275 
276  // g(p)
277  if (!is_null(outArgs.get_g(0))) {
278  Scalar g_val = ST::zero();
279  for (Ordinal i = 0; i < p.subDim(); ++i) {
280  const Scalar p_ps = p[i] - ps[i];
281  g_val += diag[i] * p_ps*p_ps;
282  if (nonlinearTermFactor_ != ST::zero()) {
283  g_val += nonlinearTermFactor_ * p_ps * p_ps * p_ps;
284  }
285  }
286  Scalar global_g_val;
287  Teuchos::reduceAll<Ordinal, Scalar>(*comm_, Teuchos::REDUCE_SUM,
288  g_val, outArg(global_g_val) );
289  DetachedSpmdVectorView<Scalar>(outArgs.get_g(0))[0] =
290  as<Scalar>(0.5) * global_g_val + g_offset_;
291  }
292 
293  // DgDp[i]
294  if (!outArgs.get_DgDp(0,0).isEmpty()) {
295  const RCP<Thyra::MultiVectorBase<Scalar> > DgDp_trans_mv =
296  get_mv<Scalar>(outArgs.get_DgDp(0,0), "DgDp^T", MEB::DERIV_TRANS_MV_BY_ROW);
297  const DetachedSpmdVectorView<Scalar> DgDp_grad(DgDp_trans_mv->col(0));
298  for (Thyra::Ordinal i = 0; i < p.subDim(); ++i) {
299  const Scalar p_ps = p[i] - ps[i];
300  Scalar DgDp_grad_i = diag[i] * p_ps;
301  if (nonlinearTermFactor_ != ST::zero()) {
302  DgDp_grad_i += as<Scalar>(1.5) * nonlinearTermFactor_ * p_ps * p_ps;
303  }
304  DgDp_grad[i] = DgDp_grad_i / s_bar[i];
305 
306  }
307  }
308 
309 }
310 
311 
312 } // namespace Thyra
313 
314 
315 #endif // OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
Create an explicit detached mutable (non-const) view of all of the local elements on this process of ...
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...
const RCP< const VectorBase< Scalar > > getSolutionVector() const
Get the solution vector ps .
Derivative< Scalar > get_DgDp(int j, int l) const
Precondition: supports(OUT_ARG_DgDp,j,l)==true.
void setSolutionVector(const RCP< const VectorBase< Scalar > > &ps)
Set the solution vector ps .
DiagonalQuadraticResponseOnlyModelEvaluator(const int localDim, const RCP< const Teuchos::Comm< Ordinal > > &comm=Teuchos::null)
RCP< const VectorBase< Scalar > > get_p(int l) const
Get p(l) where 0 <= l && l < this->Np().
void setNonlinearTermFactor(const Scalar &nonlinearTermFactor)
Set nonlinear term factory.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Abstract interface for finite-dimensional dense vectors.
Evaluation< VectorBase< Scalar > > get_g(int j) const
Precondition: supports(OUT_ARG_g)==true..
void setDiagonalVector(const RCP< const VectorBase< Scalar > > &diag)
Set the diagonal vector diag.
Base subclass for ModelEvaluator that defines some basic types.
virtual void setScalarProd(const RCP< const ScalarProdBase< Scalar > > &scalarProd)
Set a different scalar product.
void setDiagonalBarVector(const RCP< const VectorBase< Scalar > > &diag_bar)
Set the diagonal vector diag_bar.
Create an explicit detached non-mutable (const) view of all of the local elements on this process of ...
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...