ROL
ROL_Constraint_SimOpt.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_CONSTRAINT_SIMOPT_H
11 #define ROL_CONSTRAINT_SIMOPT_H
12 
13 #include "ROL_Constraint.hpp"
14 #include "ROL_Vector_SimOpt.hpp"
15 #include "ROL_Types.hpp"
16 #include <iostream>
17 
18 namespace ROL {
19 
20 template <class Real>
21 class Constraint_SimOpt;
22 
23 }
24 
25 #include "ROL_NonlinearLeastSquaresObjective.hpp"
26 #include "ROL_SimConstraint.hpp"
27 #include "ROL_Objective_FSsolver.hpp"
28 #include "ROL_TypeU_TrustRegionAlgorithm.hpp"
29 #include "ROL_TypeE_AugmentedLagrangianAlgorithm.hpp"
30 
71 namespace ROL {
72 
73 template <class Real>
74 class Constraint_SimOpt : public Constraint<Real> {
75 private:
76  // Additional vector storage for solve
77  Ptr<Vector<Real>> unew_;
78  Ptr<Vector<Real>> jv_;
79 
80  // Default parameters for solve (backtracking Newton)
81  const Real DEFAULT_atol_;
82  const Real DEFAULT_rtol_;
83  const Real DEFAULT_stol_;
84  const Real DEFAULT_factor_;
85  const Real DEFAULT_decr_;
86  const int DEFAULT_maxit_;
87  const bool DEFAULT_print_;
88  const bool DEFAULT_zero_;
89  const int DEFAULT_solverType_;
90 
91  // User-set parameters for solve (backtracking Newton)
92 
93 protected:
94  Real atol_;
95  Real rtol_;
96  Real stol_;
97  Real factor_;
98  Real decr_;
99  int maxit_;
100  bool print_;
101  bool zero_;
102  int solverType_;
103 
104  // Flag to initialize vector storage in solve
105  bool firstSolve_;
106 
107 public:
108  Constraint_SimOpt()
109  : Constraint<Real>(),
110  unew_(nullPtr), jv_(nullPtr),
111  DEFAULT_atol_(1.e-4*std::sqrt(ROL_EPSILON<Real>())),
112  DEFAULT_rtol_(1.e0),
113  DEFAULT_stol_(std::sqrt(ROL_EPSILON<Real>())),
114  DEFAULT_factor_(0.5),
115  DEFAULT_decr_(1.e-4),
116  DEFAULT_maxit_(500),
117  DEFAULT_print_(false),
118  DEFAULT_zero_(false),
119  DEFAULT_solverType_(0),
120  atol_(DEFAULT_atol_), rtol_(DEFAULT_rtol_), stol_(DEFAULT_stol_), factor_(DEFAULT_factor_),
121  decr_(DEFAULT_decr_), maxit_(DEFAULT_maxit_), print_(DEFAULT_print_), zero_(DEFAULT_zero_),
122  solverType_(DEFAULT_solverType_), firstSolve_(true) {}
123 
129  virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) {
130  update_1(u,flag,iter);
131  update_2(z,flag,iter);
132  }
133  virtual void update( const Vector<Real> &u, const Vector<Real> &z, UpdateType type, int iter = -1 ) {
134  update_1(u,type,iter);
135  update_2(z,type,iter);
136  }
137 
143  virtual void update_1( const Vector<Real> &u, bool flag = true, int iter = -1 ) {}
144  virtual void update_1( const Vector<Real> &u, UpdateType type, int iter = -1 ) {}
145 
151  virtual void update_2( const Vector<Real> &z, bool flag = true, int iter = -1 ) {}
152  virtual void update_2( const Vector<Real> &z, UpdateType type, int iter = -1 ) {}
153 
160  virtual void solve_update( const Vector<Real> &u, const Vector<Real> &z, UpdateType type, int iter = -1) {}
161 
175  virtual void value(Vector<Real> &c,
176  const Vector<Real> &u,
177  const Vector<Real> &z,
178  Real &tol) = 0;
179 
191  virtual void solve(Vector<Real> &c,
192  Vector<Real> &u,
193  const Vector<Real> &z,
194  Real &tol) {
195  if ( zero_ ) u.zero();
196  Ptr<std::ostream> stream = makeStreamPtr(std::cout, print_);
197  solve_update(u,z,UpdateType::Initial,0);
198  value(c,u,z,tol);
199  Real cnorm = c.norm();
200  const Real ctol = std::min(atol_, rtol_*cnorm);
201  if (solverType_==0 || solverType_==3 || solverType_==4) {
202  if ( firstSolve_ ) {
203  unew_ = u.clone();
204  jv_ = u.clone();
205  firstSolve_ = false;
206  }
207  const Real one(1);
208  Real alpha(1), tmp(0);
209  int cnt = 0;
210  *stream << std::endl;
211  *stream << " Default Constraint_SimOpt::solve" << std::endl;
212  *stream << " ";
213  *stream << std::setw(6) << std::left << "iter";
214  *stream << std::setw(15) << std::left << "rnorm";
215  *stream << std::setw(15) << std::left << "alpha";
216  *stream << std::endl;
217  for (cnt = 0; cnt < maxit_; ++cnt) {
218  // Compute Newton step
219  applyInverseJacobian_1(*jv_,c,u,z,tol);
220  unew_->set(u);
221  unew_->axpy(-alpha, *jv_);
222  solve_update(*unew_,z,UpdateType::Trial);
223  value(c,*unew_,z,tol);
224  tmp = c.norm();
225  // Perform backtracking line search
226  while ( tmp > (one-decr_*alpha)*cnorm &&
227  alpha > stol_ ) {
228  alpha *= factor_;
229  unew_->set(u);
230  unew_->axpy(-alpha,*jv_);
231  solve_update(*unew_,z,UpdateType::Trial);
232  value(c,*unew_,z,tol);
233  tmp = c.norm();
234  }
235  *stream << " ";
236  *stream << std::setw(6) << std::left << cnt;
237  *stream << std::scientific << std::setprecision(6);
238  *stream << std::setw(15) << std::left << tmp;
239  *stream << std::scientific << std::setprecision(6);
240  *stream << std::setw(15) << std::left << alpha;
241  *stream << std::endl;
242  // Update iterate
243  cnorm = tmp;
244  u.set(*unew_);
245  solve_update(u,z,UpdateType::Accept,cnt);
246  if (cnorm <= ctol) break; // = covers the case of identically zero residual
247  alpha = one;
248  }
249  }
250  if (solverType_==1 || (solverType_==3 && cnorm > ctol)) {
251  Ptr<Constraint<Real>> con = makePtr<SimConstraint<Real>>(makePtrFromRef(*this),makePtrFromRef(z),true);
252  Ptr<Objective<Real>> obj = makePtr<NonlinearLeastSquaresObjective<Real>>(con,u,c,true);
253  ParameterList parlist;
254  parlist.sublist("General").set("Output Level",1);
255  parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100);
256  parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG");
257  parlist.sublist("Status Test").set("Gradient Tolerance",ctol);
258  parlist.sublist("Status Test").set("Step Tolerance",stol_);
259  parlist.sublist("Status Test").set("Iteration Limit",maxit_);
260  Ptr<TypeU::Algorithm<Real>> algo = makePtr<TypeU::TrustRegionAlgorithm<Real>>(parlist);
261  algo->run(u,*obj,*stream);
262  value(c,u,z,tol);
263  }
264  if (solverType_==2 || (solverType_==4 && cnorm > ctol)) {
265  Ptr<Constraint<Real>> con = makePtr<SimConstraint<Real>>(makePtrFromRef(*this),makePtrFromRef(z),true);
266  Ptr<Objective<Real>> obj = makePtr<Objective_FSsolver<Real>>();
267  Ptr<Vector<Real>> l = c.dual().clone();
268  ParameterList parlist;
269  parlist.sublist("General").set("Output Level",1);
270  parlist.sublist("Status Test").set("Constraint Tolerance",ctol);
271  parlist.sublist("Status Test").set("Step Tolerance",stol_);
272  parlist.sublist("Status Test").set("Iteration Limit",maxit_);
273  Ptr<TypeE::Algorithm<Real>> algo = makePtr<TypeE::AugmentedLagrangianAlgorithm<Real>>(parlist);
274  algo->run(u,*obj,*con,*l,*stream);
275  value(c,u,z,tol);
276  }
277  if (solverType_ > 4 || solverType_ < 0) {
278  ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
279  ">>> ERROR (ROL:Constraint_SimOpt:solve): Invalid solver type!");
280  }
281  }
282 
303  virtual void setSolveParameters(ParameterList &parlist) {
304  ParameterList & list = parlist.sublist("SimOpt").sublist("Solve");
305  atol_ = list.get("Absolute Residual Tolerance", DEFAULT_atol_);
306  rtol_ = list.get("Relative Residual Tolerance", DEFAULT_rtol_);
307  maxit_ = list.get("Iteration Limit", DEFAULT_maxit_);
308  decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_);
309  stol_ = list.get("Step Tolerance", DEFAULT_stol_);
310  factor_ = list.get("Backtracking Factor", DEFAULT_factor_);
311  print_ = list.get("Output Iteration History", DEFAULT_print_);
312  zero_ = list.get("Zero Initial Guess", DEFAULT_zero_);
313  solverType_ = list.get("Solver Type", DEFAULT_solverType_);
314  }
315 
331  virtual void applyJacobian_1(Vector<Real> &jv,
332  const Vector<Real> &v,
333  const Vector<Real> &u,
334  const Vector<Real> &z,
335  Real &tol) {
336  Real ctol = std::sqrt(ROL_EPSILON<Real>());
337  // Compute step length
338  Real h = tol;
339  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
340  h = std::max(1.0,u.norm()/v.norm())*tol;
341  }
342  // Update state vector to u + hv
343  Ptr<Vector<Real>> unew = u.clone();
344  unew->set(u);
345  unew->axpy(h,v);
346  // Compute new constraint value
347  update(*unew,z,UpdateType::Temp);
348  value(jv,*unew,z,ctol);
349  // Compute current constraint value
350  Ptr<Vector<Real>> cold = jv.clone();
351  update(u,z,UpdateType::Temp);
352  value(*cold,u,z,ctol);
353  // Compute Newton quotient
354  jv.axpy(-1.0,*cold);
355  jv.scale(1.0/h);
356  }
357 
358 
374  virtual void applyJacobian_2(Vector<Real> &jv,
375  const Vector<Real> &v,
376  const Vector<Real> &u,
377  const Vector<Real> &z,
378  Real &tol) {
379  Real ctol = std::sqrt(ROL_EPSILON<Real>());
380  // Compute step length
381  Real h = tol;
382  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
383  h = std::max(1.0,u.norm()/v.norm())*tol;
384  }
385  // Update state vector to u + hv
386  Ptr<Vector<Real>> znew = z.clone();
387  znew->set(z);
388  znew->axpy(h,v);
389  // Compute new constraint value
390  update(u,*znew,UpdateType::Temp);
391  value(jv,u,*znew,ctol);
392  // Compute current constraint value
393  Ptr<Vector<Real>> cold = jv.clone();
394  update(u,z,UpdateType::Temp);
395  value(*cold,u,z,ctol);
396  // Compute Newton quotient
397  jv.axpy(-1.0,*cold);
398  jv.scale(1.0/h);
399  }
400 
416  virtual void applyInverseJacobian_1(Vector<Real> &ijv,
417  const Vector<Real> &v,
418  const Vector<Real> &u,
419  const Vector<Real> &z,
420  Real &tol) {
421  ROL_TEST_FOR_EXCEPTION(true, std::logic_error,
422  "The method applyInverseJacobian_1 is used but not implemented!\n");
423  }
424 
440  virtual void applyAdjointJacobian_1(Vector<Real> &ajv,
441  const Vector<Real> &v,
442  const Vector<Real> &u,
443  const Vector<Real> &z,
444  Real &tol) {
445  applyAdjointJacobian_1(ajv, v, u, z, v.dual(), tol);
446  }
447 
448 
466  virtual void applyAdjointJacobian_1(Vector<Real> &ajv,
467  const Vector<Real> &v,
468  const Vector<Real> &u,
469  const Vector<Real> &z,
470  const Vector<Real> &dualv,
471  Real &tol) {
472  Real ctol = std::sqrt(ROL_EPSILON<Real>());
473  Real h = tol;
474  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
475  h = std::max(1.0,u.norm()/v.norm())*tol;
476  }
477  Ptr<Vector<Real>> cold = dualv.clone();
478  Ptr<Vector<Real>> cnew = dualv.clone();
479  update(u,z,UpdateType::Temp);
480  value(*cold,u,z,ctol);
481  Ptr<Vector<Real>> unew = u.clone();
482  ajv.zero();
483  for (int i = 0; i < u.dimension(); i++) {
484  unew->set(u);
485  unew->axpy(h,*(u.basis(i)));
486  update(*unew,z,UpdateType::Temp);
487  value(*cnew,*unew,z,ctol);
488  cnew->axpy(-1.0,*cold);
489  cnew->scale(1.0/h);
490  ajv.axpy(cnew->dot(v),*((u.dual()).basis(i)));
491  }
492  update(u,z,UpdateType::Temp);
493  }
494 
495 
511  virtual void applyAdjointJacobian_2(Vector<Real> &ajv,
512  const Vector<Real> &v,
513  const Vector<Real> &u,
514  const Vector<Real> &z,
515  Real &tol) {
516  applyAdjointJacobian_2(ajv, v, u, z, v.dual(), tol);
517  }
518 
519 
537  virtual void applyAdjointJacobian_2(Vector<Real> &ajv,
538  const Vector<Real> &v,
539  const Vector<Real> &u,
540  const Vector<Real> &z,
541  const Vector<Real> &dualv,
542  Real &tol) {
543  Real ctol = std::sqrt(ROL_EPSILON<Real>());
544  Real h = tol;
545  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
546  h = std::max(1.0,u.norm()/v.norm())*tol;
547  }
548  Ptr<Vector<Real>> cold = dualv.clone();
549  Ptr<Vector<Real>> cnew = dualv.clone();
550  update(u,z,UpdateType::Temp);
551  value(*cold,u,z,ctol);
552  Ptr<Vector<Real>> znew = z.clone();
553  ajv.zero();
554  for (int i = 0; i < z.dimension(); i++) {
555  znew->set(z);
556  znew->axpy(h,*(z.basis(i)));
557  update(u,*znew,UpdateType::Temp);
558  value(*cnew,u,*znew,ctol);
559  cnew->axpy(-1.0,*cold);
560  cnew->scale(1.0/h);
561  ajv.axpy(cnew->dot(v),*((z.dual()).basis(i)));
562  }
563  update(u,z,UpdateType::Temp);
564  }
565 
581  virtual void applyInverseAdjointJacobian_1(Vector<Real> &iajv,
582  const Vector<Real> &v,
583  const Vector<Real> &u,
584  const Vector<Real> &z,
585  Real &tol) {
586  ROL_TEST_FOR_EXCEPTION(true, std::logic_error,
587  "The method applyInverseAdjointJacobian_1 is used but not implemented!\n");
588  };
589 
607  virtual void applyAdjointHessian_11(Vector<Real> &ahwv,
608  const Vector<Real> &w,
609  const Vector<Real> &v,
610  const Vector<Real> &u,
611  const Vector<Real> &z,
612  Real &tol) {
613  Real jtol = std::sqrt(ROL_EPSILON<Real>());
614  // Compute step size
615  Real h = tol;
616  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
617  h = std::max(1.0,u.norm()/v.norm())*tol;
618  }
619  // Evaluate Jacobian at new state
620  Ptr<Vector<Real>> unew = u.clone();
621  unew->set(u);
622  unew->axpy(h,v);
623  update(*unew,z,UpdateType::Temp);
624  applyAdjointJacobian_1(ahwv,w,*unew,z,jtol);
625  // Evaluate Jacobian at old state
626  Ptr<Vector<Real>> jv = ahwv.clone();
627  update(u,z,UpdateType::Temp);
628  applyAdjointJacobian_1(*jv,w,u,z,jtol);
629  // Compute Newton quotient
630  ahwv.axpy(-1.0,*jv);
631  ahwv.scale(1.0/h);
632  }
633 
634 
652  virtual void applyAdjointHessian_12(Vector<Real> &ahwv,
653  const Vector<Real> &w,
654  const Vector<Real> &v,
655  const Vector<Real> &u,
656  const Vector<Real> &z,
657  Real &tol) {
658  Real jtol = std::sqrt(ROL_EPSILON<Real>());
659  // Compute step size
660  Real h = tol;
661  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
662  h = std::max(1.0,u.norm()/v.norm())*tol;
663  }
664  // Evaluate Jacobian at new state
665  Ptr<Vector<Real>> unew = u.clone();
666  unew->set(u);
667  unew->axpy(h,v);
668  update(*unew,z,UpdateType::Temp);
669  applyAdjointJacobian_2(ahwv,w,*unew,z,jtol);
670  // Evaluate Jacobian at old state
671  Ptr<Vector<Real>> jv = ahwv.clone();
672  update(u,z,UpdateType::Temp);
673  applyAdjointJacobian_2(*jv,w,u,z,jtol);
674  // Compute Newton quotient
675  ahwv.axpy(-1.0,*jv);
676  ahwv.scale(1.0/h);
677  }
678 
679 
697  virtual void applyAdjointHessian_21(Vector<Real> &ahwv,
698  const Vector<Real> &w,
699  const Vector<Real> &v,
700  const Vector<Real> &u,
701  const Vector<Real> &z,
702  Real &tol) {
703  Real jtol = std::sqrt(ROL_EPSILON<Real>());
704  // Compute step size
705  Real h = tol;
706  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
707  h = std::max(1.0,u.norm()/v.norm())*tol;
708  }
709  // Evaluate Jacobian at new control
710  Ptr<Vector<Real>> znew = z.clone();
711  znew->set(z);
712  znew->axpy(h,v);
713  update(u,*znew,UpdateType::Temp);
714  applyAdjointJacobian_1(ahwv,w,u,*znew,jtol);
715  // Evaluate Jacobian at old control
716  Ptr<Vector<Real>> jv = ahwv.clone();
717  update(u,z,UpdateType::Temp);
718  applyAdjointJacobian_1(*jv,w,u,z,jtol);
719  // Compute Newton quotient
720  ahwv.axpy(-1.0,*jv);
721  ahwv.scale(1.0/h);
722  }
723 
741  virtual void applyAdjointHessian_22(Vector<Real> &ahwv,
742  const Vector<Real> &w,
743  const Vector<Real> &v,
744  const Vector<Real> &u,
745  const Vector<Real> &z,
746  Real &tol) {
747  Real jtol = std::sqrt(ROL_EPSILON<Real>());
748  // Compute step size
749  Real h = tol;
750  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
751  h = std::max(1.0,u.norm()/v.norm())*tol;
752  }
753  // Evaluate Jacobian at new control
754  Ptr<Vector<Real>> znew = z.clone();
755  znew->set(z);
756  znew->axpy(h,v);
757  update(u,*znew,UpdateType::Temp);
758  applyAdjointJacobian_2(ahwv,w,u,*znew,jtol);
759  // Evaluate Jacobian at old control
760  Ptr<Vector<Real>> jv = ahwv.clone();
761  update(u,z,UpdateType::Temp);
762  applyAdjointJacobian_2(*jv,w,u,z,jtol);
763  // Compute Newton quotient
764  ahwv.axpy(-1.0,*jv);
765  ahwv.scale(1.0/h);
766 }
767 
806  virtual std::vector<Real> solveAugmentedSystem(Vector<Real> &v1,
807  Vector<Real> &v2,
808  const Vector<Real> &b1,
809  const Vector<Real> &b2,
810  const Vector<Real> &x,
811  Real &tol) {
812  return Constraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,x,tol);
813  }
814 
815 
834  virtual void applyPreconditioner(Vector<Real> &pv,
835  const Vector<Real> &v,
836  const Vector<Real> &x,
837  const Vector<Real> &g,
838  Real &tol) {
839  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(x);
840  Ptr<Vector<Real>> ijv = (xs.get_1())->clone();
841 
842  try {
843  applyInverseJacobian_1(*ijv, v, *(xs.get_1()), *(xs.get_2()), tol);
844  }
845  catch (const std::logic_error &e) {
846  Constraint<Real>::applyPreconditioner(pv, v, x, g, tol);
847  return;
848  }
849 
850  const Vector_SimOpt<Real> &gs = dynamic_cast<const Vector_SimOpt<Real>&>(g);
851  Ptr<Vector<Real>> ijv_dual = (gs.get_1())->clone();
852  ijv_dual->set(ijv->dual());
853 
854  try {
855  applyInverseAdjointJacobian_1(pv, *ijv_dual, *(xs.get_1()), *(xs.get_2()), tol);
856  }
857  catch (const std::logic_error &e) {
858  Constraint<Real>::applyPreconditioner(pv, v, x, g, tol);
859  return;
860  }
861 
862  }
863 
869  virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
870  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
871  dynamic_cast<const Vector<Real>&>(x));
872  update(*(xs.get_1()),*(xs.get_2()),flag,iter);
873  }
874  virtual void update( const Vector<Real> &x, UpdateType type, int iter = -1 ) {
875  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
876  dynamic_cast<const Vector<Real>&>(x));
877  update(*(xs.get_1()),*(xs.get_2()),type,iter);
878  }
879 
880  virtual void value(Vector<Real> &c,
881  const Vector<Real> &x,
882  Real &tol) {
883  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
884  dynamic_cast<const Vector<Real>&>(x));
885  value(c,*(xs.get_1()),*(xs.get_2()),tol);
886  }
887 
888 
889  virtual void applyJacobian(Vector<Real> &jv,
890  const Vector<Real> &v,
891  const Vector<Real> &x,
892  Real &tol) {
893  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
894  dynamic_cast<const Vector<Real>&>(x));
895  const Vector_SimOpt<Real> &vs = dynamic_cast<const Vector_SimOpt<Real>&>(
896  dynamic_cast<const Vector<Real>&>(v));
897  applyJacobian_1(jv,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
898  Ptr<Vector<Real>> jv2 = jv.clone();
899  applyJacobian_2(*jv2,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
900  jv.plus(*jv2);
901  }
902 
903 
904  using Constraint<Real>::applyAdjointJacobian;
905  virtual void applyAdjointJacobian(Vector<Real> &ajv,
906  const Vector<Real> &v,
907  const Vector<Real> &x,
908  Real &tol) {
909  Vector_SimOpt<Real> &ajvs = dynamic_cast<Vector_SimOpt<Real>&>(
910  dynamic_cast<Vector<Real>&>(ajv));
911  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
912  dynamic_cast<const Vector<Real>&>(x));
913  Ptr<Vector<Real>> ajv1 = (ajvs.get_1())->clone();
914  applyAdjointJacobian_1(*ajv1,v,*(xs.get_1()),*(xs.get_2()),tol);
915  ajvs.set_1(*ajv1);
916  Ptr<Vector<Real>> ajv2 = (ajvs.get_2())->clone();
917  applyAdjointJacobian_2(*ajv2,v,*(xs.get_1()),*(xs.get_2()),tol);
918  ajvs.set_2(*ajv2);
919  }
920 
921 
922  virtual void applyAdjointHessian(Vector<Real> &ahwv,
923  const Vector<Real> &w,
924  const Vector<Real> &v,
925  const Vector<Real> &x,
926  Real &tol) {
927  Vector_SimOpt<Real> &ahwvs = dynamic_cast<Vector_SimOpt<Real>&>(
928  dynamic_cast<Vector<Real>&>(ahwv));
929  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
930  dynamic_cast<const Vector<Real>&>(x));
931  const Vector_SimOpt<Real> &vs = dynamic_cast<const Vector_SimOpt<Real>&>(
932  dynamic_cast<const Vector<Real>&>(v));
933  // Block-row 1
934  Ptr<Vector<Real>> C11 = (ahwvs.get_1())->clone();
935  Ptr<Vector<Real>> C21 = (ahwvs.get_1())->clone();
936  applyAdjointHessian_11(*C11,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
937  applyAdjointHessian_21(*C21,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
938  C11->plus(*C21);
939  ahwvs.set_1(*C11);
940  // Block-row 2
941  Ptr<Vector<Real>> C12 = (ahwvs.get_2())->clone();
942  Ptr<Vector<Real>> C22 = (ahwvs.get_2())->clone();
943  applyAdjointHessian_12(*C12,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
944  applyAdjointHessian_22(*C22,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
945  C22->plus(*C12);
946  ahwvs.set_2(*C22);
947  }
948 
949 
950 
951  virtual Real checkSolve(const Vector<Real> &u,
952  const Vector<Real> &z,
953  const Vector<Real> &c,
954  const bool printToStream = true,
955  std::ostream & outStream = std::cout) {
956  // Solve constraint for u.
957  Real tol = ROL_EPSILON<Real>();
958  Ptr<Vector<Real>> r = c.clone();
959  Ptr<Vector<Real>> s = u.clone();
960  solve(*r,*s,z,tol);
961  // Evaluate constraint residual at (u,z).
962  Ptr<Vector<Real>> cs = c.clone();
963  update(*s,z,UpdateType::Temp);
964  value(*cs,*s,z,tol);
965  // Output norm of residual.
966  Real rnorm = r->norm();
967  Real cnorm = cs->norm();
968  if ( printToStream ) {
969  std::stringstream hist;
970  hist << std::scientific << std::setprecision(8);
971  hist << "\nTest SimOpt solve at feasible (u,z):\n";
972  hist << " Solver Residual = " << rnorm << "\n";
973  hist << " ||c(u,z)|| = " << cnorm << "\n";
974  outStream << hist.str();
975  }
976  return cnorm;
977  }
978 
979 
992  virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w,
993  const Vector<Real> &v,
994  const Vector<Real> &u,
995  const Vector<Real> &z,
996  const bool printToStream = true,
997  std::ostream & outStream = std::cout) {
998  return checkAdjointConsistencyJacobian_1(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
999  }
1000 
1001 
1017  virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w,
1018  const Vector<Real> &v,
1019  const Vector<Real> &u,
1020  const Vector<Real> &z,
1021  const Vector<Real> &dualw,
1022  const Vector<Real> &dualv,
1023  const bool printToStream = true,
1024  std::ostream & outStream = std::cout) {
1025  Real tol = ROL_EPSILON<Real>();
1026  Ptr<Vector<Real>> Jv = dualw.clone();
1027  update(u,z,UpdateType::Temp);
1028  applyJacobian_1(*Jv,v,u,z,tol);
1029  //Real wJv = w.dot(Jv->dual());
1030  Real wJv = w.apply(*Jv);
1031  Ptr<Vector<Real>> Jw = dualv.clone();
1032  update(u,z,UpdateType::Temp);
1033  applyAdjointJacobian_1(*Jw,w,u,z,tol);
1034  //Real vJw = v.dot(Jw->dual());
1035  Real vJw = v.apply(*Jw);
1036  Real diff = std::abs(wJv-vJw);
1037  if ( printToStream ) {
1038  std::stringstream hist;
1039  hist << std::scientific << std::setprecision(8);
1040  hist << "\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1041  << diff << "\n";
1042  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1043  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1044  outStream << hist.str();
1045  }
1046  return diff;
1047  }
1048 
1049 
1062  virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w,
1063  const Vector<Real> &v,
1064  const Vector<Real> &u,
1065  const Vector<Real> &z,
1066  const bool printToStream = true,
1067  std::ostream & outStream = std::cout) {
1068  return checkAdjointConsistencyJacobian_2(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
1069  }
1070 
1086  virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w,
1087  const Vector<Real> &v,
1088  const Vector<Real> &u,
1089  const Vector<Real> &z,
1090  const Vector<Real> &dualw,
1091  const Vector<Real> &dualv,
1092  const bool printToStream = true,
1093  std::ostream & outStream = std::cout) {
1094  Real tol = ROL_EPSILON<Real>();
1095  Ptr<Vector<Real>> Jv = dualw.clone();
1096  update(u,z,UpdateType::Temp);
1097  applyJacobian_2(*Jv,v,u,z,tol);
1098  //Real wJv = w.dot(Jv->dual());
1099  Real wJv = w.apply(*Jv);
1100  Ptr<Vector<Real>> Jw = dualv.clone();
1101  update(u,z,UpdateType::Temp);
1102  applyAdjointJacobian_2(*Jw,w,u,z,tol);
1103  //Real vJw = v.dot(Jw->dual());
1104  Real vJw = v.apply(*Jw);
1105  Real diff = std::abs(wJv-vJw);
1106  if ( printToStream ) {
1107  std::stringstream hist;
1108  hist << std::scientific << std::setprecision(8);
1109  hist << "\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1110  << diff << "\n";
1111  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1112  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1113  outStream << hist.str();
1114  }
1115  return diff;
1116  }
1117 
1118  virtual Real checkInverseJacobian_1(const Vector<Real> &jv,
1119  const Vector<Real> &v,
1120  const Vector<Real> &u,
1121  const Vector<Real> &z,
1122  const bool printToStream = true,
1123  std::ostream & outStream = std::cout) {
1124  Real tol = ROL_EPSILON<Real>();
1125  Ptr<Vector<Real>> Jv = jv.clone();
1126  update(u,z,UpdateType::Temp);
1127  applyJacobian_1(*Jv,v,u,z,tol);
1128  Ptr<Vector<Real>> iJJv = u.clone();
1129  //update(u,z); // Does this update do anything?
1130  applyInverseJacobian_1(*iJJv,*Jv,u,z,tol);
1131  Ptr<Vector<Real>> diff = v.clone();
1132  diff->set(v);
1133  diff->axpy(-1.0,*iJJv);
1134  Real dnorm = diff->norm();
1135  Real vnorm = v.norm();
1136  if ( printToStream ) {
1137  std::stringstream hist;
1138  hist << std::scientific << std::setprecision(8);
1139  hist << "\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = "
1140  << dnorm << "\n";
1141  hist << " ||v|| = " << vnorm << "\n";
1142  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1143  outStream << hist.str();
1144  }
1145  return dnorm;
1146  }
1147 
1148  virtual Real checkInverseAdjointJacobian_1(const Vector<Real> &jv,
1149  const Vector<Real> &v,
1150  const Vector<Real> &u,
1151  const Vector<Real> &z,
1152  const bool printToStream = true,
1153  std::ostream & outStream = std::cout) {
1154  Real tol = ROL_EPSILON<Real>();
1155  Ptr<Vector<Real>> Jv = jv.clone();
1156  update(u,z,UpdateType::Temp);
1157  applyAdjointJacobian_1(*Jv,v,u,z,tol);
1158  Ptr<Vector<Real>> iJJv = v.clone();
1159  //update(u,z);
1160  applyInverseAdjointJacobian_1(*iJJv,*Jv,u,z,tol);
1161  Ptr<Vector<Real>> diff = v.clone();
1162  diff->set(v);
1163  diff->axpy(-1.0,*iJJv);
1164  Real dnorm = diff->norm();
1165  Real vnorm = v.norm();
1166  if ( printToStream ) {
1167  std::stringstream hist;
1168  hist << std::scientific << std::setprecision(8);
1169  hist << "\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = "
1170  << dnorm << "\n";
1171  hist << " ||v|| = " << vnorm << "\n";
1172  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1173  outStream << hist.str();
1174  }
1175  return dnorm;
1176  }
1177 
1178 
1179 
1180  std::vector<std::vector<Real>> checkApplyJacobian_1(const Vector<Real> &u,
1181  const Vector<Real> &z,
1182  const Vector<Real> &v,
1183  const Vector<Real> &jv,
1184  const bool printToStream = true,
1185  std::ostream & outStream = std::cout,
1186  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1187  const int order = 1) {
1188  std::vector<Real> steps(numSteps);
1189  for(int i=0;i<numSteps;++i) {
1190  steps[i] = pow(10,-i);
1191  }
1192 
1193  return checkApplyJacobian_1(u,z,v,jv,steps,printToStream,outStream,order);
1194  }
1195 
1196 
1197 
1198 
1199  std::vector<std::vector<Real>> checkApplyJacobian_1(const Vector<Real> &u,
1200  const Vector<Real> &z,
1201  const Vector<Real> &v,
1202  const Vector<Real> &jv,
1203  const std::vector<Real> &steps,
1204  const bool printToStream = true,
1205  std::ostream & outStream = std::cout,
1206  const int order = 1) {
1207 
1208  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1209  "Error: finite difference order must be 1,2,3, or 4" );
1210 
1211  Real one(1.0);
1212 
1213  using Finite_Difference_Arrays::shifts;
1214  using Finite_Difference_Arrays::weights;
1215 
1216  Real tol = std::sqrt(ROL_EPSILON<Real>());
1217 
1218  int numSteps = steps.size();
1219  int numVals = 4;
1220  std::vector<Real> tmp(numVals);
1221  std::vector<std::vector<Real>> jvCheck(numSteps, tmp);
1222 
1223  // Save the format state of the original outStream.
1224  nullstream oldFormatState;
1225  oldFormatState.copyfmt(outStream);
1226 
1227  // Compute constraint value at x.
1228  Ptr<Vector<Real>> c = jv.clone();
1229  this->update(u,z,UpdateType::Temp);
1230  this->value(*c, u, z, tol);
1231 
1232  // Compute (Jacobian at x) times (vector v).
1233  Ptr<Vector<Real>> Jv = jv.clone();
1234  this->applyJacobian_1(*Jv, v, u, z, tol);
1235  Real normJv = Jv->norm();
1236 
1237  // Temporary vectors.
1238  Ptr<Vector<Real>> cdif = jv.clone();
1239  Ptr<Vector<Real>> cnew = jv.clone();
1240  Ptr<Vector<Real>> unew = u.clone();
1241 
1242  for (int i=0; i<numSteps; i++) {
1243 
1244  Real eta = steps[i];
1245 
1246  unew->set(u);
1247 
1248  cdif->set(*c);
1249  cdif->scale(weights[order-1][0]);
1250 
1251  for(int j=0; j<order; ++j) {
1252 
1253  unew->axpy(eta*shifts[order-1][j], v);
1254 
1255  if( weights[order-1][j+1] != 0 ) {
1256  this->update(*unew,z,UpdateType::Temp);
1257  this->value(*cnew,*unew,z,tol);
1258  cdif->axpy(weights[order-1][j+1],*cnew);
1259  }
1260 
1261  }
1262 
1263  cdif->scale(one/eta);
1264 
1265  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1266  jvCheck[i][0] = eta;
1267  jvCheck[i][1] = normJv;
1268  jvCheck[i][2] = cdif->norm();
1269  cdif->axpy(-one, *Jv);
1270  jvCheck[i][3] = cdif->norm();
1271 
1272  if (printToStream) {
1273  std::stringstream hist;
1274  if (i==0) {
1275  hist << std::right
1276  << std::setw(20) << "Step size"
1277  << std::setw(20) << "norm(Jac*vec)"
1278  << std::setw(20) << "norm(FD approx)"
1279  << std::setw(20) << "norm(abs error)"
1280  << "\n"
1281  << std::setw(20) << "---------"
1282  << std::setw(20) << "-------------"
1283  << std::setw(20) << "---------------"
1284  << std::setw(20) << "---------------"
1285  << "\n";
1286  }
1287  hist << std::scientific << std::setprecision(11) << std::right
1288  << std::setw(20) << jvCheck[i][0]
1289  << std::setw(20) << jvCheck[i][1]
1290  << std::setw(20) << jvCheck[i][2]
1291  << std::setw(20) << jvCheck[i][3]
1292  << "\n";
1293  outStream << hist.str();
1294  }
1295 
1296  }
1297 
1298  // Reset format state of outStream.
1299  outStream.copyfmt(oldFormatState);
1300 
1301  return jvCheck;
1302  } // checkApplyJacobian
1303 
1304 
1305  std::vector<std::vector<Real>> checkApplyJacobian_2(const Vector<Real> &u,
1306  const Vector<Real> &z,
1307  const Vector<Real> &v,
1308  const Vector<Real> &jv,
1309  const bool printToStream = true,
1310  std::ostream & outStream = std::cout,
1311  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1312  const int order = 1) {
1313  std::vector<Real> steps(numSteps);
1314  for(int i=0;i<numSteps;++i) {
1315  steps[i] = pow(10,-i);
1316  }
1317 
1318  return checkApplyJacobian_2(u,z,v,jv,steps,printToStream,outStream,order);
1319  }
1320 
1321 
1322 
1323 
1324  std::vector<std::vector<Real>> checkApplyJacobian_2(const Vector<Real> &u,
1325  const Vector<Real> &z,
1326  const Vector<Real> &v,
1327  const Vector<Real> &jv,
1328  const std::vector<Real> &steps,
1329  const bool printToStream = true,
1330  std::ostream & outStream = std::cout,
1331  const int order = 1) {
1332 
1333  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1334  "Error: finite difference order must be 1,2,3, or 4" );
1335 
1336  Real one(1.0);
1337 
1338  using Finite_Difference_Arrays::shifts;
1339  using Finite_Difference_Arrays::weights;
1340 
1341  Real tol = std::sqrt(ROL_EPSILON<Real>());
1342 
1343  int numSteps = steps.size();
1344  int numVals = 4;
1345  std::vector<Real> tmp(numVals);
1346  std::vector<std::vector<Real>> jvCheck(numSteps, tmp);
1347 
1348  // Save the format state of the original outStream.
1349  nullstream oldFormatState;
1350  oldFormatState.copyfmt(outStream);
1351 
1352  // Compute constraint value at x.
1353  Ptr<Vector<Real>> c = jv.clone();
1354  this->update(u,z,UpdateType::Temp);
1355  this->value(*c, u, z, tol);
1356 
1357  // Compute (Jacobian at x) times (vector v).
1358  Ptr<Vector<Real>> Jv = jv.clone();
1359  this->applyJacobian_2(*Jv, v, u, z, tol);
1360  Real normJv = Jv->norm();
1361 
1362  // Temporary vectors.
1363  Ptr<Vector<Real>> cdif = jv.clone();
1364  Ptr<Vector<Real>> cnew = jv.clone();
1365  Ptr<Vector<Real>> znew = z.clone();
1366 
1367  for (int i=0; i<numSteps; i++) {
1368 
1369  Real eta = steps[i];
1370 
1371  znew->set(z);
1372 
1373  cdif->set(*c);
1374  cdif->scale(weights[order-1][0]);
1375 
1376  for(int j=0; j<order; ++j) {
1377 
1378  znew->axpy(eta*shifts[order-1][j], v);
1379 
1380  if( weights[order-1][j+1] != 0 ) {
1381  this->update(u,*znew,UpdateType::Temp);
1382  this->value(*cnew,u,*znew,tol);
1383  cdif->axpy(weights[order-1][j+1],*cnew);
1384  }
1385 
1386  }
1387 
1388  cdif->scale(one/eta);
1389 
1390  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1391  jvCheck[i][0] = eta;
1392  jvCheck[i][1] = normJv;
1393  jvCheck[i][2] = cdif->norm();
1394  cdif->axpy(-one, *Jv);
1395  jvCheck[i][3] = cdif->norm();
1396 
1397  if (printToStream) {
1398  std::stringstream hist;
1399  if (i==0) {
1400  hist << std::right
1401  << std::setw(20) << "Step size"
1402  << std::setw(20) << "norm(Jac*vec)"
1403  << std::setw(20) << "norm(FD approx)"
1404  << std::setw(20) << "norm(abs error)"
1405  << "\n"
1406  << std::setw(20) << "---------"
1407  << std::setw(20) << "-------------"
1408  << std::setw(20) << "---------------"
1409  << std::setw(20) << "---------------"
1410  << "\n";
1411  }
1412  hist << std::scientific << std::setprecision(11) << std::right
1413  << std::setw(20) << jvCheck[i][0]
1414  << std::setw(20) << jvCheck[i][1]
1415  << std::setw(20) << jvCheck[i][2]
1416  << std::setw(20) << jvCheck[i][3]
1417  << "\n";
1418  outStream << hist.str();
1419  }
1420 
1421  }
1422 
1423  // Reset format state of outStream.
1424  outStream.copyfmt(oldFormatState);
1425 
1426  return jvCheck;
1427  } // checkApplyJacobian
1428 
1429 
1430 
1431  std::vector<std::vector<Real>> checkApplyAdjointHessian_11(const Vector<Real> &u,
1432  const Vector<Real> &z,
1433  const Vector<Real> &p,
1434  const Vector<Real> &v,
1435  const Vector<Real> &hv,
1436  const bool printToStream = true,
1437  std::ostream & outStream = std::cout,
1438  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1439  const int order = 1 ) {
1440  std::vector<Real> steps(numSteps);
1441  for(int i=0;i<numSteps;++i) {
1442  steps[i] = pow(10,-i);
1443  }
1444 
1445  return checkApplyAdjointHessian_11(u,z,p,v,hv,steps,printToStream,outStream,order);
1446 
1447  }
1448 
1449  std::vector<std::vector<Real>> checkApplyAdjointHessian_11(const Vector<Real> &u,
1450  const Vector<Real> &z,
1451  const Vector<Real> &p,
1452  const Vector<Real> &v,
1453  const Vector<Real> &hv,
1454  const std::vector<Real> &steps,
1455  const bool printToStream = true,
1456  std::ostream & outStream = std::cout,
1457  const int order = 1 ) {
1458  using Finite_Difference_Arrays::shifts;
1459  using Finite_Difference_Arrays::weights;
1460 
1461  Real one(1.0);
1462 
1463  Real tol = std::sqrt(ROL_EPSILON<Real>());
1464 
1465  int numSteps = steps.size();
1466  int numVals = 4;
1467  std::vector<Real> tmp(numVals);
1468  std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1469 
1470  // Temporary vectors.
1471  Ptr<Vector<Real>> AJdif = hv.clone();
1472  Ptr<Vector<Real>> AJp = hv.clone();
1473  Ptr<Vector<Real>> AHpv = hv.clone();
1474  Ptr<Vector<Real>> AJnew = hv.clone();
1475  Ptr<Vector<Real>> unew = u.clone();
1476 
1477  // Save the format state of the original outStream.
1478  nullstream oldFormatState;
1479  oldFormatState.copyfmt(outStream);
1480 
1481  // Apply adjoint Jacobian to p.
1482  this->update(u,z,UpdateType::Temp);
1483  this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1484 
1485  // Apply adjoint Hessian at (u,z), in direction v, to p.
1486  this->applyAdjointHessian_11(*AHpv, p, v, u, z, tol);
1487  Real normAHpv = AHpv->norm();
1488 
1489  for (int i=0; i<numSteps; i++) {
1490 
1491  Real eta = steps[i];
1492 
1493  // Apply adjoint Jacobian to p at (u+eta*v,z).
1494  unew->set(u);
1495 
1496  AJdif->set(*AJp);
1497  AJdif->scale(weights[order-1][0]);
1498 
1499  for(int j=0; j<order; ++j) {
1500 
1501  unew->axpy(eta*shifts[order-1][j],v);
1502 
1503  if( weights[order-1][j+1] != 0 ) {
1504  this->update(*unew,z,UpdateType::Temp);
1505  this->applyAdjointJacobian_1(*AJnew, p, *unew, z, tol);
1506  AJdif->axpy(weights[order-1][j+1],*AJnew);
1507  }
1508  }
1509 
1510  AJdif->scale(one/eta);
1511 
1512  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1513  ahpvCheck[i][0] = eta;
1514  ahpvCheck[i][1] = normAHpv;
1515  ahpvCheck[i][2] = AJdif->norm();
1516  AJdif->axpy(-one, *AHpv);
1517  ahpvCheck[i][3] = AJdif->norm();
1518 
1519  if (printToStream) {
1520  std::stringstream hist;
1521  if (i==0) {
1522  hist << std::right
1523  << std::setw(20) << "Step size"
1524  << std::setw(20) << "norm(adj(H)(u,v))"
1525  << std::setw(20) << "norm(FD approx)"
1526  << std::setw(20) << "norm(abs error)"
1527  << "\n"
1528  << std::setw(20) << "---------"
1529  << std::setw(20) << "-----------------"
1530  << std::setw(20) << "---------------"
1531  << std::setw(20) << "---------------"
1532  << "\n";
1533  }
1534  hist << std::scientific << std::setprecision(11) << std::right
1535  << std::setw(20) << ahpvCheck[i][0]
1536  << std::setw(20) << ahpvCheck[i][1]
1537  << std::setw(20) << ahpvCheck[i][2]
1538  << std::setw(20) << ahpvCheck[i][3]
1539  << "\n";
1540  outStream << hist.str();
1541  }
1542 
1543  }
1544 
1545  // Reset format state of outStream.
1546  outStream.copyfmt(oldFormatState);
1547 
1548  return ahpvCheck;
1549  } // checkApplyAdjointHessian_11
1550 
1554  std::vector<std::vector<Real>> checkApplyAdjointHessian_21(const Vector<Real> &u,
1555  const Vector<Real> &z,
1556  const Vector<Real> &p,
1557  const Vector<Real> &v,
1558  const Vector<Real> &hv,
1559  const bool printToStream = true,
1560  std::ostream & outStream = std::cout,
1561  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1562  const int order = 1 ) {
1563  std::vector<Real> steps(numSteps);
1564  for(int i=0;i<numSteps;++i) {
1565  steps[i] = pow(10,-i);
1566  }
1567 
1568  return checkApplyAdjointHessian_21(u,z,p,v,hv,steps,printToStream,outStream,order);
1569 
1570  }
1571 
1575  std::vector<std::vector<Real>> checkApplyAdjointHessian_21(const Vector<Real> &u,
1576  const Vector<Real> &z,
1577  const Vector<Real> &p,
1578  const Vector<Real> &v,
1579  const Vector<Real> &hv,
1580  const std::vector<Real> &steps,
1581  const bool printToStream = true,
1582  std::ostream & outStream = std::cout,
1583  const int order = 1 ) {
1584  using Finite_Difference_Arrays::shifts;
1585  using Finite_Difference_Arrays::weights;
1586 
1587  Real one(1.0);
1588 
1589  Real tol = std::sqrt(ROL_EPSILON<Real>());
1590 
1591  int numSteps = steps.size();
1592  int numVals = 4;
1593  std::vector<Real> tmp(numVals);
1594  std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1595 
1596  // Temporary vectors.
1597  Ptr<Vector<Real>> AJdif = hv.clone();
1598  Ptr<Vector<Real>> AJp = hv.clone();
1599  Ptr<Vector<Real>> AHpv = hv.clone();
1600  Ptr<Vector<Real>> AJnew = hv.clone();
1601  Ptr<Vector<Real>> znew = z.clone();
1602 
1603  // Save the format state of the original outStream.
1604  nullstream oldFormatState;
1605  oldFormatState.copyfmt(outStream);
1606 
1607  // Apply adjoint Jacobian to p.
1608  this->update(u,z,UpdateType::Temp);
1609  this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1610 
1611  // Apply adjoint Hessian at (u,z), in direction v, to p.
1612  this->applyAdjointHessian_21(*AHpv, p, v, u, z, tol);
1613  Real normAHpv = AHpv->norm();
1614 
1615  for (int i=0; i<numSteps; i++) {
1616 
1617  Real eta = steps[i];
1618 
1619  // Apply adjoint Jacobian to p at (u,z+eta*v).
1620  znew->set(z);
1621 
1622  AJdif->set(*AJp);
1623  AJdif->scale(weights[order-1][0]);
1624 
1625  for(int j=0; j<order; ++j) {
1626 
1627  znew->axpy(eta*shifts[order-1][j],v);
1628 
1629  if( weights[order-1][j+1] != 0 ) {
1630  this->update(u,*znew,UpdateType::Temp);
1631  this->applyAdjointJacobian_1(*AJnew, p, u, *znew, tol);
1632  AJdif->axpy(weights[order-1][j+1],*AJnew);
1633  }
1634  }
1635 
1636  AJdif->scale(one/eta);
1637 
1638  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1639  ahpvCheck[i][0] = eta;
1640  ahpvCheck[i][1] = normAHpv;
1641  ahpvCheck[i][2] = AJdif->norm();
1642  AJdif->axpy(-one, *AHpv);
1643  ahpvCheck[i][3] = AJdif->norm();
1644 
1645  if (printToStream) {
1646  std::stringstream hist;
1647  if (i==0) {
1648  hist << std::right
1649  << std::setw(20) << "Step size"
1650  << std::setw(20) << "norm(adj(H)(u,v))"
1651  << std::setw(20) << "norm(FD approx)"
1652  << std::setw(20) << "norm(abs error)"
1653  << "\n"
1654  << std::setw(20) << "---------"
1655  << std::setw(20) << "-----------------"
1656  << std::setw(20) << "---------------"
1657  << std::setw(20) << "---------------"
1658  << "\n";
1659  }
1660  hist << std::scientific << std::setprecision(11) << std::right
1661  << std::setw(20) << ahpvCheck[i][0]
1662  << std::setw(20) << ahpvCheck[i][1]
1663  << std::setw(20) << ahpvCheck[i][2]
1664  << std::setw(20) << ahpvCheck[i][3]
1665  << "\n";
1666  outStream << hist.str();
1667  }
1668 
1669  }
1670 
1671  // Reset format state of outStream.
1672  outStream.copyfmt(oldFormatState);
1673 
1674  return ahpvCheck;
1675  } // checkApplyAdjointHessian_21
1676 
1680  std::vector<std::vector<Real>> checkApplyAdjointHessian_12(const Vector<Real> &u,
1681  const Vector<Real> &z,
1682  const Vector<Real> &p,
1683  const Vector<Real> &v,
1684  const Vector<Real> &hv,
1685  const bool printToStream = true,
1686  std::ostream & outStream = std::cout,
1687  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1688  const int order = 1 ) {
1689  std::vector<Real> steps(numSteps);
1690  for(int i=0;i<numSteps;++i) {
1691  steps[i] = pow(10,-i);
1692  }
1693 
1694  return checkApplyAdjointHessian_12(u,z,p,v,hv,steps,printToStream,outStream,order);
1695 
1696  }
1697 
1698 
1699  std::vector<std::vector<Real>> checkApplyAdjointHessian_12(const Vector<Real> &u,
1700  const Vector<Real> &z,
1701  const Vector<Real> &p,
1702  const Vector<Real> &v,
1703  const Vector<Real> &hv,
1704  const std::vector<Real> &steps,
1705  const bool printToStream = true,
1706  std::ostream & outStream = std::cout,
1707  const int order = 1 ) {
1708  using Finite_Difference_Arrays::shifts;
1709  using Finite_Difference_Arrays::weights;
1710 
1711  Real one(1.0);
1712 
1713  Real tol = std::sqrt(ROL_EPSILON<Real>());
1714 
1715  int numSteps = steps.size();
1716  int numVals = 4;
1717  std::vector<Real> tmp(numVals);
1718  std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1719 
1720  // Temporary vectors.
1721  Ptr<Vector<Real>> AJdif = hv.clone();
1722  Ptr<Vector<Real>> AJp = hv.clone();
1723  Ptr<Vector<Real>> AHpv = hv.clone();
1724  Ptr<Vector<Real>> AJnew = hv.clone();
1725  Ptr<Vector<Real>> unew = u.clone();
1726 
1727  // Save the format state of the original outStream.
1728  nullstream oldFormatState;
1729  oldFormatState.copyfmt(outStream);
1730 
1731  // Apply adjoint Jacobian to p.
1732  this->update(u,z,UpdateType::Temp);
1733  this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1734 
1735  // Apply adjoint Hessian at (u,z), in direction v, to p.
1736  this->applyAdjointHessian_12(*AHpv, p, v, u, z, tol);
1737  Real normAHpv = AHpv->norm();
1738 
1739  for (int i=0; i<numSteps; i++) {
1740 
1741  Real eta = steps[i];
1742 
1743  // Apply adjoint Jacobian to p at (u+eta*v,z).
1744  unew->set(u);
1745 
1746  AJdif->set(*AJp);
1747  AJdif->scale(weights[order-1][0]);
1748 
1749  for(int j=0; j<order; ++j) {
1750 
1751  unew->axpy(eta*shifts[order-1][j],v);
1752 
1753  if( weights[order-1][j+1] != 0 ) {
1754  this->update(*unew,z,UpdateType::Temp);
1755  this->applyAdjointJacobian_2(*AJnew, p, *unew, z, tol);
1756  AJdif->axpy(weights[order-1][j+1],*AJnew);
1757  }
1758  }
1759 
1760  AJdif->scale(one/eta);
1761 
1762  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1763  ahpvCheck[i][0] = eta;
1764  ahpvCheck[i][1] = normAHpv;
1765  ahpvCheck[i][2] = AJdif->norm();
1766  AJdif->axpy(-one, *AHpv);
1767  ahpvCheck[i][3] = AJdif->norm();
1768 
1769  if (printToStream) {
1770  std::stringstream hist;
1771  if (i==0) {
1772  hist << std::right
1773  << std::setw(20) << "Step size"
1774  << std::setw(20) << "norm(adj(H)(u,v))"
1775  << std::setw(20) << "norm(FD approx)"
1776  << std::setw(20) << "norm(abs error)"
1777  << "\n"
1778  << std::setw(20) << "---------"
1779  << std::setw(20) << "-----------------"
1780  << std::setw(20) << "---------------"
1781  << std::setw(20) << "---------------"
1782  << "\n";
1783  }
1784  hist << std::scientific << std::setprecision(11) << std::right
1785  << std::setw(20) << ahpvCheck[i][0]
1786  << std::setw(20) << ahpvCheck[i][1]
1787  << std::setw(20) << ahpvCheck[i][2]
1788  << std::setw(20) << ahpvCheck[i][3]
1789  << "\n";
1790  outStream << hist.str();
1791  }
1792 
1793  }
1794 
1795  // Reset format state of outStream.
1796  outStream.copyfmt(oldFormatState);
1797 
1798  return ahpvCheck;
1799  } // checkApplyAdjointHessian_12
1800 
1801  std::vector<std::vector<Real>> checkApplyAdjointHessian_22(const Vector<Real> &u,
1802  const Vector<Real> &z,
1803  const Vector<Real> &p,
1804  const Vector<Real> &v,
1805  const Vector<Real> &hv,
1806  const bool printToStream = true,
1807  std::ostream & outStream = std::cout,
1808  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1809  const int order = 1 ) {
1810  std::vector<Real> steps(numSteps);
1811  for(int i=0;i<numSteps;++i) {
1812  steps[i] = pow(10,-i);
1813  }
1814 
1815  return checkApplyAdjointHessian_22(u,z,p,v,hv,steps,printToStream,outStream,order);
1816 
1817  }
1818 
1819  std::vector<std::vector<Real>> checkApplyAdjointHessian_22(const Vector<Real> &u,
1820  const Vector<Real> &z,
1821  const Vector<Real> &p,
1822  const Vector<Real> &v,
1823  const Vector<Real> &hv,
1824  const std::vector<Real> &steps,
1825  const bool printToStream = true,
1826  std::ostream & outStream = std::cout,
1827  const int order = 1 ) {
1828  using Finite_Difference_Arrays::shifts;
1829  using Finite_Difference_Arrays::weights;
1830 
1831  Real one(1.0);
1832 
1833  Real tol = std::sqrt(ROL_EPSILON<Real>());
1834 
1835  int numSteps = steps.size();
1836  int numVals = 4;
1837  std::vector<Real> tmp(numVals);
1838  std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1839 
1840  // Temporary vectors.
1841  Ptr<Vector<Real>> AJdif = hv.clone();
1842  Ptr<Vector<Real>> AJp = hv.clone();
1843  Ptr<Vector<Real>> AHpv = hv.clone();
1844  Ptr<Vector<Real>> AJnew = hv.clone();
1845  Ptr<Vector<Real>> znew = z.clone();
1846 
1847  // Save the format state of the original outStream.
1848  nullstream oldFormatState;
1849  oldFormatState.copyfmt(outStream);
1850 
1851  // Apply adjoint Jacobian to p.
1852  this->update(u,z,UpdateType::Temp);
1853  this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1854 
1855  // Apply adjoint Hessian at (u,z), in direction v, to p.
1856  this->applyAdjointHessian_22(*AHpv, p, v, u, z, tol);
1857  Real normAHpv = AHpv->norm();
1858 
1859  for (int i=0; i<numSteps; i++) {
1860 
1861  Real eta = steps[i];
1862 
1863  // Apply adjoint Jacobian to p at (u,z+eta*v).
1864  znew->set(z);
1865 
1866  AJdif->set(*AJp);
1867  AJdif->scale(weights[order-1][0]);
1868 
1869  for(int j=0; j<order; ++j) {
1870 
1871  znew->axpy(eta*shifts[order-1][j],v);
1872 
1873  if( weights[order-1][j+1] != 0 ) {
1874  this->update(u,*znew,UpdateType::Temp);
1875  this->applyAdjointJacobian_2(*AJnew, p, u, *znew, tol);
1876  AJdif->axpy(weights[order-1][j+1],*AJnew);
1877  }
1878  }
1879 
1880  AJdif->scale(one/eta);
1881 
1882  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1883  ahpvCheck[i][0] = eta;
1884  ahpvCheck[i][1] = normAHpv;
1885  ahpvCheck[i][2] = AJdif->norm();
1886  AJdif->axpy(-one, *AHpv);
1887  ahpvCheck[i][3] = AJdif->norm();
1888 
1889  if (printToStream) {
1890  std::stringstream hist;
1891  if (i==0) {
1892  hist << std::right
1893  << std::setw(20) << "Step size"
1894  << std::setw(20) << "norm(adj(H)(u,v))"
1895  << std::setw(20) << "norm(FD approx)"
1896  << std::setw(20) << "norm(abs error)"
1897  << "\n"
1898  << std::setw(20) << "---------"
1899  << std::setw(20) << "-----------------"
1900  << std::setw(20) << "---------------"
1901  << std::setw(20) << "---------------"
1902  << "\n";
1903  }
1904  hist << std::scientific << std::setprecision(11) << std::right
1905  << std::setw(20) << ahpvCheck[i][0]
1906  << std::setw(20) << ahpvCheck[i][1]
1907  << std::setw(20) << ahpvCheck[i][2]
1908  << std::setw(20) << ahpvCheck[i][3]
1909  << "\n";
1910  outStream << hist.str();
1911  }
1912 
1913  }
1914 
1915  // Reset format state of outStream.
1916  outStream.copyfmt(oldFormatState);
1917 
1918  return ahpvCheck;
1919  } // checkApplyAdjointHessian_22
1920 
1921 }; // class Constraint_SimOpt
1922 
1923 } // namespace ROL
1924 
1925 #endif
ROL::Constraint_SimOpt::zero_
bool zero_
Definition: ROL_Constraint_SimOpt.hpp:101
ROL::Constraint_SimOpt::DEFAULT_stol_
const Real DEFAULT_stol_
Definition: ROL_Constraint_SimOpt.hpp:83
ROL::Constraint_SimOpt::value
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
Definition: ROL_Constraint_SimOpt.hpp:880
ROL::Vector::dual
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
ROL::Constraint::applyPreconditioner
virtual void applyPreconditioner(Vector< Real > &pv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g, Real &tol)
Apply a constraint preconditioner at , , to vector . Ideally, this preconditioner satisfies the follo...
Definition: ROL_Constraint.hpp:249
ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Constraint_SimOpt::applyJacobian_2
virtual void applyJacobian_2(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
Definition: ROL_Constraint_SimOpt.hpp:374
ROL::Constraint_SimOpt::applyAdjointHessian
virtual void applyAdjointHessian(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ...
Definition: ROL_Constraint_SimOpt.hpp:922
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Constraint_SimOpt::DEFAULT_decr_
const Real DEFAULT_decr_
Definition: ROL_Constraint_SimOpt.hpp:85
ROL::Vector_SimOpt::get_2
ROL::Ptr< const Vector< Real > > get_2() const
Definition: ROL_Vector_SimOpt.hpp:144
ROL::Vector::dimension
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:162
ROL::Vector::apply
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: ROL_Vector.hpp:204
ROL::Constraint_SimOpt::checkApplyAdjointHessian_12
std::vector< std::vector< Real > > checkApplyAdjointHessian_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
, , , ,
Definition: ROL_Constraint_SimOpt.hpp:1680
ROL::Constraint_SimOpt::Constraint_SimOpt
Constraint_SimOpt()
Definition: ROL_Constraint_SimOpt.hpp:108
ROL::Vector::basis
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:148
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Finite_Difference_Arrays::weights
const double weights[4][5]
Definition: ROL_Types.hpp:834
ROL::Constraint_SimOpt::solve_update
virtual void solve_update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-1)
Update SimOpt constraint during solve (disconnected from optimization updates).
Definition: ROL_Constraint_SimOpt.hpp:160
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
ROL_Objective_FSsolver.hpp
ROL::Constraint_SimOpt::applyJacobian_1
virtual void applyJacobian_1(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
Definition: ROL_Constraint_SimOpt.hpp:331
ROL::Constraint_SimOpt::checkApplyAdjointHessian_12
std::vector< std::vector< Real > > checkApplyAdjointHessian_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1699
ROL::Constraint_SimOpt::checkSolve
virtual Real checkSolve(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &c, const bool printToStream=true, std::ostream &outStream=std::cout)
Definition: ROL_Constraint_SimOpt.hpp:951
ROL::Constraint_SimOpt::update
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update constraint function.
Definition: ROL_Constraint_SimOpt.hpp:874
ROL::Vector_SimOpt
Defines the linear algebra or vector space interface for simulation-based optimization.
Definition: ROL_Vector_SimOpt.hpp:23
ROL::Constraint_SimOpt::update_1
virtual void update_1(const Vector< Real > &u, UpdateType type, int iter=-1)
Definition: ROL_Constraint_SimOpt.hpp:144
ROL::Constraint_SimOpt::applyAdjointHessian_11
virtual void applyAdjointHessian_11(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
Definition: ROL_Constraint_SimOpt.hpp:607
ROL::Constraint_SimOpt::checkAdjointConsistencyJacobian_2
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualw, const Vector< Real > &dualv, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the secondary interface, for use with dual spaces where the user does not define the dual() operation.
Definition: ROL_Constraint_SimOpt.hpp:1086
ROL::Constraint_SimOpt::maxit_
int maxit_
Definition: ROL_Constraint_SimOpt.hpp:99
ROL::Constraint_SimOpt::update
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable is ...
Definition: ROL_Constraint_SimOpt.hpp:869
ROL::Constraint_SimOpt::solverType_
int solverType_
Definition: ROL_Constraint_SimOpt.hpp:102
ROL::Constraint_SimOpt::update_1
virtual void update_1(const Vector< Real > &u, bool flag=true, int iter=-1)
Update constraint functions with respect to Sim variable. x is the optimization variable, flag = true if optimization variable is changed, iter is the outer algorithm iterations count.
Definition: ROL_Constraint_SimOpt.hpp:143
ROL::Constraint_SimOpt::applyAdjointHessian_22
virtual void applyAdjointHessian_22(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
Definition: ROL_Constraint_SimOpt.hpp:741
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL::Finite_Difference_Arrays::shifts
const int shifts[4][4]
Definition: ROL_Types.hpp:827
ROL::Vector_SimOpt::set_1
void set_1(const Vector< Real > &vec)
Definition: ROL_Vector_SimOpt.hpp:156
ROL::Constraint_SimOpt::checkApplyJacobian_1
std::vector< std::vector< Real > > checkApplyJacobian_1(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1180
ROL::Constraint_SimOpt::applyJacobian
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
Definition: ROL_Constraint_SimOpt.hpp:889
ROL::Constraint_SimOpt::checkAdjointConsistencyJacobian_2
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
Definition: ROL_Constraint_SimOpt.hpp:1062
ROL::Constraint_SimOpt::update
virtual void update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-1)
Definition: ROL_Constraint_SimOpt.hpp:133
ROL::Constraint_SimOpt::applyPreconditioner
virtual void applyPreconditioner(Vector< Real > &pv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g, Real &tol)
Apply a constraint preconditioner at , , to vector . In general, this preconditioner satisfies the fo...
Definition: ROL_Constraint_SimOpt.hpp:834
ROL::Constraint_SimOpt::rtol_
Real rtol_
Definition: ROL_Constraint_SimOpt.hpp:95
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL::Constraint_SimOpt::factor_
Real factor_
Definition: ROL_Constraint_SimOpt.hpp:97
ROL::Constraint_SimOpt::update
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable is ...
Definition: ROL_Constraint_SimOpt.hpp:129
ROL_SimConstraint.hpp
ROL::details::makeStreamPtr
Ptr< ostream > makeStreamPtr(ostream &os, bool noSuppressOutput=true)
Definition: ROL_Stream.hpp:41
ROL::Constraint_SimOpt::checkApplyJacobian_2
std::vector< std::vector< Real > > checkApplyJacobian_2(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1305
ROL::Constraint_SimOpt::jv_
Ptr< Vector< Real > > jv_
Definition: ROL_Constraint_SimOpt.hpp:78
ROL_NonlinearLeastSquaresObjective.hpp
ROL::Constraint_SimOpt::checkInverseJacobian_1
virtual Real checkInverseJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Definition: ROL_Constraint_SimOpt.hpp:1118
ROL::Constraint_SimOpt::applyInverseJacobian_1
virtual void applyInverseJacobian_1(Vector< Real > &ijv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the inverse partial constraint Jacobian at , , to the vector .
Definition: ROL_Constraint_SimOpt.hpp:416
ROL::Constraint_SimOpt::DEFAULT_solverType_
const int DEFAULT_solverType_
Definition: ROL_Constraint_SimOpt.hpp:89
ROL::Constraint_SimOpt::applyAdjointHessian_12
virtual void applyAdjointHessian_12(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
Definition: ROL_Constraint_SimOpt.hpp:652
ROL::Constraint_SimOpt::solveAugmentedSystem
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system where , , , , is an identity operator, and is a zero operator.
Definition: ROL_Constraint_SimOpt.hpp:806
ROL::Constraint_SimOpt::applyAdjointJacobian_1
virtual void applyAdjointJacobian_1(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
Definition: ROL_Constraint_SimOpt.hpp:440
ROL::Constraint_SimOpt::value
virtual void value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol)=0
Evaluate the constraint operator at .
ROL::Constraint_SimOpt::DEFAULT_atol_
const Real DEFAULT_atol_
Definition: ROL_Constraint_SimOpt.hpp:81
ROL::Constraint_SimOpt::applyAdjointJacobian_2
virtual void applyAdjointJacobian_2(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualv, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the secondary interfa...
Definition: ROL_Constraint_SimOpt.hpp:537
ROL_Constraint.hpp
ROL::Constraint_SimOpt::firstSolve_
bool firstSolve_
Definition: ROL_Constraint_SimOpt.hpp:105
ROL::Constraint_SimOpt::checkApplyAdjointHessian_11
std::vector< std::vector< Real > > checkApplyAdjointHessian_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1449
ROL::Constraint_SimOpt::checkApplyJacobian_2
std::vector< std::vector< Real > > checkApplyJacobian_2(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1324
ROL::Constraint_SimOpt::DEFAULT_rtol_
const Real DEFAULT_rtol_
Definition: ROL_Constraint_SimOpt.hpp:82
ROL::Constraint_SimOpt::stol_
Real stol_
Definition: ROL_Constraint_SimOpt.hpp:96
ROL::Constraint_SimOpt::applyAdjointJacobian_1
virtual void applyAdjointJacobian_1(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualv, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the secondary int...
Definition: ROL_Constraint_SimOpt.hpp:466
ROL::Constraint_SimOpt::applyAdjointJacobian_2
virtual void applyAdjointJacobian_2(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
Definition: ROL_Constraint_SimOpt.hpp:511
ROL_TypeE_AugmentedLagrangianAlgorithm.hpp
ROL::Constraint_SimOpt::DEFAULT_factor_
const Real DEFAULT_factor_
Definition: ROL_Constraint_SimOpt.hpp:84
ROL::Constraint_SimOpt::decr_
Real decr_
Definition: ROL_Constraint_SimOpt.hpp:98
ROL_NUM_CHECKDERIV_STEPS
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
Definition: ROL_Types.hpp:40
ROL::Constraint_SimOpt::applyAdjointHessian_21
virtual void applyAdjointHessian_21(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
Definition: ROL_Constraint_SimOpt.hpp:697
ROL::Constraint_SimOpt::print_
bool print_
Definition: ROL_Constraint_SimOpt.hpp:100
ROL::Constraint_SimOpt::DEFAULT_maxit_
const int DEFAULT_maxit_
Definition: ROL_Constraint_SimOpt.hpp:86
ROL::details::nullstream
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
ROL::Constraint_SimOpt::solve
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
Definition: ROL_Constraint_SimOpt.hpp:191
ROL::Constraint::solveAugmentedSystem
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system where , , , , is an identity or Riesz operator...
Definition: ROL_ConstraintDef.hpp:161
ROL::Constraint_SimOpt::update_2
virtual void update_2(const Vector< Real > &z, UpdateType type, int iter=-1)
Definition: ROL_Constraint_SimOpt.hpp:152
ROL::Constraint_SimOpt::checkApplyJacobian_1
std::vector< std::vector< Real > > checkApplyJacobian_1(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1199
ROL::Constraint_SimOpt::checkAdjointConsistencyJacobian_1
virtual Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualw, const Vector< Real > &dualv, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the secondary interface, for use with dual spaces where the user does not define the dual() operation.
Definition: ROL_Constraint_SimOpt.hpp:1017
ROL::Constraint_SimOpt::setSolveParameters
virtual void setSolveParameters(ParameterList &parlist)
Set solve parameters.
Definition: ROL_Constraint_SimOpt.hpp:303
ROL::Constraint_SimOpt::checkAdjointConsistencyJacobian_1
virtual Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
Definition: ROL_Constraint_SimOpt.hpp:992
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
ROL::Constraint_SimOpt::update_2
virtual void update_2(const Vector< Real > &z, bool flag=true, int iter=-1)
Update constraint functions with respect to Opt variable. x is the optimization variable, flag = true if optimization variable is changed, iter is the outer algorithm iterations count.
Definition: ROL_Constraint_SimOpt.hpp:151
ROL::Constraint_SimOpt::unew_
Ptr< Vector< Real > > unew_
Definition: ROL_Constraint_SimOpt.hpp:77
ROL::Constraint_SimOpt::checkApplyAdjointHessian_11
std::vector< std::vector< Real > > checkApplyAdjointHessian_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1431
ROL::Constraint_SimOpt::atol_
Real atol_
Definition: ROL_Constraint_SimOpt.hpp:94
ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:57
ROL::Constraint_SimOpt
Defines the constraint operator interface for simulation-based optimization.
Definition: ROL_Constraint_SimOpt.hpp:21
ROL::Constraint_SimOpt::applyAdjointJacobian
virtual void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
Definition: ROL_Constraint_SimOpt.hpp:905
ROL::Constraint_SimOpt::checkApplyAdjointHessian_22
std::vector< std::vector< Real > > checkApplyAdjointHessian_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1801
ROL::Constraint_SimOpt::checkApplyAdjointHessian_21
std::vector< std::vector< Real > > checkApplyAdjointHessian_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
, , , ,
Definition: ROL_Constraint_SimOpt.hpp:1554
ROL::Constraint_SimOpt::applyInverseAdjointJacobian_1
virtual void applyInverseAdjointJacobian_1(Vector< Real > &iajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector ...
Definition: ROL_Constraint_SimOpt.hpp:581
ROL::Vector_SimOpt::set_2
void set_2(const Vector< Real > &vec)
Definition: ROL_Vector_SimOpt.hpp:160
ROL::Constraint_SimOpt::checkApplyAdjointHessian_21
std::vector< std::vector< Real > > checkApplyAdjointHessian_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
, , , ,
Definition: ROL_Constraint_SimOpt.hpp:1575
ROL::Constraint_SimOpt::DEFAULT_zero_
const bool DEFAULT_zero_
Definition: ROL_Constraint_SimOpt.hpp:88
ROL::Constraint
Defines the general constraint operator interface.
Definition: ROL_Constraint.hpp:52
ROL_TypeU_TrustRegionAlgorithm.hpp
ROL::Constraint_SimOpt::DEFAULT_print_
const bool DEFAULT_print_
Definition: ROL_Constraint_SimOpt.hpp:87
ROL::Constraint_SimOpt::checkApplyAdjointHessian_22
std::vector< std::vector< Real > > checkApplyAdjointHessian_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition: ROL_Constraint_SimOpt.hpp:1819
ROL::Constraint_SimOpt::checkInverseAdjointJacobian_1
virtual Real checkInverseAdjointJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Definition: ROL_Constraint_SimOpt.hpp:1148
ROL::Vector_SimOpt::dual
const Vector< Real > & dual(void) const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector_SimOpt.hpp:73
ROL::Vector_SimOpt::get_1
ROL::Ptr< const Vector< Real > > get_1() const
Definition: ROL_Vector_SimOpt.hpp:140
ROL_Vector_SimOpt.hpp
Generated by   doxygen 1.8.5