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>
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_rtol_;
103  const Real DEFAULT_stol_;
104  const Real DEFAULT_factor_;
105  const Real DEFAULT_decr_;
106  const int DEFAULT_maxit_;
107  const bool DEFAULT_print_;
108 
109  // User-set parameters for solve (backtracking Newton)
110  Real rtol_;
111  Real stol_;
112  Real factor_;
113  Real decr_;
114  int maxit_;
115  bool print_;
116 
117  // Flag to initialize vector storage in solve
119 
120 public:
122  : EqualityConstraint<Real>(),
123  unew_(Teuchos::null), jv_(Teuchos::null),
124  DEFAULT_rtol_(1.e-4*std::sqrt(ROL_EPSILON)),
125  DEFAULT_stol_(std::sqrt(ROL_EPSILON)),
126  DEFAULT_factor_(0.5),
127  DEFAULT_decr_(1.e-4),
128  DEFAULT_maxit_(500),
129  DEFAULT_print_(false),
130  rtol_(DEFAULT_rtol_), stol_(DEFAULT_stol_), factor_(DEFAULT_factor_),
131  decr_(DEFAULT_decr_), maxit_(DEFAULT_maxit_), print_(DEFAULT_print_),
132  firstSolve_(true) {}
133 
139  virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) {}
140 
141 
155  virtual void value(Vector<Real> &c,
156  const Vector<Real> &u,
157  const Vector<Real> &z,
158  Real &tol) = 0;
159 
171  virtual void solve(Vector<Real> &c,
172  Vector<Real> &u,
173  const Vector<Real> &z,
174  Real &tol) {
175  if ( firstSolve_ ) {
176  unew_ = u.clone();
177  jv_ = u.clone();
178  firstSolve_ = false;
179  }
180  update(u,z);
181  value(c,u,z,tol);
182  Real cnorm = c.norm(), alpha = 1.0, tmp = 0.0;
183  int cnt = 0;
184  if ( print_ ) {
185  std::cout << "\n Default EqualityConstraint_SimOpt::solve\n";
186  std::cout << " ";
187  std::cout << std::setw(6) << std::left << "iter";
188  std::cout << std::setw(15) << std::left << "rnorm";
189  std::cout << std::setw(15) << std::left << "alpha";
190  std::cout << "\n";
191  }
192  while ( cnorm > rtol_ && cnt < maxit_) {
193  // Compute Newton step
194  applyInverseJacobian_1(*jv_,c,u,z,tol);
195  unew_->set(u);
196  unew_->axpy(-alpha, *jv_);
197  update(*unew_,z);
198  value(c,*unew_,z,tol);
199  tmp = c.norm();
200  // Perform backtracking line search
201  while ( tmp > (1.0-decr_*alpha)*cnorm &&
202  alpha > stol_ ) {
203  alpha *= factor_;
204  unew_->set(u);
205  unew_->axpy(-alpha,*jv_);
206  update(*unew_,z);
207  value(c,*unew_,z,tol);
208  tmp = c.norm();
209  }
210  if ( print_ ) {
211  std::cout << " ";
212  std::cout << std::setw(6) << std::left << cnt;
213  std::cout << std::scientific << std::setprecision(6);
214  std::cout << std::setw(15) << std::left << tmp;
215  std::cout << std::scientific << std::setprecision(6);
216  std::cout << std::setw(15) << std::left << alpha;
217  std::cout << "\n";
218  }
219  // Update iterate
220  cnorm = tmp;
221  u.set(*unew_);
222  alpha = 1.0;
223  cnt++;
224  }
225  }
226 
245  virtual void setSolveParameters(Teuchos::ParameterList &parlist) {
246  Teuchos::ParameterList & list = parlist.sublist("SimOpt").sublist("Solve");
247  rtol_ = list.get("Residual Tolerance", DEFAULT_rtol_);
248  maxit_ = list.get("Iteration Limit", DEFAULT_maxit_);
249  decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_);
250  stol_ = list.get("Step Tolerance", DEFAULT_stol_);
251  factor_ = list.get("Backtracking Factor", DEFAULT_factor_);
252  print_ = list.get("Output Iteration History", DEFAULT_print_);
253  }
254 
270  virtual void applyJacobian_1(Vector<Real> &jv,
271  const Vector<Real> &v,
272  const Vector<Real> &u,
273  const Vector<Real> &z,
274  Real &tol) {
275  Real ctol = std::sqrt(ROL_EPSILON);
276  // Compute step length
277  Real h = tol;
278  if (v.norm() > std::sqrt(ROL_EPSILON)) {
279  h = std::max(1.0,u.norm()/v.norm())*tol;
280  }
281  // Update state vector to u + hv
282  Teuchos::RCP<Vector<Real> > unew = u.clone();
283  unew->set(u);
284  unew->axpy(h,v);
285  // Compute new constraint value
286  update(*unew,z);
287  value(jv,*unew,z,ctol);
288  // Compute current constraint value
289  Teuchos::RCP<Vector<Real> > cold = jv.clone();
290  update(u,z);
291  value(*cold,u,z,ctol);
292  // Compute Newton quotient
293  jv.axpy(-1.0,*cold);
294  jv.scale(1.0/h);
295  }
296 
297 
313  virtual void applyJacobian_2(Vector<Real> &jv,
314  const Vector<Real> &v,
315  const Vector<Real> &u,
316  const Vector<Real> &z,
317  Real &tol) {
318  Real ctol = std::sqrt(ROL_EPSILON);
319  // Compute step length
320  Real h = tol;
321  if (v.norm() > std::sqrt(ROL_EPSILON)) {
322  h = std::max(1.0,u.norm()/v.norm())*tol;
323  }
324  // Update state vector to u + hv
325  Teuchos::RCP<Vector<Real> > znew = z.clone();
326  znew->set(z);
327  znew->axpy(h,v);
328  // Compute new constraint value
329  update(u,*znew);
330  value(jv,u,*znew,ctol);
331  // Compute current constraint value
332  Teuchos::RCP<Vector<Real> > cold = jv.clone();
333  update(u,z);
334  value(*cold,u,z,ctol);
335  // Compute Newton quotient
336  jv.axpy(-1.0,*cold);
337  jv.scale(1.0/h);
338  }
339 
356  const Vector<Real> &v,
357  const Vector<Real> &u,
358  const Vector<Real> &z,
359  Real &tol) {
360  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
361  "The method applyInverseJacobian_1 is used but not implemented!\n");
362  }
363 
380  const Vector<Real> &v,
381  const Vector<Real> &u,
382  const Vector<Real> &z,
383  Real &tol) {
384  applyAdjointJacobian_1(ajv, v, u, z, v.dual(), tol);
385  }
386 
387 
406  const Vector<Real> &v,
407  const Vector<Real> &u,
408  const Vector<Real> &z,
409  const Vector<Real> &dualv,
410  Real &tol) {
411  Real ctol = std::sqrt(ROL_EPSILON);
412  Real h = tol;
413  if (v.norm() > std::sqrt(ROL_EPSILON)) {
414  h = std::max(1.0,u.norm()/v.norm())*tol;
415  }
416  Teuchos::RCP<Vector<Real> > cold = dualv.clone();
417  Teuchos::RCP<Vector<Real> > cnew = dualv.clone();
418  update(u,z);
419  value(*cold,u,z,ctol);
420  Teuchos::RCP<Vector<Real> > unew = u.clone();
421  ajv.zero();
422  for (int i = 0; i < u.dimension(); i++) {
423  unew->set(u);
424  unew->axpy(h,*(u.basis(i)));
425  update(*unew,z);
426  value(*cnew,*unew,z,ctol);
427  cnew->axpy(-1.0,*cold);
428  cnew->scale(1.0/h);
429  ajv.axpy(cnew->dot(v),*((u.dual()).basis(i)));
430  }
431  update(u,z);
432  }
433 
434 
451  const Vector<Real> &v,
452  const Vector<Real> &u,
453  const Vector<Real> &z,
454  Real &tol) {
455  applyAdjointJacobian_2(ajv, v, u, z, v.dual(), tol);
456  }
457 
458 
477  const Vector<Real> &v,
478  const Vector<Real> &u,
479  const Vector<Real> &z,
480  const Vector<Real> &dualv,
481  Real &tol) {
482  Real ctol = std::sqrt(ROL_EPSILON);
483  Real h = tol;
484  if (v.norm() > std::sqrt(ROL_EPSILON)) {
485  h = std::max(1.0,u.norm()/v.norm())*tol;
486  }
487  Teuchos::RCP<Vector<Real> > cold = dualv.clone();
488  Teuchos::RCP<Vector<Real> > cnew = dualv.clone();
489  update(u,z);
490  value(*cold,u,z,ctol);
491  Teuchos::RCP<Vector<Real> > znew = z.clone();
492  ajv.zero();
493  for (int i = 0; i < z.dimension(); i++) {
494  znew->set(z);
495  znew->axpy(h,*(z.basis(i)));
496  update(u,*znew);
497  value(*cnew,u,*znew,ctol);
498  cnew->axpy(-1.0,*cold);
499  cnew->scale(1.0/h);
500  ajv.axpy(cnew->dot(v),*((z.dual()).basis(i)));
501  }
502  update(u,z);
503  }
504 
521  const Vector<Real> &v,
522  const Vector<Real> &u,
523  const Vector<Real> &z,
524  Real &tol) {
525  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
526  "The method applyInverseAdjointJacobian_1 is used but not implemented!\n");
527  };
528 
546  const Vector<Real> &w,
547  const Vector<Real> &v,
548  const Vector<Real> &u,
549  const Vector<Real> &z,
550  Real &tol) {
551  Real jtol = std::sqrt(ROL_EPSILON);
552  // Compute step size
553  Real h = tol;
554  if (v.norm() > std::sqrt(ROL_EPSILON)) {
555  h = std::max(1.0,u.norm()/v.norm())*tol;
556  }
557  // Evaluate Jacobian at new state
558  Teuchos::RCP<Vector<Real> > unew = u.clone();
559  unew->set(u);
560  unew->axpy(h,v);
561  update(*unew,z);
562  applyAdjointJacobian_1(ahwv,w,*unew,z,jtol);
563  // Evaluate Jacobian at old state
564  Teuchos::RCP<Vector<Real> > jv = ahwv.clone();
565  update(u,z);
566  applyAdjointJacobian_1(*jv,w,u,z,jtol);
567  // Compute Newton quotient
568  ahwv.axpy(-1.0,*jv);
569  ahwv.scale(1.0/h);
570  }
571 
572 
590  const Vector<Real> &w,
591  const Vector<Real> &v,
592  const Vector<Real> &u,
593  const Vector<Real> &z,
594  Real &tol) {
595  Real jtol = std::sqrt(ROL_EPSILON);
596  // Compute step size
597  Real h = tol;
598  if (v.norm() > std::sqrt(ROL_EPSILON)) {
599  h = std::max(1.0,u.norm()/v.norm())*tol;
600  }
601  // Evaluate Jacobian at new state
602  Teuchos::RCP<Vector<Real> > unew = u.clone();
603  unew->set(u);
604  unew->axpy(h,v);
605  update(*unew,z);
606  applyAdjointJacobian_2(ahwv,w,*unew,z,jtol);
607  // Evaluate Jacobian at old state
608  Teuchos::RCP<Vector<Real> > jv = ahwv.clone();
609  update(u,z);
610  applyAdjointJacobian_2(*jv,w,u,z,jtol);
611  // Compute Newton quotient
612  ahwv.axpy(-1.0,*jv);
613  ahwv.scale(1.0/h);
614  }
615 
616 
634  const Vector<Real> &w,
635  const Vector<Real> &v,
636  const Vector<Real> &u,
637  const Vector<Real> &z,
638  Real &tol) {
639  Real jtol = std::sqrt(ROL_EPSILON);
640  // Compute step size
641  Real h = tol;
642  if (v.norm() > std::sqrt(ROL_EPSILON)) {
643  h = std::max(1.0,u.norm()/v.norm())*tol;
644  }
645  // Evaluate Jacobian at new control
646  Teuchos::RCP<Vector<Real> > znew = z.clone();
647  znew->set(z);
648  znew->axpy(h,v);
649  update(u,*znew);
650  applyAdjointJacobian_1(ahwv,w,u,*znew,jtol);
651  // Evaluate Jacobian at old control
652  Teuchos::RCP<Vector<Real> > jv = ahwv.clone();
653  update(u,z);
654  applyAdjointJacobian_1(*jv,w,u,z,jtol);
655  // Compute Newton quotient
656  ahwv.axpy(-1.0,*jv);
657  ahwv.scale(1.0/h);
658  }
659 
677  const Vector<Real> &w,
678  const Vector<Real> &v,
679  const Vector<Real> &u,
680  const Vector<Real> &z,
681  Real &tol) {
682  Real jtol = std::sqrt(ROL_EPSILON);
683  // Compute step size
684  Real h = tol;
685  if (v.norm() > std::sqrt(ROL_EPSILON)) {
686  h = std::max(1.0,u.norm()/v.norm())*tol;
687  }
688  // Evaluate Jacobian at new control
689  Teuchos::RCP<Vector<Real> > znew = z.clone();
690  znew->set(z);
691  znew->axpy(h,v);
692  update(u,*znew);
693  applyAdjointJacobian_2(ahwv,w,u,*znew,jtol);
694  // Evaluate Jacobian at old control
695  Teuchos::RCP<Vector<Real> > jv = ahwv.clone();
696  update(u,z);
697  applyAdjointJacobian_2(*jv,w,u,z,jtol);
698  // Compute Newton quotient
699  ahwv.axpy(-1.0,*jv);
700  ahwv.scale(1.0/h);
701 }
702 
741  virtual std::vector<Real> solveAugmentedSystem(Vector<Real> &v1,
742  Vector<Real> &v2,
743  const Vector<Real> &b1,
744  const Vector<Real> &b2,
745  const Vector<Real> &x,
746  Real &tol) {
747  return EqualityConstraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,x,tol);
748  }
749 
750 
770  const Vector<Real> &v,
771  const Vector<Real> &x,
772  const Vector<Real> &g,
773  Real &tol) {
774  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(x);
775  Teuchos::RCP<ROL::Vector<Real> > ijv = (xs.get_1())->clone();
776 
777  try {
778  applyInverseJacobian_1(*ijv, v, *(xs.get_1()), *(xs.get_2()), tol);
779  }
780  catch (const std::logic_error &e) {
782  return;
783  }
784 
785  const Vector_SimOpt<Real> &gs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(g);
786  Teuchos::RCP<ROL::Vector<Real> > ijv_dual = (gs.get_1())->clone();
787  ijv_dual->set(ijv->dual());
788 
789  try {
790  applyInverseAdjointJacobian_1(pv, *ijv_dual, *(xs.get_1()), *(xs.get_2()), tol);
791  }
792  catch (const std::logic_error &e) {
794  return;
795  }
796 
797  }
798 
799 //
800 // EqualityConstraint_SimOpt(void) : EqualityConstraint<Real>() {}
801 
807  virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
808  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
809  Teuchos::dyn_cast<const Vector<Real> >(x));
810  update(*(xs.get_1()),*(xs.get_2()),flag,iter);
811  }
812 
815  virtual bool isFeasible( const Vector<Real> &v ) { return true; }
816 
817  virtual void value(Vector<Real> &c,
818  const Vector<Real> &x,
819  Real &tol) {
820  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
821  Teuchos::dyn_cast<const Vector<Real> >(x));
822  value(c,*(xs.get_1()),*(xs.get_2()),tol);
823  }
824 
825 
826  virtual void applyJacobian(Vector<Real> &jv,
827  const Vector<Real> &v,
828  const Vector<Real> &x,
829  Real &tol) {
830  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
831  Teuchos::dyn_cast<const Vector<Real> >(x));
832  const Vector_SimOpt<Real> &vs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
833  Teuchos::dyn_cast<const Vector<Real> >(v));
834  applyJacobian_1(jv,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
835  Teuchos::RCP<Vector<Real> > jv2 = jv.clone();
836  applyJacobian_2(*jv2,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
837  jv.plus(*jv2);
838  }
839 
840 
842  const Vector<Real> &v,
843  const Vector<Real> &x,
844  Real &tol) {
845  Vector_SimOpt<Real> &ajvs = Teuchos::dyn_cast<Vector_SimOpt<Real> >(
846  Teuchos::dyn_cast<Vector<Real> >(ajv));
847  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
848  Teuchos::dyn_cast<const Vector<Real> >(x));
849  Teuchos::RCP<Vector<Real> > ajv1 = (ajvs.get_1())->clone();
850  applyAdjointJacobian_1(*ajv1,v,*(xs.get_1()),*(xs.get_2()),tol);
851  ajvs.set_1(*ajv1);
852  Teuchos::RCP<Vector<Real> > ajv2 = (ajvs.get_2())->clone();
853  applyAdjointJacobian_2(*ajv2,v,*(xs.get_1()),*(xs.get_2()),tol);
854  ajvs.set_2(*ajv2);
855  }
856 
857 
858  virtual void applyAdjointHessian(Vector<Real> &ahwv,
859  const Vector<Real> &w,
860  const Vector<Real> &v,
861  const Vector<Real> &x,
862  Real &tol) {
863  Vector_SimOpt<Real> &ahwvs = Teuchos::dyn_cast<Vector_SimOpt<Real> >(
864  Teuchos::dyn_cast<Vector<Real> >(ahwv));
865  const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
866  Teuchos::dyn_cast<const Vector<Real> >(x));
867  const Vector_SimOpt<Real> &vs = Teuchos::dyn_cast<const Vector_SimOpt<Real> >(
868  Teuchos::dyn_cast<const Vector<Real> >(v));
869  // Block-row 1
870  Teuchos::RCP<Vector<Real> > C11 = (ahwvs.get_1())->clone();
871  Teuchos::RCP<Vector<Real> > C21 = (ahwvs.get_1())->clone();
872  applyAdjointHessian_11(*C11,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
873  applyAdjointHessian_21(*C21,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
874  C11->plus(*C21);
875  ahwvs.set_1(*C11);
876  // Block-row 2
877  Teuchos::RCP<Vector<Real> > C12 = (ahwvs.get_2())->clone();
878  Teuchos::RCP<Vector<Real> > C22 = (ahwvs.get_2())->clone();
879  applyAdjointHessian_12(*C12,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
880  applyAdjointHessian_22(*C22,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
881  C22->plus(*C12);
882  ahwvs.set_2(*C22);
883  }
884 
885 
886 
887  virtual Real checkSolve(const ROL::Vector<Real> &u,
888  const ROL::Vector<Real> &z,
889  const ROL::Vector<Real> &c,
890  const bool printToStream = true,
891  std::ostream & outStream = std::cout) {
892  // Solve equality constraint for u.
893  Real tol = ROL_EPSILON;
894  Teuchos::RCP<ROL::Vector<Real> > r = c.clone();
895  Teuchos::RCP<ROL::Vector<Real> > s = u.clone();
896  solve(*r,*s,z,tol);
897  // Evaluate equality constraint residual at (u,z).
898  Teuchos::RCP<ROL::Vector<Real> > cs = c.clone();
899  update(*s,z);
900  value(*cs,*s,z,tol);
901  // Output norm of residual.
902  Real rnorm = r->norm();
903  Real cnorm = cs->norm();
904  if ( printToStream ) {
905  std::stringstream hist;
906  hist << std::scientific << std::setprecision(8);
907  hist << "\nTest SimOpt solve at feasible (u,z):\n";
908  hist << " Solver Residual = " << rnorm << "\n";
909  hist << " ||c(u,z)|| = " << cnorm << "\n";
910  outStream << hist.str();
911  }
912  return cnorm;
913  }
914 
915 
929  const Vector<Real> &v,
930  const Vector<Real> &u,
931  const Vector<Real> &z,
932  const bool printToStream = true,
933  std::ostream & outStream = std::cout) {
934  return checkAdjointConsistencyJacobian_1(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
935  }
936 
937 
954  const Vector<Real> &v,
955  const Vector<Real> &u,
956  const Vector<Real> &z,
957  const Vector<Real> &dualw,
958  const Vector<Real> &dualv,
959  const bool printToStream = true,
960  std::ostream & outStream = std::cout) {
961  Real tol = ROL_EPSILON;
962  Teuchos::RCP<Vector<Real> > Jv = dualw.clone();
963  update(u,z);
964  applyJacobian_1(*Jv,v,u,z,tol);
965  Real wJv = w.dot(Jv->dual());
966  Teuchos::RCP<Vector<Real> > Jw = dualv.clone();
967  update(u,z);
968  applyAdjointJacobian_1(*Jw,w,u,z,tol);
969  Real vJw = v.dot(Jw->dual());
970  Real diff = std::abs(wJv-vJw);
971  if ( printToStream ) {
972  std::stringstream hist;
973  hist << std::scientific << std::setprecision(8);
974  hist << "\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
975  << diff << "\n";
976  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
977  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW) << "\n";
978  outStream << hist.str();
979  }
980  return diff;
981  }
982 
983 
997  const Vector<Real> &v,
998  const Vector<Real> &u,
999  const Vector<Real> &z,
1000  const bool printToStream = true,
1001  std::ostream & outStream = std::cout) {
1002  return checkAdjointConsistencyJacobian_2(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
1003  }
1004 
1021  const Vector<Real> &v,
1022  const Vector<Real> &u,
1023  const Vector<Real> &z,
1024  const Vector<Real> &dualw,
1025  const Vector<Real> &dualv,
1026  const bool printToStream = true,
1027  std::ostream & outStream = std::cout) {
1028  Real tol = ROL_EPSILON;
1029  Teuchos::RCP<Vector<Real> > Jv = dualw.clone();
1030  update(u,z);
1031  applyJacobian_2(*Jv,v,u,z,tol);
1032  Real wJv = w.dot(Jv->dual());
1033  Teuchos::RCP<Vector<Real> > Jw = dualv.clone();
1034  update(u,z);
1035  applyAdjointJacobian_2(*Jw,w,u,z,tol);
1036  Real vJw = v.dot(Jw->dual());
1037  Real diff = std::abs(wJv-vJw);
1038  if ( printToStream ) {
1039  std::stringstream hist;
1040  hist << std::scientific << std::setprecision(8);
1041  hist << "\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1042  << diff << "\n";
1043  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1044  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW) << "\n";
1045  outStream << hist.str();
1046  }
1047  return diff;
1048  }
1049 
1050  virtual Real checkInverseJacobian_1(const Vector<Real> &jv,
1051  const Vector<Real> &v,
1052  const Vector<Real> &u,
1053  const Vector<Real> &z,
1054  const bool printToStream = true,
1055  std::ostream & outStream = std::cout) {
1056  Real tol = ROL_EPSILON;
1057  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1058  update(u,z);
1059  applyJacobian_1(*Jv,v,u,z,tol);
1060  Teuchos::RCP<Vector<Real> > iJJv = u.clone();
1061  update(u,z);
1062  applyInverseJacobian_1(*iJJv,*Jv,u,z,tol);
1063  Teuchos::RCP<Vector<Real> > diff = v.clone();
1064  diff->set(v);
1065  diff->axpy(-1.0,*iJJv);
1066  Real dnorm = diff->norm();
1067  Real vnorm = v.norm();
1068  if ( printToStream ) {
1069  std::stringstream hist;
1070  hist << std::scientific << std::setprecision(8);
1071  hist << "\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = "
1072  << dnorm << "\n";
1073  hist << " ||v|| = " << vnorm << "\n";
1074  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW) << "\n";
1075  outStream << hist.str();
1076  }
1077  return dnorm;
1078  }
1079 
1081  const Vector<Real> &v,
1082  const Vector<Real> &u,
1083  const Vector<Real> &z,
1084  const bool printToStream = true,
1085  std::ostream & outStream = std::cout) {
1086  Real tol = ROL_EPSILON;
1087  Teuchos::RCP<Vector<Real> > Jv = jv.clone();
1088  update(u,z);
1089  applyAdjointJacobian_1(*Jv,v,u,z,tol);
1090  Teuchos::RCP<Vector<Real> > iJJv = v.clone();
1091  update(u,z);
1092  applyInverseAdjointJacobian_1(*iJJv,*Jv,u,z,tol);
1093  Teuchos::RCP<Vector<Real> > diff = v.clone();
1094  diff->set(v);
1095  diff->axpy(-1.0,*iJJv);
1096  Real dnorm = diff->norm();
1097  Real vnorm = v.norm();
1098  if ( printToStream ) {
1099  std::stringstream hist;
1100  hist << std::scientific << std::setprecision(8);
1101  hist << "\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = "
1102  << dnorm << "\n";
1103  hist << " ||v|| = " << vnorm << "\n";
1104  hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW) << "\n";
1105  outStream << hist.str();
1106  }
1107  return dnorm;
1108  }
1109 
1110 }; // class EqualityConstraint_SimOpt
1111 
1112 } // namespace ROL
1113 
1114 #endif
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 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
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 .
Teuchos::RCP< const Vector< Real > > get_1() const
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:183
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.
virtual void plus(const Vector &x)=0
Compute , where .
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 .
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 ...
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
Contains definitions of custom data types in ROL.
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
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 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 is ...
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...
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 adjoint of the partial constraint Hessian at , , to vector in direction ...
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
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 .
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.
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 applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
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 is ...
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 Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:172
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.
Teuchos::RCP< const Vector< Real > > get_2() const
void set_2(const Vector< Real > &vec)
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.
static const double ROL_UNDERFLOW
Platform-dependent minimum double.
Definition: ROL_Types.hpp:133
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118