ROL
ROL_ReducedDynamicObjective.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Rapid Optimization Library (ROL) Package
4//
5// Copyright 2014 NTESS and the ROL contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef ROL_REDUCEDDYNAMICOBJECTIVE_HPP
11#define ROL_REDUCEDDYNAMICOBJECTIVE_HPP
12
13#include "ROL_Ptr.hpp"
14#include "ROL_Sketch.hpp"
15#include "ROL_Objective.hpp"
16#include "ROL_DynamicObjective.hpp"
17#include "ROL_DynamicConstraint.hpp"
18
19#include <fstream>
20
43namespace ROL {
44
45template<typename Real>
46class ReducedDynamicObjective : public Objective<Real> {
47 using size_type = typename std::vector<Real>::size_type;
48private:
49 // Problem data.
50 const Ptr<DynamicObjective<Real>> obj_;
51 const Ptr<DynamicConstraint<Real>> con_;
52 const Ptr<Vector<Real>> u0_;
53 const std::vector<TimeStamp<Real>> timeStamp_;
54 const size_type Nt_;
55 // General sketch information.
56 const bool useSketch_;
57 // State sketch information.
58 size_type rankState_;
59 Ptr<Sketch<Real>> stateSketch_;
60 Ptr<Sketch<Real>> stateSketchCache_;
61 // Adjoint sketch information.
62 size_type rankAdjoint_;
63 Ptr<Sketch<Real>> adjointSketch_;
64 // State sensitivity sketch information.
65 size_type rankStateSens_;
66 Ptr<Sketch<Real>> stateSensSketch_;
67 // Inexactness information.
68 const bool useInexact_;
69 const size_type updateFactor_;
70 const size_type maxRank_;
71 const bool syncHessRank_;
72 const bool sumRankUpdate_;
73 const bool useDefaultRankUpdate_;
74 const Real a_, b_;
75 const Real maxTol_;
76 // Vector storage for intermediate computations.
77 std::vector<Ptr<Vector<Real>>> uhist_; // State history.
78 std::vector<Ptr<Vector<Real>>> lhist_; // Adjoint history.
79 std::vector<Ptr<Vector<Real>>> whist_; // State sensitivity history.
80 std::vector<Ptr<Vector<Real>>> phist_; // Adjoint sensitivity history.
81 Ptr<Vector<Real>> cprimal_; // Primal constraint vector.
82 Ptr<Vector<Real>> crhs_; // Primal constraint vector.
83 Ptr<Vector<Real>> udual_; // Dual state vector.
84 Ptr<Vector<Real>> rhs_; // Dual state vector.
85 Ptr<Vector<Real>> zdual_; // Dual control vector.
86 Real val_, valCache_; // Objective function value.
87 // Flags.
88 bool isValueComputed_; // Whether objective has been evaluated.
89 bool isStateComputed_; // Whether state has been solved.
90 bool isAdjointComputed_; // Whether adjoint has been solved.
91 bool isValueCached_; // Whether objective has been evaluated.
92 bool isStateCached_; // Whether state has been solved.
93 bool useHessian_; // Whether to use Hessian or FD.
94 bool useSymHess_; // Whether to use symmetric Hessian approximation.
95 // Output information
96 Ptr<std::ostream> stream_;
97 const bool print_;
98 const int freq_;
99
100 PartitionedVector<Real>& partition ( Vector<Real>& x ) const {
101 return static_cast<PartitionedVector<Real>&>(x);
102 }
103
104 const PartitionedVector<Real>& partition ( const Vector<Real>& x ) const {
105 return static_cast<const PartitionedVector<Real>&>(x);
106 }
107
108 void throwError(const int err, const std::string &sfunc,
109 const std::string &func, const int line) const {
110 if (err != 0) {
111 std::stringstream out;
112 out << ">>> ROL::ReducedDynamicObjective::" << func << " (Line " << line << "): ";
113 if (err == 1)
114 out << "Reconstruction has already been called!";
115 else if (err == 2)
116 out << "Input column index exceeds total number of columns!";
117 else if (err == 3)
118 out << "Reconstruction failed to compute domain-space basis!";
119 else if (err == 4)
120 out << "Reconstruction failed to compute range-space basis!";
121 else if (err == 5)
122 out << "Reconstruction failed to generate core matrix!";
123 throw Exception::NotImplemented(out.str());
124 }
125 }
126
127public:
128 ReducedDynamicObjective(const Ptr<DynamicObjective<Real>> &obj,
129 const Ptr<DynamicConstraint<Real>> &con,
130 const Ptr<Vector<Real>> &u0,
131 const Ptr<Vector<Real>> &zvec,
132 const Ptr<Vector<Real>> &cvec,
133 const std::vector<TimeStamp<Real>> &timeStamp,
134 ROL::ParameterList &pl,
135 const Ptr<std::ostream> &stream = nullPtr)
136 : obj_ ( obj ), // Dynamic objective function.
137 con_ ( con ), // Dynamic constraint function.
138 u0_ ( u0 ), // Initial condition.
139 timeStamp_ ( timeStamp ), // Vector of time stamps.
140 Nt_ ( timeStamp.size() ), // Number of time intervals.
141 useSketch_ ( pl.get("Use Sketching", false) ), // Use state sketch if true.
142 rankState_ ( pl.get("State Rank", 10) ), // Rank of state sketch.
143 stateSketch_ ( nullPtr ), // State sketch object.
144 rankAdjoint_ ( pl.get("Adjoint Rank", 10) ), // Rank of adjoint sketch.
145 adjointSketch_ ( nullPtr ), // Adjoint sketch object.
146 rankStateSens_ ( pl.get("State Sensitivity Rank", 10) ), // Rank of state sensitivity sketch.
147 stateSensSketch_ ( nullPtr ), // State sensitivity sketch object.
148 useInexact_ ( pl.get("Adaptive Rank", false) ), // Update rank adaptively.
149 updateFactor_ ( pl.get("Rank Update Factor", 2) ), // Rank update factor.
150 maxRank_ ( pl.get("Maximum Rank", 100) ), // Maximum rank.
151 syncHessRank_ ( pl.get("Sync Hessian Rank", useInexact_) ), // Sync rank for Hessian storage.
152 sumRankUpdate_ ( pl.get("Additive Rank Update", true) ), // Use additive rank update, otherwise use multiplicative
153 useDefaultRankUpdate_ ( pl.get("Use Basic Rank Update", true) ), // Use basic additive/multiplicative rank update
154 a_ ( pl.get("Log Rank Update Slope", 1.0) ), // Slope of semilogy tail energy
155 b_ ( pl.get("Log Rank Update Shift", 1.0) ), // Shift of semilogy tail energy
156 maxTol_ ( pl.get("Maximum Tolerance", ROL_INF<Real>()) ), // Maximum rank update tolerance
157 isValueComputed_ ( false ), // Flag indicating whether value has been computed.
158 isStateComputed_ ( false ), // Flag indicating whether state has been computed.
159 isAdjointComputed_ ( false ), // Flag indicating whether adjoint has been computed.
160 isValueCached_ ( false ), // Flag indicating whether value has been computed.
161 isStateCached_ ( false ), // Flag indicating whether state has been computed.
162 useHessian_ ( pl.get("Use Hessian", true) ), // Flag indicating whether to use the Hessian.
163 useSymHess_ ( pl.get("Use Only Sketched Sensitivity", true) ), // Flag indicating whether to use symmetric sketched Hessian.
164 stream_ ( stream ), // Output stream to print sketch information.
165 print_ ( pl.get("Print Optimization Vector", false) ), // Print control vector to file
166 freq_ ( pl.get("Output Frequency", 5) ) { // Print frequency for control vector
167 uhist_.clear(); lhist_.clear(); whist_.clear(); phist_.clear();
168 if (useSketch_) { // Only maintain a sketch of the state time history
169 Real orthTol = pl.get("Orthogonality Tolerance", 1e2*ROL_EPSILON<Real>());
170 int orthIt = pl.get("Reorthogonalization Iterations", 5);
171 bool trunc = pl.get("Truncate Approximation", false);
172 if (syncHessRank_) {
173 rankAdjoint_ = rankState_;
174 rankStateSens_ = rankState_;
175 }
176 unsigned dseed = pl.get("State Domain Seed",0);
177 unsigned rseed = pl.get("State Range Seed",0);
178 stateSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,
179 rankState_,orthTol,orthIt,trunc,dseed,rseed);
180 stateSketchCache_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,
181 rankState_,orthTol,orthIt,trunc,dseed,rseed);
182 uhist_.push_back(u0_->clone());
183 uhist_.push_back(u0_->clone());
184 lhist_.push_back(cvec->dual().clone());
185 dseed = pl.get("Adjoint Domain Seed",0);
186 rseed = pl.get("Adjoint Range Seed",0);
187 adjointSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,rankAdjoint_,
188 orthTol,orthIt,trunc,dseed,rseed);
189 if (useHessian_) {
190 dseed = pl.get("State Sensitivity Domain Seed",0);
191 rseed = pl.get("State Sensitivity Range Seed",0);
192 stateSensSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,
193 rankStateSens_,orthTol,orthIt,trunc,dseed,rseed);
194 whist_.push_back(u0_->clone());
195 whist_.push_back(u0_->clone());
196 phist_.push_back(cvec->dual().clone());
197 }
198 }
199 else { // Store entire state time history
200 for (size_type i = 0; i < Nt_; ++i) {
201 uhist_.push_back(u0_->clone());
202 lhist_.push_back(cvec->dual().clone());
203 if (useHessian_) {
204 whist_.push_back(u0_->clone());
205 phist_.push_back(cvec->dual().clone());
206 }
207 }
208 }
209 cprimal_ = cvec->clone();
210 crhs_ = cvec->clone();
211 udual_ = u0_->dual().clone();
212 rhs_ = u0_->dual().clone();
213 zdual_ = zvec->dual().clone();
214 }
215
216 Ptr<Vector<Real>> makeDynamicVector(const Vector<Real> &x) const {
217 return ROL::PartitionedVector<Real>::create(x, Nt_);
218 }
219
220 size_type getStateRank() const { return rankState_; }
221 size_type getAdjointRank() const { return rankAdjoint_; }
222 size_type getStateSensitivityRank() const { return rankStateSens_; }
223
224 void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
225 if (useSketch_) {
226 stateSketch_->reset(true);
227 if (flag == true) {
228 adjointSketch_->reset(true);
229 if (useHessian_) stateSensSketch_->reset(true);
230 }
231 }
232 for (size_type i = 0; i < uhist_.size(); ++i) uhist_[i]->zero();
233 val_ = static_cast<Real>(0);
234 isValueComputed_ = false;
235 isStateComputed_ = false;
236 if (flag == true) isAdjointComputed_ = false;
237
238 if (iter >= 0 && iter%freq_==0 && print_) {
239 std::stringstream name;
240 name << "optvector." << iter << ".txt";
241 std::ofstream file;
242 file.open(name.str());
243 x.print(file);
244 file.close();
245 }
246 }
247
248 void update(const Vector<Real> &x, UpdateType type, int iter = -1) {
249 // This just resets all storage independent of update type.
250 // If sketching is not in use, we could be smarter with the
251 // storage. When sketching is in use, you may want to reset
252 // everything and re-randomize the sketch object. It is
253 // unclear to me (Drew) if there is a benefit to saving the
254 // previous computations in this case.
255 if (useSketch_) {
256 switch(type) {
257 case UpdateType::Initial:
258 {
259 for (size_type i = 0; i < uhist_.size(); ++i) uhist_[i]->zero();
260 val_ = static_cast<Real>(0);
261 isValueComputed_ = false;
262 isStateComputed_ = false;
263 isAdjointComputed_ = false;
264 stateSketch_->reset(true);
265 stateSketchCache_->reset(true);
266 adjointSketch_->reset(true);
267 if (useHessian_) stateSensSketch_->reset(true);
268 break;
269 }
270 case UpdateType::Trial:
271 {
272 for (size_type i = 0; i < uhist_.size(); ++i) uhist_[i]->zero();
273 valCache_ = val_;
274 isValueCached_ = isValueComputed_;
275 isStateCached_ = isStateComputed_;
276 val_ = static_cast<Real>(0);
277 isValueComputed_ = false;
278 isStateComputed_ = false;
279 auto tmp = stateSketch_;
280 stateSketch_ = stateSketchCache_;
281 stateSketchCache_ = tmp;
282 stateSketch_->reset(true);
283 break;
284 }
285 case UpdateType::Accept:
286 {
287 isAdjointComputed_ = false;
288 adjointSketch_->reset(true);
289 if (useHessian_) stateSensSketch_->reset(true);
290 break;
291 }
292 case UpdateType::Revert:
293 {
294 for (size_type i = 0; i < uhist_.size(); ++i) uhist_[i]->zero();
295 val_ = valCache_;
296 isValueComputed_ = isValueCached_;
297 isStateComputed_ = isStateCached_;
298 auto tmp = stateSketchCache_;
299 stateSketchCache_ = stateSketch_;
300 stateSketch_ = tmp;
301 if (useHessian_) stateSensSketch_->reset(true);
302 break;
303 }
304 case UpdateType::Temp:
305 {
306 for (size_type i = 0; i < uhist_.size(); ++i) uhist_[i]->zero();
307 valCache_ = val_;
308 isValueCached_ = isValueComputed_;
309 isStateCached_ = isStateComputed_;
310 val_ = static_cast<Real>(0);
311 isValueComputed_ = false;
312 isStateComputed_ = false;
313 isAdjointComputed_ = false;
314 auto tmp = stateSketch_;
315 stateSketch_ = stateSketchCache_;
316 stateSketchCache_ = tmp;
317 stateSketch_->reset(true);
318 adjointSketch_->reset(true);
319 if (useHessian_) stateSensSketch_->reset(true);
320 break;
321 }
322 }
323 }
324 else {
325 switch(type) {
326 case UpdateType::Initial:
327 case UpdateType::Accept:
328 case UpdateType::Revert:
329 case UpdateType::Trial:
330 case UpdateType::Temp:
331 {
332 for (size_type i = 0; i < uhist_.size(); ++i) uhist_[i]->zero();
333 val_ = static_cast<Real>(0);
334 isValueComputed_ = false;
335 isStateComputed_ = false;
336 isAdjointComputed_ = false;
337 break;
338 }
339 }
340 }
341
342 if (iter >= 0 && iter%freq_==0 && print_) {
343 std::stringstream name;
344 name << "optvector." << iter << ".txt";
345 std::ofstream file;
346 file.open(name.str());
347 x.print(file);
348 file.close();
349 }
350 }
351
352 Real value( const Vector<Real> &x, Real &tol ) {
353 if (!isValueComputed_) {
354 int eflag(0);
355 const Real one(1);
356 const PartitionedVector<Real> &xp = partition(x);
357 // Set initial condition
358 uhist_[0]->set(*u0_);
359 if (useSketch_) {
360 stateSketch_->update();
361 uhist_[1]->set(*u0_);
362 }
363 // Run time stepper
364 Real valk(0);
365 size_type index;
366 for (size_type k = 1; k < Nt_; ++k) {
367 index = (useSketch_ ? 1 : k);
368 if (!useSketch_) uhist_[index]->set(*uhist_[index-1]);
369 // Update dynamic constraint
370 con_->update_uo(*uhist_[index-1], timeStamp_[k]);
371 con_->update_z(*xp.get(k), timeStamp_[k]);
372 // Solve state on current time interval
373 con_->solve(*cprimal_, *uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
374 // Update dynamic objective
375 obj_->update(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
376 // Compute objective function value on current time interval
377 valk = obj_->value(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
378 // Update total objective function value
379 val_ += valk;
380 // Sketch state
381 if (useSketch_) {
382 eflag = stateSketch_->advance(one,*uhist_[1],static_cast<int>(k)-1,one);
383 throwError(eflag,"advance","value",315);
384 uhist_[0]->set(*uhist_[1]);
385 }
386 }
387 isValueComputed_ = true;
388 isStateComputed_ = true;
389 }
390 return val_;
391 }
392
393 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
394 int eflag(0);
395 PartitionedVector<Real> &gp = partition(g);
396 gp.get(0)->zero(); // zero for the nonexistant zeroth control interval
397 const PartitionedVector<Real> &xp = partition(x);
398 const Real one(1);
399 size_type uindex = (useSketch_ ? 1 : Nt_-1);
400 size_type lindex = (useSketch_ ? 0 : Nt_-1);
401 // Must first compute the value
402 solveState(x);
403 if (useSketch_ && useInexact_) {
404 if (stream_ != nullPtr) {
405 *stream_ << std::string(80,'=') << std::endl;
406 *stream_ << " ROL::ReducedDynamicObjective::gradient" << std::endl;
407 }
408 tol = updateSketch(x,tol);
409 if (stream_ != nullPtr) {
410 *stream_ << " State Rank for Gradient Computation: " << rankState_ << std::endl;
411 *stream_ << " Residual Norm: " << tol << std::endl;
412 *stream_ << std::string(80,'=') << std::endl;
413 }
414 }
415 // Recover state from sketch
416 if (useSketch_) {
417 uhist_[1]->set(*uhist_[0]);
418 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
419 throwError(eflag,"reconstruct","gradient",351);
420 if (isAdjointComputed_) {
421 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(Nt_)-2);
422 throwError(eflag,"reconstruct","gradient",354);
423 }
424 }
425 // Update dynamic constraint and objective
426 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
427 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
428 // Solve for terminal condition
429 if (!isAdjointComputed_) {
430 setTerminalCondition(*lhist_[lindex],
431 *uhist_[uindex-1], *uhist_[uindex],
432 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
433 if (useSketch_) {
434 eflag = adjointSketch_->advance(one,*lhist_[0],static_cast<int>(Nt_)-2,one);
435 throwError(eflag,"advance","gradient",367);
436 }
437 }
438 // Update gradient on terminal interval
439 updateGradient(*gp.get(Nt_-1), *lhist_[lindex],
440 *uhist_[uindex-1], *uhist_[uindex],
441 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
442 // Run reverse time stepper
443 for (size_type k = Nt_-2; k > 0; --k) {
444 if (!isAdjointComputed_) {
445 // Compute k+1 component of rhs
446 computeAdjointRHS(*rhs_, *lhist_[lindex],
447 *uhist_[uindex-1], *uhist_[uindex],
448 *xp.get(k+1), timeStamp_[k+1]);
449 }
450 uindex = (useSketch_ ? 1 : k);
451 lindex = (useSketch_ ? 0 : k);
452 // Recover state from sketch
453 if (useSketch_) {
454 uhist_[1]->set(*uhist_[0]);
455 if (k==1) {
456 uhist_[0]->set(*u0_);
457 }
458 else {
459 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
460 throwError(eflag,"reconstruct","gradient",392);
461 }
462 if (isAdjointComputed_) {
463 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
464 throwError(eflag,"reconstruct","gradient",396);
465 }
466 }
467 // Update dynamic constraint and objective
468 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
469 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
470 // Solve for adjoint on interval k
471 if (!isAdjointComputed_) {
472 advanceAdjoint(*lhist_[lindex], *rhs_,
473 *uhist_[uindex-1], *uhist_[uindex],
474 *xp.get(k), timeStamp_[k]);
475 if (useSketch_) {
476 eflag = adjointSketch_->advance(one,*lhist_[0],static_cast<int>(k)-1,one);
477 throwError(eflag,"advance","gradient",367);
478 }
479 }
480 // Update gradient on interval k
481 updateGradient(*gp.get(k), *lhist_[lindex],
482 *uhist_[uindex-1], *uhist_[uindex],
483 *xp.get(k), timeStamp_[k]);
484 }
485 isAdjointComputed_ = true;
486 }
487
488 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
489 int eflag(0);
490 if (useHessian_) {
491 // Must first solve the state and adjoint equations
492 solveState(x);
493 solveAdjoint(x);
494 // Now compute Hessian
495 const Real one(1);
496 const PartitionedVector<Real> &xp = partition(x);
497 const PartitionedVector<Real> &vp = partition(v);
498 PartitionedVector<Real> &hvp = partition(hv);
499 hvp.get(0)->zero(); // zero for the nonexistant zeroth control interval
500 // Compute state sensitivity
501 whist_[0]->zero();
502 if (useSketch_) {
503 stateSensSketch_->reset(false);
504 //stateSensSketch_->advance(one,*whist_[0],0,one);
505 uhist_[0]->set(*u0_);
506 }
507 size_type uindex, lindex;
508 for (size_type k = 1; k < Nt_; ++k) {
509 uindex = (useSketch_ ? 1 : k);
510 lindex = (useSketch_ ? 0 : k);
511 // Reconstruct sketched state
512 if (useSketch_) {
513 eflag = stateSketch_->reconstruct(*uhist_[1],static_cast<int>(k)-1);
514 throwError(eflag,"reconstruct","hessVec",446);
515 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
516 throwError(eflag,"reconstruct","hessVec",448);
517 }
518 // Update dynamic constraint and objective
519 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
520 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
521 // Advance state sensitivity on current time interval
522 advanceStateSens(*whist_[uindex],
523 *vp.get(k), *whist_[uindex-1],
524 *uhist_[uindex-1], *uhist_[uindex],
525 *xp.get(k), timeStamp_[k]);
526 // Set control only Hessian
527 computeControlHessLag(*hvp.get(k),
528 *vp.get(k), *lhist_[lindex],
529 *uhist_[uindex-1], *uhist_[uindex],
530 *xp.get(k), timeStamp_[k]);
531 if (!useSymHess_) {
532 // Add mixed derivative Hessian
533 addMixedHessLag(*hvp.get(k), *lhist_[lindex],
534 *whist_[uindex-1], *whist_[uindex],
535 *uhist_[uindex-1], *uhist_[uindex],
536 *xp.get(k), timeStamp_[k]);
537 }
538 // Sketch state sensitivity
539 if (useSketch_) {
540 eflag = stateSensSketch_->advance(one,*whist_[1],static_cast<int>(k)-1,one);
541 throwError(eflag,"advance","hessVec",473);
542 whist_[0]->set(*whist_[1]);
543 uhist_[0]->set(*uhist_[1]);
544 }
545 }
546
547 // Compute adjoint sensitivity
548 uindex = (useSketch_ ? 1 : Nt_-1);
549 lindex = (useSketch_ ? 0 : Nt_-1);
550 if (useSketch_) {
551 // Recover terminal state
552 uhist_[1]->set(*uhist_[0]);
553 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
554 throwError(eflag,"reconstruct","hessVec",486);
555 // Recover terminal adjoint
556 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(Nt_)-2);
557 throwError(eflag,"reconstruct","hessVec",489);
558 // Recover terminal state sensitivity
559 whist_[1]->set(*whist_[0]);
560 eflag = stateSensSketch_->reconstruct(*whist_[0],static_cast<int>(Nt_)-3);
561 throwError(eflag,"reconstruct","hessVec",493);
562 }
563 // Update dynamic constraint and objective
564 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
565 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
566 // Solve for terminal condition
567 setTerminalConditionHess(*phist_[lindex],
568 *vp.get(Nt_-1), *lhist_[lindex],
569 *whist_[uindex-1], *whist_[uindex],
570 *uhist_[uindex-1], *uhist_[uindex],
571 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
572 if (useSymHess_) {
573 // Add mixed derivative Hessian
574 addMixedHessLag(*hvp.get(Nt_-1), *lhist_[lindex],
575 *whist_[uindex-1], *whist_[uindex],
576 *uhist_[uindex-1], *uhist_[uindex],
577 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
578 }
579 // Add adjoint sensitivity to Hessian
580 addAdjointSens(*hvp.get(Nt_-1), *phist_[lindex],
581 *uhist_[uindex-1], *uhist_[uindex],
582 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
583 // Run reverse time stepper
584 for (size_type k = Nt_-2; k > 0; --k) {
585 // Compute new components of rhs
586 computeNewStateHessLag(*rhs_, *lhist_[lindex],
587 *whist_[uindex-1], *whist_[uindex],
588 *uhist_[uindex-1], *uhist_[uindex],
589 *xp.get(k+1), timeStamp_[k+1],
590 false);
591 computeNewMixedHessLag(*rhs_,
592 *vp.get(k+1), *lhist_[lindex],
593 *uhist_[uindex-1], *uhist_[uindex],
594 *xp.get(k+1), timeStamp_[k+1],
595 true);
596 computeNewStateJacobian(*rhs_, *phist_[lindex],
597 *uhist_[uindex-1], *uhist_[uindex],
598 *xp.get(k+1), timeStamp_[k+1],
599 true);
600 // Recover state, adjoint and state sensitivity from sketch
601 if (useSketch_) {
602 uhist_[1]->set(*uhist_[0]);
603 whist_[1]->set(*whist_[0]);
604 if (k==1) {
605 uhist_[0]->set(*u0_);
606 whist_[0]->zero();
607 }
608 else {
609 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
610 throwError(eflag,"reconstruct","hessVec",542);
611 eflag = stateSensSketch_->reconstruct(*whist_[0],static_cast<int>(k)-2);
612 throwError(eflag,"reconstruct","hessVec",544);
613 }
614 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
615 throwError(eflag,"reconstruct","hessVec",547);
616 }
617 uindex = (useSketch_ ? 1 : k);
618 lindex = (useSketch_ ? 0 : k);
619 // Update dynamic constraint and objective
620 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
621 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
622 // Compute old components of rhs
623 computeOldStateHessLag(*rhs_, *lhist_[lindex],
624 *whist_[uindex-1], *whist_[uindex],
625 *uhist_[uindex-1], *uhist_[uindex],
626 *xp.get(k), timeStamp_[k],
627 true);
628 computeOldMixedHessLag(*rhs_,
629 *vp.get(k), *lhist_[lindex],
630 *uhist_[uindex-1], *uhist_[uindex],
631 *xp.get(k), timeStamp_[k],
632 true);
633 // Solve for adjoint on interval k
634 advanceAdjointSens(*phist_[lindex], *rhs_,
635 *uhist_[uindex-1], *uhist_[uindex],
636 *xp.get(k), timeStamp_[k]);
637 if (useSymHess_) {
638 // Add mixed derivative Hessian
639 addMixedHessLag(*hvp.get(k), *lhist_[lindex],
640 *whist_[uindex-1], *whist_[uindex],
641 *uhist_[uindex-1], *uhist_[uindex],
642 *xp.get(k), timeStamp_[k]);
643 }
644 // Add adjoint sensitivity to Hessian
645 addAdjointSens(*hvp.get(k), *phist_[lindex],
646 *uhist_[uindex-1], *uhist_[uindex],
647 *xp.get(k), timeStamp_[k]);
648 }
649 }
650 else {
651 Objective<Real>::hessVec(hv,v,x,tol);
652 }
653 }
654
655private:
656 /***************************************************************************/
657 /************ Method to solve state equation *******************************/
658 /***************************************************************************/
659 Real solveState(const Vector<Real> &x) {
660 Real cnorm(0);
661 if (!isStateComputed_) {
662 int eflag(0);
663 const Real one(1);
664 const PartitionedVector<Real> &xp = partition(x);
665 // Set initial condition
666 uhist_[0]->set(*u0_);
667 if (useSketch_) { // Set initial guess for solve
668 stateSketch_->update();
669 uhist_[1]->set(*u0_);
670 }
671 // Run time stepper
672 size_type index;
673 for (size_type k = 1; k < Nt_; ++k) {
674 index = (useSketch_ ? 1 : k);
675 // Set initial guess for solve
676 if (!useSketch_) uhist_[index]->set(*uhist_[index-1]);
677 // Solve state on current time interval
678 con_->update_uo(*uhist_[index-1], timeStamp_[k]);
679 con_->update_z(*xp.get(k), timeStamp_[k]);
680 con_->solve(*cprimal_, *uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
681 cnorm = std::max(cnorm,cprimal_->norm());
682 // Sketch state
683 if (useSketch_) {
684 eflag = stateSketch_->advance(one, *uhist_[1], static_cast<int>(k)-1, one);
685 throwError(eflag,"advance","solveState",574);
686 uhist_[0]->set(*uhist_[1]);
687 }
688 }
689 isStateComputed_ = true;
690 }
691 return cnorm;
692 }
693
694 Real updateSketch(const Vector<Real> &x, const Real tol) {
695 int eflag(0);
696 const PartitionedVector<Real> &xp = partition(x);
697 Real err(0), err2(0), serr(0), cdot(0), dt(0); //, cnorm(0)
698 Real tol0 = std::min(maxTol_,tol);
699 bool flag = true;
700 while (flag) {
701 err = static_cast<Real>(0);
702 err2 = err;
703 uhist_[0]->set(*u0_);
704 for (size_type k = 1; k < Nt_; ++k) {
705 eflag = stateSketch_->reconstruct(*uhist_[1],static_cast<int>(k)-1);
706 throwError(eflag,"reconstruct","updateSketch",592);
707 con_->update(*uhist_[0], *uhist_[1], *xp.get(k), timeStamp_[k]);
708 con_->value(*cprimal_, *uhist_[0], *uhist_[1], *xp.get(k), timeStamp_[k]);
709 /**** Linf norm for residual error ****/
710 //cnorm = cprimal_->norm();
711 //err = (cnorm > err ? cnorm : err);
712 /**** L2 norm for residual error ****/
713 cdot = cprimal_->dot(*cprimal_);
714 dt = timeStamp_[k].t[1]-timeStamp_[k].t[0];
715 err2 += cdot / dt;
716 err = std::sqrt(err2);
717 if (err > tol0) break;
718 uhist_[0]->set(*uhist_[1]);
719 }
720 if (stream_ != nullPtr) {
721 *stream_ << " *** State Rank: " << rankState_ << std::endl;
722 *stream_ << " *** Required Tolerance: " << tol0 << std::endl;
723 *stream_ << " *** Residual Norm: " << err << std::endl;
724 }
725 if (err > tol0) {
726 rankState_ = (sumRankUpdate_ ? rankState_ + updateFactor_ : rankState_ * updateFactor_);
727 if (!useDefaultRankUpdate_)
728 rankState_ = std::max(rankState_,static_cast<size_type>(std::ceil((b_-std::log(tol0))/a_)));
729 //Real a(0.1838), b(3.1451); // Navier-Stokes
730 //Real a(2.6125), b(2.4841); // Semilinear
731 //rankState_ = std::max(rankState_+2,static_cast<size_t>(std::ceil((b-std::log(tol0))/a)));
732 //rankState_ *= updateFactor_; // Perhaps there is a better update strategy
733 rankState_ = (maxRank_ < rankState_ ? maxRank_ : rankState_);
734 stateSketch_->setRank(rankState_);
735 if (syncHessRank_) {
736 rankAdjoint_ = rankState_;
737 rankStateSens_ = rankState_;
738 adjointSketch_->setRank(rankAdjoint_);
739 stateSensSketch_->setRank(rankStateSens_);
740 }
741 isStateComputed_ = false;
742 isAdjointComputed_ = false;
743 serr = solveState(x);
744 if (stream_ != nullPtr)
745 *stream_ << " *** Maximum Solver Error: " << serr << std::endl;
746 }
747 else {
748 flag = false;
749 break;
750 }
751 }
752 return err;
753 }
754
755 /***************************************************************************/
756 /************ Methods to solve adjoint equation ****************************/
757 /***************************************************************************/
758 void setTerminalCondition(Vector<Real> &l,
759 const Vector<Real> &uold, const Vector<Real> &unew,
760 const Vector<Real> &z, const TimeStamp<Real> &ts) {
761 obj_->gradient_un(*udual_, uold, unew, z, ts);
762 con_->applyInverseAdjointJacobian_un(l, *udual_, uold, unew, z, ts);
763 }
764
765 void computeAdjointRHS(Vector<Real> &rhs, const Vector<Real> &l,
766 const Vector<Real> &uold, const Vector<Real> &unew,
767 const Vector<Real> &z, const TimeStamp<Real> &ts) {
768 const Real one(1);
769 obj_->gradient_uo(rhs, uold, unew, z, ts);
770 con_->applyAdjointJacobian_uo(*udual_, l, uold, unew, z, ts);
771 rhs.axpy(-one,*udual_);
772 }
773
774 void advanceAdjoint(Vector<Real> &l, Vector<Real> &rhs,
775 const Vector<Real> &uold, const Vector<Real> &unew,
776 const Vector<Real> &z, const TimeStamp<Real> &ts) {
777 obj_->gradient_un(*udual_, uold, unew, z, ts);
778 rhs.plus(*udual_);
779 con_->applyInverseAdjointJacobian_un(l, rhs, uold, unew, z, ts);
780 }
781
782 void solveAdjoint(const Vector<Real> &x) {
783 int eflag(0);
784 if (!isAdjointComputed_) {
785 const PartitionedVector<Real> &xp = partition(x);
786 const Real one(1);
787 size_type uindex = (useSketch_ ? 1 : Nt_-1);
788 size_type lindex = (useSketch_ ? 0 : Nt_-1);
789 // Must first compute solve the state equation
790 solveState(x);
791 // Recover state from sketch
792 if (useSketch_) {
793 uhist_[1]->set(*uhist_[0]);
794 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
795 throwError(eflag,"reconstruct","solveAdjoint",672);
796 }
797 // Update dynamic constraint and objective
798 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
799 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
800 // Solve for terminal condition
801 setTerminalCondition(*lhist_[lindex],
802 *uhist_[uindex-1], *uhist_[uindex],
803 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
804 if (useSketch_) {
805 eflag = adjointSketch_->advance(one,*lhist_[lindex],static_cast<int>(Nt_)-2,one);
806 throwError(eflag,"advance","solveAdjoint",683);
807 }
808 // Run reverse time stepper
809 for (size_type k = Nt_-2; k > 0; --k) {
810 // Compute k+1 component of rhs
811 computeAdjointRHS(*rhs_, *lhist_[lindex],
812 *uhist_[uindex-1], *uhist_[uindex],
813 *xp.get(k+1), timeStamp_[k+1]);
814 // Recover state from sketch
815 if (useSketch_) {
816 uhist_[1]->set(*uhist_[0]);
817 if (k==1) {
818 uhist_[0]->set(*u0_);
819 }
820 else {
821 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
822 throwError(eflag,"reconstruct","solveAdjoint",699);
823 }
824 }
825 uindex = (useSketch_ ? 1 : k);
826 lindex = (useSketch_ ? 0 : k);
827 // Update dynamic constraint and objective
828 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
829 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
830 // Solve for adjoint on interval k
831 advanceAdjoint(*lhist_[lindex], *rhs_,
832 *uhist_[uindex-1], *uhist_[uindex],
833 *xp.get(k), timeStamp_[k]);
834 if (useSketch_) {
835 eflag = adjointSketch_->advance(one,*lhist_[lindex],static_cast<int>(k)-1,one);
836 throwError(eflag,"advance","solveAdjoint",713);
837 }
838 }
839 isAdjointComputed_ = true;
840 }
841 }
842
843 /***************************************************************************/
844 /************ Method for gradient computation ******************************/
845 /***************************************************************************/
846 void updateGradient(Vector<Real> &g, const Vector<Real> &l,
847 const Vector<Real> &uold, const Vector<Real> &unew,
848 const Vector<Real> &z, const TimeStamp<Real> &ts) {
849 const Real one(1);
850 obj_->gradient_z(g, uold, unew, z, ts);
851 con_->applyAdjointJacobian_z(*zdual_, l, uold, unew, z, ts);
852 g.axpy(-one,*zdual_);
853 }
854
855 /***************************************************************************/
856 /************ Method to solve state sensitivity equation *******************/
857 /***************************************************************************/
858 void advanceStateSens(Vector<Real> &wnew, const Vector<Real> &v,
859 const Vector<Real> &wold, const Vector<Real> &uold,
860 const Vector<Real> &unew, const Vector<Real> &z,
861 const TimeStamp<Real> &ts) {
862 const Real one(1);
863 con_->applyJacobian_z(*crhs_, v, uold, unew, z, ts);
864 con_->applyJacobian_uo(*cprimal_, wold, uold, unew, z, ts);
865 crhs_->axpy(-one, *cprimal_);
866 con_->applyInverseJacobian_un(wnew, *crhs_, uold, unew, z, ts);
867 }
868
869 /***************************************************************************/
870 /************ Methods to solve adjoint sensitivity equation ****************/
871 /***************************************************************************/
872 void setTerminalConditionHess(Vector<Real> &p,
873 const Vector<Real> &v, const Vector<Real> &l,
874 const Vector<Real> &wold, const Vector<Real> &wnew,
875 const Vector<Real> &uold, const Vector<Real> &unew,
876 const Vector<Real> &z, const TimeStamp<Real> &ts) {
877 const Real one(1);
878 // Mixed derivative rhs term
879 con_->applyAdjointHessian_z_un(*rhs_, l, v, uold, unew, z, ts);
880 obj_->hessVec_un_z(*udual_, v, uold, unew, z, ts);
881 rhs_->axpy(-one, *udual_);
882 // State derivative rhs term
883 con_->applyAdjointHessian_un_un(*udual_, l, wnew, uold, unew, z, ts);
884 rhs_->axpy(-one, *udual_);
885 obj_->hessVec_un_un(*udual_, wnew, uold, unew, z, ts);
886 rhs_->plus(*udual_);
887 con_->applyAdjointHessian_uo_un(*udual_, l, wold, uold, unew, z, ts);
888 rhs_->axpy(-one, *udual_);
889 obj_->hessVec_un_uo(*udual_, wold, uold, unew, z, ts);
890 rhs_->plus(*udual_);
891 // Invert adjoint Jacobian
892 con_->applyInverseAdjointJacobian_un(p, *rhs_, uold, unew, z, ts);
893 }
894
895 void computeOldStateHessLag(Vector<Real> &Hv, const Vector<Real> &l,
896 const Vector<Real> &wold, const Vector<Real> &wnew,
897 const Vector<Real> &uold, const Vector<Real> &unew,
898 const Vector<Real> &z, const TimeStamp<Real> &ts,
899 const bool sumInto = false) {
900 const Real one(1);
901 // Compute new/old Hessian of Lagrangian
902 if (!sumInto) {
903 obj_->hessVec_un_uo(Hv, wold, uold, unew, z, ts);
904 }
905 else {
906 obj_->hessVec_un_uo(*udual_, wold, uold, unew, z, ts);
907 Hv.plus(*udual_);
908 }
909 con_->applyAdjointHessian_uo_un(*udual_, l, wold, uold, unew, z, ts);
910 Hv.axpy(-one,*udual_);
911 // Compute new/new Hessian of Lagrangian
912 obj_->hessVec_un_un(*udual_, wnew, uold, unew, z, ts);
913 Hv.plus(*udual_);
914 con_->applyAdjointHessian_un_un(*udual_, l, wnew, uold, unew, z, ts);
915 Hv.axpy(-one,*udual_);
916 }
917
918 void computeOldMixedHessLag(Vector<Real> &Hv,
919 const Vector<Real> &v, const Vector<Real> &l,
920 const Vector<Real> &uold, const Vector<Real> &unew,
921 const Vector<Real> &z, const TimeStamp<Real> &ts,
922 const bool sumInto = false) {
923 const Real one(1);
924 // Compute new/old Hessian of Lagrangian
925 if (!sumInto) {
926 con_->applyAdjointHessian_z_un(Hv, l, v, uold, unew, z, ts);
927 }
928 else {
929 con_->applyAdjointHessian_z_un(*udual_, l, v, uold, unew, z, ts);
930 Hv.plus(*udual_);
931 }
932 obj_->hessVec_un_z(*udual_, v, uold, unew, z, ts);
933 Hv.axpy(-one, *udual_);
934 }
935
936 void computeNewStateJacobian(Vector<Real> &Hv, const Vector<Real> &p,
937 const Vector<Real> &uold, const Vector<Real> &unew,
938 const Vector<Real> &z, const TimeStamp<Real> &ts,
939 const bool sumInto = false) {
940 const Real one(1);
941 if (!sumInto) {
942 con_->applyAdjointJacobian_uo(Hv, p, uold, unew, z, ts);
943 Hv.scale(-one);
944 }
945 else {
946 con_->applyAdjointJacobian_uo(*udual_, p, uold, unew, z, ts);
947 Hv.axpy(-one, *udual_);
948 }
949 }
950
951 void computeNewStateHessLag(Vector<Real> &Hv, const Vector<Real> &l,
952 const Vector<Real> &wold, const Vector<Real> &wnew,
953 const Vector<Real> &uold, const Vector<Real> &unew,
954 const Vector<Real> &z, const TimeStamp<Real> &ts,
955 const bool sumInto = false) {
956 const Real one(1);
957 // Compute old/new Hessian of Lagrangian
958 if (!sumInto) {
959 obj_->hessVec_uo_un(Hv, wnew, uold, unew, z, ts);
960 }
961 else {
962 obj_->hessVec_uo_un(*udual_, wnew, uold, unew, z, ts);
963 Hv.plus(*udual_);
964 }
965 con_->applyAdjointHessian_un_uo(*udual_, l, wnew, uold, unew, z, ts);
966 Hv.axpy(-one,*udual_);
967 // Compute old/old Hessian of Lagrangian
968 obj_->hessVec_uo_uo(*udual_, wold, uold, unew, z, ts);
969 Hv.plus(*udual_);
970 con_->applyAdjointHessian_uo_uo(*udual_, l, wold, uold, unew, z, ts);
971 Hv.axpy(-one,*udual_);
972 }
973
974 void computeNewMixedHessLag(Vector<Real> &Hv,
975 const Vector<Real> &v, const Vector<Real> &l,
976 const Vector<Real> &uold, const Vector<Real> &unew,
977 const Vector<Real> &z, const TimeStamp<Real> &ts,
978 const bool sumInto = false) {
979 const Real one(1);
980 // Compute new/old Hessian of Lagrangian
981 if (!sumInto) {
982 con_->applyAdjointHessian_z_uo(Hv, l, v, uold, unew, z, ts);
983 }
984 else {
985 con_->applyAdjointHessian_z_uo(*udual_, l, v, uold, unew, z, ts);
986 Hv.plus(*udual_);
987 }
988 obj_->hessVec_uo_z(*udual_, v, uold, unew, z, ts);
989 Hv.axpy(-one,*udual_);
990 }
991
992 void advanceAdjointSens(Vector<Real> &p, Vector<Real> &rhs,
993 const Vector<Real> &uold, const Vector<Real> &unew,
994 const Vector<Real> &z, const TimeStamp<Real> &ts) {
995 // Solve adjoint sensitivity on current time interval
996 con_->applyInverseAdjointJacobian_un(p, rhs, uold, unew, z, ts);
997 }
998
999 /***************************************************************************/
1000 /************ Method for Hessian-times-a-vector computation ****************/
1001 /***************************************************************************/
1002 void computeControlHessLag(Vector<Real> &Hv,
1003 const Vector<Real> &v, const Vector<Real> &l,
1004 const Vector<Real> &uold, const Vector<Real> &unew,
1005 const Vector<Real> &z, const TimeStamp<Real> &ts) {
1006 const Real one(1);
1007 // Compute Hessian of Lagrangian
1008 obj_->hessVec_z_z(Hv, v, uold, unew, z, ts);
1009 con_->applyAdjointHessian_z_z(*zdual_, l, v, uold, unew, z, ts);
1010 Hv.axpy(-one, *zdual_);
1011 }
1012
1013 void addMixedHessLag(Vector<Real> &Hv, const Vector<Real> &l,
1014 const Vector<Real> &wold, const Vector<Real> &wnew,
1015 const Vector<Real> &uold, const Vector<Real> &unew,
1016 const Vector<Real> &z, const TimeStamp<Real> &ts) {
1017 const Real one(1);
1018 // Compute Hessian of Lagrangian on previous time interval
1019 obj_->hessVec_z_uo(*zdual_, wold, uold, unew, z, ts);
1020 Hv.axpy(-one, *zdual_);
1021 con_->applyAdjointHessian_uo_z(*zdual_, l, wold, uold, unew, z, ts);
1022 Hv.plus(*zdual_);
1023 // Compute Hessian of Lagrangian on current time interval
1024 obj_->hessVec_z_un(*zdual_, wnew, uold, unew, z, ts);
1025 Hv.axpy(-one, *zdual_);
1026 con_->applyAdjointHessian_un_z(*zdual_, l, wnew, uold, unew, z, ts);
1027 Hv.plus(*zdual_);
1028 }
1029
1030 void addAdjointSens(Vector<Real> &Hv, const Vector<Real> &p,
1031 const Vector<Real> &uold, const Vector<Real> &unew,
1032 const Vector<Real> &z, const TimeStamp<Real> &ts) {
1033 con_->applyAdjointJacobian_z(*zdual_, p, uold, unew, z, ts);
1034 Hv.plus(*zdual_);
1035 }
1036};
1037
1038} // namespace ROL
1039
1040#endif // ROL_REDUCEDDYNAMICOBJECTIVE_HPP
ROL_DynamicConstraint.hpp
ROL_DynamicObjective.hpp
ROL_Objective.hpp
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition ROL_Objective_SerialSimOpt.hpp:77
ROL_Sketch.hpp
ROL::DynamicConstraint
Defines the time-dependent constraint operator interface for simulation-based optimization.
Definition ROL_DynamicConstraint.hpp:53
ROL::DynamicObjective
Defines the time-dependent objective function interface for simulation-based optimization....
Definition ROL_DynamicObjective.hpp:38
ROL::Objective
Provides the interface to evaluate objective functions.
Definition ROL_Objective.hpp:44
ROL::Objective::hessVec
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Definition ROL_ObjectiveDef.hpp:61
ROL::PartitionedVector
Defines the linear algebra of vector space on a generic partitioned vector.
Definition ROL_PartitionedVector.hpp:26
ROL::PartitionedVector::get
ROL::Ptr< const Vector< Real > > get(size_type i) const
Definition ROL_PartitionedVector.hpp:221
ROL::PartitionedVector::create
static Ptr< PartitionedVector > create(std::initializer_list< Vp > vs)
Definition ROL_PartitionedVector.hpp:244
ROL::ReducedDynamicObjective
Defines the reduced time-dependent objective function interface for simulation-based optimization.
Definition ROL_ReducedDynamicObjective.hpp:46
ROL::ReducedDynamicObjective::hessVec
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Definition ROL_ReducedDynamicObjective.hpp:488
ROL::ReducedDynamicObjective::updateSketch
Real updateSketch(const Vector< Real > &x, const Real tol)
Definition ROL_ReducedDynamicObjective.hpp:694
ROL::ReducedDynamicObjective::partition
const PartitionedVector< Real > & partition(const Vector< Real > &x) const
Definition ROL_ReducedDynamicObjective.hpp:104
ROL::ReducedDynamicObjective::makeDynamicVector
Ptr< Vector< Real > > makeDynamicVector(const Vector< Real > &x) const
Definition ROL_ReducedDynamicObjective.hpp:216
ROL::ReducedDynamicObjective::addMixedHessLag
void addMixedHessLag(Vector< Real > &Hv, const Vector< Real > &l, const Vector< Real > &wold, const Vector< Real > &wnew, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:1013
ROL::ReducedDynamicObjective::advanceStateSens
void advanceStateSens(Vector< Real > &wnew, const Vector< Real > &v, const Vector< Real > &wold, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:858
ROL::ReducedDynamicObjective::computeNewMixedHessLag
void computeNewMixedHessLag(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
Definition ROL_ReducedDynamicObjective.hpp:974
ROL::ReducedDynamicObjective::computeNewStateJacobian
void computeNewStateJacobian(Vector< Real > &Hv, const Vector< Real > &p, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
Definition ROL_ReducedDynamicObjective.hpp:936
ROL::ReducedDynamicObjective::size_type
typename std::vector< Real >::size_type size_type
Definition ROL_ReducedDynamicObjective.hpp:47
ROL::ReducedDynamicObjective::addAdjointSens
void addAdjointSens(Vector< Real > &Hv, const Vector< Real > &p, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:1030
ROL::ReducedDynamicObjective::setTerminalConditionHess
void setTerminalConditionHess(Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &l, const Vector< Real > &wold, const Vector< Real > &wnew, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:872
ROL::ReducedDynamicObjective::computeOldMixedHessLag
void computeOldMixedHessLag(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
Definition ROL_ReducedDynamicObjective.hpp:918
ROL::ReducedDynamicObjective::uhist_
std::vector< Ptr< Vector< Real > > > uhist_
Definition ROL_ReducedDynamicObjective.hpp:77
ROL::ReducedDynamicObjective::phist_
std::vector< Ptr< Vector< Real > > > phist_
Definition ROL_ReducedDynamicObjective.hpp:80
ROL::ReducedDynamicObjective::solveState
Real solveState(const Vector< Real > &x)
Definition ROL_ReducedDynamicObjective.hpp:659
ROL::ReducedDynamicObjective::computeNewStateHessLag
void computeNewStateHessLag(Vector< Real > &Hv, const Vector< Real > &l, const Vector< Real > &wold, const Vector< Real > &wnew, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
Definition ROL_ReducedDynamicObjective.hpp:951
ROL::ReducedDynamicObjective::setTerminalCondition
void setTerminalCondition(Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:758
ROL::ReducedDynamicObjective::gradient
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition ROL_ReducedDynamicObjective.hpp:393
ROL::ReducedDynamicObjective::whist_
std::vector< Ptr< Vector< Real > > > whist_
Definition ROL_ReducedDynamicObjective.hpp:79
ROL::ReducedDynamicObjective::update
void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Definition ROL_ReducedDynamicObjective.hpp:248
ROL::ReducedDynamicObjective::timeStamp_
const std::vector< TimeStamp< Real > > timeStamp_
Definition ROL_ReducedDynamicObjective.hpp:53
ROL::ReducedDynamicObjective::solveAdjoint
void solveAdjoint(const Vector< Real > &x)
Definition ROL_ReducedDynamicObjective.hpp:782
ROL::ReducedDynamicObjective::update
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Definition ROL_ReducedDynamicObjective.hpp:224
ROL::ReducedDynamicObjective::updateGradient
void updateGradient(Vector< Real > &g, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:846
ROL::ReducedDynamicObjective::value
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Definition ROL_ReducedDynamicObjective.hpp:352
ROL::ReducedDynamicObjective::advanceAdjointSens
void advanceAdjointSens(Vector< Real > &p, Vector< Real > &rhs, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:992
ROL::ReducedDynamicObjective::computeOldStateHessLag
void computeOldStateHessLag(Vector< Real > &Hv, const Vector< Real > &l, const Vector< Real > &wold, const Vector< Real > &wnew, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
Definition ROL_ReducedDynamicObjective.hpp:895
ROL::ReducedDynamicObjective::con_
const Ptr< DynamicConstraint< Real > > con_
Definition ROL_ReducedDynamicObjective.hpp:51
ROL::ReducedDynamicObjective::lhist_
std::vector< Ptr< Vector< Real > > > lhist_
Definition ROL_ReducedDynamicObjective.hpp:78
ROL::ReducedDynamicObjective::advanceAdjoint
void advanceAdjoint(Vector< Real > &l, Vector< Real > &rhs, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:774
ROL::ReducedDynamicObjective::partition
PartitionedVector< Real > & partition(Vector< Real > &x) const
Definition ROL_ReducedDynamicObjective.hpp:100
ROL::ReducedDynamicObjective::obj_
const Ptr< DynamicObjective< Real > > obj_
Definition ROL_ReducedDynamicObjective.hpp:50
ROL::ReducedDynamicObjective::computeControlHessLag
void computeControlHessLag(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:1002
ROL::ReducedDynamicObjective::ReducedDynamicObjective
ReducedDynamicObjective(const Ptr< DynamicObjective< Real > > &obj, const Ptr< DynamicConstraint< Real > > &con, const Ptr< Vector< Real > > &u0, const Ptr< Vector< Real > > &zvec, const Ptr< Vector< Real > > &cvec, const std::vector< TimeStamp< Real > > &timeStamp, ROL::ParameterList &pl, const Ptr< std::ostream > &stream=nullPtr)
Definition ROL_ReducedDynamicObjective.hpp:128
ROL::ReducedDynamicObjective::computeAdjointRHS
void computeAdjointRHS(Vector< Real > &rhs, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Definition ROL_ReducedDynamicObjective.hpp:765
ROL::ReducedDynamicObjective::throwError
void throwError(const int err, const std::string &sfunc, const std::string &func, const int line) const
Definition ROL_ReducedDynamicObjective.hpp:108
ROL::Vector
Defines the linear algebra or vector space interface.
Definition ROL_Vector.hpp:50
ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Vector::print
virtual void print(std::ostream &outStream) const
Definition ROL_Vector.hpp:227
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition ROL_Vector.hpp:119
ROL::ROL_INF
Real ROL_INF(void)
Definition ROL_Types.hpp:71
ROL::TimeStamp
Contains local time step information.
Definition ROL_TimeStamp.hpp:35
Generated by doxygen 1.9.8