ROL
ROL_Constraint_SimOpt.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_CONSTRAINT_SIMOPT_H
45 #define ROL_CONSTRAINT_SIMOPT_H
46 
47 #include "ROL_Constraint.hpp"
48 #include "ROL_Vector_SimOpt.hpp"
49 #include "ROL_Types.hpp"
50 #include <iostream>
51 
52 #include "ROL_NonlinearLeastSquaresObjective_SimOpt.hpp"
53 #include "ROL_Constraint_State.hpp"
54 #include "ROL_Objective_FSsolver.hpp"
55 #include "ROL_Algorithm.hpp"
56 #include "ROL_TrustRegionStep.hpp"
57 #include "ROL_CompositeStep.hpp"
58 #include "ROL_ConstraintStatusTest.hpp"
59 
100 namespace ROL {
101 
102 template <class Real>
103 class Constraint_SimOpt : public Constraint<Real> {
104 private:
105  // Additional vector storage for solve
106  Ptr<Vector<Real>> unew_;
107  Ptr<Vector<Real>> jv_;
108 
109  // Default parameters for solve (backtracking Newton)
110  const Real DEFAULT_atol_;
111  const Real DEFAULT_rtol_;
112  const Real DEFAULT_stol_;
113  const Real DEFAULT_factor_;
114  const Real DEFAULT_decr_;
115  const int DEFAULT_maxit_;
116  const bool DEFAULT_print_;
117  const bool DEFAULT_zero_;
118  const int DEFAULT_solverType_;
119 
120  // User-set parameters for solve (backtracking Newton)
121 
122 protected:
123  Real atol_;
124  Real rtol_;
125  Real stol_;
126  Real factor_;
127  Real decr_;
128  int maxit_;
129  bool print_;
130  bool zero_;
131  int solverType_;
132 
133  // Flag to initialize vector storage in solve
134  bool firstSolve_;
135 
136 public:
137  Constraint_SimOpt()
138  : Constraint<Real>(),
139  unew_(nullPtr), jv_(nullPtr),
140  DEFAULT_atol_(1.e-4*std::sqrt(ROL_EPSILON<Real>())),
141  DEFAULT_rtol_(1.e0),
142  DEFAULT_stol_(std::sqrt(ROL_EPSILON<Real>())),
143  DEFAULT_factor_(0.5),
144  DEFAULT_decr_(1.e-4),
145  DEFAULT_maxit_(500),
146  DEFAULT_print_(false),
147  DEFAULT_zero_(false),
148  DEFAULT_solverType_(0),
149  atol_(DEFAULT_atol_), rtol_(DEFAULT_rtol_), stol_(DEFAULT_stol_), factor_(DEFAULT_factor_),
150  decr_(DEFAULT_decr_), maxit_(DEFAULT_maxit_), print_(DEFAULT_print_), zero_(DEFAULT_zero_),
151  solverType_(DEFAULT_solverType_), firstSolve_(true) {}
152 
158  virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) {
159  update_1(u,flag,iter);
160  update_2(z,flag,iter);
161  }
162 
168  virtual void update_1( const Vector<Real> &u, bool flag = true, int iter = -1 ) {}
169 
175  virtual void update_2( const Vector<Real> &z, bool flag = true, int iter = -1 ) {}
176 
177 
191  virtual void value(Vector<Real> &c,
192  const Vector<Real> &u,
193  const Vector<Real> &z,
194  Real &tol) = 0;
195 
207  virtual void solve(Vector<Real> &c,
208  Vector<Real> &u,
209  const Vector<Real> &z,
210  Real &tol) {
211  if ( zero_ ) {
212  u.zero();
213  }
214  update(u,z,true);
215  value(c,u,z,tol);
216  Real cnorm = c.norm();
217  const Real ctol = std::min(atol_, rtol_*cnorm);
218  if (solverType_==0 || solverType_==3 || solverType_==4) {
219  if ( firstSolve_ ) {
220  unew_ = u.clone();
221  jv_ = u.clone();
222  firstSolve_ = false;
223  }
224  Real alpha(1), tmp(0);
225  int cnt = 0;
226  if ( print_ ) {
227  std::cout << "\n Default Constraint_SimOpt::solve\n";
228  std::cout << " ";
229  std::cout << std::setw(6) << std::left << "iter";
230  std::cout << std::setw(15) << std::left << "rnorm";
231  std::cout << std::setw(15) << std::left << "alpha";
232  std::cout << "\n";
233  }
234  for (cnt = 0; cnt < maxit_; ++cnt) {
235  // Compute Newton step
236  applyInverseJacobian_1(*jv_,c,u,z,tol);
237  unew_->set(u);
238  unew_->axpy(-alpha, *jv_);
239  update_1(*unew_);
240  //update(*unew_,z);
241  value(c,*unew_,z,tol);
242  tmp = c.norm();
243  // Perform backtracking line search
244  while ( tmp > (1.0-decr_*alpha)*cnorm &&
245  alpha > stol_ ) {
246  alpha *= factor_;
247  unew_->set(u);
248  unew_->axpy(-alpha,*jv_);
249  update_1(*unew_);
250  //update(*unew_,z);
251  value(c,*unew_,z,tol);
252  tmp = c.norm();
253  }
254  if ( print_ ) {
255  std::cout << " ";
256  std::cout << std::setw(6) << std::left << cnt;
257  std::cout << std::scientific << std::setprecision(6);
258  std::cout << std::setw(15) << std::left << tmp;
259  std::cout << std::scientific << std::setprecision(6);
260  std::cout << std::setw(15) << std::left << alpha;
261  std::cout << "\n";
262  }
263  // Update iterate
264  cnorm = tmp;
265  u.set(*unew_);
266  if (cnorm <= ctol) { // = covers the case of identically zero residual
267  break;
268  }
269  update(u,z,true);
270  alpha = 1.0;
271  }
272  }
273  if (solverType_==1 || (solverType_==3 && cnorm > ctol)) {
274  Ptr<Constraint_SimOpt<Real>> con = makePtrFromRef(*this);
275  Ptr<Objective<Real>> obj = makePtr<NonlinearLeastSquaresObjective_SimOpt<Real>>(con,u,z,c,true);
276  ParameterList parlist;
277  parlist.sublist("Status Test").set("Gradient Tolerance",ctol);
278  parlist.sublist("Status Test").set("Step Tolerance",stol_);
279  parlist.sublist("Status Test").set("Iteration Limit",maxit_);
280  parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG");
281  parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100);
282  Ptr<Step<Real>> step = makePtr<TrustRegionStep<Real>>(parlist);
283  Ptr<StatusTest<Real>> status = makePtr<StatusTest<Real>>(parlist);
284  Ptr<Algorithm<Real>> algo = makePtr<Algorithm<Real>>(step,status,false);
285  algo->run(u,*obj,print_);
286  value(c,u,z,tol);
287  }
288  if (solverType_==2 || (solverType_==4 && cnorm > ctol)) {
289  Ptr<Constraint_SimOpt<Real>> con = makePtrFromRef(*this);
290  Ptr<const Vector<Real>> zVec = makePtrFromRef(z);
291  Ptr<Constraint<Real>> conU
292  = makePtr<Constraint_State<Real>>(con,zVec);
293  Ptr<Objective<Real>> objU
294  = makePtr<Objective_FSsolver<Real>>();
295  ParameterList parlist;
296  parlist.sublist("Status Test").set("Constraint Tolerance",ctol);
297  parlist.sublist("Status Test").set("Step Tolerance",stol_);
298  parlist.sublist("Status Test").set("Iteration Limit",maxit_);
299  Ptr<Step<Real>> step = makePtr<CompositeStep<Real>>(parlist);
300  Ptr<StatusTest<Real>> status = makePtr<ConstraintStatusTest<Real>>(parlist);
301  Ptr<Algorithm<Real>> algo = makePtr<Algorithm<Real>>(step,status,false);
302  Ptr<Vector<Real>> l = c.dual().clone();
303  algo->run(u,*l,*objU,*conU,print_);
304  value(c,u,z,tol);
305  }
306  if (solverType_ > 4 || solverType_ < 0) {
307  ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
308  ">>> ERROR (ROL:Constraint_SimOpt:solve): Invalid solver type!");
309  }
310  }
311 
332  virtual void setSolveParameters(ParameterList &parlist) {
333  ParameterList & list = parlist.sublist("SimOpt").sublist("Solve");
334  atol_ = list.get("Absolute Residual Tolerance", DEFAULT_atol_);
335  rtol_ = list.get("Relative Residual Tolerance", DEFAULT_rtol_);
336  maxit_ = list.get("Iteration Limit", DEFAULT_maxit_);
337  decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_);
338  stol_ = list.get("Step Tolerance", DEFAULT_stol_);
339  factor_ = list.get("Backtracking Factor", DEFAULT_factor_);
340  print_ = list.get("Output Iteration History", DEFAULT_print_);
341  zero_ = list.get("Zero Initial Guess", DEFAULT_zero_);
342  solverType_ = list.get("Solver Type", DEFAULT_solverType_);
343  }
344 
360  virtual void applyJacobian_1(Vector<Real> &jv,
361  const Vector<Real> &v,
362  const Vector<Real> &u,
363  const Vector<Real> &z,
364  Real &tol) {
365  Real ctol = std::sqrt(ROL_EPSILON<Real>());
366  // Compute step length
367  Real h = tol;
368  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
369  h = std::max(1.0,u.norm()/v.norm())*tol;
370  }
371  // Update state vector to u + hv
372  Ptr<Vector<Real>> unew = u.clone();
373  unew->set(u);
374  unew->axpy(h,v);
375  // Compute new constraint value
376  update(*unew,z);
377  value(jv,*unew,z,ctol);
378  // Compute current constraint value
379  Ptr<Vector<Real>> cold = jv.clone();
380  update(u,z);
381  value(*cold,u,z,ctol);
382  // Compute Newton quotient
383  jv.axpy(-1.0,*cold);
384  jv.scale(1.0/h);
385  }
386 
387 
403  virtual void applyJacobian_2(Vector<Real> &jv,
404  const Vector<Real> &v,
405  const Vector<Real> &u,
406  const Vector<Real> &z,
407  Real &tol) {
408  Real ctol = std::sqrt(ROL_EPSILON<Real>());
409  // Compute step length
410  Real h = tol;
411  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
412  h = std::max(1.0,u.norm()/v.norm())*tol;
413  }
414  // Update state vector to u + hv
415  Ptr<Vector<Real>> znew = z.clone();
416  znew->set(z);
417  znew->axpy(h,v);
418  // Compute new constraint value
419  update(u,*znew);
420  value(jv,u,*znew,ctol);
421  // Compute current constraint value
422  Ptr<Vector<Real>> cold = jv.clone();
423  update(u,z);
424  value(*cold,u,z,ctol);
425  // Compute Newton quotient
426  jv.axpy(-1.0,*cold);
427  jv.scale(1.0/h);
428  }
429 
445  virtual void applyInverseJacobian_1(Vector<Real> &ijv,
446  const Vector<Real> &v,
447  const Vector<Real> &u,
448  const Vector<Real> &z,
449  Real &tol) {
450  ROL_TEST_FOR_EXCEPTION(true, std::logic_error,
451  "The method applyInverseJacobian_1 is used but not implemented!\n");
452  }
453 
469  virtual void applyAdjointJacobian_1(Vector<Real> &ajv,
470  const Vector<Real> &v,
471  const Vector<Real> &u,
472  const Vector<Real> &z,
473  Real &tol) {
474  applyAdjointJacobian_1(ajv, v, u, z, v.dual(), tol);
475  }
476 
477 
495  virtual void applyAdjointJacobian_1(Vector<Real> &ajv,
496  const Vector<Real> &v,
497  const Vector<Real> &u,
498  const Vector<Real> &z,
499  const Vector<Real> &dualv,
500  Real &tol) {
501  Real ctol = std::sqrt(ROL_EPSILON<Real>());
502  Real h = tol;
503  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
504  h = std::max(1.0,u.norm()/v.norm())*tol;
505  }
506  Ptr<Vector<Real>> cold = dualv.clone();
507  Ptr<Vector<Real>> cnew = dualv.clone();
508  update(u,z);
509  value(*cold,u,z,ctol);
510  Ptr<Vector<Real>> unew = u.clone();
511  ajv.zero();
512  for (int i = 0; i < u.dimension(); i++) {
513  unew->set(u);
514  unew->axpy(h,*(u.basis(i)));
515  update(*unew,z);
516  value(*cnew,*unew,z,ctol);
517  cnew->axpy(-1.0,*cold);
518  cnew->scale(1.0/h);
519  ajv.axpy(cnew->dot(v),*((u.dual()).basis(i)));
520  }
521  update(u,z);
522  }
523 
524 
540  virtual void applyAdjointJacobian_2(Vector<Real> &ajv,
541  const Vector<Real> &v,
542  const Vector<Real> &u,
543  const Vector<Real> &z,
544  Real &tol) {
545  applyAdjointJacobian_2(ajv, v, u, z, v.dual(), tol);
546  }
547 
548 
566  virtual void applyAdjointJacobian_2(Vector<Real> &ajv,
567  const Vector<Real> &v,
568  const Vector<Real> &u,
569  const Vector<Real> &z,
570  const Vector<Real> &dualv,
571  Real &tol) {
572  Real ctol = std::sqrt(ROL_EPSILON<Real>());
573  Real h = tol;
574  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
575  h = std::max(1.0,u.norm()/v.norm())*tol;
576  }
577  Ptr<Vector<Real>> cold = dualv.clone();
578  Ptr<Vector<Real>> cnew = dualv.clone();
579  update(u,z);
580  value(*cold,u,z,ctol);
581  Ptr<Vector<Real>> znew = z.clone();
582  ajv.zero();
583  for (int i = 0; i < z.dimension(); i++) {
584  znew->set(z);
585  znew->axpy(h,*(z.basis(i)));
586  update(u,*znew);
587  value(*cnew,u,*znew,ctol);
588  cnew->axpy(-1.0,*cold);
589  cnew->scale(1.0/h);
590  ajv.axpy(cnew->dot(v),*((z.dual()).basis(i)));
591  }
592  update(u,z);
593  }
594 
610  virtual void applyInverseAdjointJacobian_1(Vector<Real> &iajv,
611  const Vector<Real> &v,
612  const Vector<Real> &u,
613  const Vector<Real> &z,
614  Real &tol) {
615  ROL_TEST_FOR_EXCEPTION(true, std::logic_error,
616  "The method applyInverseAdjointJacobian_1 is used but not implemented!\n");
617  };
618 
636  virtual void applyAdjointHessian_11(Vector<Real> &ahwv,
637  const Vector<Real> &w,
638  const Vector<Real> &v,
639  const Vector<Real> &u,
640  const Vector<Real> &z,
641  Real &tol) {
642  Real jtol = std::sqrt(ROL_EPSILON<Real>());
643  // Compute step size
644  Real h = tol;
645  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
646  h = std::max(1.0,u.norm()/v.norm())*tol;
647  }
648  // Evaluate Jacobian at new state
649  Ptr<Vector<Real>> unew = u.clone();
650  unew->set(u);
651  unew->axpy(h,v);
652  update(*unew,z);
653  applyAdjointJacobian_1(ahwv,w,*unew,z,jtol);
654  // Evaluate Jacobian at old state
655  Ptr<Vector<Real>> jv = ahwv.clone();
656  update(u,z);
657  applyAdjointJacobian_1(*jv,w,u,z,jtol);
658  // Compute Newton quotient
659  ahwv.axpy(-1.0,*jv);
660  ahwv.scale(1.0/h);
661  }
662 
663 
681  virtual void applyAdjointHessian_12(Vector<Real> &ahwv,
682  const Vector<Real> &w,
683  const Vector<Real> &v,
684  const Vector<Real> &u,
685  const Vector<Real> &z,
686  Real &tol) {
687  Real jtol = std::sqrt(ROL_EPSILON<Real>());
688  // Compute step size
689  Real h = tol;
690  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
691  h = std::max(1.0,u.norm()/v.norm())*tol;
692  }
693  // Evaluate Jacobian at new state
694  Ptr<Vector<Real>> unew = u.clone();
695  unew->set(u);
696  unew->axpy(h,v);
697  update(*unew,z);
698  applyAdjointJacobian_2(ahwv,w,*unew,z,jtol);
699  // Evaluate Jacobian at old state
700  Ptr<Vector<Real>> jv = ahwv.clone();
701  update(u,z);
702  applyAdjointJacobian_2(*jv,w,u,z,jtol);
703  // Compute Newton quotient
704  ahwv.axpy(-1.0,*jv);
705  ahwv.scale(1.0/h);
706  }
707 
708 
726  virtual void applyAdjointHessian_21(Vector<Real> &ahwv,
727  const Vector<Real> &w,
728  const Vector<Real> &v,
729  const Vector<Real> &u,
730  const Vector<Real> &z,
731  Real &tol) {
732  Real jtol = std::sqrt(ROL_EPSILON<Real>());
733  // Compute step size
734  Real h = tol;
735  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
736  h = std::max(1.0,u.norm()/v.norm())*tol;
737  }
738  // Evaluate Jacobian at new control
739  Ptr<Vector<Real>> znew = z.clone();
740  znew->set(z);
741  znew->axpy(h,v);
742  update(u,*znew);
743  applyAdjointJacobian_1(ahwv,w,u,*znew,jtol);
744  // Evaluate Jacobian at old control
745  Ptr<Vector<Real>> jv = ahwv.clone();
746  update(u,z);
747  applyAdjointJacobian_1(*jv,w,u,z,jtol);
748  // Compute Newton quotient
749  ahwv.axpy(-1.0,*jv);
750  ahwv.scale(1.0/h);
751  }
752 
770  virtual void applyAdjointHessian_22(Vector<Real> &ahwv,
771  const Vector<Real> &w,
772  const Vector<Real> &v,
773  const Vector<Real> &u,
774  const Vector<Real> &z,
775  Real &tol) {
776  Real jtol = std::sqrt(ROL_EPSILON<Real>());
777  // Compute step size
778  Real h = tol;
779  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
780  h = std::max(1.0,u.norm()/v.norm())*tol;
781  }
782  // Evaluate Jacobian at new control
783  Ptr<Vector<Real>> znew = z.clone();
784  znew->set(z);
785  znew->axpy(h,v);
786  update(u,*znew);
787  applyAdjointJacobian_2(ahwv,w,u,*znew,jtol);
788  // Evaluate Jacobian at old control
789  Ptr<Vector<Real>> jv = ahwv.clone();
790  update(u,z);
791  applyAdjointJacobian_2(*jv,w,u,z,jtol);
792  // Compute Newton quotient
793  ahwv.axpy(-1.0,*jv);
794  ahwv.scale(1.0/h);
795 }
796 
835  virtual std::vector<Real> solveAugmentedSystem(Vector<Real> &v1,
836  Vector<Real> &v2,
837  const Vector<Real> &b1,
838  const Vector<Real> &b2,
839  const Vector<Real> &x,
840  Real &tol) {
841  return Constraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,x,tol);
842  }
843 
844 
863  virtual void applyPreconditioner(Vector<Real> &pv,
864  const Vector<Real> &v,
865  const Vector<Real> &x,
866  const Vector<Real> &g,
867  Real &tol) {
868  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(x);
869  Ptr<Vector<Real>> ijv = (xs.get_1())->clone();
870 
871  try {
872  applyInverseJacobian_1(*ijv, v, *(xs.get_1()), *(xs.get_2()), tol);
873  }
874  catch (const std::logic_error &e) {
875  Constraint<Real>::applyPreconditioner(pv, v, x, g, tol);
876  return;
877  }
878 
879  const Vector_SimOpt<Real> &gs = dynamic_cast<const Vector_SimOpt<Real>&>(g);
880  Ptr<Vector<Real>> ijv_dual = (gs.get_1())->clone();
881  ijv_dual->set(ijv->dual());
882 
883  try {
884  applyInverseAdjointJacobian_1(pv, *ijv_dual, *(xs.get_1()), *(xs.get_2()), tol);
885  }
886  catch (const std::logic_error &e) {
887  Constraint<Real>::applyPreconditioner(pv, v, x, g, tol);
888  return;
889  }
890 
891  }
892 
898  virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
899  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
900  dynamic_cast<const Vector<Real>&>(x));
901  update(*(xs.get_1()),*(xs.get_2()),flag,iter);
902  }
903 
904  virtual void value(Vector<Real> &c,
905  const Vector<Real> &x,
906  Real &tol) {
907  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
908  dynamic_cast<const Vector<Real>&>(x));
909  value(c,*(xs.get_1()),*(xs.get_2()),tol);
910  }
911 
912 
913  virtual void applyJacobian(Vector<Real> &jv,
914  const Vector<Real> &v,
915  const Vector<Real> &x,
916  Real &tol) {
917  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
918  dynamic_cast<const Vector<Real>&>(x));
919  const Vector_SimOpt<Real> &vs = dynamic_cast<const Vector_SimOpt<Real>&>(
920  dynamic_cast<const Vector<Real>&>(v));
921  applyJacobian_1(jv,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
922  Ptr<Vector<Real>> jv2 = jv.clone();
923  applyJacobian_2(*jv2,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
924  jv.plus(*jv2);
925  }
926 
927 
928  using Constraint<Real>::applyAdjointJacobian;
929  virtual void applyAdjointJacobian(Vector<Real> &ajv,
930  const Vector<Real> &v,
931  const Vector<Real> &x,
932  Real &tol) {
933  Vector_SimOpt<Real> &ajvs = dynamic_cast<Vector_SimOpt<Real>&>(
934  dynamic_cast<Vector<Real>&>(ajv));
935  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
936  dynamic_cast<const Vector<Real>&>(x));
937  Ptr<Vector<Real>> ajv1 = (ajvs.get_1())->clone();
938  applyAdjointJacobian_1(*ajv1,v,*(xs.get_1()),*(xs.get_2()),tol);
939  ajvs.set_1(*ajv1);
940  Ptr<Vector<Real>> ajv2 = (ajvs.get_2())->clone();
941  applyAdjointJacobian_2(*ajv2,v,*(xs.get_1()),*(xs.get_2()),tol);
942  ajvs.set_2(*ajv2);
943  }
944 
945 
946  virtual void applyAdjointHessian(Vector<Real> &ahwv,
947  const Vector<Real> &w,
948  const Vector<Real> &v,
949  const Vector<Real> &x,
950  Real &tol) {
951  Vector_SimOpt<Real> &ahwvs = dynamic_cast<Vector_SimOpt<Real>&>(
952  dynamic_cast<Vector<Real>&>(ahwv));
953  const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
954  dynamic_cast<const Vector<Real>&>(x));
955  const Vector_SimOpt<Real> &vs = dynamic_cast<const Vector_SimOpt<Real>&>(
956  dynamic_cast<const Vector<Real>&>(v));
957  // Block-row 1
958  Ptr<Vector<Real>> C11 = (ahwvs.get_1())->clone();
959  Ptr<Vector<Real>> C21 = (ahwvs.get_1())->clone();
960  applyAdjointHessian_11(*C11,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
961  applyAdjointHessian_21(*C21,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
962  C11->plus(*C21);
963  ahwvs.set_1(*C11);
964  // Block-row 2
965  Ptr<Vector<Real>> C12 = (ahwvs.get_2())->clone();
966  Ptr<Vector<Real>> C22 = (ahwvs.get_2())->clone();
967  applyAdjointHessian_12(*C12,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
968  applyAdjointHessian_22(*C22,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
969  C22->plus(*C12);
970  ahwvs.set_2(*C22);
971  }
972 
973 
974 
975  virtual Real checkSolve(const Vector<Real> &u,
976  const Vector<Real> &z,
977  const Vector<Real> &c,
978  const bool printToStream = true,
979  std::ostream & outStream = std::cout) {
980  // Solve constraint for u.
981  Real tol = ROL_EPSILON<Real>();
982  Ptr<Vector<Real>> r = c.clone();
983  Ptr<Vector<Real>> s = u.clone();
984  solve(*r,*s,z,tol);
985  // Evaluate constraint residual at (u,z).
986  Ptr<Vector<Real>> cs = c.clone();
987  update(*s,z);
988  value(*cs,*s,z,tol);
989  // Output norm of residual.
990  Real rnorm = r->norm();
991  Real cnorm = cs->norm();
992  if ( printToStream ) {
993  std::stringstream hist;
994  hist << std::scientific << std::setprecision(8);
995  hist << "\nTest SimOpt solve at feasible (u,z):\n";
996  hist << " Solver Residual = " << rnorm << "\n";
997  hist << " ||c(u,z)|| = " << cnorm << "\n";
998  outStream << hist.str();
999  }
1000  return cnorm;
1001  }
1002 
1003 
1016  virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w,
1017  const Vector<Real> &v,
1018  const Vector<Real> &u,
1019  const Vector<Real> &z,
1020  const bool printToStream = true,
1021  std::ostream & outStream = std::cout) {
1022  return checkAdjointConsistencyJacobian_1(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
1023  }
1024 
1025 
1043  virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w,
1044  const Vector<Real> &v,
1045  const Vector<Real> &u,
1046  const Vector<Real> &z,
1047  const Vector<Real> &dualw,
1048  const Vector<Real> &dualv,
1049  const bool printToStream = true,
1050  std::ostream & outStream = std::cout) {
1051  Real tol = ROL_EPSILON<Real>();
1052  Ptr<Vector<Real>> Jv = dualw.clone();
1053  update(u,z);
1054  applyJacobian_1(*Jv,v,u,z,tol);
1055  Real wJv = w.dot(Jv->dual());
1056  Ptr<Vector<Real>> Jw = dualv.clone();
1057  update(u,z);
1058  applyAdjointJacobian_1(*Jw,w,u,z,tol);
1059  Real vJw = v.dot(Jw->dual());
1060  Real diff = std::abs(wJv-vJw);
1061  if ( printToStream ) {
1062  std::stringstream hist;
1063  hist << std::scientific << std::setprecision(8);
1064  hist << "\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1065  << diff << "\n";
1066  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1067  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1068  outStream << hist.str();
1069  }
1070  return diff;
1071  }
1072 
1073 
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 bool printToStream = true,
1091  std::ostream & outStream = std::cout) {
1092  return checkAdjointConsistencyJacobian_2(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
1093  }
1094 
1110  virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w,
1111  const Vector<Real> &v,
1112  const Vector<Real> &u,
1113  const Vector<Real> &z,
1114  const Vector<Real> &dualw,
1115  const Vector<Real> &dualv,
1116  const bool printToStream = true,
1117  std::ostream & outStream = std::cout) {
1118  Real tol = ROL_EPSILON<Real>();
1119  Ptr<Vector<Real>> Jv = dualw.clone();
1120  update(u,z);
1121  applyJacobian_2(*Jv,v,u,z,tol);
1122  Real wJv = w.dot(Jv->dual());
1123  Ptr<Vector<Real>> Jw = dualv.clone();
1124  update(u,z);
1125  applyAdjointJacobian_2(*Jw,w,u,z,tol);
1126  Real vJw = v.dot(Jw->dual());
1127  Real diff = std::abs(wJv-vJw);
1128  if ( printToStream ) {
1129  std::stringstream hist;
1130  hist << std::scientific << std::setprecision(8);
1131  hist << "\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1132  << diff << "\n";
1133  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1134  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1135  outStream << hist.str();
1136  }
1137  return diff;
1138  }
1139 
1140  virtual Real checkInverseJacobian_1(const Vector<Real> &jv,
1141  const Vector<Real> &v,
1142  const Vector<Real> &u,
1143  const Vector<Real> &z,
1144  const bool printToStream = true,
1145  std::ostream & outStream = std::cout) {
1146  Real tol = ROL_EPSILON<Real>();
1147  Ptr<Vector<Real>> Jv = jv.clone();
1148  update(u,z);
1149  applyJacobian_1(*Jv,v,u,z,tol);
1150  Ptr<Vector<Real>> iJJv = u.clone();
1151  update(u,z); // Does this update do anything?
1152  applyInverseJacobian_1(*iJJv,*Jv,u,z,tol);
1153  Ptr<Vector<Real>> diff = v.clone();
1154  diff->set(v);
1155  diff->axpy(-1.0,*iJJv);
1156  Real dnorm = diff->norm();
1157  Real vnorm = v.norm();
1158  if ( printToStream ) {
1159  std::stringstream hist;
1160  hist << std::scientific << std::setprecision(8);
1161  hist << "\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = "
1162  << dnorm << "\n";
1163  hist << " ||v|| = " << vnorm << "\n";
1164  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1165  outStream << hist.str();
1166  }
1167  return dnorm;
1168  }
1169 
1170  virtual Real checkInverseAdjointJacobian_1(const Vector<Real> &jv,
1171  const Vector<Real> &v,
1172  const Vector<Real> &u,
1173  const Vector<Real> &z,
1174  const bool printToStream = true,
1175  std::ostream & outStream = std::cout) {
1176  Real tol = ROL_EPSILON<Real>();
1177  Ptr<Vector<Real>> Jv = jv.clone();
1178  update(u,z);
1179  applyAdjointJacobian_1(*Jv,v,u,z,tol);
1180  Ptr<Vector<Real>> iJJv = v.clone();
1181  update(u,z);
1182  applyInverseAdjointJacobian_1(*iJJv,*Jv,u,z,tol);
1183  Ptr<Vector<Real>> diff = v.clone();
1184  diff->set(v);
1185  diff->axpy(-1.0,*iJJv);
1186  Real dnorm = diff->norm();
1187  Real vnorm = v.norm();
1188  if ( printToStream ) {
1189  std::stringstream hist;
1190  hist << std::scientific << std::setprecision(8);
1191  hist << "\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = "
1192  << dnorm << "\n";
1193  hist << " ||v|| = " << vnorm << "\n";
1194  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1195  outStream << hist.str();
1196  }
1197  return dnorm;
1198  }
1199 
1200 
1201 
1202  std::vector<std::vector<Real>> checkApplyJacobian_1(const Vector<Real> &u,
1203  const Vector<Real> &z,
1204  const Vector<Real> &v,
1205  const Vector<Real> &jv,
1206  const bool printToStream = true,
1207  std::ostream & outStream = std::cout,
1208  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1209  const int order = 1) {
1210  std::vector<Real> steps(numSteps);
1211  for(int i=0;i<numSteps;++i) {
1212  steps[i] = pow(10,-i);
1213  }
1214 
1215  return checkApplyJacobian_1(u,z,v,jv,steps,printToStream,outStream,order);
1216  }
1217 
1218 
1219 
1220 
1221  std::vector<std::vector<Real>> checkApplyJacobian_1(const Vector<Real> &u,
1222  const Vector<Real> &z,
1223  const Vector<Real> &v,
1224  const Vector<Real> &jv,
1225  const std::vector<Real> &steps,
1226  const bool printToStream = true,
1227  std::ostream & outStream = std::cout,
1228  const int order = 1) {
1229 
1230  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1231  "Error: finite difference order must be 1,2,3, or 4" );
1232 
1233  Real one(1.0);
1234 
1235  using Finite_Difference_Arrays::shifts;
1236  using Finite_Difference_Arrays::weights;
1237 
1238  Real tol = std::sqrt(ROL_EPSILON<Real>());
1239 
1240  int numSteps = steps.size();
1241  int numVals = 4;
1242  std::vector<Real> tmp(numVals);
1243  std::vector<std::vector<Real>> jvCheck(numSteps, tmp);
1244 
1245  // Save the format state of the original outStream.
1246  nullstream oldFormatState;
1247  oldFormatState.copyfmt(outStream);
1248 
1249  // Compute constraint value at x.
1250  Ptr<Vector<Real>> c = jv.clone();
1251  this->update(u,z);
1252  this->value(*c, u, z, tol);
1253 
1254  // Compute (Jacobian at x) times (vector v).
1255  Ptr<Vector<Real>> Jv = jv.clone();
1256  this->applyJacobian_1(*Jv, v, u, z, tol);
1257  Real normJv = Jv->norm();
1258 
1259  // Temporary vectors.
1260  Ptr<Vector<Real>> cdif = jv.clone();
1261  Ptr<Vector<Real>> cnew = jv.clone();
1262  Ptr<Vector<Real>> unew = u.clone();
1263 
1264  for (int i=0; i<numSteps; i++) {
1265 
1266  Real eta = steps[i];
1267 
1268  unew->set(u);
1269 
1270  cdif->set(*c);
1271  cdif->scale(weights[order-1][0]);
1272 
1273  for(int j=0; j<order; ++j) {
1274 
1275  unew->axpy(eta*shifts[order-1][j], v);
1276 
1277  if( weights[order-1][j+1] != 0 ) {
1278  this->update(*unew,z);
1279  this->value(*cnew,*unew,z,tol);
1280  cdif->axpy(weights[order-1][j+1],*cnew);
1281  }
1282 
1283  }
1284 
1285  cdif->scale(one/eta);
1286 
1287  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1288  jvCheck[i][0] = eta;
1289  jvCheck[i][1] = normJv;
1290  jvCheck[i][2] = cdif->norm();
1291  cdif->axpy(-one, *Jv);
1292  jvCheck[i][3] = cdif->norm();
1293 
1294  if (printToStream) {
1295  std::stringstream hist;
1296  if (i==0) {
1297  hist << std::right
1298  << std::setw(20) << "Step size"
1299  << std::setw(20) << "norm(Jac*vec)"
1300  << std::setw(20) << "norm(FD approx)"
1301  << std::setw(20) << "norm(abs error)"
1302  << "\n"
1303  << std::setw(20) << "---------"
1304  << std::setw(20) << "-------------"
1305  << std::setw(20) << "---------------"
1306  << std::setw(20) << "---------------"
1307  << "\n";
1308  }
1309  hist << std::scientific << std::setprecision(11) << std::right
1310  << std::setw(20) << jvCheck[i][0]
1311  << std::setw(20) << jvCheck[i][1]
1312  << std::setw(20) << jvCheck[i][2]
1313  << std::setw(20) << jvCheck[i][3]
1314  << "\n";
1315  outStream << hist.str();
1316  }
1317 
1318  }
1319 
1320  // Reset format state of outStream.
1321  outStream.copyfmt(oldFormatState);
1322 
1323  return jvCheck;
1324  } // checkApplyJacobian
1325 
1326 
1327  std::vector<std::vector<Real>> checkApplyJacobian_2(const Vector<Real> &u,
1328  const Vector<Real> &z,
1329  const Vector<Real> &v,
1330  const Vector<Real> &jv,
1331  const bool printToStream = true,
1332  std::ostream & outStream = std::cout,
1333  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1334  const int order = 1) {
1335  std::vector<Real> steps(numSteps);
1336  for(int i=0;i<numSteps;++i) {
1337  steps[i] = pow(10,-i);
1338  }
1339 
1340  return checkApplyJacobian_2(u,z,v,jv,steps,printToStream,outStream,order);
1341  }
1342 
1343 
1344 
1345 
1346  std::vector<std::vector<Real>> checkApplyJacobian_2(const Vector<Real> &u,
1347  const Vector<Real> &z,
1348  const Vector<Real> &v,
1349  const Vector<Real> &jv,
1350  const std::vector<Real> &steps,
1351  const bool printToStream = true,
1352  std::ostream & outStream = std::cout,
1353  const int order = 1) {
1354 
1355  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1356  "Error: finite difference order must be 1,2,3, or 4" );
1357 
1358  Real one(1.0);
1359 
1360  using Finite_Difference_Arrays::shifts;
1361  using Finite_Difference_Arrays::weights;
1362 
1363  Real tol = std::sqrt(ROL_EPSILON<Real>());
1364 
1365  int numSteps = steps.size();
1366  int numVals = 4;
1367  std::vector<Real> tmp(numVals);
1368  std::vector<std::vector<Real>> jvCheck(numSteps, tmp);
1369 
1370  // Save the format state of the original outStream.
1371  nullstream oldFormatState;
1372  oldFormatState.copyfmt(outStream);
1373 
1374  // Compute constraint value at x.
1375  Ptr<Vector<Real>> c = jv.clone();
1376  this->update(u,z);
1377  this->value(*c, u, z, tol);
1378 
1379  // Compute (Jacobian at x) times (vector v).
1380  Ptr<Vector<Real>> Jv = jv.clone();
1381  this->applyJacobian_2(*Jv, v, u, z, tol);
1382  Real normJv = Jv->norm();
1383 
1384  // Temporary vectors.
1385  Ptr<Vector<Real>> cdif = jv.clone();
1386  Ptr<Vector<Real>> cnew = jv.clone();
1387  Ptr<Vector<Real>> znew = z.clone();
1388 
1389  for (int i=0; i<numSteps; i++) {
1390 
1391  Real eta = steps[i];
1392 
1393  znew->set(z);
1394 
1395  cdif->set(*c);
1396  cdif->scale(weights[order-1][0]);
1397 
1398  for(int j=0; j<order; ++j) {
1399 
1400  znew->axpy(eta*shifts[order-1][j], v);
1401 
1402  if( weights[order-1][j+1] != 0 ) {
1403  this->update(u,*znew);
1404  this->value(*cnew,u,*znew,tol);
1405  cdif->axpy(weights[order-1][j+1],*cnew);
1406  }
1407 
1408  }
1409 
1410  cdif->scale(one/eta);
1411 
1412  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1413  jvCheck[i][0] = eta;
1414  jvCheck[i][1] = normJv;
1415  jvCheck[i][2] = cdif->norm();
1416  cdif->axpy(-one, *Jv);
1417  jvCheck[i][3] = cdif->norm();
1418 
1419  if (printToStream) {
1420  std::stringstream hist;
1421  if (i==0) {
1422  hist << std::right
1423  << std::setw(20) << "Step size"
1424  << std::setw(20) << "norm(Jac*vec)"
1425  << std::setw(20) << "norm(FD approx)"
1426  << std::setw(20) << "norm(abs error)"
1427  << "\n"
1428  << std::setw(20) << "---------"
1429  << std::setw(20) << "-------------"
1430  << std::setw(20) << "---------------"
1431  << std::setw(20) << "---------------"
1432  << "\n";
1433  }
1434  hist << std::scientific << std::setprecision(11) << std::right
1435  << std::setw(20) << jvCheck[i][0]
1436  << std::setw(20) << jvCheck[i][1]
1437  << std::setw(20) << jvCheck[i][2]
1438  << std::setw(20) << jvCheck[i][3]
1439  << "\n";
1440  outStream << hist.str();
1441  }
1442 
1443  }
1444 
1445  // Reset format state of outStream.
1446  outStream.copyfmt(oldFormatState);
1447 
1448  return jvCheck;
1449  } // checkApplyJacobian
1450 
1451 
1452 
1453  std::vector<std::vector<Real>> checkApplyAdjointHessian_11(const Vector<Real> &u,
1454  const Vector<Real> &z,
1455  const Vector<Real> &p,
1456  const Vector<Real> &v,
1457  const Vector<Real> &hv,
1458  const bool printToStream = true,
1459  std::ostream & outStream = std::cout,
1460  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1461  const int order = 1 ) {
1462  std::vector<Real> steps(numSteps);
1463  for(int i=0;i<numSteps;++i) {
1464  steps[i] = pow(10,-i);
1465  }
1466 
1467  return checkApplyAdjointHessian_11(u,z,p,v,hv,steps,printToStream,outStream,order);
1468 
1469  }
1470 
1471  std::vector<std::vector<Real>> checkApplyAdjointHessian_11(const Vector<Real> &u,
1472  const Vector<Real> &z,
1473  const Vector<Real> &p,
1474  const Vector<Real> &v,
1475  const Vector<Real> &hv,
1476  const std::vector<Real> &steps,
1477  const bool printToStream = true,
1478  std::ostream & outStream = std::cout,
1479  const int order = 1 ) {
1480  using Finite_Difference_Arrays::shifts;
1481  using Finite_Difference_Arrays::weights;
1482 
1483  Real one(1.0);
1484 
1485  Real tol = std::sqrt(ROL_EPSILON<Real>());
1486 
1487  int numSteps = steps.size();
1488  int numVals = 4;
1489  std::vector<Real> tmp(numVals);
1490  std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1491 
1492  // Temporary vectors.
1493  Ptr<Vector<Real>> AJdif = hv.clone();
1494  Ptr<Vector<Real>> AJp = hv.clone();
1495  Ptr<Vector<Real>> AHpv = hv.clone();
1496  Ptr<Vector<Real>> AJnew = hv.clone();
1497  Ptr<Vector<Real>> unew = u.clone();
1498 
1499  // Save the format state of the original outStream.
1500  nullstream oldFormatState;
1501  oldFormatState.copyfmt(outStream);
1502 
1503  // Apply adjoint Jacobian to p.
1504  this->update(u,z);
1505  this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1506 
1507  // Apply adjoint Hessian at (u,z), in direction v, to p.
1508  this->applyAdjointHessian_11(*AHpv, p, v, u, z, tol);
1509  Real normAHpv = AHpv->norm();
1510 
1511  for (int i=0; i<numSteps; i++) {
1512 
1513  Real eta = steps[i];
1514 
1515  // Apply adjoint Jacobian to p at (u+eta*v,z).
1516  unew->set(u);
1517 
1518  AJdif->set(*AJp);
1519  AJdif->scale(weights[order-1][0]);
1520 
1521  for(int j=0; j<order; ++j) {
1522 
1523  unew->axpy(eta*shifts[order-1][j],v);
1524 
1525  if( weights[order-1][j+1] != 0 ) {
1526  this->update(*unew,z);
1527  this->applyAdjointJacobian_1(*AJnew, p, *unew, z, tol);
1528  AJdif->axpy(weights[order-1][j+1],*AJnew);
1529  }
1530  }
1531 
1532  AJdif->scale(one/eta);
1533 
1534  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1535  ahpvCheck[i][0] = eta;
1536  ahpvCheck[i][1] = normAHpv;
1537  ahpvCheck[i][2] = AJdif->norm();
1538  AJdif->axpy(-one, *AHpv);
1539  ahpvCheck[i][3] = AJdif->norm();
1540 
1541  if (printToStream) {
1542  std::stringstream hist;
1543  if (i==0) {
1544  hist << std::right
1545  << std::setw(20) << "Step size"
1546  << std::setw(20) << "norm(adj(H)(u,v))"
1547  << std::setw(20) << "norm(FD approx)"
1548  << std::setw(20) << "norm(abs error)"
1549  << "\n"
1550  << std::setw(20) << "---------"
1551  << std::setw(20) << "-----------------"
1552  << std::setw(20) << "---------------"
1553  << std::setw(20) << "---------------"
1554  << "\n";
1555  }
1556  hist << std::scientific << std::setprecision(11) << std::right
1557  << std::setw(20) << ahpvCheck[i][0]
1558  << std::setw(20) << ahpvCheck[i][1]
1559  << std::setw(20) << ahpvCheck[i][2]
1560  << std::setw(20) << ahpvCheck[i][3]
1561  << "\n";
1562  outStream << hist.str();
1563  }
1564 
1565  }
1566 
1567  // Reset format state of outStream.
1568  outStream.copyfmt(oldFormatState);
1569 
1570  return ahpvCheck;
1571  } // checkApplyAdjointHessian_11
1572 
1576  std::vector<std::vector<Real>> checkApplyAdjointHessian_21(const Vector<Real> &u,
1577  const Vector<Real> &z,
1578  const Vector<Real> &p,
1579  const Vector<Real> &v,
1580  const Vector<Real> &hv,
1581  const bool printToStream = true,
1582  std::ostream & outStream = std::cout,
1583  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1584  const int order = 1 ) {
1585  std::vector<Real> steps(numSteps);
1586  for(int i=0;i<numSteps;++i) {
1587  steps[i] = pow(10,-i);
1588  }
1589 
1590  return checkApplyAdjointHessian_21(u,z,p,v,hv,steps,printToStream,outStream,order);
1591 
1592  }
1593 
1597  std::vector<std::vector<Real>> checkApplyAdjointHessian_21(const Vector<Real> &u,
1598  const Vector<Real> &z,
1599  const Vector<Real> &p,
1600  const Vector<Real> &v,
1601  const Vector<Real> &hv,
1602  const std::vector<Real> &steps,
1603  const bool printToStream = true,
1604  std::ostream & outStream = std::cout,
1605  const int order = 1 ) {
1606  using Finite_Difference_Arrays::shifts;
1607  using Finite_Difference_Arrays::weights;
1608 
1609  Real one(1.0);
1610 
1611  Real tol = std::sqrt(ROL_EPSILON<Real>());
1612 
1613  int numSteps = steps.size();
1614  int numVals = 4;
1615  std::vector<Real> tmp(numVals);
1616  std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1617 
1618  // Temporary vectors.
1619  Ptr<Vector<Real>> AJdif = hv.clone();
1620  Ptr<Vector<Real>> AJp = hv.clone();
1621  Ptr<Vector<Real>> AHpv = hv.clone();
1622  Ptr<Vector<Real>> AJnew = hv.clone();
1623  Ptr<Vector<Real>> znew = z.clone();
1624 
1625  // Save the format state of the original outStream.
1626  nullstream oldFormatState;
1627  oldFormatState.copyfmt(outStream);
1628 
1629  // Apply adjoint Jacobian to p.
1630  this->update(u,z);
1631  this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1632 
1633  // Apply adjoint Hessian at (u,z), in direction v, to p.
1634  this->applyAdjointHessian_21(*AHpv, p, v, u, z, tol);
1635  Real normAHpv = AHpv->norm();
1636 
1637  for (int i=0; i<numSteps; i++) {
1638 
1639  Real eta = steps[i];
1640 
1641  // Apply adjoint Jacobian to p at (u,z+eta*v).
1642  znew->set(z);
1643 
1644  AJdif->set(*AJp);
1645  AJdif->scale(weights[order-1][0]);
1646 
1647  for(int j=0; j<order; ++j) {
1648 
1649  znew->axpy(eta*shifts[order-1][j],v);
1650 
1651  if( weights[order-1][j+1] != 0 ) {
1652  this->update(u,*znew);
1653  this->applyAdjointJacobian_1(*AJnew, p, u, *znew, tol);
1654  AJdif->axpy(weights[order-1][j+1],*AJnew);
1655  }
1656  }
1657 
1658  AJdif->scale(one/eta);
1659 
1660  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1661  ahpvCheck[i][0] = eta;
1662  ahpvCheck[i][1] = normAHpv;
1663  ahpvCheck[i][2] = AJdif->norm();
1664  AJdif->axpy(-one, *AHpv);
1665  ahpvCheck[i][3] = AJdif->norm();
1666 
1667  if (printToStream) {
1668  std::stringstream hist;
1669  if (i==0) {
1670  hist << std::right
1671  << std::setw(20) << "Step size"
1672  << std::setw(20) << "norm(adj(H)(u,v))"
1673  << std::setw(20) << "norm(FD approx)"
1674  << std::setw(20) << "norm(abs error)"
1675  << "\n"
1676  << std::setw(20) << "---------"
1677  << std::setw(20) << "-----------------"
1678  << std::setw(20) << "---------------"
1679  << std::setw(20) << "---------------"
1680  << "\n";
1681  }
1682  hist << std::scientific << std::setprecision(11) << std::right
1683  << std::setw(20) << ahpvCheck[i][0]
1684  << std::setw(20) << ahpvCheck[i][1]
1685  << std::setw(20) << ahpvCheck[i][2]
1686  << std::setw(20) << ahpvCheck[i][3]
1687  << "\n";
1688  outStream << hist.str();
1689  }
1690 
1691  }
1692 
1693  // Reset format state of outStream.
1694  outStream.copyfmt(oldFormatState);
1695 
1696  return ahpvCheck;
1697  } // checkApplyAdjointHessian_21
1698 
1702  std::vector<std::vector<Real>> checkApplyAdjointHessian_12(const Vector<Real> &u,
1703  const Vector<Real> &z,
1704  const Vector<Real> &p,
1705  const Vector<Real> &v,
1706  const Vector<Real> &hv,
1707  const bool printToStream = true,
1708  std::ostream & outStream = std::cout,
1709  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1710  const int order = 1 ) {
1711  std::vector<Real> steps(numSteps);
1712  for(int i=0;i<numSteps;++i) {
1713  steps[i] = pow(10,-i);
1714  }
1715 
1716  return checkApplyAdjointHessian_12(u,z,p,v,hv,steps,printToStream,outStream,order);
1717 
1718  }
1719 
1720 
1721  std::vector<std::vector<Real>> checkApplyAdjointHessian_12(const Vector<Real> &u,
1722  const Vector<Real> &z,
1723  const Vector<Real> &p,
1724  const Vector<Real> &v,
1725  const Vector<Real> &hv,
1726  const std::vector<Real> &steps,
1727  const bool printToStream = true,
1728  std::ostream & outStream = std::cout,
1729  const int order = 1 ) {
1730  using Finite_Difference_Arrays::shifts;
1731  using Finite_Difference_Arrays::weights;
1732 
1733  Real one(1.0);
1734 
1735  Real tol = std::sqrt(ROL_EPSILON<Real>());
1736 
1737  int numSteps = steps.size();
1738  int numVals = 4;
1739  std::vector<Real> tmp(numVals);
1740  std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1741 
1742  // Temporary vectors.
1743  Ptr<Vector<Real>> AJdif = hv.clone();
1744  Ptr<Vector<Real>> AJp = hv.clone();
1745  Ptr<Vector<Real>> AHpv = hv.clone();
1746  Ptr<Vector<Real>> AJnew = hv.clone();
1747  Ptr<Vector<Real>> unew = u.clone();
1748 
1749  // Save the format state of the original outStream.
1750  nullstream oldFormatState;
1751  oldFormatState.copyfmt(outStream);
1752 
1753  // Apply adjoint Jacobian to p.
1754  this->update(u,z);
1755  this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1756 
1757  // Apply adjoint Hessian at (u,z), in direction v, to p.
1758  this->applyAdjointHessian_12(*AHpv, p, v, u, z, tol);
1759  Real normAHpv = AHpv->norm();
1760 
1761  for (int i=0; i<numSteps; i++) {
1762 
1763  Real eta = steps[i];
1764 
1765  // Apply adjoint Jacobian to p at (u+eta*v,z).
1766  unew->set(u);
1767 
1768  AJdif->set(*AJp);
1769  AJdif->scale(weights[order-1][0]);
1770 
1771  for(int j=0; j<order; ++j) {
1772 
1773  unew->axpy(eta*shifts[order-1][j],v);
1774 
1775  if( weights[order-1][j+1] != 0 ) {
1776  this->update(*unew,z);
1777  this->applyAdjointJacobian_2(*AJnew, p, *unew, z, tol);
1778  AJdif->axpy(weights[order-1][j+1],*AJnew);
1779  }
1780  }
1781 
1782  AJdif->scale(one/eta);
1783 
1784  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1785  ahpvCheck[i][0] = eta;
1786  ahpvCheck[i][1] = normAHpv;
1787  ahpvCheck[i][2] = AJdif->norm();
1788  AJdif->axpy(-one, *AHpv);
1789  ahpvCheck[i][3] = AJdif->norm();
1790 
1791  if (printToStream) {
1792  std::stringstream hist;
1793  if (i==0) {
1794  hist << std::right
1795  << std::setw(20) << "Step size"
1796  << std::setw(20) << "norm(adj(H)(u,v))"
1797  << std::setw(20) << "norm(FD approx)"
1798  << std::setw(20) << "norm(abs error)"
1799  << "\n"
1800  << std::setw(20) << "---------"
1801  << std::setw(20) << "-----------------"
1802  << std::setw(20) << "---------------"
1803  << std::setw(20) << "---------------"
1804  << "\n";
1805  }
1806  hist << std::scientific << std::setprecision(11) << std::right
1807  << std::setw(20) << ahpvCheck[i][0]
1808  << std::setw(20) << ahpvCheck[i][1]
1809  << std::setw(20) << ahpvCheck[i][2]
1810  << std::setw(20) << ahpvCheck[i][3]
1811  << "\n";
1812  outStream << hist.str();
1813  }
1814 
1815  }
1816 
1817  // Reset format state of outStream.
1818  outStream.copyfmt(oldFormatState);
1819 
1820  return ahpvCheck;
1821  } // checkApplyAdjointHessian_12
1822 
1823  std::vector<std::vector<Real>> checkApplyAdjointHessian_22(const Vector<Real> &u,
1824  const Vector<Real> &z,
1825  const Vector<Real> &p,
1826  const Vector<Real> &v,
1827  const Vector<Real> &hv,
1828  const bool printToStream = true,
1829  std::ostream & outStream = std::cout,
1830  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1831  const int order = 1 ) {
1832  std::vector<Real> steps(numSteps);
1833  for(int i=0;i<numSteps;++i) {
1834  steps[i] = pow(10,-i);
1835  }
1836 
1837  return checkApplyAdjointHessian_22(u,z,p,v,hv,steps,printToStream,outStream,order);
1838 
1839  }
1840 
1841  std::vector<std::vector<Real>> checkApplyAdjointHessian_22(const Vector<Real> &u,
1842  const Vector<Real> &z,
1843  const Vector<Real> &p,
1844  const Vector<Real> &v,
1845  const Vector<Real> &hv,
1846  const std::vector<Real> &steps,
1847  const bool printToStream = true,
1848  std::ostream & outStream = std::cout,
1849  const int order = 1 ) {
1850  using Finite_Difference_Arrays::shifts;
1851  using Finite_Difference_Arrays::weights;
1852 
1853  Real one(1.0);
1854 
1855  Real tol = std::sqrt(ROL_EPSILON<Real>());
1856 
1857  int numSteps = steps.size();
1858  int numVals = 4;
1859  std::vector<Real> tmp(numVals);
1860  std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1861 
1862  // Temporary vectors.
1863  Ptr<Vector<Real>> AJdif = hv.clone();
1864  Ptr<Vector<Real>> AJp = hv.clone();
1865  Ptr<Vector<Real>> AHpv = hv.clone();
1866  Ptr<Vector<Real>> AJnew = hv.clone();
1867  Ptr<Vector<Real>> znew = z.clone();
1868 
1869  // Save the format state of the original outStream.
1870  nullstream oldFormatState;
1871  oldFormatState.copyfmt(outStream);
1872 
1873  // Apply adjoint Jacobian to p.
1874  this->update(u,z);
1875  this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1876 
1877  // Apply adjoint Hessian at (u,z), in direction v, to p.
1878  this->applyAdjointHessian_22(*AHpv, p, v, u, z, tol);
1879  Real normAHpv = AHpv->norm();
1880 
1881  for (int i=0; i<numSteps; i++) {
1882 
1883  Real eta = steps[i];
1884 
1885  // Apply adjoint Jacobian to p at (u,z+eta*v).
1886  znew->set(z);
1887 
1888  AJdif->set(*AJp);
1889  AJdif->scale(weights[order-1][0]);
1890 
1891  for(int j=0; j<order; ++j) {
1892 
1893  znew->axpy(eta*shifts[order-1][j],v);
1894 
1895  if( weights[order-1][j+1] != 0 ) {
1896  this->update(u,*znew);
1897  this->applyAdjointJacobian_2(*AJnew, p, u, *znew, tol);
1898  AJdif->axpy(weights[order-1][j+1],*AJnew);
1899  }
1900  }
1901 
1902  AJdif->scale(one/eta);
1903 
1904  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1905  ahpvCheck[i][0] = eta;
1906  ahpvCheck[i][1] = normAHpv;
1907  ahpvCheck[i][2] = AJdif->norm();
1908  AJdif->axpy(-one, *AHpv);
1909  ahpvCheck[i][3] = AJdif->norm();
1910 
1911  if (printToStream) {
1912  std::stringstream hist;
1913  if (i==0) {
1914  hist << std::right
1915  << std::setw(20) << "Step size"
1916  << std::setw(20) << "norm(adj(H)(u,v))"
1917  << std::setw(20) << "norm(FD approx)"
1918  << std::setw(20) << "norm(abs error)"
1919  << "\n"
1920  << std::setw(20) << "---------"
1921  << std::setw(20) << "-----------------"
1922  << std::setw(20) << "---------------"
1923  << std::setw(20) << "---------------"
1924  << "\n";
1925  }
1926  hist << std::scientific << std::setprecision(11) << std::right
1927  << std::setw(20) << ahpvCheck[i][0]
1928  << std::setw(20) << ahpvCheck[i][1]
1929  << std::setw(20) << ahpvCheck[i][2]
1930  << std::setw(20) << ahpvCheck[i][3]
1931  << "\n";
1932  outStream << hist.str();
1933  }
1934 
1935  }
1936 
1937  // Reset format state of outStream.
1938  outStream.copyfmt(oldFormatState);
1939 
1940  return ahpvCheck;
1941  } // checkApplyAdjointHessian_22
1942 
1943 }; // class Constraint_SimOpt
1944 
1945 } // namespace ROL
1946 
1947 #endif
ROL_Algorithm.hpp
ROL_CompositeStep.hpp
ROL_ConstraintStatusTest.hpp
ROL_Constraint.hpp
ROL_Constraint_State.hpp
ROL_NonlinearLeastSquaresObjective_SimOpt.hpp
ROL_Objective_FSsolver.hpp
ROL_TrustRegionStep.hpp
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL_NUM_CHECKDERIV_STEPS
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
Definition: ROL_Types.hpp:74
ROL_Vector_SimOpt.hpp
ROL::Constraint_SimOpt
Defines the constraint operator interface for simulation-based optimization.
Definition: ROL_Constraint_SimOpt.hpp:103
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 i...
Definition: ROL_Constraint_SimOpt.hpp:158
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:540
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:1140
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:1016
ROL::Constraint_SimOpt::firstSolve_
bool firstSolve_
Definition: ROL_Constraint_SimOpt.hpp:134
ROL::Constraint_SimOpt::unew_
Ptr< Vector< Real > > unew_
Definition: ROL_Constraint_SimOpt.hpp:106
ROL::Constraint_SimOpt::decr_
Real decr_
Definition: ROL_Constraint_SimOpt.hpp:127
ROL::Constraint_SimOpt::rtol_
Real rtol_
Definition: ROL_Constraint_SimOpt.hpp:124
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:1576
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:445
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
Definition: ROL_Constraint_SimOpt.hpp:835
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,...
Definition: ROL_Constraint_SimOpt.hpp:1110
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:610
ROL::Constraint_SimOpt::DEFAULT_decr_
const Real DEFAULT_decr_
Definition: ROL_Constraint_SimOpt.hpp:114
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:403
ROL::Constraint_SimOpt::zero_
bool zero_
Definition: ROL_Constraint_SimOpt.hpp:130
ROL::Constraint_SimOpt::maxit_
int maxit_
Definition: ROL_Constraint_SimOpt.hpp:128
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:1841
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:681
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:1471
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:1346
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:636
ROL::Constraint_SimOpt::atol_
Real atol_
Definition: ROL_Constraint_SimOpt.hpp:123
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,...
Definition: ROL_Constraint_SimOpt.hpp:168
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,...
Definition: ROL_Constraint_SimOpt.hpp:175
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:770
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::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:495
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:566
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,...
Definition: ROL_Constraint_SimOpt.hpp:1043
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:469
ROL::Constraint_SimOpt::solverType_
int solverType_
Definition: ROL_Constraint_SimOpt.hpp:131
ROL::Constraint_SimOpt::DEFAULT_rtol_
const Real DEFAULT_rtol_
Definition: ROL_Constraint_SimOpt.hpp:111
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:929
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:1086
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:726
ROL::Constraint_SimOpt::DEFAULT_solverType_
const int DEFAULT_solverType_
Definition: ROL_Constraint_SimOpt.hpp:118
ROL::Constraint_SimOpt::setSolveParameters
virtual void setSolveParameters(ParameterList &parlist)
Set solve parameters.
Definition: ROL_Constraint_SimOpt.hpp:332
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:1170
ROL::Constraint_SimOpt::DEFAULT_factor_
const Real DEFAULT_factor_
Definition: ROL_Constraint_SimOpt.hpp:113
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:1823
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:1221
ROL::Constraint_SimOpt::DEFAULT_atol_
const Real DEFAULT_atol_
Definition: ROL_Constraint_SimOpt.hpp:110
ROL::Constraint_SimOpt::DEFAULT_zero_
const bool DEFAULT_zero_
Definition: ROL_Constraint_SimOpt.hpp:117
ROL::Constraint_SimOpt::print_
bool print_
Definition: ROL_Constraint_SimOpt.hpp:129
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:360
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:1702
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:1202
ROL::Constraint_SimOpt::DEFAULT_print_
const bool DEFAULT_print_
Definition: ROL_Constraint_SimOpt.hpp:116
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 i...
Definition: ROL_Constraint_SimOpt.hpp:898
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:946
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:1721
ROL::Constraint_SimOpt::stol_
Real stol_
Definition: ROL_Constraint_SimOpt.hpp:125
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:1327
ROL::Constraint_SimOpt::jv_
Ptr< Vector< Real > > jv_
Definition: ROL_Constraint_SimOpt.hpp:107
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:207
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:913
ROL::Constraint_SimOpt::DEFAULT_maxit_
const int DEFAULT_maxit_
Definition: ROL_Constraint_SimOpt.hpp:115
ROL::Constraint_SimOpt::factor_
Real factor_
Definition: ROL_Constraint_SimOpt.hpp:126
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:863
ROL::Constraint_SimOpt::DEFAULT_stol_
const Real DEFAULT_stol_
Definition: ROL_Constraint_SimOpt.hpp:112
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:975
ROL::Constraint_SimOpt::Constraint_SimOpt
Constraint_SimOpt()
Definition: ROL_Constraint_SimOpt.hpp:137
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:1453
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:904
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:1597
ROL::Constraint
Defines the general constraint operator interface.
Definition: ROL_Constraint.hpp:85
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:269
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
Definition: ROL_ConstraintDef.hpp:195
ROL::Vector_SimOpt
Defines the linear algebra or vector space interface for simulation-based optimization.
Definition: ROL_Vector_SimOpt.hpp:57
ROL::Vector_SimOpt::get_1
ROL::Ptr< const Vector< Real > > get_1() const
Definition: ROL_Vector_SimOpt.hpp:168
ROL::Vector_SimOpt::get_2
ROL::Ptr< const Vector< Real > > get_2() const
Definition: ROL_Vector_SimOpt.hpp:172
ROL::Vector_SimOpt::set_1
void set_1(const Vector< Real > &vec)
Definition: ROL_Vector_SimOpt.hpp:184
ROL::Vector_SimOpt::set_2
void set_2(const Vector< Real > &vec)
Definition: ROL_Vector_SimOpt.hpp:188
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:107
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
ROL::Vector::basis
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:182
ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
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:226
ROL::Vector::dimension
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:196
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
ROL::Vector::dot
virtual Real dot(const Vector &x) const =0
Compute where .
ROL::Finite_Difference_Arrays::shifts
const int shifts[4][4]
Definition: ROL_Types.hpp:859
ROL::Finite_Difference_Arrays::weights
const double weights[4][5]
Definition: ROL_Types.hpp:866
ROL::details::nullstream
basic_nullstream< char, char_traits< char > > nullstream
Definition: ROL_Stream.hpp:72
ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91
Generated by doxygen 1.9.1