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