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