ROL
ROL_EqualityConstraint_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_EQUALITY_CONSTRAINT_SIMOPT_H
45 #define ROL_EQUALITY_CONSTRAINT_SIMOPT_H
46 
47 #include "ROL_EqualityConstraint.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_EqualityConstraint_State.hpp"
54 #include "ROL_Objective_FSsolver.hpp"
55 #include "ROL_Algorithm.hpp"
56 
97 namespace ROL {
98 
99 template <class Real>
100 class EqualityConstraint_SimOpt : public virtual EqualityConstraint<Real> {
101 private:
102  // Additional vector storage for solve
103  Teuchos::RCP<Vector<Real> > unew_;
104  Teuchos::RCP<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  EqualityConstraint_SimOpt()
133  : EqualityConstraint<Real>(),
134  unew_(Teuchos::null), jv_(Teuchos::null),
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 EqualityConstraint_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  Teuchos::RCP<EqualityConstraint_SimOpt<Real> > con = Teuchos::rcp(this,false);
270  Teuchos::RCP<Objective<Real> > obj = Teuchos::rcp(new NonlinearLeastSquaresObjective_SimOpt<Real>(con,u,z,c,true));
271  Teuchos::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  Teuchos::RCP<Algorithm<Real> > algo = Teuchos::rcp(new 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  Teuchos::RCP<EqualityConstraint_SimOpt<Real> > con = Teuchos::rcp(this,false);
283  Teuchos::RCP<const Vector<Real> > zVec = Teuchos::rcpFromRef(z);
284  Teuchos::RCP<EqualityConstraint<Real> > conU
285  = Teuchos::rcp(new EqualityConstraint_State<Real>(con,zVec));
286  Teuchos::RCP<Objective<Real> > objU
287  = Teuchos::rcp(new Objective_FSsolver<Real>());
288  Teuchos::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  Teuchos::RCP<Algorithm<Real> > algo = Teuchos::rcp(new Algorithm<Real>("Composite Step",parlist,false));
293  Teuchos::RCP<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  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
299  ">>> ERROR (ROL:EqualityConstraint_SimOpt:solve): Invalid solver type!");
300  }
301  }
302 
323  virtual void setSolveParameters(Teuchos::ParameterList &parlist) {
324  Teuchos::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  Teuchos::RCP<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  Teuchos::RCP<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  Teuchos::RCP<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  Teuchos::RCP<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  TEUCHOS_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  Teuchos::RCP<Vector<Real> > cold = dualv.clone();
498  Teuchos::RCP<Vector<Real> > cnew = dualv.clone();
499  update(u,z);
500  value(*cold,u,z,ctol);
501  Teuchos::RCP<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  Teuchos::RCP<Vector<Real> > cold = dualv.clone();
569  Teuchos::RCP<Vector<Real> > cnew = dualv.clone();
570  update(u,z);
571  value(*cold,u,z,ctol);
572  Teuchos::RCP<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  TEUCHOS_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  Teuchos::RCP<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  Teuchos::RCP<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  Teuchos::RCP<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  Teuchos::RCP<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  Teuchos::RCP<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  Teuchos::RCP<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  Teuchos::RCP<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  Teuchos::RCP<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 EqualityConstraint<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 = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(x);
860  Teuchos::RCP<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  EqualityConstraint<Real>::applyPreconditioner(pv, v, x, g, tol);
867  return;
868  }
869 
870  const Vector_SimOpt<Real> &gs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(g);
871  Teuchos::RCP<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  EqualityConstraint<Real>::applyPreconditioner(pv, v, x, g, tol);
879  return;
880  }
881 
882  }
883 
884 //
885 // EqualityConstraint_SimOpt(void) : EqualityConstraint<Real>() {}
886 
892  virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
893  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
894  Teuchos::dyn_cast<const Vector<Real> >(x));
895  update(*(xs.get_1()),*(xs.get_2()),flag,iter);
896  }
897 
900  virtual bool isFeasible( const Vector<Real> &v ) { return true; }
901 
902  virtual void value(Vector<Real> &c,
903  const Vector<Real> &x,
904  Real &tol) {
905  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
906  Teuchos::dyn_cast<const Vector<Real> >(x));
907  value(c,*(xs.get_1()),*(xs.get_2()),tol);
908  }
909 
910 
911  virtual void applyJacobian(Vector<Real> &jv,
912  const Vector<Real> &v,
913  const Vector<Real> &x,
914  Real &tol) {
915  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
916  Teuchos::dyn_cast<const Vector<Real> >(x));
917  const Vector_SimOpt<Real> &vs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
918  Teuchos::dyn_cast<const Vector<Real> >(v));
919  applyJacobian_1(jv,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
920  Teuchos::RCP<Vector<Real> > jv2 = jv.clone();
921  applyJacobian_2(*jv2,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
922  jv.plus(*jv2);
923  }
924 
925 
926  virtual void applyAdjointJacobian(Vector<Real> &ajv,
927  const Vector<Real> &v,
928  const Vector<Real> &x,
929  Real &tol) {
930  Vector_SimOpt<Real> &ajvs = Teuchos::dyn_cast<Vector_SimOpt<Real> >(
931  Teuchos::dyn_cast<Vector<Real> >(ajv));
932  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
933  Teuchos::dyn_cast<const Vector<Real> >(x));
934  Teuchos::RCP<Vector<Real> > ajv1 = (ajvs.get_1())->clone();
935  applyAdjointJacobian_1(*ajv1,v,*(xs.get_1()),*(xs.get_2()),tol);
936  ajvs.set_1(*ajv1);
937  Teuchos::RCP<Vector<Real> > ajv2 = (ajvs.get_2())->clone();
938  applyAdjointJacobian_2(*ajv2,v,*(xs.get_1()),*(xs.get_2()),tol);
939  ajvs.set_2(*ajv2);
940  }
941 
942 
943  virtual void applyAdjointHessian(Vector<Real> &ahwv,
944  const Vector<Real> &w,
945  const Vector<Real> &v,
946  const Vector<Real> &x,
947  Real &tol) {
948  Vector_SimOpt<Real> &ahwvs = Teuchos::dyn_cast<Vector_SimOpt<Real> >(
949  Teuchos::dyn_cast<Vector<Real> >(ahwv));
950  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
951  Teuchos::dyn_cast<const Vector<Real> >(x));
952  const Vector_SimOpt<Real> &vs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
953  Teuchos::dyn_cast<const Vector<Real> >(v));
954  // Block-row 1
955  Teuchos::RCP<Vector<Real> > C11 = (ahwvs.get_1())->clone();
956  Teuchos::RCP<Vector<Real> > C21 = (ahwvs.get_1())->clone();
957  applyAdjointHessian_11(*C11,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
958  applyAdjointHessian_21(*C21,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
959  C11->plus(*C21);
960  ahwvs.set_1(*C11);
961  // Block-row 2
962  Teuchos::RCP<Vector<Real> > C12 = (ahwvs.get_2())->clone();
963  Teuchos::RCP<Vector<Real> > C22 = (ahwvs.get_2())->clone();
964  applyAdjointHessian_12(*C12,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
965  applyAdjointHessian_22(*C22,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
966  C22->plus(*C12);
967  ahwvs.set_2(*C22);
968  }
969 
970 
971 
972  virtual Real checkSolve(const ROL::Vector<Real> &u,
973  const ROL::Vector<Real> &z,
974  const ROL::Vector<Real> &c,
975  const bool printToStream = true,
976  std::ostream & outStream = std::cout) {
977  // Solve equality constraint for u.
978  Real tol = ROL_EPSILON<Real>();
979  Teuchos::RCP<ROL::Vector<Real> > r = c.clone();
980  Teuchos::RCP<ROL::Vector<Real> > s = u.clone();
981  solve(*r,*s,z,tol);
982  // Evaluate equality constraint residual at (u,z).
983  Teuchos::RCP<ROL::Vector<Real> > cs = c.clone();
984  update(*s,z);
985  value(*cs,*s,z,tol);
986  // Output norm of residual.
987  Real rnorm = r->norm();
988  Real cnorm = cs->norm();
989  if ( printToStream ) {
990  std::stringstream hist;
991  hist << std::scientific << std::setprecision(8);
992  hist << "\nTest SimOpt solve at feasible (u,z):\n";
993  hist << " Solver Residual = " << rnorm << "\n";
994  hist << " ||c(u,z)|| = " << cnorm << "\n";
995  outStream << hist.str();
996  }
997  return cnorm;
998  }
999 
1000 
1013  virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w,
1014  const Vector<Real> &v,
1015  const Vector<Real> &u,
1016  const Vector<Real> &z,
1017  const bool printToStream = true,
1018  std::ostream & outStream = std::cout) {
1019  return checkAdjointConsistencyJacobian_1(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
1020  }
1021 
1022 
1038  virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w,
1039  const Vector<Real> &v,
1040  const Vector<Real> &u,
1041  const Vector<Real> &z,
1042  const Vector<Real> &dualw,
1043  const Vector<Real> &dualv,
1044  const bool printToStream = true,
1045  std::ostream & outStream = std::cout) {
1046  Real tol = ROL_EPSILON<Real>();
1047  Teuchos::RCP<Vector<Real> > Jv = dualw.clone();
1048  update(u,z);
1049  applyJacobian_1(*Jv,v,u,z,tol);
1050  Real wJv = w.dot(Jv->dual());
1051  Teuchos::RCP<Vector<Real> > Jw = dualv.clone();
1052  update(u,z);
1053  applyAdjointJacobian_1(*Jw,w,u,z,tol);
1054  Real vJw = v.dot(Jw->dual());
1055  Real diff = std::abs(wJv-vJw);
1056  if ( printToStream ) {
1057  std::stringstream hist;
1058  hist << std::scientific << std::setprecision(8);
1059  hist << "\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1060  << diff << "\n";
1061  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1062  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1063  outStream << hist.str();
1064  }
1065  return diff;
1066  }
1067 
1068 
1081  virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w,
1082  const Vector<Real> &v,
1083  const Vector<Real> &u,
1084  const Vector<Real> &z,
1085  const bool printToStream = true,
1086  std::ostream & outStream = std::cout) {
1087  return checkAdjointConsistencyJacobian_2(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
1088  }
1089 
1105  virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w,
1106  const Vector<Real> &v,
1107  const Vector<Real> &u,
1108  const Vector<Real> &z,
1109  const Vector<Real> &dualw,
1110  const Vector<Real> &dualv,
1111  const bool printToStream = true,
1112  std::ostream & outStream = std::cout) {
1113  Real tol = ROL_EPSILON<Real>();
1114  Teuchos::RCP<Vector<Real> > Jv = dualw.clone();
1115  update(u,z);
1116  applyJacobian_2(*Jv,v,u,z,tol);
1117  Real wJv = w.dot(Jv->dual());
1118  Teuchos::RCP<Vector<Real> > Jw = dualv.clone();
1119  update(u,z);
1120  applyAdjointJacobian_2(*Jw,w,u,z,tol);
1121  Real vJw = v.dot(Jw->dual());
1122  Real diff = std::abs(wJv-vJw);
1123  if ( printToStream ) {
1124  std::stringstream hist;
1125  hist << std::scientific << std::setprecision(8);
1126  hist << "\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1127  << diff << "\n";
1128  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1129  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1130  outStream << hist.str();
1131  }
1132  return diff;
1133  }
1134 
1135  virtual Real checkInverseJacobian_1(const Vector<Real> &jv,
1136  const Vector<Real> &v,
1137  const Vector<Real> &u,
1138  const Vector<Real> &z,
1139  const bool printToStream = true,
1140  std::ostream & outStream = std::cout) {
1141  Real tol = ROL_EPSILON<Real>();
1142  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1143  update(u,z);
1144  applyJacobian_1(*Jv,v,u,z,tol);
1145  Teuchos::RCP<Vector<Real> > iJJv = u.clone();
1146  update(u,z);
1147  applyInverseJacobian_1(*iJJv,*Jv,u,z,tol);
1148  Teuchos::RCP<Vector<Real> > diff = v.clone();
1149  diff->set(v);
1150  diff->axpy(-1.0,*iJJv);
1151  Real dnorm = diff->norm();
1152  Real vnorm = v.norm();
1153  if ( printToStream ) {
1154  std::stringstream hist;
1155  hist << std::scientific << std::setprecision(8);
1156  hist << "\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = "
1157  << dnorm << "\n";
1158  hist << " ||v|| = " << vnorm << "\n";
1159  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1160  outStream << hist.str();
1161  }
1162  return dnorm;
1163  }
1164 
1165  virtual Real checkInverseAdjointJacobian_1(const Vector<Real> &jv,
1166  const Vector<Real> &v,
1167  const Vector<Real> &u,
1168  const Vector<Real> &z,
1169  const bool printToStream = true,
1170  std::ostream & outStream = std::cout) {
1171  Real tol = ROL_EPSILON<Real>();
1172  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1173  update(u,z);
1174  applyAdjointJacobian_1(*Jv,v,u,z,tol);
1175  Teuchos::RCP<Vector<Real> > iJJv = v.clone();
1176  update(u,z);
1177  applyInverseAdjointJacobian_1(*iJJv,*Jv,u,z,tol);
1178  Teuchos::RCP<Vector<Real> > diff = v.clone();
1179  diff->set(v);
1180  diff->axpy(-1.0,*iJJv);
1181  Real dnorm = diff->norm();
1182  Real vnorm = v.norm();
1183  if ( printToStream ) {
1184  std::stringstream hist;
1185  hist << std::scientific << std::setprecision(8);
1186  hist << "\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = "
1187  << dnorm << "\n";
1188  hist << " ||v|| = " << vnorm << "\n";
1189  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1190  outStream << hist.str();
1191  }
1192  return dnorm;
1193  }
1194 
1195 
1196 
1197  std::vector<std::vector<Real> > checkApplyJacobian_1(const Vector<Real> &u,
1198  const Vector<Real> &z,
1199  const Vector<Real> &v,
1200  const Vector<Real> &jv,
1201  const bool printToStream = true,
1202  std::ostream & outStream = std::cout,
1203  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1204  const int order = 1) {
1205  std::vector<Real> steps(numSteps);
1206  for(int i=0;i<numSteps;++i) {
1207  steps[i] = pow(10,-i);
1208  }
1209 
1210  return checkApplyJacobian_1(u,z,v,jv,steps,printToStream,outStream,order);
1211  }
1212 
1213 
1214 
1215 
1216  std::vector<std::vector<Real> > checkApplyJacobian_1(const Vector<Real> &u,
1217  const Vector<Real> &z,
1218  const Vector<Real> &v,
1219  const Vector<Real> &jv,
1220  const std::vector<Real> &steps,
1221  const bool printToStream = true,
1222  std::ostream & outStream = std::cout,
1223  const int order = 1) {
1224 
1225  TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1226  "Error: finite difference order must be 1,2,3, or 4" );
1227 
1228  Real one(1.0);
1229 
1230  using Finite_Difference_Arrays::shifts;
1231  using Finite_Difference_Arrays::weights;
1232 
1233  Real tol = std::sqrt(ROL_EPSILON<Real>());
1234 
1235  int numSteps = steps.size();
1236  int numVals = 4;
1237  std::vector<Real> tmp(numVals);
1238  std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
1239 
1240  // Save the format state of the original outStream.
1241  Teuchos::oblackholestream oldFormatState;
1242  oldFormatState.copyfmt(outStream);
1243 
1244  // Compute constraint value at x.
1245  Teuchos::RCP<Vector<Real> > c = jv.clone();
1246  this->update(u,z);
1247  this->value(*c, u, z, tol);
1248 
1249  // Compute (Jacobian at x) times (vector v).
1250  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1251  this->applyJacobian_1(*Jv, v, u, z, tol);
1252  Real normJv = Jv->norm();
1253 
1254  // Temporary vectors.
1255  Teuchos::RCP<Vector<Real> > cdif = jv.clone();
1256  Teuchos::RCP<Vector<Real> > cnew = jv.clone();
1257  Teuchos::RCP<Vector<Real> > unew = u.clone();
1258 
1259  for (int i=0; i<numSteps; i++) {
1260 
1261  Real eta = steps[i];
1262 
1263  unew->set(u);
1264 
1265  cdif->set(*c);
1266  cdif->scale(weights[order-1][0]);
1267 
1268  for(int j=0; j<order; ++j) {
1269 
1270  unew->axpy(eta*shifts[order-1][j], v);
1271 
1272  if( weights[order-1][j+1] != 0 ) {
1273  this->update(*unew,z);
1274  this->value(*cnew,*unew,z,tol);
1275  cdif->axpy(weights[order-1][j+1],*cnew);
1276  }
1277 
1278  }
1279 
1280  cdif->scale(one/eta);
1281 
1282  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1283  jvCheck[i][0] = eta;
1284  jvCheck[i][1] = normJv;
1285  jvCheck[i][2] = cdif->norm();
1286  cdif->axpy(-one, *Jv);
1287  jvCheck[i][3] = cdif->norm();
1288 
1289  if (printToStream) {
1290  std::stringstream hist;
1291  if (i==0) {
1292  hist << std::right
1293  << std::setw(20) << "Step size"
1294  << std::setw(20) << "norm(Jac*vec)"
1295  << std::setw(20) << "norm(FD approx)"
1296  << std::setw(20) << "norm(abs error)"
1297  << "\n"
1298  << std::setw(20) << "---------"
1299  << std::setw(20) << "-------------"
1300  << std::setw(20) << "---------------"
1301  << std::setw(20) << "---------------"
1302  << "\n";
1303  }
1304  hist << std::scientific << std::setprecision(11) << std::right
1305  << std::setw(20) << jvCheck[i][0]
1306  << std::setw(20) << jvCheck[i][1]
1307  << std::setw(20) << jvCheck[i][2]
1308  << std::setw(20) << jvCheck[i][3]
1309  << "\n";
1310  outStream << hist.str();
1311  }
1312 
1313  }
1314 
1315  // Reset format state of outStream.
1316  outStream.copyfmt(oldFormatState);
1317 
1318  return jvCheck;
1319  } // checkApplyJacobian
1320 
1321 
1322  std::vector<std::vector<Real> > checkApplyJacobian_2(const Vector<Real> &u,
1323  const Vector<Real> &z,
1324  const Vector<Real> &v,
1325  const Vector<Real> &jv,
1326  const bool printToStream = true,
1327  std::ostream & outStream = std::cout,
1328  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1329  const int order = 1) {
1330  std::vector<Real> steps(numSteps);
1331  for(int i=0;i<numSteps;++i) {
1332  steps[i] = pow(10,-i);
1333  }
1334 
1335  return checkApplyJacobian_2(u,z,v,jv,steps,printToStream,outStream,order);
1336  }
1337 
1338 
1339 
1340 
1341  std::vector<std::vector<Real> > checkApplyJacobian_2(const Vector<Real> &u,
1342  const Vector<Real> &z,
1343  const Vector<Real> &v,
1344  const Vector<Real> &jv,
1345  const std::vector<Real> &steps,
1346  const bool printToStream = true,
1347  std::ostream & outStream = std::cout,
1348  const int order = 1) {
1349 
1350  TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1351  "Error: finite difference order must be 1,2,3, or 4" );
1352 
1353  Real one(1.0);
1354 
1355  using Finite_Difference_Arrays::shifts;
1356  using Finite_Difference_Arrays::weights;
1357 
1358  Real tol = std::sqrt(ROL_EPSILON<Real>());
1359 
1360  int numSteps = steps.size();
1361  int numVals = 4;
1362  std::vector<Real> tmp(numVals);
1363  std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
1364 
1365  // Save the format state of the original outStream.
1366  Teuchos::oblackholestream oldFormatState;
1367  oldFormatState.copyfmt(outStream);
1368 
1369  // Compute constraint value at x.
1370  Teuchos::RCP<Vector<Real> > c = jv.clone();
1371  this->update(u,z);
1372  this->value(*c, u, z, tol);
1373 
1374  // Compute (Jacobian at x) times (vector v).
1375  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1376  this->applyJacobian_2(*Jv, v, u, z, tol);
1377  Real normJv = Jv->norm();
1378 
1379  // Temporary vectors.
1380  Teuchos::RCP<Vector<Real> > cdif = jv.clone();
1381  Teuchos::RCP<Vector<Real> > cnew = jv.clone();
1382  Teuchos::RCP<Vector<Real> > znew = z.clone();
1383 
1384  for (int i=0; i<numSteps; i++) {
1385 
1386  Real eta = steps[i];
1387 
1388  znew->set(z);
1389 
1390  cdif->set(*c);
1391  cdif->scale(weights[order-1][0]);
1392 
1393  for(int j=0; j<order; ++j) {
1394 
1395  znew->axpy(eta*shifts[order-1][j], v);
1396 
1397  if( weights[order-1][j+1] != 0 ) {
1398  this->update(u,*znew);
1399  this->value(*cnew,u,*znew,tol);
1400  cdif->axpy(weights[order-1][j+1],*cnew);
1401  }
1402 
1403  }
1404 
1405  cdif->scale(one/eta);
1406 
1407  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1408  jvCheck[i][0] = eta;
1409  jvCheck[i][1] = normJv;
1410  jvCheck[i][2] = cdif->norm();
1411  cdif->axpy(-one, *Jv);
1412  jvCheck[i][3] = cdif->norm();
1413 
1414  if (printToStream) {
1415  std::stringstream hist;
1416  if (i==0) {
1417  hist << std::right
1418  << std::setw(20) << "Step size"
1419  << std::setw(20) << "norm(Jac*vec)"
1420  << std::setw(20) << "norm(FD approx)"
1421  << std::setw(20) << "norm(abs error)"
1422  << "\n"
1423  << std::setw(20) << "---------"
1424  << std::setw(20) << "-------------"
1425  << std::setw(20) << "---------------"
1426  << std::setw(20) << "---------------"
1427  << "\n";
1428  }
1429  hist << std::scientific << std::setprecision(11) << std::right
1430  << std::setw(20) << jvCheck[i][0]
1431  << std::setw(20) << jvCheck[i][1]
1432  << std::setw(20) << jvCheck[i][2]
1433  << std::setw(20) << jvCheck[i][3]
1434  << "\n";
1435  outStream << hist.str();
1436  }
1437 
1438  }
1439 
1440  // Reset format state of outStream.
1441  outStream.copyfmt(oldFormatState);
1442 
1443  return jvCheck;
1444  } // checkApplyJacobian
1445 
1446 
1447 
1448  std::vector<std::vector<Real> > checkApplyAdjointHessian_11(const Vector<Real> &u,
1449  const Vector<Real> &z,
1450  const Vector<Real> &p,
1451  const Vector<Real> &v,
1452  const Vector<Real> &hv,
1453  const bool printToStream = true,
1454  std::ostream & outStream = std::cout,
1455  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1456  const int order = 1 ) {
1457  std::vector<Real> steps(numSteps);
1458  for(int i=0;i<numSteps;++i) {
1459  steps[i] = pow(10,-i);
1460  }
1461 
1462  return checkApplyAdjointHessian_11(u,z,p,v,hv,steps,printToStream,outStream,order);
1463 
1464  }
1465 
1466  std::vector<std::vector<Real> > checkApplyAdjointHessian_11(const Vector<Real> &u,
1467  const Vector<Real> &z,
1468  const Vector<Real> &p,
1469  const Vector<Real> &v,
1470  const Vector<Real> &hv,
1471  const std::vector<Real> &steps,
1472  const bool printToStream = true,
1473  std::ostream & outStream = std::cout,
1474  const int order = 1 ) {
1475  using Finite_Difference_Arrays::shifts;
1476  using Finite_Difference_Arrays::weights;
1477 
1478  Real one(1.0);
1479 
1480  Real tol = std::sqrt(ROL_EPSILON<Real>());
1481 
1482  int numSteps = steps.size();
1483  int numVals = 4;
1484  std::vector<Real> tmp(numVals);
1485  std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1486 
1487  // Temporary vectors.
1488  Teuchos::RCP<Vector<Real> > AJdif = hv.clone();
1489  Teuchos::RCP<Vector<Real> > AJp = hv.clone();
1490  Teuchos::RCP<Vector<Real> > AHpv = hv.clone();
1491  Teuchos::RCP<Vector<Real> > AJnew = hv.clone();
1492  Teuchos::RCP<Vector<Real> > unew = u.clone();
1493 
1494  // Save the format state of the original outStream.
1495  Teuchos::oblackholestream oldFormatState;
1496  oldFormatState.copyfmt(outStream);
1497 
1498  // Apply adjoint Jacobian to p.
1499  this->update(u,z);
1500  this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1501 
1502  // Apply adjoint Hessian at (u,z), in direction v, to p.
1503  this->applyAdjointHessian_11(*AHpv, p, v, u, z, tol);
1504  Real normAHpv = AHpv->norm();
1505 
1506  for (int i=0; i<numSteps; i++) {
1507 
1508  Real eta = steps[i];
1509 
1510  // Apply adjoint Jacobian to p at (u+eta*v,z).
1511  unew->set(u);
1512 
1513  AJdif->set(*AJp);
1514  AJdif->scale(weights[order-1][0]);
1515 
1516  for(int j=0; j<order; ++j) {
1517 
1518  unew->axpy(eta*shifts[order-1][j],v);
1519 
1520  if( weights[order-1][j+1] != 0 ) {
1521  this->update(*unew,z);
1522  this->applyAdjointJacobian_1(*AJnew, p, *unew, z, tol);
1523  AJdif->axpy(weights[order-1][j+1],*AJnew);
1524  }
1525  }
1526 
1527  AJdif->scale(one/eta);
1528 
1529  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1530  ahpvCheck[i][0] = eta;
1531  ahpvCheck[i][1] = normAHpv;
1532  ahpvCheck[i][2] = AJdif->norm();
1533  AJdif->axpy(-one, *AHpv);
1534  ahpvCheck[i][3] = AJdif->norm();
1535 
1536  if (printToStream) {
1537  std::stringstream hist;
1538  if (i==0) {
1539  hist << std::right
1540  << std::setw(20) << "Step size"
1541  << std::setw(20) << "norm(adj(H)(u,v))"
1542  << std::setw(20) << "norm(FD approx)"
1543  << std::setw(20) << "norm(abs error)"
1544  << "\n"
1545  << std::setw(20) << "---------"
1546  << std::setw(20) << "-----------------"
1547  << std::setw(20) << "---------------"
1548  << std::setw(20) << "---------------"
1549  << "\n";
1550  }
1551  hist << std::scientific << std::setprecision(11) << std::right
1552  << std::setw(20) << ahpvCheck[i][0]
1553  << std::setw(20) << ahpvCheck[i][1]
1554  << std::setw(20) << ahpvCheck[i][2]
1555  << std::setw(20) << ahpvCheck[i][3]
1556  << "\n";
1557  outStream << hist.str();
1558  }
1559 
1560  }
1561 
1562  // Reset format state of outStream.
1563  outStream.copyfmt(oldFormatState);
1564 
1565  return ahpvCheck;
1566  } // checkApplyAdjointHessian_11
1567 
1568  std::vector<std::vector<Real> > checkApplyAdjointHessian_21(const Vector<Real> &u,
1569  const Vector<Real> &z,
1570  const Vector<Real> &p,
1571  const Vector<Real> &v,
1572  const Vector<Real> &hv,
1573  const bool printToStream = true,
1574  std::ostream & outStream = std::cout,
1575  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1576  const int order = 1 ) {
1577  std::vector<Real> steps(numSteps);
1578  for(int i=0;i<numSteps;++i) {
1579  steps[i] = pow(10,-i);
1580  }
1581 
1582  return checkApplyAdjointHessian_21(u,z,p,v,hv,steps,printToStream,outStream,order);
1583 
1584  }
1585 
1586  std::vector<std::vector<Real> > checkApplyAdjointHessian_21(const Vector<Real> &u,
1587  const Vector<Real> &z,
1588  const Vector<Real> &p,
1589  const Vector<Real> &v,
1590  const Vector<Real> &hv,
1591  const std::vector<Real> &steps,
1592  const bool printToStream = true,
1593  std::ostream & outStream = std::cout,
1594  const int order = 1 ) {
1595  using Finite_Difference_Arrays::shifts;
1596  using Finite_Difference_Arrays::weights;
1597 
1598  Real one(1.0);
1599 
1600  Real tol = std::sqrt(ROL_EPSILON<Real>());
1601 
1602  int numSteps = steps.size();
1603  int numVals = 4;
1604  std::vector<Real> tmp(numVals);
1605  std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1606 
1607  // Temporary vectors.
1608  Teuchos::RCP<Vector<Real> > AJdif = hv.clone();
1609  Teuchos::RCP<Vector<Real> > AJp = hv.clone();
1610  Teuchos::RCP<Vector<Real> > AHpv = hv.clone();
1611  Teuchos::RCP<Vector<Real> > AJnew = hv.clone();
1612  Teuchos::RCP<Vector<Real> > znew = z.clone();
1613 
1614  // Save the format state of the original outStream.
1615  Teuchos::oblackholestream oldFormatState;
1616  oldFormatState.copyfmt(outStream);
1617 
1618  // Apply adjoint Jacobian to p.
1619  this->update(u,z);
1620  this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1621 
1622  // Apply adjoint Hessian at (u,z), in direction v, to p.
1623  this->applyAdjointHessian_21(*AHpv, p, v, u, z, tol);
1624  Real normAHpv = AHpv->norm();
1625 
1626  for (int i=0; i<numSteps; i++) {
1627 
1628  Real eta = steps[i];
1629 
1630  // Apply adjoint Jacobian to p at (u,z+eta*v).
1631  znew->set(z);
1632 
1633  AJdif->set(*AJp);
1634  AJdif->scale(weights[order-1][0]);
1635 
1636  for(int j=0; j<order; ++j) {
1637 
1638  znew->axpy(eta*shifts[order-1][j],v);
1639 
1640  if( weights[order-1][j+1] != 0 ) {
1641  this->update(u,*znew);
1642  this->applyAdjointJacobian_1(*AJnew, p, u, *znew, tol);
1643  AJdif->axpy(weights[order-1][j+1],*AJnew);
1644  }
1645  }
1646 
1647  AJdif->scale(one/eta);
1648 
1649  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1650  ahpvCheck[i][0] = eta;
1651  ahpvCheck[i][1] = normAHpv;
1652  ahpvCheck[i][2] = AJdif->norm();
1653  AJdif->axpy(-one, *AHpv);
1654  ahpvCheck[i][3] = AJdif->norm();
1655 
1656  if (printToStream) {
1657  std::stringstream hist;
1658  if (i==0) {
1659  hist << std::right
1660  << std::setw(20) << "Step size"
1661  << std::setw(20) << "norm(adj(H)(u,v))"
1662  << std::setw(20) << "norm(FD approx)"
1663  << std::setw(20) << "norm(abs error)"
1664  << "\n"
1665  << std::setw(20) << "---------"
1666  << std::setw(20) << "-----------------"
1667  << std::setw(20) << "---------------"
1668  << std::setw(20) << "---------------"
1669  << "\n";
1670  }
1671  hist << std::scientific << std::setprecision(11) << std::right
1672  << std::setw(20) << ahpvCheck[i][0]
1673  << std::setw(20) << ahpvCheck[i][1]
1674  << std::setw(20) << ahpvCheck[i][2]
1675  << std::setw(20) << ahpvCheck[i][3]
1676  << "\n";
1677  outStream << hist.str();
1678  }
1679 
1680  }
1681 
1682  // Reset format state of outStream.
1683  outStream.copyfmt(oldFormatState);
1684 
1685  return ahpvCheck;
1686  } // checkApplyAdjointHessian_21
1687 
1688  std::vector<std::vector<Real> > checkApplyAdjointHessian_12(const Vector<Real> &u,
1689  const Vector<Real> &z,
1690  const Vector<Real> &p,
1691  const Vector<Real> &v,
1692  const Vector<Real> &hv,
1693  const bool printToStream = true,
1694  std::ostream & outStream = std::cout,
1695  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1696  const int order = 1 ) {
1697  std::vector<Real> steps(numSteps);
1698  for(int i=0;i<numSteps;++i) {
1699  steps[i] = pow(10,-i);
1700  }
1701 
1702  return checkApplyAdjointHessian_12(u,z,p,v,hv,steps,printToStream,outStream,order);
1703 
1704  }
1705 
1706  std::vector<std::vector<Real> > checkApplyAdjointHessian_12(const Vector<Real> &u,
1707  const Vector<Real> &z,
1708  const Vector<Real> &p,
1709  const Vector<Real> &v,
1710  const Vector<Real> &hv,
1711  const std::vector<Real> &steps,
1712  const bool printToStream = true,
1713  std::ostream & outStream = std::cout,
1714  const int order = 1 ) {
1715  using Finite_Difference_Arrays::shifts;
1716  using Finite_Difference_Arrays::weights;
1717 
1718  Real one(1.0);
1719 
1720  Real tol = std::sqrt(ROL_EPSILON<Real>());
1721 
1722  int numSteps = steps.size();
1723  int numVals = 4;
1724  std::vector<Real> tmp(numVals);
1725  std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1726 
1727  // Temporary vectors.
1728  Teuchos::RCP<Vector<Real> > AJdif = hv.clone();
1729  Teuchos::RCP<Vector<Real> > AJp = hv.clone();
1730  Teuchos::RCP<Vector<Real> > AHpv = hv.clone();
1731  Teuchos::RCP<Vector<Real> > AJnew = hv.clone();
1732  Teuchos::RCP<Vector<Real> > unew = u.clone();
1733 
1734  // Save the format state of the original outStream.
1735  Teuchos::oblackholestream oldFormatState;
1736  oldFormatState.copyfmt(outStream);
1737 
1738  // Apply adjoint Jacobian to p.
1739  this->update(u,z);
1740  this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1741 
1742  // Apply adjoint Hessian at (u,z), in direction v, to p.
1743  this->applyAdjointHessian_12(*AHpv, p, v, u, z, tol);
1744  Real normAHpv = AHpv->norm();
1745 
1746  for (int i=0; i<numSteps; i++) {
1747 
1748  Real eta = steps[i];
1749 
1750  // Apply adjoint Jacobian to p at (u+eta*v,z).
1751  unew->set(u);
1752 
1753  AJdif->set(*AJp);
1754  AJdif->scale(weights[order-1][0]);
1755 
1756  for(int j=0; j<order; ++j) {
1757 
1758  unew->axpy(eta*shifts[order-1][j],v);
1759 
1760  if( weights[order-1][j+1] != 0 ) {
1761  this->update(*unew,z);
1762  this->applyAdjointJacobian_2(*AJnew, p, *unew, z, tol);
1763  AJdif->axpy(weights[order-1][j+1],*AJnew);
1764  }
1765  }
1766 
1767  AJdif->scale(one/eta);
1768 
1769  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1770  ahpvCheck[i][0] = eta;
1771  ahpvCheck[i][1] = normAHpv;
1772  ahpvCheck[i][2] = AJdif->norm();
1773  AJdif->axpy(-one, *AHpv);
1774  ahpvCheck[i][3] = AJdif->norm();
1775 
1776  if (printToStream) {
1777  std::stringstream hist;
1778  if (i==0) {
1779  hist << std::right
1780  << std::setw(20) << "Step size"
1781  << std::setw(20) << "norm(adj(H)(u,v))"
1782  << std::setw(20) << "norm(FD approx)"
1783  << std::setw(20) << "norm(abs error)"
1784  << "\n"
1785  << std::setw(20) << "---------"
1786  << std::setw(20) << "-----------------"
1787  << std::setw(20) << "---------------"
1788  << std::setw(20) << "---------------"
1789  << "\n";
1790  }
1791  hist << std::scientific << std::setprecision(11) << std::right
1792  << std::setw(20) << ahpvCheck[i][0]
1793  << std::setw(20) << ahpvCheck[i][1]
1794  << std::setw(20) << ahpvCheck[i][2]
1795  << std::setw(20) << ahpvCheck[i][3]
1796  << "\n";
1797  outStream << hist.str();
1798  }
1799 
1800  }
1801 
1802  // Reset format state of outStream.
1803  outStream.copyfmt(oldFormatState);
1804 
1805  return ahpvCheck;
1806  } // checkApplyAdjointHessian_12
1807 
1808  std::vector<std::vector<Real> > checkApplyAdjointHessian_22(const Vector<Real> &u,
1809  const Vector<Real> &z,
1810  const Vector<Real> &p,
1811  const Vector<Real> &v,
1812  const Vector<Real> &hv,
1813  const bool printToStream = true,
1814  std::ostream & outStream = std::cout,
1815  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1816  const int order = 1 ) {
1817  std::vector<Real> steps(numSteps);
1818  for(int i=0;i<numSteps;++i) {
1819  steps[i] = pow(10,-i);
1820  }
1821 
1822  return checkApplyAdjointHessian_22(u,z,p,v,hv,steps,printToStream,outStream,order);
1823 
1824  }
1825 
1826  std::vector<std::vector<Real> > checkApplyAdjointHessian_22(const Vector<Real> &u,
1827  const Vector<Real> &z,
1828  const Vector<Real> &p,
1829  const Vector<Real> &v,
1830  const Vector<Real> &hv,
1831  const std::vector<Real> &steps,
1832  const bool printToStream = true,
1833  std::ostream & outStream = std::cout,
1834  const int order = 1 ) {
1835  using Finite_Difference_Arrays::shifts;
1836  using Finite_Difference_Arrays::weights;
1837 
1838  Real one(1.0);
1839 
1840  Real tol = std::sqrt(ROL_EPSILON<Real>());
1841 
1842  int numSteps = steps.size();
1843  int numVals = 4;
1844  std::vector<Real> tmp(numVals);
1845  std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1846 
1847  // Temporary vectors.
1848  Teuchos::RCP<Vector<Real> > AJdif = hv.clone();
1849  Teuchos::RCP<Vector<Real> > AJp = hv.clone();
1850  Teuchos::RCP<Vector<Real> > AHpv = hv.clone();
1851  Teuchos::RCP<Vector<Real> > AJnew = hv.clone();
1852  Teuchos::RCP<Vector<Real> > znew = z.clone();
1853 
1854  // Save the format state of the original outStream.
1855  Teuchos::oblackholestream oldFormatState;
1856  oldFormatState.copyfmt(outStream);
1857 
1858  // Apply adjoint Jacobian to p.
1859  this->update(u,z);
1860  this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1861 
1862  // Apply adjoint Hessian at (u,z), in direction v, to p.
1863  this->applyAdjointHessian_22(*AHpv, p, v, u, z, tol);
1864  Real normAHpv = AHpv->norm();
1865 
1866  for (int i=0; i<numSteps; i++) {
1867 
1868  Real eta = steps[i];
1869 
1870  // Apply adjoint Jacobian to p at (u,z+eta*v).
1871  znew->set(z);
1872 
1873  AJdif->set(*AJp);
1874  AJdif->scale(weights[order-1][0]);
1875 
1876  for(int j=0; j<order; ++j) {
1877 
1878  znew->axpy(eta*shifts[order-1][j],v);
1879 
1880  if( weights[order-1][j+1] != 0 ) {
1881  this->update(u,*znew);
1882  this->applyAdjointJacobian_2(*AJnew, p, u, *znew, tol);
1883  AJdif->axpy(weights[order-1][j+1],*AJnew);
1884  }
1885  }
1886 
1887  AJdif->scale(one/eta);
1888 
1889  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1890  ahpvCheck[i][0] = eta;
1891  ahpvCheck[i][1] = normAHpv;
1892  ahpvCheck[i][2] = AJdif->norm();
1893  AJdif->axpy(-one, *AHpv);
1894  ahpvCheck[i][3] = AJdif->norm();
1895 
1896  if (printToStream) {
1897  std::stringstream hist;
1898  if (i==0) {
1899  hist << std::right
1900  << std::setw(20) << "Step size"
1901  << std::setw(20) << "norm(adj(H)(u,v))"
1902  << std::setw(20) << "norm(FD approx)"
1903  << std::setw(20) << "norm(abs error)"
1904  << "\n"
1905  << std::setw(20) << "---------"
1906  << std::setw(20) << "-----------------"
1907  << std::setw(20) << "---------------"
1908  << std::setw(20) << "---------------"
1909  << "\n";
1910  }
1911  hist << std::scientific << std::setprecision(11) << std::right
1912  << std::setw(20) << ahpvCheck[i][0]
1913  << std::setw(20) << ahpvCheck[i][1]
1914  << std::setw(20) << ahpvCheck[i][2]
1915  << std::setw(20) << ahpvCheck[i][3]
1916  << "\n";
1917  outStream << hist.str();
1918  }
1919 
1920  }
1921 
1922  // Reset format state of outStream.
1923  outStream.copyfmt(oldFormatState);
1924 
1925  return ahpvCheck;
1926  } // checkApplyAdjointHessian_22
1927 
1928 }; // class EqualityConstraint_SimOpt
1929 
1930 } // namespace ROL
1931 
1932 #endif
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1688
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1135
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:557
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:215
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:460
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:531
ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Vector_SimOpt::get_1
Teuchos::RCP< const Vector< Real > > get_1() const
Definition: ROL_Vector_SimOpt.hpp:158
ROL::Vector::dimension
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:185
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:394
ROL::EqualityConstraint_SimOpt::setSolveParameters
virtual void setSolveParameters(Teuchos::ParameterList &parlist)
Set solve parameters.
Definition: ROL_EqualityConstraint_SimOpt.hpp:323
ROL_Algorithm.hpp
ROL::EqualityConstraint_SimOpt::DEFAULT_atol_
const Real DEFAULT_atol_
Definition: ROL_EqualityConstraint_SimOpt.hpp:107
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1448
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Finite_Difference_Arrays::weights
const double weights[4][5]
Definition: ROL_Types.hpp:1160
ROL::EqualityConstraint_SimOpt::solve
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
Definition: ROL_EqualityConstraint_SimOpt.hpp:202
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:145
ROL_Objective_FSsolver.hpp
ROL::Vector_SimOpt
Defines the linear algebra or vector space interface for simulation-based optimization.
Definition: ROL_Vector_SimOpt.hpp:57
ROL::EqualityConstraint_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::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1341
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:601
ROL::EqualityConstraint_SimOpt::DEFAULT_maxit_
const int DEFAULT_maxit_
Definition: ROL_EqualityConstraint_SimOpt.hpp:112
ROL::Objective_FSsolver
Definition: ROL_Objective_FSsolver.hpp:52
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:627
ROL::EqualityConstraint_SimOpt::atol_
Real atol_
Definition: ROL_EqualityConstraint_SimOpt.hpp:118
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:761
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1808
ROL::EqualityConstraint_SimOpt::DEFAULT_factor_
const Real DEFAULT_factor_
Definition: ROL_EqualityConstraint_SimOpt.hpp:110
ROL::Finite_Difference_Arrays::shifts
const int shifts[4][4]
Definition: ROL_Types.hpp:1153
ROL::Vector_SimOpt::set_1
void set_1(const Vector< Real > &vec)
Definition: ROL_Vector_SimOpt.hpp:174
ROL::EqualityConstraint_SimOpt::maxit_
int maxit_
Definition: ROL_EqualityConstraint_SimOpt.hpp:123
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:717
ROL_EqualityConstraint_State.hpp
ROL_EqualityConstraint.hpp
ROL::Vector::clone
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::EqualityConstraint_SimOpt::firstSolve_
bool firstSolve_
Definition: ROL_EqualityConstraint_SimOpt.hpp:129
ROL::NonlinearLeastSquaresObjective_SimOpt
Provides the interface to evaluate nonlinear least squares objective functions.
Definition: ROL_NonlinearLeastSquaresObjective_SimOpt.hpp:76
ROL::EqualityConstraint_SimOpt::solverType_
int solverType_
Definition: ROL_EqualityConstraint_SimOpt.hpp:126
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:943
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:159
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:76
ROL::EqualityConstraint_SimOpt::jv_
Teuchos::RCP< Vector< Real > > jv_
Definition: ROL_EqualityConstraint_SimOpt.hpp:104
ROL::EqualityConstraint_SimOpt::DEFAULT_print_
const bool DEFAULT_print_
Definition: ROL_EqualityConstraint_SimOpt.hpp:113
ROL::EqualityConstraint_SimOpt
Defines the equality constraint operator interface for simulation-based optimization.
Definition: ROL_EqualityConstraint_SimOpt.hpp:100
ROL::Vector::dot
virtual Real dot(const Vector &x) const =0
Compute where .
ROL::EqualityConstraint_SimOpt::DEFAULT_rtol_
const Real DEFAULT_rtol_
Definition: ROL_EqualityConstraint_SimOpt.hpp:108
ROL::EqualityConstraint_SimOpt::unew_
Teuchos::RCP< Vector< Real > > unew_
Definition: ROL_EqualityConstraint_SimOpt.hpp:103
ROL::EqualityConstraint_SimOpt::value
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
Definition: ROL_EqualityConstraint_SimOpt.hpp:902
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:854
ROL::EqualityConstraint_SimOpt::update
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable is ...
Definition: ROL_EqualityConstraint_SimOpt.hpp:153
ROL::EqualityConstraint::solveAugmentedSystem
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system where , , , , is an identity or Riesz operator...
Definition: ROL_EqualityConstraintDef.hpp:200
ROL::EqualityConstraint
Defines the equality constraint operator interface.
Definition: ROL_EqualityConstraint.hpp:88
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:972
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:436
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:672
ROL::EqualityConstraint_SimOpt::rtol_
Real rtol_
Definition: ROL_EqualityConstraint_SimOpt.hpp:119
ROL::EqualityConstraint_SimOpt::EqualityConstraint_SimOpt
EqualityConstraint_SimOpt()
Definition: ROL_EqualityConstraint_SimOpt.hpp:132
ROL::EqualityConstraint_SimOpt::DEFAULT_stol_
const Real DEFAULT_stol_
Definition: ROL_EqualityConstraint_SimOpt.hpp:109
ROL::EqualityConstraint_SimOpt::decr_
Real decr_
Definition: ROL_EqualityConstraint_SimOpt.hpp:122
ROL::EqualityConstraint_SimOpt::isFeasible
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
Definition: ROL_EqualityConstraint_SimOpt.hpp:900
ROL::EqualityConstraint_SimOpt::zero_
bool zero_
Definition: ROL_EqualityConstraint_SimOpt.hpp:125
ROL::EqualityConstraint_SimOpt::print_
bool print_
Definition: ROL_EqualityConstraint_SimOpt.hpp:124
ROL::Algorithm
Provides an interface to run optimization algorithms.
Definition: ROL_Algorithm.hpp:72
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:351
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:926
ROL::EqualityConstraint_SimOpt::checkAdjointConsistencyJacobian_2
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualw, const Vector< Real > &dualv, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the secondary interface, for use with dual spaces where the user does not define the dual() operation.
Definition: ROL_EqualityConstraint_SimOpt.hpp:1105
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1322
ROL::EqualityConstraint_SimOpt::DEFAULT_decr_
const Real DEFAULT_decr_
Definition: ROL_EqualityConstraint_SimOpt.hpp:111
ROL::EqualityConstraint_SimOpt::DEFAULT_zero_
const bool DEFAULT_zero_
Definition: ROL_EqualityConstraint_SimOpt.hpp:114
ROL::EqualityConstraint::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_EqualityConstraint.hpp:264
ROL::EqualityConstraint_SimOpt::update_1
virtual void update_1(const Vector< Real > &u, bool flag=true, int iter=-1)
Update constraint functions with respect to Sim variable. x is the optimization variable, flag = true if optimization variable is changed, iter is the outer algorithm iterations count.
Definition: ROL_EqualityConstraint_SimOpt.hpp:163
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1586
ROL::EqualityConstraint_SimOpt::stol_
Real stol_
Definition: ROL_EqualityConstraint_SimOpt.hpp:120
ROL_NUM_CHECKDERIV_STEPS
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
Definition: ROL_Types.hpp:73
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1826
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:911
ROL::EqualityConstraint_SimOpt::update_2
virtual void update_2(const Vector< Real > &z, bool flag=true, int iter=-1)
Update constraint functions with respect to Opt variable. x is the optimization variable, flag = true if optimization variable is changed, iter is the outer algorithm iterations count.
Definition: ROL_EqualityConstraint_SimOpt.hpp:170
ROL::EqualityConstraint_SimOpt::checkAdjointConsistencyJacobian_1
virtual Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualw, const Vector< Real > &dualv, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the secondary interface, for use with dual spaces where the user does not define the dual() operation.
Definition: ROL_EqualityConstraint_SimOpt.hpp:1038
ROL::EqualityConstraint_SimOpt::update
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable is ...
Definition: ROL_EqualityConstraint_SimOpt.hpp:892
ROL::EqualityConstraint_SimOpt::factor_
Real factor_
Definition: ROL_EqualityConstraint_SimOpt.hpp:121
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1013
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:198
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:486
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1165
ROL::Vector::basis
virtual Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:174
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
ROL::EqualityConstraint_SimOpt::solveAugmentedSystem
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system where , , , , is an identity operator, and is a zero operator.
Definition: ROL_EqualityConstraint_SimOpt.hpp:826
ROL::EqualityConstraint_SimOpt::DEFAULT_solverType_
const int DEFAULT_solverType_
Definition: ROL_EqualityConstraint_SimOpt.hpp:115
ROL_NonlinearLeastSquaresObjective_SimOpt.hpp
ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:139
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1216
ROL::Vector_SimOpt::get_2
Teuchos::RCP< const Vector< Real > > get_2() const
Definition: ROL_Vector_SimOpt.hpp:162
ROL::EqualityConstraint_State
Definition: ROL_EqualityConstraint_State.hpp:52
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1466
ROL::Vector_SimOpt::set_2
void set_2(const Vector< Real > &vec)
Definition: ROL_Vector_SimOpt.hpp:178
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1706
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1197
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1081
ROL::EqualityConstraint_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_EqualityConstraint_SimOpt.hpp:1568
ROL_Vector_SimOpt.hpp
Generated by   doxygen 1.8.6