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