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 
43 namespace ROL {
44 
45 template<typename Real>
46 class ReducedDynamicObjective : public Objective<Real> {
47  using size_type = typename std::vector<Real>::size_type;
48 private:
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 
127 public:
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 
655 private:
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::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::Objective
Provides the interface to evaluate objective functions.
Definition: ROL_Objective.hpp:44
ROL::PartitionedVector::create
static Ptr< PartitionedVector > create(std::initializer_list< Vp > vs)
Definition: ROL_PartitionedVector.hpp:244
ROL::ReducedDynamicObjective::update
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Definition: ROL_ReducedDynamicObjective.hpp:224
size_type
typename PV< Real >::size_type size_type
Definition: ROL_Objective_SerialSimOpt.hpp:90
ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
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::partition
const PartitionedVector< Real > & partition(const Vector< Real > &x) const
Definition: ROL_ReducedDynamicObjective.hpp:104
ROL::ReducedDynamicObjective::update
void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Definition: ROL_ReducedDynamicObjective.hpp:248
ROL::PartitionedVector::get
ROL::Ptr< const Vector< Real > > get(size_type i) const
Definition: ROL_PartitionedVector.hpp:221
ROL::ReducedDynamicObjective
Defines the reduced time-dependent objective function interface for simulation-based optimization...
Definition: ROL_ReducedDynamicObjective.hpp:46
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::ReducedDynamicObjective::makeDynamicVector
Ptr< Vector< Real > > makeDynamicVector(const Vector< Real > &x) const
Definition: ROL_ReducedDynamicObjective.hpp:216
ROL::Vector::print
virtual void print(std::ostream &outStream) const
Definition: ROL_Vector.hpp:227
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
ROL::DynamicConstraint
Defines the time-dependent constraint operator interface for simulation-based optimization.
Definition: ROL_DynamicConstraint.hpp:24
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::PartitionedVector
Defines the linear algebra of vector space on a generic partitioned vector.
Definition: ROL_PartitionedVector.hpp:26
ROL::ROL_INF
Real ROL_INF(void)
Definition: ROL_Types.hpp:71
ROL_DynamicConstraint.hpp
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::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::solveState
Real solveState(const Vector< Real > &x)
Definition: ROL_ReducedDynamicObjective.hpp:659
ROL::ReducedDynamicObjective::timeStamp_
const std::vector< TimeStamp< Real > > timeStamp_
Definition: ROL_ReducedDynamicObjective.hpp:53
ROL_DynamicObjective.hpp
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL::DynamicObjective
Defines the time-dependent objective function interface for simulation-based optimization. Computes time-local contributions of value, gradient, Hessian-vector product etc to a larger composite objective defined over the simulation time. In contrast to other objective classes Objective_TimeSimOpt has a default implementation of value which returns zero, as time-dependent simulation based optimization problems may have an objective value which depends only on the final state of the system.
Definition: ROL_DynamicObjective.hpp:38
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::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
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:111
ROL::ReducedDynamicObjective::value
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Definition: ROL_ReducedDynamicObjective.hpp:352
ROL::ReducedDynamicObjective::partition
PartitionedVector< Real > & partition(Vector< Real > &x) const
Definition: ROL_ReducedDynamicObjective.hpp:100
ROL_Sketch.hpp
ROL::ReducedDynamicObjective::updateSketch
Real updateSketch(const Vector< Real > &x, const Real tol)
Definition: ROL_ReducedDynamicObjective.hpp:694
ROL::ReducedDynamicObjective::size_type
typename std::vector< Real >::size_type size_type
Definition: ROL_ReducedDynamicObjective.hpp:47
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::obj_
const Ptr< DynamicObjective< Real > > obj_
Definition: ROL_ReducedDynamicObjective.hpp:50
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::ReducedDynamicObjective::con_
const Ptr< DynamicConstraint< Real > > con_
Definition: ROL_ReducedDynamicObjective.hpp:51
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_Objective.hpp
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::whist_
std::vector< Ptr< Vector< Real > > > whist_
Definition: ROL_ReducedDynamicObjective.hpp:79
ROL::ReducedDynamicObjective::solveAdjoint
void solveAdjoint(const Vector< Real > &x)
Definition: ROL_ReducedDynamicObjective.hpp:782
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::uhist_
std::vector< Ptr< Vector< Real > > > uhist_
Definition: ROL_ReducedDynamicObjective.hpp:77
ROL::ReducedDynamicObjective::lhist_
std::vector< Ptr< Vector< Real > > > lhist_
Definition: ROL_ReducedDynamicObjective.hpp:78
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::phist_
std::vector< Ptr< Vector< Real > > > phist_
Definition: ROL_ReducedDynamicObjective.hpp:80
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::gradient
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ReducedDynamicObjective.hpp:393
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::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::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::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::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
Generated by   doxygen 1.8.5