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 
48 #include "ROL_Vector_SimOpt.hpp"
49 #include "ROL_Types.hpp"
50 #include <iostream>
51 
92 namespace ROL {
93 
94 template <class Real>
95 class EqualityConstraint_SimOpt : public virtual EqualityConstraint<Real> {
96 private:
97  // Additional vector storage for solve
98  Teuchos::RCP<Vector<Real> > unew_;
99  Teuchos::RCP<Vector<Real> > jv_;
100 
101  // Default parameters for solve (backtracking Newton)
102  const Real DEFAULT_atol_;
103  const Real DEFAULT_rtol_;
104  const Real DEFAULT_stol_;
105  const Real DEFAULT_factor_;
106  const Real DEFAULT_decr_;
107  const int DEFAULT_maxit_;
108  const bool DEFAULT_print_;
109 
110  // User-set parameters for solve (backtracking Newton)
111  Real atol_;
112  Real rtol_;
113  Real stol_;
114  Real factor_;
115  Real decr_;
116  int maxit_;
117  bool print_;
118 
119  // Flag to initialize vector storage in solve
121 
122 public:
124  : EqualityConstraint<Real>(),
125  unew_(Teuchos::null), jv_(Teuchos::null),
126  DEFAULT_atol_(1.e-4*std::sqrt(ROL_EPSILON<Real>())),
127  DEFAULT_rtol_(1.e0),
128  DEFAULT_stol_(std::sqrt(ROL_EPSILON<Real>())),
129  DEFAULT_factor_(0.5),
130  DEFAULT_decr_(1.e-4),
131  DEFAULT_maxit_(500),
132  DEFAULT_print_(false),
135  firstSolve_(true) {}
136 
142  virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) {
143  update_1(u,flag,iter);
144  update_2(z,flag,iter);
145  }
146 
152  virtual void update_1( const Vector<Real> &u, bool flag = true, int iter = -1 ) {}
153 
159  virtual void update_2( const Vector<Real> &z, bool flag = true, int iter = -1 ) {}
160 
161 
175  virtual void value(Vector<Real> &c,
176  const Vector<Real> &u,
177  const Vector<Real> &z,
178  Real &tol) = 0;
179 
191  virtual void solve(Vector<Real> &c,
192  Vector<Real> &u,
193  const Vector<Real> &z,
194  Real &tol) {
195  if ( firstSolve_ ) {
196  unew_ = u.clone();
197  jv_ = u.clone();
198  firstSolve_ = false;
199  }
200  update(u,z,true);
201  value(c,u,z,tol);
202  Real cnorm = c.norm(), alpha = 1.0, tmp = 0.0;
203  const Real ctol = std::min(atol_, rtol_*cnorm);
204  int cnt = 0;
205  if ( print_ ) {
206  std::cout << "\n Default EqualityConstraint_SimOpt::solve\n";
207  std::cout << " ";
208  std::cout << std::setw(6) << std::left << "iter";
209  std::cout << std::setw(15) << std::left << "rnorm";
210  std::cout << std::setw(15) << std::left << "alpha";
211  std::cout << "\n";
212  }
213  for (cnt = 0; cnt < maxit_; ++cnt) {
214  // Compute Newton step
215  applyInverseJacobian_1(*jv_,c,u,z,tol);
216  unew_->set(u);
217  unew_->axpy(-alpha, *jv_);
218  update(*unew_,z);
219  value(c,*unew_,z,tol);
220  tmp = c.norm();
221  // Perform backtracking line search
222  while ( tmp > (1.0-decr_*alpha)*cnorm &&
223  alpha > stol_ ) {
224  alpha *= factor_;
225  unew_->set(u);
226  unew_->axpy(-alpha,*jv_);
227  update(*unew_,z);
228  value(c,*unew_,z,tol);
229  tmp = c.norm();
230  }
231  if ( print_ ) {
232  std::cout << " ";
233  std::cout << std::setw(6) << std::left << cnt;
234  std::cout << std::scientific << std::setprecision(6);
235  std::cout << std::setw(15) << std::left << tmp;
236  std::cout << std::scientific << std::setprecision(6);
237  std::cout << std::setw(15) << std::left << alpha;
238  std::cout << "\n";
239  }
240  // Update iterate
241  cnorm = tmp;
242  u.set(*unew_);
243  if (cnorm < ctol) {
244  break;
245  }
246  update(u,z,true);
247  alpha = 1.0;
248  }
249  }
250 
269  virtual void setSolveParameters(Teuchos::ParameterList &parlist) {
270  Teuchos::ParameterList & list = parlist.sublist("SimOpt").sublist("Solve");
271  atol_ = list.get("Absolute Residual Tolerance", DEFAULT_atol_);
272  rtol_ = list.get("Relative Residual Tolerance", DEFAULT_rtol_);
273  maxit_ = list.get("Iteration Limit", DEFAULT_maxit_);
274  decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_);
275  stol_ = list.get("Step Tolerance", DEFAULT_stol_);
276  factor_ = list.get("Backtracking Factor", DEFAULT_factor_);
277  print_ = list.get("Output Iteration History", DEFAULT_print_);
278  }
279 
295  virtual void applyJacobian_1(Vector<Real> &jv,
296  const Vector<Real> &v,
297  const Vector<Real> &u,
298  const Vector<Real> &z,
299  Real &tol) {
300  Real ctol = std::sqrt(ROL_EPSILON<Real>());
301  // Compute step length
302  Real h = tol;
303  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
304  h = std::max(1.0,u.norm()/v.norm())*tol;
305  }
306  // Update state vector to u + hv
307  Teuchos::RCP<Vector<Real> > unew = u.clone();
308  unew->set(u);
309  unew->axpy(h,v);
310  // Compute new constraint value
311  update(*unew,z);
312  value(jv,*unew,z,ctol);
313  // Compute current constraint value
314  Teuchos::RCP<Vector<Real> > cold = jv.clone();
315  update(u,z);
316  value(*cold,u,z,ctol);
317  // Compute Newton quotient
318  jv.axpy(-1.0,*cold);
319  jv.scale(1.0/h);
320  }
321 
322 
338  virtual void applyJacobian_2(Vector<Real> &jv,
339  const Vector<Real> &v,
340  const Vector<Real> &u,
341  const Vector<Real> &z,
342  Real &tol) {
343  Real ctol = std::sqrt(ROL_EPSILON<Real>());
344  // Compute step length
345  Real h = tol;
346  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
347  h = std::max(1.0,u.norm()/v.norm())*tol;
348  }
349  // Update state vector to u + hv
350  Teuchos::RCP<Vector<Real> > znew = z.clone();
351  znew->set(z);
352  znew->axpy(h,v);
353  // Compute new constraint value
354  update(u,*znew);
355  value(jv,u,*znew,ctol);
356  // Compute current constraint value
357  Teuchos::RCP<Vector<Real> > cold = jv.clone();
358  update(u,z);
359  value(*cold,u,z,ctol);
360  // Compute Newton quotient
361  jv.axpy(-1.0,*cold);
362  jv.scale(1.0/h);
363  }
364 
381  const Vector<Real> &v,
382  const Vector<Real> &u,
383  const Vector<Real> &z,
384  Real &tol) {
385  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
386  "The method applyInverseJacobian_1 is used but not implemented!\n");
387  }
388 
405  const Vector<Real> &v,
406  const Vector<Real> &u,
407  const Vector<Real> &z,
408  Real &tol) {
409  applyAdjointJacobian_1(ajv, v, u, z, v.dual(), tol);
410  }
411 
412 
431  const Vector<Real> &v,
432  const Vector<Real> &u,
433  const Vector<Real> &z,
434  const Vector<Real> &dualv,
435  Real &tol) {
436  Real ctol = std::sqrt(ROL_EPSILON<Real>());
437  Real h = tol;
438  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
439  h = std::max(1.0,u.norm()/v.norm())*tol;
440  }
441  Teuchos::RCP<Vector<Real> > cold = dualv.clone();
442  Teuchos::RCP<Vector<Real> > cnew = dualv.clone();
443  update(u,z);
444  value(*cold,u,z,ctol);
445  Teuchos::RCP<Vector<Real> > unew = u.clone();
446  ajv.zero();
447  for (int i = 0; i < u.dimension(); i++) {
448  unew->set(u);
449  unew->axpy(h,*(u.basis(i)));
450  update(*unew,z);
451  value(*cnew,*unew,z,ctol);
452  cnew->axpy(-1.0,*cold);
453  cnew->scale(1.0/h);
454  ajv.axpy(cnew->dot(v),*((u.dual()).basis(i)));
455  }
456  update(u,z);
457  }
458 
459 
476  const Vector<Real> &v,
477  const Vector<Real> &u,
478  const Vector<Real> &z,
479  Real &tol) {
480  applyAdjointJacobian_2(ajv, v, u, z, v.dual(), tol);
481  }
482 
483 
502  const Vector<Real> &v,
503  const Vector<Real> &u,
504  const Vector<Real> &z,
505  const Vector<Real> &dualv,
506  Real &tol) {
507  Real ctol = std::sqrt(ROL_EPSILON<Real>());
508  Real h = tol;
509  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
510  h = std::max(1.0,u.norm()/v.norm())*tol;
511  }
512  Teuchos::RCP<Vector<Real> > cold = dualv.clone();
513  Teuchos::RCP<Vector<Real> > cnew = dualv.clone();
514  update(u,z);
515  value(*cold,u,z,ctol);
516  Teuchos::RCP<Vector<Real> > znew = z.clone();
517  ajv.zero();
518  for (int i = 0; i < z.dimension(); i++) {
519  znew->set(z);
520  znew->axpy(h,*(z.basis(i)));
521  update(u,*znew);
522  value(*cnew,u,*znew,ctol);
523  cnew->axpy(-1.0,*cold);
524  cnew->scale(1.0/h);
525  ajv.axpy(cnew->dot(v),*((z.dual()).basis(i)));
526  }
527  update(u,z);
528  }
529 
546  const Vector<Real> &v,
547  const Vector<Real> &u,
548  const Vector<Real> &z,
549  Real &tol) {
550  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
551  "The method applyInverseAdjointJacobian_1 is used but not implemented!\n");
552  };
553 
572  const Vector<Real> &w,
573  const Vector<Real> &v,
574  const Vector<Real> &u,
575  const Vector<Real> &z,
576  Real &tol) {
577  Real jtol = std::sqrt(ROL_EPSILON<Real>());
578  // Compute step size
579  Real h = tol;
580  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
581  h = std::max(1.0,u.norm()/v.norm())*tol;
582  }
583  // Evaluate Jacobian at new state
584  Teuchos::RCP<Vector<Real> > unew = u.clone();
585  unew->set(u);
586  unew->axpy(h,v);
587  update(*unew,z);
588  applyAdjointJacobian_1(ahwv,w,*unew,z,jtol);
589  // Evaluate Jacobian at old state
590  Teuchos::RCP<Vector<Real> > jv = ahwv.clone();
591  update(u,z);
592  applyAdjointJacobian_1(*jv,w,u,z,jtol);
593  // Compute Newton quotient
594  ahwv.axpy(-1.0,*jv);
595  ahwv.scale(1.0/h);
596  }
597 
598 
617  const Vector<Real> &w,
618  const Vector<Real> &v,
619  const Vector<Real> &u,
620  const Vector<Real> &z,
621  Real &tol) {
622  Real jtol = std::sqrt(ROL_EPSILON<Real>());
623  // Compute step size
624  Real h = tol;
625  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
626  h = std::max(1.0,u.norm()/v.norm())*tol;
627  }
628  // Evaluate Jacobian at new state
629  Teuchos::RCP<Vector<Real> > unew = u.clone();
630  unew->set(u);
631  unew->axpy(h,v);
632  update(*unew,z);
633  applyAdjointJacobian_2(ahwv,w,*unew,z,jtol);
634  // Evaluate Jacobian at old state
635  Teuchos::RCP<Vector<Real> > jv = ahwv.clone();
636  update(u,z);
637  applyAdjointJacobian_2(*jv,w,u,z,jtol);
638  // Compute Newton quotient
639  ahwv.axpy(-1.0,*jv);
640  ahwv.scale(1.0/h);
641  }
642 
643 
662  const Vector<Real> &w,
663  const Vector<Real> &v,
664  const Vector<Real> &u,
665  const Vector<Real> &z,
666  Real &tol) {
667  Real jtol = std::sqrt(ROL_EPSILON<Real>());
668  // Compute step size
669  Real h = tol;
670  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
671  h = std::max(1.0,u.norm()/v.norm())*tol;
672  }
673  // Evaluate Jacobian at new control
674  Teuchos::RCP<Vector<Real> > znew = z.clone();
675  znew->set(z);
676  znew->axpy(h,v);
677  update(u,*znew);
678  applyAdjointJacobian_1(ahwv,w,u,*znew,jtol);
679  // Evaluate Jacobian at old control
680  Teuchos::RCP<Vector<Real> > jv = ahwv.clone();
681  update(u,z);
682  applyAdjointJacobian_1(*jv,w,u,z,jtol);
683  // Compute Newton quotient
684  ahwv.axpy(-1.0,*jv);
685  ahwv.scale(1.0/h);
686  }
687 
706  const Vector<Real> &w,
707  const Vector<Real> &v,
708  const Vector<Real> &u,
709  const Vector<Real> &z,
710  Real &tol) {
711  Real jtol = std::sqrt(ROL_EPSILON<Real>());
712  // Compute step size
713  Real h = tol;
714  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
715  h = std::max(1.0,u.norm()/v.norm())*tol;
716  }
717  // Evaluate Jacobian at new control
718  Teuchos::RCP<Vector<Real> > znew = z.clone();
719  znew->set(z);
720  znew->axpy(h,v);
721  update(u,*znew);
722  applyAdjointJacobian_2(ahwv,w,u,*znew,jtol);
723  // Evaluate Jacobian at old control
724  Teuchos::RCP<Vector<Real> > jv = ahwv.clone();
725  update(u,z);
726  applyAdjointJacobian_2(*jv,w,u,z,jtol);
727  // Compute Newton quotient
728  ahwv.axpy(-1.0,*jv);
729  ahwv.scale(1.0/h);
730 }
731 
770  virtual std::vector<Real> solveAugmentedSystem(Vector<Real> &v1,
771  Vector<Real> &v2,
772  const Vector<Real> &b1,
773  const Vector<Real> &b2,
774  const Vector<Real> &x,
775  Real &tol) {
776  return EqualityConstraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,x,tol);
777  }
778 
779 
799  const Vector<Real> &v,
800  const Vector<Real> &x,
801  const Vector<Real> &g,
802  Real &tol) {
803  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(x);
804  Teuchos::RCP<ROL::Vector<Real> > ijv = (xs.get_1())->clone();
805 
806  try {
807  applyInverseJacobian_1(*ijv, v, *(xs.get_1()), *(xs.get_2()), tol);
808  }
809  catch (const std::logic_error &e) {
811  return;
812  }
813 
814  const Vector_SimOpt<Real> &gs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(g);
815  Teuchos::RCP<ROL::Vector<Real> > ijv_dual = (gs.get_1())->clone();
816  ijv_dual->set(ijv->dual());
817 
818  try {
819  applyInverseAdjointJacobian_1(pv, *ijv_dual, *(xs.get_1()), *(xs.get_2()), tol);
820  }
821  catch (const std::logic_error &e) {
823  return;
824  }
825 
826  }
827 
828 //
829 // EqualityConstraint_SimOpt(void) : EqualityConstraint<Real>() {}
830 
836  virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
837  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
838  Teuchos::dyn_cast<const Vector<Real> >(x));
839  update(*(xs.get_1()),*(xs.get_2()),flag,iter);
840  }
841 
844  virtual bool isFeasible( const Vector<Real> &v ) { return true; }
845 
846  virtual void value(Vector<Real> &c,
847  const Vector<Real> &x,
848  Real &tol) {
849  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
850  Teuchos::dyn_cast<const Vector<Real> >(x));
851  value(c,*(xs.get_1()),*(xs.get_2()),tol);
852  }
853 
854 
855  virtual void applyJacobian(Vector<Real> &jv,
856  const Vector<Real> &v,
857  const Vector<Real> &x,
858  Real &tol) {
859  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
860  Teuchos::dyn_cast<const Vector<Real> >(x));
861  const Vector_SimOpt<Real> &vs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
862  Teuchos::dyn_cast<const Vector<Real> >(v));
863  applyJacobian_1(jv,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
864  Teuchos::RCP<Vector<Real> > jv2 = jv.clone();
865  applyJacobian_2(*jv2,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
866  jv.plus(*jv2);
867  }
868 
869 
871  const Vector<Real> &v,
872  const Vector<Real> &x,
873  Real &tol) {
874  Vector_SimOpt<Real> &ajvs = Teuchos::dyn_cast<Vector_SimOpt<Real> >(
875  Teuchos::dyn_cast<Vector<Real> >(ajv));
876  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
877  Teuchos::dyn_cast<const Vector<Real> >(x));
878  Teuchos::RCP<Vector<Real> > ajv1 = (ajvs.get_1())->clone();
879  applyAdjointJacobian_1(*ajv1,v,*(xs.get_1()),*(xs.get_2()),tol);
880  ajvs.set_1(*ajv1);
881  Teuchos::RCP<Vector<Real> > ajv2 = (ajvs.get_2())->clone();
882  applyAdjointJacobian_2(*ajv2,v,*(xs.get_1()),*(xs.get_2()),tol);
883  ajvs.set_2(*ajv2);
884  }
885 
886 
887  virtual void applyAdjointHessian(Vector<Real> &ahwv,
888  const Vector<Real> &w,
889  const Vector<Real> &v,
890  const Vector<Real> &x,
891  Real &tol) {
892  Vector_SimOpt<Real> &ahwvs = Teuchos::dyn_cast<Vector_SimOpt<Real> >(
893  Teuchos::dyn_cast<Vector<Real> >(ahwv));
894  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
895  Teuchos::dyn_cast<const Vector<Real> >(x));
896  const Vector_SimOpt<Real> &vs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
897  Teuchos::dyn_cast<const Vector<Real> >(v));
898  // Block-row 1
899  Teuchos::RCP<Vector<Real> > C11 = (ahwvs.get_1())->clone();
900  Teuchos::RCP<Vector<Real> > C21 = (ahwvs.get_1())->clone();
901  applyAdjointHessian_11(*C11,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
902  applyAdjointHessian_21(*C21,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
903  C11->plus(*C21);
904  ahwvs.set_1(*C11);
905  // Block-row 2
906  Teuchos::RCP<Vector<Real> > C12 = (ahwvs.get_2())->clone();
907  Teuchos::RCP<Vector<Real> > C22 = (ahwvs.get_2())->clone();
908  applyAdjointHessian_12(*C12,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
909  applyAdjointHessian_22(*C22,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
910  C22->plus(*C12);
911  ahwvs.set_2(*C22);
912  }
913 
914 
915 
916  virtual Real checkSolve(const ROL::Vector<Real> &u,
917  const ROL::Vector<Real> &z,
918  const ROL::Vector<Real> &c,
919  const bool printToStream = true,
920  std::ostream & outStream = std::cout) {
921  // Solve equality constraint for u.
922  Real tol = ROL_EPSILON<Real>();
923  Teuchos::RCP<ROL::Vector<Real> > r = c.clone();
924  Teuchos::RCP<ROL::Vector<Real> > s = u.clone();
925  solve(*r,*s,z,tol);
926  // Evaluate equality constraint residual at (u,z).
927  Teuchos::RCP<ROL::Vector<Real> > cs = c.clone();
928  update(*s,z);
929  value(*cs,*s,z,tol);
930  // Output norm of residual.
931  Real rnorm = r->norm();
932  Real cnorm = cs->norm();
933  if ( printToStream ) {
934  std::stringstream hist;
935  hist << std::scientific << std::setprecision(8);
936  hist << "\nTest SimOpt solve at feasible (u,z):\n";
937  hist << " Solver Residual = " << rnorm << "\n";
938  hist << " ||c(u,z)|| = " << cnorm << "\n";
939  outStream << hist.str();
940  }
941  return cnorm;
942  }
943 
944 
958  const Vector<Real> &v,
959  const Vector<Real> &u,
960  const Vector<Real> &z,
961  const bool printToStream = true,
962  std::ostream & outStream = std::cout) {
963  return checkAdjointConsistencyJacobian_1(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
964  }
965 
966 
983  const Vector<Real> &v,
984  const Vector<Real> &u,
985  const Vector<Real> &z,
986  const Vector<Real> &dualw,
987  const Vector<Real> &dualv,
988  const bool printToStream = true,
989  std::ostream & outStream = std::cout) {
990  Real tol = ROL_EPSILON<Real>();
991  Teuchos::RCP<Vector<Real> > Jv = dualw.clone();
992  update(u,z);
993  applyJacobian_1(*Jv,v,u,z,tol);
994  Real wJv = w.dot(Jv->dual());
995  Teuchos::RCP<Vector<Real> > Jw = dualv.clone();
996  update(u,z);
997  applyAdjointJacobian_1(*Jw,w,u,z,tol);
998  Real vJw = v.dot(Jw->dual());
999  Real diff = std::abs(wJv-vJw);
1000  if ( printToStream ) {
1001  std::stringstream hist;
1002  hist << std::scientific << std::setprecision(8);
1003  hist << "\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1004  << diff << "\n";
1005  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1006  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1007  outStream << hist.str();
1008  }
1009  return diff;
1010  }
1011 
1012 
1026  const Vector<Real> &v,
1027  const Vector<Real> &u,
1028  const Vector<Real> &z,
1029  const bool printToStream = true,
1030  std::ostream & outStream = std::cout) {
1031  return checkAdjointConsistencyJacobian_2(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
1032  }
1033 
1050  const Vector<Real> &v,
1051  const Vector<Real> &u,
1052  const Vector<Real> &z,
1053  const Vector<Real> &dualw,
1054  const Vector<Real> &dualv,
1055  const bool printToStream = true,
1056  std::ostream & outStream = std::cout) {
1057  Real tol = ROL_EPSILON<Real>();
1058  Teuchos::RCP<Vector<Real> > Jv = dualw.clone();
1059  update(u,z);
1060  applyJacobian_2(*Jv,v,u,z,tol);
1061  Real wJv = w.dot(Jv->dual());
1062  Teuchos::RCP<Vector<Real> > Jw = dualv.clone();
1063  update(u,z);
1064  applyAdjointJacobian_2(*Jw,w,u,z,tol);
1065  Real vJw = v.dot(Jw->dual());
1066  Real diff = std::abs(wJv-vJw);
1067  if ( printToStream ) {
1068  std::stringstream hist;
1069  hist << std::scientific << std::setprecision(8);
1070  hist << "\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1071  << diff << "\n";
1072  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1073  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1074  outStream << hist.str();
1075  }
1076  return diff;
1077  }
1078 
1079  virtual Real checkInverseJacobian_1(const Vector<Real> &jv,
1080  const Vector<Real> &v,
1081  const Vector<Real> &u,
1082  const Vector<Real> &z,
1083  const bool printToStream = true,
1084  std::ostream & outStream = std::cout) {
1085  Real tol = ROL_EPSILON<Real>();
1086  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1087  update(u,z);
1088  applyJacobian_1(*Jv,v,u,z,tol);
1089  Teuchos::RCP<Vector<Real> > iJJv = u.clone();
1090  update(u,z);
1091  applyInverseJacobian_1(*iJJv,*Jv,u,z,tol);
1092  Teuchos::RCP<Vector<Real> > diff = v.clone();
1093  diff->set(v);
1094  diff->axpy(-1.0,*iJJv);
1095  Real dnorm = diff->norm();
1096  Real vnorm = v.norm();
1097  if ( printToStream ) {
1098  std::stringstream hist;
1099  hist << std::scientific << std::setprecision(8);
1100  hist << "\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = "
1101  << dnorm << "\n";
1102  hist << " ||v|| = " << vnorm << "\n";
1103  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1104  outStream << hist.str();
1105  }
1106  return dnorm;
1107  }
1108 
1110  const Vector<Real> &v,
1111  const Vector<Real> &u,
1112  const Vector<Real> &z,
1113  const bool printToStream = true,
1114  std::ostream & outStream = std::cout) {
1115  Real tol = ROL_EPSILON<Real>();
1116  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1117  update(u,z);
1118  applyAdjointJacobian_1(*Jv,v,u,z,tol);
1119  Teuchos::RCP<Vector<Real> > iJJv = v.clone();
1120  update(u,z);
1121  applyInverseAdjointJacobian_1(*iJJv,*Jv,u,z,tol);
1122  Teuchos::RCP<Vector<Real> > diff = v.clone();
1123  diff->set(v);
1124  diff->axpy(-1.0,*iJJv);
1125  Real dnorm = diff->norm();
1126  Real vnorm = v.norm();
1127  if ( printToStream ) {
1128  std::stringstream hist;
1129  hist << std::scientific << std::setprecision(8);
1130  hist << "\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = "
1131  << dnorm << "\n";
1132  hist << " ||v|| = " << vnorm << "\n";
1133  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1134  outStream << hist.str();
1135  }
1136  return dnorm;
1137  }
1138 
1139 
1140 
1141  std::vector<std::vector<Real> > checkApplyJacobian_1(const Vector<Real> &u,
1142  const Vector<Real> &z,
1143  const Vector<Real> &v,
1144  const Vector<Real> &jv,
1145  const bool printToStream = true,
1146  std::ostream & outStream = std::cout,
1147  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1148  const int order = 1) {
1149  std::vector<Real> steps(numSteps);
1150  for(int i=0;i<numSteps;++i) {
1151  steps[i] = pow(10,-i);
1152  }
1153 
1154  return checkApplyJacobian_1(u,z,v,jv,steps,printToStream,outStream,order);
1155  }
1156 
1157 
1158 
1159 
1160  std::vector<std::vector<Real> > checkApplyJacobian_1(const Vector<Real> &u,
1161  const Vector<Real> &z,
1162  const Vector<Real> &v,
1163  const Vector<Real> &jv,
1164  const std::vector<Real> &steps,
1165  const bool printToStream = true,
1166  std::ostream & outStream = std::cout,
1167  const int order = 1) {
1168 
1169  TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1170  "Error: finite difference order must be 1,2,3, or 4" );
1171 
1172  Real one(1.0);
1173 
1176 
1177  Real tol = std::sqrt(ROL_EPSILON<Real>());
1178 
1179  int numSteps = steps.size();
1180  int numVals = 4;
1181  std::vector<Real> tmp(numVals);
1182  std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
1183 
1184  // Save the format state of the original outStream.
1185  Teuchos::oblackholestream oldFormatState;
1186  oldFormatState.copyfmt(outStream);
1187 
1188  // Compute constraint value at x.
1189  Teuchos::RCP<Vector<Real> > c = jv.clone();
1190  this->update(u,z);
1191  this->value(*c, u, z, tol);
1192 
1193  // Compute (Jacobian at x) times (vector v).
1194  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1195  this->applyJacobian_1(*Jv, v, u, z, tol);
1196  Real normJv = Jv->norm();
1197 
1198  // Temporary vectors.
1199  Teuchos::RCP<Vector<Real> > cdif = jv.clone();
1200  Teuchos::RCP<Vector<Real> > cnew = jv.clone();
1201  Teuchos::RCP<Vector<Real> > unew = u.clone();
1202 
1203  for (int i=0; i<numSteps; i++) {
1204 
1205  Real eta = steps[i];
1206 
1207  unew->set(u);
1208 
1209  cdif->set(*c);
1210  cdif->scale(weights[order-1][0]);
1211 
1212  for(int j=0; j<order; ++j) {
1213 
1214  unew->axpy(eta*shifts[order-1][j], v);
1215 
1216  if( weights[order-1][j+1] != 0 ) {
1217  this->update(*unew,z);
1218  this->value(*cnew,*unew,z,tol);
1219  cdif->axpy(weights[order-1][j+1],*cnew);
1220  }
1221 
1222  }
1223 
1224  cdif->scale(one/eta);
1225 
1226  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1227  jvCheck[i][0] = eta;
1228  jvCheck[i][1] = normJv;
1229  jvCheck[i][2] = cdif->norm();
1230  cdif->axpy(-one, *Jv);
1231  jvCheck[i][3] = cdif->norm();
1232 
1233  if (printToStream) {
1234  std::stringstream hist;
1235  if (i==0) {
1236  hist << std::right
1237  << std::setw(20) << "Step size"
1238  << std::setw(20) << "norm(Jac*vec)"
1239  << std::setw(20) << "norm(FD approx)"
1240  << std::setw(20) << "norm(abs error)"
1241  << "\n"
1242  << std::setw(20) << "---------"
1243  << std::setw(20) << "-------------"
1244  << std::setw(20) << "---------------"
1245  << std::setw(20) << "---------------"
1246  << "\n";
1247  }
1248  hist << std::scientific << std::setprecision(11) << std::right
1249  << std::setw(20) << jvCheck[i][0]
1250  << std::setw(20) << jvCheck[i][1]
1251  << std::setw(20) << jvCheck[i][2]
1252  << std::setw(20) << jvCheck[i][3]
1253  << "\n";
1254  outStream << hist.str();
1255  }
1256 
1257  }
1258 
1259  // Reset format state of outStream.
1260  outStream.copyfmt(oldFormatState);
1261 
1262  return jvCheck;
1263  } // checkApplyJacobian
1264 
1265 
1266  std::vector<std::vector<Real> > checkApplyJacobian_2(const Vector<Real> &u,
1267  const Vector<Real> &z,
1268  const Vector<Real> &v,
1269  const Vector<Real> &jv,
1270  const bool printToStream = true,
1271  std::ostream & outStream = std::cout,
1272  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1273  const int order = 1) {
1274  std::vector<Real> steps(numSteps);
1275  for(int i=0;i<numSteps;++i) {
1276  steps[i] = pow(10,-i);
1277  }
1278 
1279  return checkApplyJacobian_2(u,z,v,jv,steps,printToStream,outStream,order);
1280  }
1281 
1282 
1283 
1284 
1285  std::vector<std::vector<Real> > checkApplyJacobian_2(const Vector<Real> &u,
1286  const Vector<Real> &z,
1287  const Vector<Real> &v,
1288  const Vector<Real> &jv,
1289  const std::vector<Real> &steps,
1290  const bool printToStream = true,
1291  std::ostream & outStream = std::cout,
1292  const int order = 1) {
1293 
1294  TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1295  "Error: finite difference order must be 1,2,3, or 4" );
1296 
1297  Real one(1.0);
1298 
1301 
1302  Real tol = std::sqrt(ROL_EPSILON<Real>());
1303 
1304  int numSteps = steps.size();
1305  int numVals = 4;
1306  std::vector<Real> tmp(numVals);
1307  std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
1308 
1309  // Save the format state of the original outStream.
1310  Teuchos::oblackholestream oldFormatState;
1311  oldFormatState.copyfmt(outStream);
1312 
1313  // Compute constraint value at x.
1314  Teuchos::RCP<Vector<Real> > c = jv.clone();
1315  this->update(u,z);
1316  this->value(*c, u, z, tol);
1317 
1318  // Compute (Jacobian at x) times (vector v).
1319  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1320  this->applyJacobian_2(*Jv, v, u, z, tol);
1321  Real normJv = Jv->norm();
1322 
1323  // Temporary vectors.
1324  Teuchos::RCP<Vector<Real> > cdif = jv.clone();
1325  Teuchos::RCP<Vector<Real> > cnew = jv.clone();
1326  Teuchos::RCP<Vector<Real> > znew = z.clone();
1327 
1328  for (int i=0; i<numSteps; i++) {
1329 
1330  Real eta = steps[i];
1331 
1332  znew->set(z);
1333 
1334  cdif->set(*c);
1335  cdif->scale(weights[order-1][0]);
1336 
1337  for(int j=0; j<order; ++j) {
1338 
1339  znew->axpy(eta*shifts[order-1][j], v);
1340 
1341  if( weights[order-1][j+1] != 0 ) {
1342  this->update(u,*znew);
1343  this->value(*cnew,u,*znew,tol);
1344  cdif->axpy(weights[order-1][j+1],*cnew);
1345  }
1346 
1347  }
1348 
1349  cdif->scale(one/eta);
1350 
1351  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1352  jvCheck[i][0] = eta;
1353  jvCheck[i][1] = normJv;
1354  jvCheck[i][2] = cdif->norm();
1355  cdif->axpy(-one, *Jv);
1356  jvCheck[i][3] = cdif->norm();
1357 
1358  if (printToStream) {
1359  std::stringstream hist;
1360  if (i==0) {
1361  hist << std::right
1362  << std::setw(20) << "Step size"
1363  << std::setw(20) << "norm(Jac*vec)"
1364  << std::setw(20) << "norm(FD approx)"
1365  << std::setw(20) << "norm(abs error)"
1366  << "\n"
1367  << std::setw(20) << "---------"
1368  << std::setw(20) << "-------------"
1369  << std::setw(20) << "---------------"
1370  << std::setw(20) << "---------------"
1371  << "\n";
1372  }
1373  hist << std::scientific << std::setprecision(11) << std::right
1374  << std::setw(20) << jvCheck[i][0]
1375  << std::setw(20) << jvCheck[i][1]
1376  << std::setw(20) << jvCheck[i][2]
1377  << std::setw(20) << jvCheck[i][3]
1378  << "\n";
1379  outStream << hist.str();
1380  }
1381 
1382  }
1383 
1384  // Reset format state of outStream.
1385  outStream.copyfmt(oldFormatState);
1386 
1387  return jvCheck;
1388  } // checkApplyJacobian
1389 
1390 
1391 
1392  std::vector<std::vector<Real> > checkApplyAdjointHessian_11(const Vector<Real> &u,
1393  const Vector<Real> &z,
1394  const Vector<Real> &p,
1395  const Vector<Real> &v,
1396  const Vector<Real> &hv,
1397  const bool printToStream = true,
1398  std::ostream & outStream = std::cout,
1399  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1400  const int order = 1 ) {
1401  std::vector<Real> steps(numSteps);
1402  for(int i=0;i<numSteps;++i) {
1403  steps[i] = pow(10,-i);
1404  }
1405 
1406  return checkApplyAdjointHessian_11(u,z,p,v,hv,steps,printToStream,outStream,order);
1407 
1408  }
1409 
1410  std::vector<std::vector<Real> > checkApplyAdjointHessian_11(const Vector<Real> &u,
1411  const Vector<Real> &z,
1412  const Vector<Real> &p,
1413  const Vector<Real> &v,
1414  const Vector<Real> &hv,
1415  const std::vector<Real> &steps,
1416  const bool printToStream = true,
1417  std::ostream & outStream = std::cout,
1418  const int order = 1 ) {
1421 
1422  Real one(1.0);
1423 
1424  Real tol = std::sqrt(ROL_EPSILON<Real>());
1425 
1426  int numSteps = steps.size();
1427  int numVals = 4;
1428  std::vector<Real> tmp(numVals);
1429  std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1430 
1431  // Temporary vectors.
1432  Teuchos::RCP<Vector<Real> > AJdif = hv.clone();
1433  Teuchos::RCP<Vector<Real> > AJp = hv.clone();
1434  Teuchos::RCP<Vector<Real> > AHpv = hv.clone();
1435  Teuchos::RCP<Vector<Real> > AJnew = hv.clone();
1436  Teuchos::RCP<Vector<Real> > unew = u.clone();
1437 
1438  // Save the format state of the original outStream.
1439  Teuchos::oblackholestream oldFormatState;
1440  oldFormatState.copyfmt(outStream);
1441 
1442  // Apply adjoint Jacobian to p.
1443  this->update(u,z);
1444  this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1445 
1446  // Apply adjoint Hessian at (u,z), in direction v, to p.
1447  this->applyAdjointHessian_11(*AHpv, p, v, u, z, tol);
1448  Real normAHpv = AHpv->norm();
1449 
1450  for (int i=0; i<numSteps; i++) {
1451 
1452  Real eta = steps[i];
1453 
1454  // Apply adjoint Jacobian to p at (u+eta*v,z).
1455  unew->set(u);
1456 
1457  AJdif->set(*AJp);
1458  AJdif->scale(weights[order-1][0]);
1459 
1460  for(int j=0; j<order; ++j) {
1461 
1462  unew->axpy(eta*shifts[order-1][j],v);
1463 
1464  if( weights[order-1][j+1] != 0 ) {
1465  this->update(*unew,z);
1466  this->applyAdjointJacobian_1(*AJnew, p, *unew, z, tol);
1467  AJdif->axpy(weights[order-1][j+1],*AJnew);
1468  }
1469  }
1470 
1471  AJdif->scale(one/eta);
1472 
1473  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1474  ahpvCheck[i][0] = eta;
1475  ahpvCheck[i][1] = normAHpv;
1476  ahpvCheck[i][2] = AJdif->norm();
1477  AJdif->axpy(-one, *AHpv);
1478  ahpvCheck[i][3] = AJdif->norm();
1479 
1480  if (printToStream) {
1481  std::stringstream hist;
1482  if (i==0) {
1483  hist << std::right
1484  << std::setw(20) << "Step size"
1485  << std::setw(20) << "norm(adj(H)(u,v))"
1486  << std::setw(20) << "norm(FD approx)"
1487  << std::setw(20) << "norm(abs error)"
1488  << "\n"
1489  << std::setw(20) << "---------"
1490  << std::setw(20) << "-----------------"
1491  << std::setw(20) << "---------------"
1492  << std::setw(20) << "---------------"
1493  << "\n";
1494  }
1495  hist << std::scientific << std::setprecision(11) << std::right
1496  << std::setw(20) << ahpvCheck[i][0]
1497  << std::setw(20) << ahpvCheck[i][1]
1498  << std::setw(20) << ahpvCheck[i][2]
1499  << std::setw(20) << ahpvCheck[i][3]
1500  << "\n";
1501  outStream << hist.str();
1502  }
1503 
1504  }
1505 
1506  // Reset format state of outStream.
1507  outStream.copyfmt(oldFormatState);
1508 
1509  return ahpvCheck;
1510  } // checkApplyAdjointHessian_11
1511 
1512  std::vector<std::vector<Real> > checkApplyAdjointHessian_21(const Vector<Real> &u,
1513  const Vector<Real> &z,
1514  const Vector<Real> &p,
1515  const Vector<Real> &v,
1516  const Vector<Real> &hv,
1517  const bool printToStream = true,
1518  std::ostream & outStream = std::cout,
1519  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1520  const int order = 1 ) {
1521  std::vector<Real> steps(numSteps);
1522  for(int i=0;i<numSteps;++i) {
1523  steps[i] = pow(10,-i);
1524  }
1525 
1526  return checkApplyAdjointHessian_21(u,z,p,v,hv,steps,printToStream,outStream,order);
1527 
1528  }
1529 
1530  std::vector<std::vector<Real> > checkApplyAdjointHessian_21(const Vector<Real> &u,
1531  const Vector<Real> &z,
1532  const Vector<Real> &p,
1533  const Vector<Real> &v,
1534  const Vector<Real> &hv,
1535  const std::vector<Real> &steps,
1536  const bool printToStream = true,
1537  std::ostream & outStream = std::cout,
1538  const int order = 1 ) {
1541 
1542  Real one(1.0);
1543 
1544  Real tol = std::sqrt(ROL_EPSILON<Real>());
1545 
1546  int numSteps = steps.size();
1547  int numVals = 4;
1548  std::vector<Real> tmp(numVals);
1549  std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1550 
1551  // Temporary vectors.
1552  Teuchos::RCP<Vector<Real> > AJdif = hv.clone();
1553  Teuchos::RCP<Vector<Real> > AJp = hv.clone();
1554  Teuchos::RCP<Vector<Real> > AHpv = hv.clone();
1555  Teuchos::RCP<Vector<Real> > AJnew = hv.clone();
1556  Teuchos::RCP<Vector<Real> > znew = z.clone();
1557 
1558  // Save the format state of the original outStream.
1559  Teuchos::oblackholestream oldFormatState;
1560  oldFormatState.copyfmt(outStream);
1561 
1562  // Apply adjoint Jacobian to p.
1563  this->update(u,z);
1564  this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1565 
1566  // Apply adjoint Hessian at (u,z), in direction v, to p.
1567  this->applyAdjointHessian_21(*AHpv, p, v, u, z, tol);
1568  Real normAHpv = AHpv->norm();
1569 
1570  for (int i=0; i<numSteps; i++) {
1571 
1572  Real eta = steps[i];
1573 
1574  // Apply adjoint Jacobian to p at (u,z+eta*v).
1575  znew->set(z);
1576 
1577  AJdif->set(*AJp);
1578  AJdif->scale(weights[order-1][0]);
1579 
1580  for(int j=0; j<order; ++j) {
1581 
1582  znew->axpy(eta*shifts[order-1][j],v);
1583 
1584  if( weights[order-1][j+1] != 0 ) {
1585  this->update(u,*znew);
1586  this->applyAdjointJacobian_1(*AJnew, p, u, *znew, tol);
1587  AJdif->axpy(weights[order-1][j+1],*AJnew);
1588  }
1589  }
1590 
1591  AJdif->scale(one/eta);
1592 
1593  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1594  ahpvCheck[i][0] = eta;
1595  ahpvCheck[i][1] = normAHpv;
1596  ahpvCheck[i][2] = AJdif->norm();
1597  AJdif->axpy(-one, *AHpv);
1598  ahpvCheck[i][3] = AJdif->norm();
1599 
1600  if (printToStream) {
1601  std::stringstream hist;
1602  if (i==0) {
1603  hist << std::right
1604  << std::setw(20) << "Step size"
1605  << std::setw(20) << "norm(adj(H)(u,v))"
1606  << std::setw(20) << "norm(FD approx)"
1607  << std::setw(20) << "norm(abs error)"
1608  << "\n"
1609  << std::setw(20) << "---------"
1610  << std::setw(20) << "-----------------"
1611  << std::setw(20) << "---------------"
1612  << std::setw(20) << "---------------"
1613  << "\n";
1614  }
1615  hist << std::scientific << std::setprecision(11) << std::right
1616  << std::setw(20) << ahpvCheck[i][0]
1617  << std::setw(20) << ahpvCheck[i][1]
1618  << std::setw(20) << ahpvCheck[i][2]
1619  << std::setw(20) << ahpvCheck[i][3]
1620  << "\n";
1621  outStream << hist.str();
1622  }
1623 
1624  }
1625 
1626  // Reset format state of outStream.
1627  outStream.copyfmt(oldFormatState);
1628 
1629  return ahpvCheck;
1630  } // checkApplyAdjointHessian_21
1631 
1632  std::vector<std::vector<Real> > checkApplyAdjointHessian_12(const Vector<Real> &u,
1633  const Vector<Real> &z,
1634  const Vector<Real> &p,
1635  const Vector<Real> &v,
1636  const Vector<Real> &hv,
1637  const bool printToStream = true,
1638  std::ostream & outStream = std::cout,
1639  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1640  const int order = 1 ) {
1641  std::vector<Real> steps(numSteps);
1642  for(int i=0;i<numSteps;++i) {
1643  steps[i] = pow(10,-i);
1644  }
1645 
1646  return checkApplyAdjointHessian_12(u,z,p,v,hv,steps,printToStream,outStream,order);
1647 
1648  }
1649 
1650  std::vector<std::vector<Real> > checkApplyAdjointHessian_12(const Vector<Real> &u,
1651  const Vector<Real> &z,
1652  const Vector<Real> &p,
1653  const Vector<Real> &v,
1654  const Vector<Real> &hv,
1655  const std::vector<Real> &steps,
1656  const bool printToStream = true,
1657  std::ostream & outStream = std::cout,
1658  const int order = 1 ) {
1661 
1662  Real one(1.0);
1663 
1664  Real tol = std::sqrt(ROL_EPSILON<Real>());
1665 
1666  int numSteps = steps.size();
1667  int numVals = 4;
1668  std::vector<Real> tmp(numVals);
1669  std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1670 
1671  // Temporary vectors.
1672  Teuchos::RCP<Vector<Real> > AJdif = hv.clone();
1673  Teuchos::RCP<Vector<Real> > AJp = hv.clone();
1674  Teuchos::RCP<Vector<Real> > AHpv = hv.clone();
1675  Teuchos::RCP<Vector<Real> > AJnew = hv.clone();
1676  Teuchos::RCP<Vector<Real> > unew = u.clone();
1677 
1678  // Save the format state of the original outStream.
1679  Teuchos::oblackholestream oldFormatState;
1680  oldFormatState.copyfmt(outStream);
1681 
1682  // Apply adjoint Jacobian to p.
1683  this->update(u,z);
1684  this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1685 
1686  // Apply adjoint Hessian at (u,z), in direction v, to p.
1687  this->applyAdjointHessian_12(*AHpv, p, v, u, z, tol);
1688  Real normAHpv = AHpv->norm();
1689 
1690  for (int i=0; i<numSteps; i++) {
1691 
1692  Real eta = steps[i];
1693 
1694  // Apply adjoint Jacobian to p at (u+eta*v,z).
1695  unew->set(u);
1696 
1697  AJdif->set(*AJp);
1698  AJdif->scale(weights[order-1][0]);
1699 
1700  for(int j=0; j<order; ++j) {
1701 
1702  unew->axpy(eta*shifts[order-1][j],v);
1703 
1704  if( weights[order-1][j+1] != 0 ) {
1705  this->update(*unew,z);
1706  this->applyAdjointJacobian_2(*AJnew, p, *unew, z, tol);
1707  AJdif->axpy(weights[order-1][j+1],*AJnew);
1708  }
1709  }
1710 
1711  AJdif->scale(one/eta);
1712 
1713  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1714  ahpvCheck[i][0] = eta;
1715  ahpvCheck[i][1] = normAHpv;
1716  ahpvCheck[i][2] = AJdif->norm();
1717  AJdif->axpy(-one, *AHpv);
1718  ahpvCheck[i][3] = AJdif->norm();
1719 
1720  if (printToStream) {
1721  std::stringstream hist;
1722  if (i==0) {
1723  hist << std::right
1724  << std::setw(20) << "Step size"
1725  << std::setw(20) << "norm(adj(H)(u,v))"
1726  << std::setw(20) << "norm(FD approx)"
1727  << std::setw(20) << "norm(abs error)"
1728  << "\n"
1729  << std::setw(20) << "---------"
1730  << std::setw(20) << "-----------------"
1731  << std::setw(20) << "---------------"
1732  << std::setw(20) << "---------------"
1733  << "\n";
1734  }
1735  hist << std::scientific << std::setprecision(11) << std::right
1736  << std::setw(20) << ahpvCheck[i][0]
1737  << std::setw(20) << ahpvCheck[i][1]
1738  << std::setw(20) << ahpvCheck[i][2]
1739  << std::setw(20) << ahpvCheck[i][3]
1740  << "\n";
1741  outStream << hist.str();
1742  }
1743 
1744  }
1745 
1746  // Reset format state of outStream.
1747  outStream.copyfmt(oldFormatState);
1748 
1749  return ahpvCheck;
1750  } // checkApplyAdjointHessian_12
1751 
1752  std::vector<std::vector<Real> > checkApplyAdjointHessian_22(const Vector<Real> &u,
1753  const Vector<Real> &z,
1754  const Vector<Real> &p,
1755  const Vector<Real> &v,
1756  const Vector<Real> &hv,
1757  const bool printToStream = true,
1758  std::ostream & outStream = std::cout,
1759  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1760  const int order = 1 ) {
1761  std::vector<Real> steps(numSteps);
1762  for(int i=0;i<numSteps;++i) {
1763  steps[i] = pow(10,-i);
1764  }
1765 
1766  return checkApplyAdjointHessian_22(u,z,p,v,hv,steps,printToStream,outStream,order);
1767 
1768  }
1769 
1770  std::vector<std::vector<Real> > checkApplyAdjointHessian_22(const Vector<Real> &u,
1771  const Vector<Real> &z,
1772  const Vector<Real> &p,
1773  const Vector<Real> &v,
1774  const Vector<Real> &hv,
1775  const std::vector<Real> &steps,
1776  const bool printToStream = true,
1777  std::ostream & outStream = std::cout,
1778  const int order = 1 ) {
1781 
1782  Real one(1.0);
1783 
1784  Real tol = std::sqrt(ROL_EPSILON<Real>());
1785 
1786  int numSteps = steps.size();
1787  int numVals = 4;
1788  std::vector<Real> tmp(numVals);
1789  std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1790 
1791  // Temporary vectors.
1792  Teuchos::RCP<Vector<Real> > AJdif = hv.clone();
1793  Teuchos::RCP<Vector<Real> > AJp = hv.clone();
1794  Teuchos::RCP<Vector<Real> > AHpv = hv.clone();
1795  Teuchos::RCP<Vector<Real> > AJnew = hv.clone();
1796  Teuchos::RCP<Vector<Real> > znew = z.clone();
1797 
1798  // Save the format state of the original outStream.
1799  Teuchos::oblackholestream oldFormatState;
1800  oldFormatState.copyfmt(outStream);
1801 
1802  // Apply adjoint Jacobian to p.
1803  this->update(u,z);
1804  this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1805 
1806  // Apply adjoint Hessian at (u,z), in direction v, to p.
1807  this->applyAdjointHessian_22(*AHpv, p, v, u, z, tol);
1808  Real normAHpv = AHpv->norm();
1809 
1810  for (int i=0; i<numSteps; i++) {
1811 
1812  Real eta = steps[i];
1813 
1814  // Apply adjoint Jacobian to p at (u,z+eta*v).
1815  znew->set(z);
1816 
1817  AJdif->set(*AJp);
1818  AJdif->scale(weights[order-1][0]);
1819 
1820  for(int j=0; j<order; ++j) {
1821 
1822  znew->axpy(eta*shifts[order-1][j],v);
1823 
1824  if( weights[order-1][j+1] != 0 ) {
1825  this->update(u,*znew);
1826  this->applyAdjointJacobian_2(*AJnew, p, u, *znew, tol);
1827  AJdif->axpy(weights[order-1][j+1],*AJnew);
1828  }
1829  }
1830 
1831  AJdif->scale(one/eta);
1832 
1833  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1834  ahpvCheck[i][0] = eta;
1835  ahpvCheck[i][1] = normAHpv;
1836  ahpvCheck[i][2] = AJdif->norm();
1837  AJdif->axpy(-one, *AHpv);
1838  ahpvCheck[i][3] = AJdif->norm();
1839 
1840  if (printToStream) {
1841  std::stringstream hist;
1842  if (i==0) {
1843  hist << std::right
1844  << std::setw(20) << "Step size"
1845  << std::setw(20) << "norm(adj(H)(u,v))"
1846  << std::setw(20) << "norm(FD approx)"
1847  << std::setw(20) << "norm(abs error)"
1848  << "\n"
1849  << std::setw(20) << "---------"
1850  << std::setw(20) << "-----------------"
1851  << std::setw(20) << "---------------"
1852  << std::setw(20) << "---------------"
1853  << "\n";
1854  }
1855  hist << std::scientific << std::setprecision(11) << std::right
1856  << std::setw(20) << ahpvCheck[i][0]
1857  << std::setw(20) << ahpvCheck[i][1]
1858  << std::setw(20) << ahpvCheck[i][2]
1859  << std::setw(20) << ahpvCheck[i][3]
1860  << "\n";
1861  outStream << hist.str();
1862  }
1863 
1864  }
1865 
1866  // Reset format state of outStream.
1867  outStream.copyfmt(oldFormatState);
1868 
1869  return ahpvCheck;
1870  } // checkApplyAdjointHessian_22
1871 
1872 }; // class EqualityConstraint_SimOpt
1873 
1874 } // namespace ROL
1875 
1876 #endif
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)
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)
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...
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...
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...
virtual void scale(const Real alpha)=0
Compute where .
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 .
virtual void setSolveParameters(Teuchos::ParameterList &parlist)
Set solve parameters.
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)
virtual void plus(const Vector &x)=0
Compute , where .
const double weights[4][5]
Definition: ROL_Types.hpp:1105
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
Defines the linear algebra or vector space interface for simulation-based optimization.
virtual void value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol)=0
Evaluate the constraint operator at .
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)
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 ...
Teuchos::RCP< const Vector< Real > > get_2() const
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 ...
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 ...
Contains definitions of custom data types in ROL.
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)
void set_1(const Vector< Real > &vec)
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...
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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 ...
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Defines the equality constraint operator interface for simulation-based optimization.
virtual Real dot(const Vector &x) const =0
Compute where .
virtual Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:172
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
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...
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...
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...
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:213
Defines the equality constraint operator interface.
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)
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 .
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...
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:183
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 .
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 .
Teuchos::RCP< const Vector< Real > > get_1() const
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.
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)
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...
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...
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)
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
Definition: ROL_Types.hpp:73
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)
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
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.
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.
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...
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.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
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...
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)
virtual Real norm() const =0
Returns where .
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.
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:139
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)
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)
void set_2(const Vector< Real > &vec)
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)
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)
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.
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)