ROL
ROL_Constraint_SerialSimOpt.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) 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 lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44
45#pragma once
46#ifndef ROL_CONSTRAINT_SERIALSIMOPT_HPP
47#define ROL_CONSTRAINT_SERIALSIMOPT_HPP
48
49namespace ROL {
50
76
77
78template<typename Real>
80
81 using V = Vector<Real>;
84
86
87private:
88
89 const Ptr<Con> con_; // Constraint for single time step
90 const Ptr<V> ui_; // Initial condition
91 Ptr<V> u0_; // Zero sim vector on time step
92 Pre<V> z0_; // Zero opt vector on time step
93
94 VectorWorkspace<Real> workspace_; // Memory management
95
96protected:
97
98 PV& partition( V& x ) const { return static_cast<PV&>(x); }
99
100 const V& partition( const V& x ) const { return static_cast<const PV&>(x); }
101
102
103public:
104
105 Constraint_SerialSimOpt( const Ptr<Con>& con, const V& ui, const V& zi ) :
106 Constraint_SimOpt<Real>(), con_(con),
107 ui_( workspace_.copy( ui ) ),
108 u0_( workspace_.clone( ui ) ),
109 z0_( workspace_.clone( ui ) ) {
110 u0_->zero(); z0_->zero();
111 }
112
113// virtual void update_1( const V& u, const V& z, bool flag = true, int iter = -1 ) override {
114//
115// }
116//
117// virtual void update_2( const V& u, const V& z, bool flag = true, int iter = -1 ) override {
118//
119// }
120
121 virtual void value( V& c, const V& u, const V& z, Real& tol ) override {
122
123 auto& cp = partition(c);
124 auto& up = partition(u);
125 auto& zp = partition(z);
126
127 // First interval value uses initial condition
128 con_->value( *(cp.get(0)), *ui_, *(zp.get(0)), tol );
129
130 for( size_type k=1, k<cp.numVectors(), ++k ) {
131 con_->value( *(cp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
132 }
133 }
134
135 virtual void solve( V& c, V& u, const V& z, Real& tol ) override {
136
137 auto& cp = partition(c);
138 auto& up = partition(u);
139 auto& zp = partition(z);
140
141 // First interval value uses initial condition
142 con_->solve( *(cp.get(0)), *ui_, *(zp.get(0)), tol );
143
144 for( size_type k=1, k<cp.numVectors(), ++k ) {
145 con_->solve( *(cp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
146 }
147 }
148
149 virtual void applyJacobian_1( V& jv, const V& v, const V& u, const V& z, Real& tol ) override {
150
151 auto& jvp = partition(jv); auto& vp = partition(v);
152 auto& up = partition(u); auto& zp = partition(z);
153
154 auto tmp = workspace_.clone( up.get(0) );
155
156 con_->applyJacobian_1_new( *(jvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)), tol );
157
158 for( size_type k=1; k<jvp.numVectors(); ++jvp ) {
159 con_->applyJacobian_1_new( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
160 con_->applyJacobian_1_old( *tmp, *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
161 jvp.plus(tmp);
162 }
163 } // applyJacobian_1
164
165 virtual void applyJacobian_2( V& jv, const V& v, const V& u, const V& z, Real& tol) override {
166
167 auto& jvp = partition(jv); auto& vp = partition(v);
168 auto& up = partition(u); auto& zp = partition(z);
169
170 for( size_type k=0; k<jvp.numVectors(); ++k ) {
171 con_->applyJacobian_2( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
172 }
173 } // applyJacobian_2
174
175
176 virtual void applyInverseJacobian_1( V& ijv, const V& v, const V& u,
177 const V& z, Real& tol ) override {
178
179 auto& ijvp = partition(ijv); auto& vp = partition(v);
180 auto& up = partition(u); auto& zp = partition(z);
181
182 // First step
183 con_->applyInverseJacobian_1_new( *(ijvp.get(0)), *(vp.get(0), *u0_, *(up.get(0)), *(zp.get(0)), tol );
184
185 auto tmp = workspace.clone( vp.get(0) );
186
187 for( size_type k=1; k<ijvp.numVectors(); ++k ) {
188
189 con_->applyJacobian_1_old( *tmp, *(ijvp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
190 tmp->scale(-1.0);
191 tmp->plus( *(vp.get(k) );
192
193 con_->applyInverseJacobian_1_new( *(ijvp.get(k)), *tmp, *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
194 }
195 } // applyInverseJacobian_1
196
197
198
199 virtual void applyAdjointJacobian_1( V& ajv, const V& v, const V& u,
200 const V& z, Real& tol) override {
201
202 auto& ajvp = partition(ajv); auto& vp = partition(v);
203 auto& up = partition(u); auto& zp = partition(z);
204
205 auto tmp = workspace.clone( ajvp.get(0) );
206
207 size_type N = ajvp.numVectors()-1;
208
209 // First step
210 con_->applyAdjointJacobian_1_new( *(ajvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)) );
211 con_->applyAdjointJacobian_1_old( *(ajvp.get(0)), *(vp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)) );
212
213 for( size_type k=1; k<N; ++k ) {
214 con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
215 con_->applyAdjointJacobian_1_old( *tmp, *(vp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
216 ajvp.plus(*tmp);
217 }
218
219 // Last step
220 con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
221
222 } // applyAdjointJacobian_1
223
224
225
226 virtual void applyAdjointJacobian_1( V& ajv, const V& v, const V& u, const V &z,
227 const V& dualv, Real& tol) override {
228
229 auto& ajvp = partition(ajv); auto& vp = partition(dualv);
230 auto& up = partition(u); auto& zp = partition(z);
231 auto tmp = workspace.clone( ajvp.get(0) );
232
233 size_type N = ajvp.numVectors()-1;
234
235 // First step
236 con_->applyAdjointJacobian_1_new( *(ajvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)) );
237 con_->applyAdjointJacobian_1_old( *(ajvp.get(0)), *(vp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)) );
238
239 for( size_type k=1; k<N; ++k ) {
240 con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
241 con_->applyAdjointJacobian_1_old( *tmp, *(vp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
242 ajvp.plus(*tmp);
243 }
244
245 // Last step
246 con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
247
248 } // applyAdjointJacobian_1
249
250
251
252 virtual void applyAdjointJacobian_2( V& ajv, const V& v, const V& u,
253 const V& z, Real& tol) override {
254 auto& ajvp = partition(ajv); auto& vp = partition(v);
255 auto& up = partition(u); auto& zp = partition(z);
256
257 for( size_type k=0; k<jvp.numVectors(); ++k ) {
258 con_->applyAdjointJacobian_2_time( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
259 }
260 } // applyAdjointJacobian_2
261
262
263 virtual void applyAdjointJacobian_2( V& ajv, const V& v, const V& u, const V& z,
264 const V& dualv, Real& tol) override {
265
266 auto& ajvp = partition(ajv); auto& vp = partition(v);
267 auto& up = partition(u); auto& zp = partition(z);
268
269 for( size_type k=0; k<jvp.numVectors(); ++k ) {
270 con_->applyAdjointJacobian_2_time( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
271 }
272 }
273
274
275 virtual void applyInverseAdjointJacobian_1( V& iajv, const V& v, const V& u,
276 const V& z, Real& tol) override {
277
278 auto& iajvp = partition(iajv); auto& vp = partition(v);
279 auto& up = partition(u); auto& zp = partition(z);
280
281 size_type N = iajvp.numVectors()-1;
282
283 // Last Step
284 con_->applyInverseAdjointJacobian_1_new( *(iajvp.get(N)), *(vp.get(N)), *(up.get(N-1)), *(up.get(N)), *(zp.get(N)), tol );
285
286 auto tmp = workspace_.clone( vp.get(0) );
287
288 for( size_type k=N-1; k>0; --k ) {
289 con_->applyAdjointJacobian_1_old( *tmp, *(iajvp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
290
291 tmp->scale( -1.0 );
292 tmp->plus( *(vp.get(k) );
293
294 con_->applyAdjointJacobian_1_new( *(iajvp.get(k)), *tmp, *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
295 }
296
297 // "First step"
298 con_->applyAdjointJacobian_1_old( *tmp, *(iajvp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)), tol );
299
300 tmp->scale( -1.0 );
301 tmp->plus( *(vp.get(0) );
302
303 con_->applyAdjointJacobian_1_new( *(iajvp.get(0)), *tmp, *u0_, *(up.get(0)), *(zp.get(0)), tol );
304 }
305
306
307 virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v,
308 const V& u, const V& z, Real& tol) override {
309 // TODO
310 }
311
312
313 virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v,
314 const V& u, const V& z, Real& tol) override {
315 ahwv.zero();
316 }
317
318 virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v,
319 const V& u, const V& z, Real& tol) override {
320 ahwv.zero();
321 }
322
323 virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v,
324 const V& u, const V& z, Real& tol) override {
325 ahwv.zero();
326 }
327
328
329
330
331
332
333
334
335
336
337}; // class Constraint_SerialSimOpt
338
339
340
341} // namespace ROL
342
343
344#endif // ROL_CONSTRAINT_SERIALSIMOPT_HPP
345
Unifies the constraint defined on a single time step that are defined through the Constraint_TimeSimO...
virtual void applyJacobian_1(V &jv, const V &v, const V &u, const V &z, Real &tol) override
Apply the partial constraint Jacobian at , , to the vector .
Constraint_SerialSimOpt(const Ptr< Con > &con, const V &ui, const V &zi)
virtual void applyJacobian_2(V &jv, const V &v, const V &u, const V &z, Real &tol) override
Apply the partial constraint Jacobian at , , to the vector .
virtual void applyInverseJacobian_1(V &ijv, const V &v, const V &u, const V &z, Real &tol) override
Apply the inverse partial constraint Jacobian at , , to the vector .
virtual void applyAdjointJacobian_2(V &ajv, const V &v, const V &u, const V &z, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
virtual void applyAdjointJacobian_1(V &ajv, const V &v, const V &u, const V &z, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
virtual void applyInverseAdjointJacobian_1(V &iajv, const V &v, const V &u, const V &z, Real &tol) override
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector .
virtual void solve(V &c, V &u, const V &z, Real &tol) override
Given , solve for .
virtual void value(V &c, const V &u, const V &z, Real &tol) override
Evaluate the constraint operator at .
virtual void applyAdjointHessian_11(V &ahwv, const V &w, const V &v, const V &u, const V &z, Real &tol) override
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
typename PV< Real >::size_type size_type
virtual void applyAdjointJacobian_2(V &ajv, const V &v, const V &u, const V &z, const V &dualv, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the secondary interfa...
virtual void applyAdjointJacobian_1(V &ajv, const V &v, const V &u, const V &z, const V &dualv, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the secondary int...
Defines the constraint operator interface for simulation-based optimization.
Defines the time dependent constraint operator interface for simulation-based optimization.
Defines the linear algebra of vector space on a generic partitioned vector.
std::vector< PV >::size_type size_type
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167