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