ROL
ROL_Constraint_SimOpt.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Rapid Optimization Library (ROL) Package
4//
5// Copyright 2014 NTESS and the ROL contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef ROL_CONSTRAINT_SIMOPT_H
11#define ROL_CONSTRAINT_SIMOPT_H
12
13#include "ROL_Constraint.hpp"
14#include "ROL_Vector_SimOpt.hpp"
15#include "ROL_Types.hpp"
16#include <iostream>
17
18namespace ROL {
19
20template <class Real>
21class Constraint_SimOpt;
22
23}
24
25#include "ROL_NonlinearLeastSquaresObjective.hpp"
26#include "ROL_SimConstraint.hpp"
27#include "ROL_Objective_FSsolver.hpp"
28#include "ROL_TypeU_TrustRegionAlgorithm.hpp"
29#include "ROL_TypeE_CompositeStepAlgorithm.hpp"
30
71namespace ROL {
72
73template <class Real>
74class Constraint_SimOpt : public Constraint<Real> {
75private:
76 // Additional vector storage for solve
77 Ptr<Vector<Real>> unew_;
78 Ptr<Vector<Real>> jv_;
79
80 // Default parameters for solve (backtracking Newton)
81 const Real DEFAULT_atol_;
82 const Real DEFAULT_rtol_;
83 const Real DEFAULT_stol_;
84 const Real DEFAULT_factor_;
85 const Real DEFAULT_decr_;
86 const int DEFAULT_maxit_;
87 const bool DEFAULT_print_;
88 const bool DEFAULT_zero_;
89 const int DEFAULT_solverType_;
90
91 // User-set parameters for solve (backtracking Newton)
92
93protected:
94 Real atol_;
95 Real rtol_;
96 Real stol_;
97 Real factor_;
98 Real decr_;
99 int maxit_;
100 bool print_;
101 bool zero_;
102 int solverType_;
103
104 // Flag to initialize vector storage in solve
105 bool firstSolve_;
106
107public:
123
129 virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) {
130 update_1(u,flag,iter);
131 update_2(z,flag,iter);
132 }
133 virtual void update( const Vector<Real> &u, const Vector<Real> &z, UpdateType type, int iter = -1 ) {
134 update_1(u,type,iter);
135 update_2(z,type,iter);
136 }
137
143 virtual void update_1( const Vector<Real> &u, bool flag = true, int iter = -1 ) {}
144 virtual void update_1( const Vector<Real> &u, UpdateType type, int iter = -1 ) {}
145
151 virtual void update_2( const Vector<Real> &z, bool flag = true, int iter = -1 ) {}
152 virtual void update_2( const Vector<Real> &z, UpdateType type, int iter = -1 ) {}
153
160 virtual void solve_update( const Vector<Real> &u, const Vector<Real> &z, UpdateType type, int iter = -1) {}
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 ( zero_ ) u.zero();
196 Ptr<std::ostream> stream = makeStreamPtr(std::cout, print_);
197 solve_update(u,z,UpdateType::Initial,0);
198 value(c,u,z,tol);
199 Real cnorm = c.norm();
200 const Real ctol = std::min(atol_, rtol_*cnorm);
201 if (solverType_==0 || solverType_==3 || solverType_==4) {
202 if ( firstSolve_ ) {
203 unew_ = u.clone();
204 jv_ = u.clone();
205 firstSolve_ = false;
206 }
207 const Real one(1);
208 Real alpha(1), tmp(0);
209 int cnt = 0;
210 *stream << std::endl;
211 *stream << " Default Constraint_SimOpt::solve" << std::endl;
212 *stream << " ";
213 *stream << std::setw(6) << std::left << "iter";
214 *stream << std::setw(15) << std::left << "rnorm";
215 *stream << std::setw(15) << std::left << "alpha";
216 *stream << std::endl;
217 for (cnt = 0; cnt < maxit_; ++cnt) {
218 // Compute Newton step
219 applyInverseJacobian_1(*jv_,c,u,z,tol);
220 unew_->set(u);
221 unew_->axpy(-alpha, *jv_);
222 solve_update(*unew_,z,UpdateType::Trial);
223 value(c,*unew_,z,tol);
224 tmp = c.norm();
225 // Perform backtracking line search
226 while ( tmp > (one-decr_*alpha)*cnorm &&
227 alpha > stol_ ) {
228 alpha *= factor_;
229 unew_->set(u);
230 unew_->axpy(-alpha,*jv_);
231 solve_update(*unew_,z,UpdateType::Trial);
232 value(c,*unew_,z,tol);
233 tmp = c.norm();
234 }
235 *stream << " ";
236 *stream << std::setw(6) << std::left << cnt;
237 *stream << std::scientific << std::setprecision(6);
238 *stream << std::setw(15) << std::left << tmp;
239 *stream << std::scientific << std::setprecision(6);
240 *stream << std::setw(15) << std::left << alpha;
241 *stream << std::endl;
242 // Update iterate
243 cnorm = tmp;
244 u.set(*unew_);
245 solve_update(u,z,UpdateType::Accept,cnt);
246 if (cnorm <= ctol) break; // = covers the case of identically zero residual
247 alpha = one;
248 }
249 }
250 if (solverType_==1 || (solverType_==3 && cnorm > ctol)) {
251 Ptr<Constraint<Real>> con = makePtr<SimConstraint<Real>>(makePtrFromRef(*this),makePtrFromRef(z),true);
252 Ptr<Objective<Real>> obj = makePtr<NonlinearLeastSquaresObjective<Real>>(con,u,c,true);
253 ParameterList parlist;
254 parlist.sublist("General").set("Output Level",1);
255 parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100);
256 parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG");
257 parlist.sublist("Status Test").set("Gradient Tolerance",ctol);
258 parlist.sublist("Status Test").set("Step Tolerance",stol_);
259 parlist.sublist("Status Test").set("Iteration Limit",maxit_);
260 Ptr<TypeU::Algorithm<Real>> algo = makePtr<TypeU::TrustRegionAlgorithm<Real>>(parlist);
261 algo->run(u,*obj,*stream);
262 value(c,u,z,tol);
263 }
264 if (solverType_==2 || (solverType_==4 && cnorm > ctol)) {
265 Ptr<Constraint<Real>> con = makePtr<SimConstraint<Real>>(makePtrFromRef(*this),makePtrFromRef(z),true);
266 Ptr<Objective<Real>> obj = makePtr<Objective_FSsolver<Real>>();
267 Ptr<Vector<Real>> l = c.dual().clone();
268 ParameterList parlist;
269 parlist.sublist("General").set("Output Level",1);
270 parlist.sublist("Status Test").set("Constraint Tolerance",ctol);
271 parlist.sublist("Status Test").set("Step Tolerance",stol_);
272 parlist.sublist("Status Test").set("Iteration Limit",maxit_);
273 Ptr<TypeE::Algorithm<Real>> algo = makePtr<TypeE::CompositeStepAlgorithm<Real>>(parlist);
274 algo->run(u,*obj,*con,*l,*stream);
275 value(c,u,z,tol);
276 }
277 if (solverType_ > 4 || solverType_ < 0) {
278 ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
279 ">>> ERROR (ROL:Constraint_SimOpt:solve): Invalid solver type!");
280 }
281 }
282
303 virtual void setSolveParameters(ParameterList &parlist) {
304 ParameterList & list = parlist.sublist("SimOpt").sublist("Solve");
305 atol_ = list.get("Absolute Residual Tolerance", DEFAULT_atol_);
306 rtol_ = list.get("Relative Residual Tolerance", DEFAULT_rtol_);
307 maxit_ = list.get("Iteration Limit", DEFAULT_maxit_);
308 decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_);
309 stol_ = list.get("Step Tolerance", DEFAULT_stol_);
310 factor_ = list.get("Backtracking Factor", DEFAULT_factor_);
311 print_ = list.get("Output Iteration History", DEFAULT_print_);
312 zero_ = list.get("Zero Initial Guess", DEFAULT_zero_);
313 solverType_ = list.get("Solver Type", DEFAULT_solverType_);
314 }
315
331 virtual void applyJacobian_1(Vector<Real> &jv,
332 const Vector<Real> &v,
333 const Vector<Real> &u,
334 const Vector<Real> &z,
335 Real &tol) {
336 Real ctol = std::sqrt(ROL_EPSILON<Real>());
337 // Compute step length
338 Real h = tol;
339 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
340 h = std::max(1.0,u.norm()/v.norm())*tol;
341 }
342 // Update state vector to u + hv
343 Ptr<Vector<Real>> unew = u.clone();
344 unew->set(u);
345 unew->axpy(h,v);
346 // Compute new constraint value
347 update(*unew,z,UpdateType::Temp);
348 value(jv,*unew,z,ctol);
349 // Compute current constraint value
350 Ptr<Vector<Real>> cold = jv.clone();
351 update(u,z,UpdateType::Temp);
352 value(*cold,u,z,ctol);
353 // Compute Newton quotient
354 jv.axpy(-1.0,*cold);
355 jv.scale(1.0/h);
356 }
357
358
374 virtual void applyJacobian_2(Vector<Real> &jv,
375 const Vector<Real> &v,
376 const Vector<Real> &u,
377 const Vector<Real> &z,
378 Real &tol) {
379 Real ctol = std::sqrt(ROL_EPSILON<Real>());
380 // Compute step length
381 Real h = tol;
382 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
383 h = std::max(1.0,u.norm()/v.norm())*tol;
384 }
385 // Update state vector to u + hv
386 Ptr<Vector<Real>> znew = z.clone();
387 znew->set(z);
388 znew->axpy(h,v);
389 // Compute new constraint value
390 update(u,*znew,UpdateType::Temp);
391 value(jv,u,*znew,ctol);
392 // Compute current constraint value
393 Ptr<Vector<Real>> cold = jv.clone();
394 update(u,z,UpdateType::Temp);
395 value(*cold,u,z,ctol);
396 // Compute Newton quotient
397 jv.axpy(-1.0,*cold);
398 jv.scale(1.0/h);
399 }
400
416 virtual void applyInverseJacobian_1(Vector<Real> &ijv,
417 const Vector<Real> &v,
418 const Vector<Real> &u,
419 const Vector<Real> &z,
420 Real &tol) {
421 ROL_TEST_FOR_EXCEPTION(true, std::logic_error,
422 "The method applyInverseJacobian_1 is used but not implemented!\n");
423 }
424
440 virtual void applyAdjointJacobian_1(Vector<Real> &ajv,
441 const Vector<Real> &v,
442 const Vector<Real> &u,
443 const Vector<Real> &z,
444 Real &tol) {
445 applyAdjointJacobian_1(ajv, v, u, z, v.dual(), tol);
446 }
447
448
466 virtual void applyAdjointJacobian_1(Vector<Real> &ajv,
467 const Vector<Real> &v,
468 const Vector<Real> &u,
469 const Vector<Real> &z,
470 const Vector<Real> &dualv,
471 Real &tol) {
472 Real ctol = std::sqrt(ROL_EPSILON<Real>());
473 Real h = tol;
474 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
475 h = std::max(1.0,u.norm()/v.norm())*tol;
476 }
477 Ptr<Vector<Real>> cold = dualv.clone();
478 Ptr<Vector<Real>> cnew = dualv.clone();
479 update(u,z,UpdateType::Temp);
480 value(*cold,u,z,ctol);
481 Ptr<Vector<Real>> unew = u.clone();
482 ajv.zero();
483 for (int i = 0; i < u.dimension(); i++) {
484 unew->set(u);
485 unew->axpy(h,*(u.basis(i)));
486 update(*unew,z,UpdateType::Temp);
487 value(*cnew,*unew,z,ctol);
488 cnew->axpy(-1.0,*cold);
489 cnew->scale(1.0/h);
490 ajv.axpy(cnew->dot(v),*((u.dual()).basis(i)));
491 }
492 update(u,z,UpdateType::Temp);
493 }
494
495
511 virtual void applyAdjointJacobian_2(Vector<Real> &ajv,
512 const Vector<Real> &v,
513 const Vector<Real> &u,
514 const Vector<Real> &z,
515 Real &tol) {
516 applyAdjointJacobian_2(ajv, v, u, z, v.dual(), tol);
517 }
518
519
537 virtual void applyAdjointJacobian_2(Vector<Real> &ajv,
538 const Vector<Real> &v,
539 const Vector<Real> &u,
540 const Vector<Real> &z,
541 const Vector<Real> &dualv,
542 Real &tol) {
543 Real ctol = std::sqrt(ROL_EPSILON<Real>());
544 Real h = tol;
545 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
546 h = std::max(1.0,u.norm()/v.norm())*tol;
547 }
548 Ptr<Vector<Real>> cold = dualv.clone();
549 Ptr<Vector<Real>> cnew = dualv.clone();
550 update(u,z,UpdateType::Temp);
551 value(*cold,u,z,ctol);
552 Ptr<Vector<Real>> znew = z.clone();
553 ajv.zero();
554 for (int i = 0; i < z.dimension(); i++) {
555 znew->set(z);
556 znew->axpy(h,*(z.basis(i)));
557 update(u,*znew,UpdateType::Temp);
558 value(*cnew,u,*znew,ctol);
559 cnew->axpy(-1.0,*cold);
560 cnew->scale(1.0/h);
561 ajv.axpy(cnew->dot(v),*((z.dual()).basis(i)));
562 }
563 update(u,z,UpdateType::Temp);
564 }
565
581 virtual void applyInverseAdjointJacobian_1(Vector<Real> &iajv,
582 const Vector<Real> &v,
583 const Vector<Real> &u,
584 const Vector<Real> &z,
585 Real &tol) {
586 ROL_TEST_FOR_EXCEPTION(true, std::logic_error,
587 "The method applyInverseAdjointJacobian_1 is used but not implemented!\n");
588 };
589
607 virtual void applyAdjointHessian_11(Vector<Real> &ahwv,
608 const Vector<Real> &w,
609 const Vector<Real> &v,
610 const Vector<Real> &u,
611 const Vector<Real> &z,
612 Real &tol) {
613 Real jtol = std::sqrt(ROL_EPSILON<Real>());
614 // Compute step size
615 Real h = tol;
616 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
617 h = std::max(1.0,u.norm()/v.norm())*tol;
618 }
619 // Evaluate Jacobian at new state
620 Ptr<Vector<Real>> unew = u.clone();
621 unew->set(u);
622 unew->axpy(h,v);
623 update(*unew,z,UpdateType::Temp);
624 applyAdjointJacobian_1(ahwv,w,*unew,z,jtol);
625 // Evaluate Jacobian at old state
626 Ptr<Vector<Real>> jv = ahwv.clone();
627 update(u,z,UpdateType::Temp);
628 applyAdjointJacobian_1(*jv,w,u,z,jtol);
629 // Compute Newton quotient
630 ahwv.axpy(-1.0,*jv);
631 ahwv.scale(1.0/h);
632 }
633
634
652 virtual void applyAdjointHessian_12(Vector<Real> &ahwv,
653 const Vector<Real> &w,
654 const Vector<Real> &v,
655 const Vector<Real> &u,
656 const Vector<Real> &z,
657 Real &tol) {
658 Real jtol = std::sqrt(ROL_EPSILON<Real>());
659 // Compute step size
660 Real h = tol;
661 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
662 h = std::max(1.0,u.norm()/v.norm())*tol;
663 }
664 // Evaluate Jacobian at new state
665 Ptr<Vector<Real>> unew = u.clone();
666 unew->set(u);
667 unew->axpy(h,v);
668 update(*unew,z,UpdateType::Temp);
669 applyAdjointJacobian_2(ahwv,w,*unew,z,jtol);
670 // Evaluate Jacobian at old state
671 Ptr<Vector<Real>> jv = ahwv.clone();
672 update(u,z,UpdateType::Temp);
673 applyAdjointJacobian_2(*jv,w,u,z,jtol);
674 // Compute Newton quotient
675 ahwv.axpy(-1.0,*jv);
676 ahwv.scale(1.0/h);
677 }
678
679
697 virtual void applyAdjointHessian_21(Vector<Real> &ahwv,
698 const Vector<Real> &w,
699 const Vector<Real> &v,
700 const Vector<Real> &u,
701 const Vector<Real> &z,
702 Real &tol) {
703 Real jtol = std::sqrt(ROL_EPSILON<Real>());
704 // Compute step size
705 Real h = tol;
706 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
707 h = std::max(1.0,u.norm()/v.norm())*tol;
708 }
709 // Evaluate Jacobian at new control
710 Ptr<Vector<Real>> znew = z.clone();
711 znew->set(z);
712 znew->axpy(h,v);
713 update(u,*znew,UpdateType::Temp);
714 applyAdjointJacobian_1(ahwv,w,u,*znew,jtol);
715 // Evaluate Jacobian at old control
716 Ptr<Vector<Real>> jv = ahwv.clone();
717 update(u,z,UpdateType::Temp);
718 applyAdjointJacobian_1(*jv,w,u,z,jtol);
719 // Compute Newton quotient
720 ahwv.axpy(-1.0,*jv);
721 ahwv.scale(1.0/h);
722 }
723
741 virtual void applyAdjointHessian_22(Vector<Real> &ahwv,
742 const Vector<Real> &w,
743 const Vector<Real> &v,
744 const Vector<Real> &u,
745 const Vector<Real> &z,
746 Real &tol) {
747 Real jtol = std::sqrt(ROL_EPSILON<Real>());
748 // Compute step size
749 Real h = tol;
750 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
751 h = std::max(1.0,u.norm()/v.norm())*tol;
752 }
753 // Evaluate Jacobian at new control
754 Ptr<Vector<Real>> znew = z.clone();
755 znew->set(z);
756 znew->axpy(h,v);
757 update(u,*znew,UpdateType::Temp);
758 applyAdjointJacobian_2(ahwv,w,u,*znew,jtol);
759 // Evaluate Jacobian at old control
760 Ptr<Vector<Real>> jv = ahwv.clone();
761 update(u,z,UpdateType::Temp);
762 applyAdjointJacobian_2(*jv,w,u,z,jtol);
763 // Compute Newton quotient
764 ahwv.axpy(-1.0,*jv);
765 ahwv.scale(1.0/h);
766}
767
806 virtual std::vector<Real> solveAugmentedSystem(Vector<Real> &v1,
807 Vector<Real> &v2,
808 const Vector<Real> &b1,
809 const Vector<Real> &b2,
810 const Vector<Real> &x,
811 Real &tol) {
812 return Constraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,x,tol);
813 }
814
815
834 virtual void applyPreconditioner(Vector<Real> &pv,
835 const Vector<Real> &v,
836 const Vector<Real> &x,
837 const Vector<Real> &g,
838 Real &tol) {
839 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(x);
840 Ptr<Vector<Real>> ijv = (xs.get_1())->clone();
841
842 try {
843 applyInverseJacobian_1(*ijv, v, *(xs.get_1()), *(xs.get_2()), tol);
844 }
845 catch (const std::logic_error &e) {
846 Constraint<Real>::applyPreconditioner(pv, v, x, g, tol);
847 return;
848 }
849
850 const Vector_SimOpt<Real> &gs = dynamic_cast<const Vector_SimOpt<Real>&>(g);
851 Ptr<Vector<Real>> ijv_dual = (gs.get_1())->clone();
852 ijv_dual->set(ijv->dual());
853
854 try {
855 applyInverseAdjointJacobian_1(pv, *ijv_dual, *(xs.get_1()), *(xs.get_2()), tol);
856 }
857 catch (const std::logic_error &e) {
858 Constraint<Real>::applyPreconditioner(pv, v, x, g, tol);
859 return;
860 }
861
862 }
863
869 virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
870 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
871 dynamic_cast<const Vector<Real>&>(x));
872 update(*(xs.get_1()),*(xs.get_2()),flag,iter);
873 }
874 virtual void update( const Vector<Real> &x, UpdateType type, int iter = -1 ) {
875 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
876 dynamic_cast<const Vector<Real>&>(x));
877 update(*(xs.get_1()),*(xs.get_2()),type,iter);
878 }
879
880 virtual void value(Vector<Real> &c,
881 const Vector<Real> &x,
882 Real &tol) {
883 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
884 dynamic_cast<const Vector<Real>&>(x));
885 value(c,*(xs.get_1()),*(xs.get_2()),tol);
886 }
887
888
889 virtual void applyJacobian(Vector<Real> &jv,
890 const Vector<Real> &v,
891 const Vector<Real> &x,
892 Real &tol) {
893 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
894 dynamic_cast<const Vector<Real>&>(x));
895 const Vector_SimOpt<Real> &vs = dynamic_cast<const Vector_SimOpt<Real>&>(
896 dynamic_cast<const Vector<Real>&>(v));
897 applyJacobian_1(jv,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
898 Ptr<Vector<Real>> jv2 = jv.clone();
899 applyJacobian_2(*jv2,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
900 jv.plus(*jv2);
901 }
902
903
904 using Constraint<Real>::applyAdjointJacobian;
905 virtual void applyAdjointJacobian(Vector<Real> &ajv,
906 const Vector<Real> &v,
907 const Vector<Real> &x,
908 Real &tol) {
909 Vector_SimOpt<Real> &ajvs = dynamic_cast<Vector_SimOpt<Real>&>(
910 dynamic_cast<Vector<Real>&>(ajv));
911 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
912 dynamic_cast<const Vector<Real>&>(x));
913 Ptr<Vector<Real>> ajv1 = (ajvs.get_1())->clone();
914 applyAdjointJacobian_1(*ajv1,v,*(xs.get_1()),*(xs.get_2()),tol);
915 ajvs.set_1(*ajv1);
916 Ptr<Vector<Real>> ajv2 = (ajvs.get_2())->clone();
917 applyAdjointJacobian_2(*ajv2,v,*(xs.get_1()),*(xs.get_2()),tol);
918 ajvs.set_2(*ajv2);
919 }
920
921
922 virtual void applyAdjointHessian(Vector<Real> &ahwv,
923 const Vector<Real> &w,
924 const Vector<Real> &v,
925 const Vector<Real> &x,
926 Real &tol) {
927 Vector_SimOpt<Real> &ahwvs = dynamic_cast<Vector_SimOpt<Real>&>(
928 dynamic_cast<Vector<Real>&>(ahwv));
929 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(
930 dynamic_cast<const Vector<Real>&>(x));
931 const Vector_SimOpt<Real> &vs = dynamic_cast<const Vector_SimOpt<Real>&>(
932 dynamic_cast<const Vector<Real>&>(v));
933 // Block-row 1
934 Ptr<Vector<Real>> C11 = (ahwvs.get_1())->clone();
935 Ptr<Vector<Real>> C21 = (ahwvs.get_1())->clone();
936 applyAdjointHessian_11(*C11,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
937 applyAdjointHessian_21(*C21,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
938 C11->plus(*C21);
939 ahwvs.set_1(*C11);
940 // Block-row 2
941 Ptr<Vector<Real>> C12 = (ahwvs.get_2())->clone();
942 Ptr<Vector<Real>> C22 = (ahwvs.get_2())->clone();
943 applyAdjointHessian_12(*C12,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
944 applyAdjointHessian_22(*C22,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
945 C22->plus(*C12);
946 ahwvs.set_2(*C22);
947 }
948
949
950
951 virtual Real checkSolve(const Vector<Real> &u,
952 const Vector<Real> &z,
953 const Vector<Real> &c,
954 const bool printToStream = true,
955 std::ostream & outStream = std::cout) {
956 // Solve constraint for u.
957 Real tol = ROL_EPSILON<Real>();
958 Ptr<Vector<Real>> r = c.clone();
959 Ptr<Vector<Real>> s = u.clone();
960 solve(*r,*s,z,tol);
961 // Evaluate constraint residual at (u,z).
962 Ptr<Vector<Real>> cs = c.clone();
963 update(*s,z,UpdateType::Temp);
964 value(*cs,*s,z,tol);
965 // Output norm of residual.
966 Real rnorm = r->norm();
967 Real cnorm = cs->norm();
968 if ( printToStream ) {
969 std::stringstream hist;
970 hist << std::scientific << std::setprecision(8);
971 hist << "\nTest SimOpt solve at feasible (u,z):\n";
972 hist << " Solver Residual = " << rnorm << "\n";
973 hist << " ||c(u,z)|| = " << cnorm << "\n";
974 outStream << hist.str();
975 }
976 return cnorm;
977 }
978
979
992 virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w,
993 const Vector<Real> &v,
994 const Vector<Real> &u,
995 const Vector<Real> &z,
996 const bool printToStream = true,
997 std::ostream & outStream = std::cout) {
998 return checkAdjointConsistencyJacobian_1(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
999 }
1000
1001
1017 virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w,
1018 const Vector<Real> &v,
1019 const Vector<Real> &u,
1020 const Vector<Real> &z,
1021 const Vector<Real> &dualw,
1022 const Vector<Real> &dualv,
1023 const bool printToStream = true,
1024 std::ostream & outStream = std::cout) {
1025 Real tol = ROL_EPSILON<Real>();
1026 Ptr<Vector<Real>> Jv = dualw.clone();
1027 update(u,z,UpdateType::Temp);
1028 applyJacobian_1(*Jv,v,u,z,tol);
1029 //Real wJv = w.dot(Jv->dual());
1030 Real wJv = w.apply(*Jv);
1031 Ptr<Vector<Real>> Jw = dualv.clone();
1032 update(u,z,UpdateType::Temp);
1033 applyAdjointJacobian_1(*Jw,w,u,z,tol);
1034 //Real vJw = v.dot(Jw->dual());
1035 Real vJw = v.apply(*Jw);
1036 Real diff = std::abs(wJv-vJw);
1037 if ( printToStream ) {
1038 std::stringstream hist;
1039 hist << std::scientific << std::setprecision(8);
1040 hist << "\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1041 << diff << "\n";
1042 hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1043 hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1044 outStream << hist.str();
1045 }
1046 return diff;
1047 }
1048
1049
1062 virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w,
1063 const Vector<Real> &v,
1064 const Vector<Real> &u,
1065 const Vector<Real> &z,
1066 const bool printToStream = true,
1067 std::ostream & outStream = std::cout) {
1068 return checkAdjointConsistencyJacobian_2(w, v, u, z, w.dual(), v.dual(), printToStream, outStream);
1069 }
1070
1086 virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w,
1087 const Vector<Real> &v,
1088 const Vector<Real> &u,
1089 const Vector<Real> &z,
1090 const Vector<Real> &dualw,
1091 const Vector<Real> &dualv,
1092 const bool printToStream = true,
1093 std::ostream & outStream = std::cout) {
1094 Real tol = ROL_EPSILON<Real>();
1095 Ptr<Vector<Real>> Jv = dualw.clone();
1096 update(u,z,UpdateType::Temp);
1097 applyJacobian_2(*Jv,v,u,z,tol);
1098 //Real wJv = w.dot(Jv->dual());
1099 Real wJv = w.apply(*Jv);
1100 Ptr<Vector<Real>> Jw = dualv.clone();
1101 update(u,z,UpdateType::Temp);
1102 applyAdjointJacobian_2(*Jw,w,u,z,tol);
1103 //Real vJw = v.dot(Jw->dual());
1104 Real vJw = v.apply(*Jw);
1105 Real diff = std::abs(wJv-vJw);
1106 if ( printToStream ) {
1107 std::stringstream hist;
1108 hist << std::scientific << std::setprecision(8);
1109 hist << "\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1110 << diff << "\n";
1111 hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
1112 hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
1113 outStream << hist.str();
1114 }
1115 return diff;
1116 }
1117
1118 virtual Real checkInverseJacobian_1(const Vector<Real> &jv,
1119 const Vector<Real> &v,
1120 const Vector<Real> &u,
1121 const Vector<Real> &z,
1122 const bool printToStream = true,
1123 std::ostream & outStream = std::cout) {
1124 Real tol = ROL_EPSILON<Real>();
1125 Ptr<Vector<Real>> Jv = jv.clone();
1126 update(u,z,UpdateType::Temp);
1127 applyJacobian_1(*Jv,v,u,z,tol);
1128 Ptr<Vector<Real>> iJJv = u.clone();
1129 //update(u,z); // Does this update do anything?
1130 applyInverseJacobian_1(*iJJv,*Jv,u,z,tol);
1131 Ptr<Vector<Real>> diff = v.clone();
1132 diff->set(v);
1133 diff->axpy(-1.0,*iJJv);
1134 Real dnorm = diff->norm();
1135 Real vnorm = v.norm();
1136 if ( printToStream ) {
1137 std::stringstream hist;
1138 hist << std::scientific << std::setprecision(8);
1139 hist << "\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = "
1140 << dnorm << "\n";
1141 hist << " ||v|| = " << vnorm << "\n";
1142 hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1143 outStream << hist.str();
1144 }
1145 return dnorm;
1146 }
1147
1148 virtual Real checkInverseAdjointJacobian_1(const Vector<Real> &jv,
1149 const Vector<Real> &v,
1150 const Vector<Real> &u,
1151 const Vector<Real> &z,
1152 const bool printToStream = true,
1153 std::ostream & outStream = std::cout) {
1154 Real tol = ROL_EPSILON<Real>();
1155 Ptr<Vector<Real>> Jv = jv.clone();
1156 update(u,z,UpdateType::Temp);
1157 applyAdjointJacobian_1(*Jv,v,u,z,tol);
1158 Ptr<Vector<Real>> iJJv = v.clone();
1159 //update(u,z);
1160 applyInverseAdjointJacobian_1(*iJJv,*Jv,u,z,tol);
1161 Ptr<Vector<Real>> diff = v.clone();
1162 diff->set(v);
1163 diff->axpy(-1.0,*iJJv);
1164 Real dnorm = diff->norm();
1165 Real vnorm = v.norm();
1166 if ( printToStream ) {
1167 std::stringstream hist;
1168 hist << std::scientific << std::setprecision(8);
1169 hist << "\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = "
1170 << dnorm << "\n";
1171 hist << " ||v|| = " << vnorm << "\n";
1172 hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n";
1173 outStream << hist.str();
1174 }
1175 return dnorm;
1176 }
1177
1178
1179
1180 std::vector<std::vector<Real>> checkApplyJacobian_1(const Vector<Real> &u,
1181 const Vector<Real> &z,
1182 const Vector<Real> &v,
1183 const Vector<Real> &jv,
1184 const bool printToStream = true,
1185 std::ostream & outStream = std::cout,
1186 const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1187 const int order = 1) {
1188 std::vector<Real> steps(numSteps);
1189 for(int i=0;i<numSteps;++i) {
1190 steps[i] = pow(10,-i);
1191 }
1192
1193 return checkApplyJacobian_1(u,z,v,jv,steps,printToStream,outStream,order);
1194 }
1195
1196
1197
1198
1199 std::vector<std::vector<Real>> checkApplyJacobian_1(const Vector<Real> &u,
1200 const Vector<Real> &z,
1201 const Vector<Real> &v,
1202 const Vector<Real> &jv,
1203 const std::vector<Real> &steps,
1204 const bool printToStream = true,
1205 std::ostream & outStream = std::cout,
1206 const int order = 1) {
1207
1208 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1209 "Error: finite difference order must be 1,2,3, or 4" );
1210
1211 Real one(1.0);
1212
1213 using Finite_Difference_Arrays::shifts;
1214 using Finite_Difference_Arrays::weights;
1215
1216 Real tol = std::sqrt(ROL_EPSILON<Real>());
1217
1218 int numSteps = steps.size();
1219 int numVals = 4;
1220 std::vector<Real> tmp(numVals);
1221 std::vector<std::vector<Real>> jvCheck(numSteps, tmp);
1222
1223 // Save the format state of the original outStream.
1224 nullstream oldFormatState;
1225 oldFormatState.copyfmt(outStream);
1226
1227 // Compute constraint value at x.
1228 Ptr<Vector<Real>> c = jv.clone();
1229 this->update(u,z,UpdateType::Temp);
1230 this->value(*c, u, z, tol);
1231
1232 // Compute (Jacobian at x) times (vector v).
1233 Ptr<Vector<Real>> Jv = jv.clone();
1234 this->applyJacobian_1(*Jv, v, u, z, tol);
1235 Real normJv = Jv->norm();
1236
1237 // Temporary vectors.
1238 Ptr<Vector<Real>> cdif = jv.clone();
1239 Ptr<Vector<Real>> cnew = jv.clone();
1240 Ptr<Vector<Real>> unew = u.clone();
1241
1242 for (int i=0; i<numSteps; i++) {
1243
1244 Real eta = steps[i];
1245
1246 unew->set(u);
1247
1248 cdif->set(*c);
1249 cdif->scale(weights[order-1][0]);
1250
1251 for(int j=0; j<order; ++j) {
1252
1253 unew->axpy(eta*shifts[order-1][j], v);
1254
1255 if( weights[order-1][j+1] != 0 ) {
1256 this->update(*unew,z,UpdateType::Temp);
1257 this->value(*cnew,*unew,z,tol);
1258 cdif->axpy(weights[order-1][j+1],*cnew);
1259 }
1260
1261 }
1262
1263 cdif->scale(one/eta);
1264
1265 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1266 jvCheck[i][0] = eta;
1267 jvCheck[i][1] = normJv;
1268 jvCheck[i][2] = cdif->norm();
1269 cdif->axpy(-one, *Jv);
1270 jvCheck[i][3] = cdif->norm();
1271
1272 if (printToStream) {
1273 std::stringstream hist;
1274 if (i==0) {
1275 hist << std::right
1276 << std::setw(20) << "Step size"
1277 << std::setw(20) << "norm(Jac*vec)"
1278 << std::setw(20) << "norm(FD approx)"
1279 << std::setw(20) << "norm(abs error)"
1280 << "\n"
1281 << std::setw(20) << "---------"
1282 << std::setw(20) << "-------------"
1283 << std::setw(20) << "---------------"
1284 << std::setw(20) << "---------------"
1285 << "\n";
1286 }
1287 hist << std::scientific << std::setprecision(11) << std::right
1288 << std::setw(20) << jvCheck[i][0]
1289 << std::setw(20) << jvCheck[i][1]
1290 << std::setw(20) << jvCheck[i][2]
1291 << std::setw(20) << jvCheck[i][3]
1292 << "\n";
1293 outStream << hist.str();
1294 }
1295
1296 }
1297
1298 // Reset format state of outStream.
1299 outStream.copyfmt(oldFormatState);
1300
1301 return jvCheck;
1302 } // checkApplyJacobian
1303
1304
1305 std::vector<std::vector<Real>> checkApplyJacobian_2(const Vector<Real> &u,
1306 const Vector<Real> &z,
1307 const Vector<Real> &v,
1308 const Vector<Real> &jv,
1309 const bool printToStream = true,
1310 std::ostream & outStream = std::cout,
1311 const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1312 const int order = 1) {
1313 std::vector<Real> steps(numSteps);
1314 for(int i=0;i<numSteps;++i) {
1315 steps[i] = pow(10,-i);
1316 }
1317
1318 return checkApplyJacobian_2(u,z,v,jv,steps,printToStream,outStream,order);
1319 }
1320
1321
1322
1323
1324 std::vector<std::vector<Real>> checkApplyJacobian_2(const Vector<Real> &u,
1325 const Vector<Real> &z,
1326 const Vector<Real> &v,
1327 const Vector<Real> &jv,
1328 const std::vector<Real> &steps,
1329 const bool printToStream = true,
1330 std::ostream & outStream = std::cout,
1331 const int order = 1) {
1332
1333 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1334 "Error: finite difference order must be 1,2,3, or 4" );
1335
1336 Real one(1.0);
1337
1338 using Finite_Difference_Arrays::shifts;
1339 using Finite_Difference_Arrays::weights;
1340
1341 Real tol = std::sqrt(ROL_EPSILON<Real>());
1342
1343 int numSteps = steps.size();
1344 int numVals = 4;
1345 std::vector<Real> tmp(numVals);
1346 std::vector<std::vector<Real>> jvCheck(numSteps, tmp);
1347
1348 // Save the format state of the original outStream.
1349 nullstream oldFormatState;
1350 oldFormatState.copyfmt(outStream);
1351
1352 // Compute constraint value at x.
1353 Ptr<Vector<Real>> c = jv.clone();
1354 this->update(u,z,UpdateType::Temp);
1355 this->value(*c, u, z, tol);
1356
1357 // Compute (Jacobian at x) times (vector v).
1358 Ptr<Vector<Real>> Jv = jv.clone();
1359 this->applyJacobian_2(*Jv, v, u, z, tol);
1360 Real normJv = Jv->norm();
1361
1362 // Temporary vectors.
1363 Ptr<Vector<Real>> cdif = jv.clone();
1364 Ptr<Vector<Real>> cnew = jv.clone();
1365 Ptr<Vector<Real>> znew = z.clone();
1366
1367 for (int i=0; i<numSteps; i++) {
1368
1369 Real eta = steps[i];
1370
1371 znew->set(z);
1372
1373 cdif->set(*c);
1374 cdif->scale(weights[order-1][0]);
1375
1376 for(int j=0; j<order; ++j) {
1377
1378 znew->axpy(eta*shifts[order-1][j], v);
1379
1380 if( weights[order-1][j+1] != 0 ) {
1381 this->update(u,*znew,UpdateType::Temp);
1382 this->value(*cnew,u,*znew,tol);
1383 cdif->axpy(weights[order-1][j+1],*cnew);
1384 }
1385
1386 }
1387
1388 cdif->scale(one/eta);
1389
1390 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1391 jvCheck[i][0] = eta;
1392 jvCheck[i][1] = normJv;
1393 jvCheck[i][2] = cdif->norm();
1394 cdif->axpy(-one, *Jv);
1395 jvCheck[i][3] = cdif->norm();
1396
1397 if (printToStream) {
1398 std::stringstream hist;
1399 if (i==0) {
1400 hist << std::right
1401 << std::setw(20) << "Step size"
1402 << std::setw(20) << "norm(Jac*vec)"
1403 << std::setw(20) << "norm(FD approx)"
1404 << std::setw(20) << "norm(abs error)"
1405 << "\n"
1406 << std::setw(20) << "---------"
1407 << std::setw(20) << "-------------"
1408 << std::setw(20) << "---------------"
1409 << std::setw(20) << "---------------"
1410 << "\n";
1411 }
1412 hist << std::scientific << std::setprecision(11) << std::right
1413 << std::setw(20) << jvCheck[i][0]
1414 << std::setw(20) << jvCheck[i][1]
1415 << std::setw(20) << jvCheck[i][2]
1416 << std::setw(20) << jvCheck[i][3]
1417 << "\n";
1418 outStream << hist.str();
1419 }
1420
1421 }
1422
1423 // Reset format state of outStream.
1424 outStream.copyfmt(oldFormatState);
1425
1426 return jvCheck;
1427 } // checkApplyJacobian
1428
1429
1430
1431 std::vector<std::vector<Real>> checkApplyAdjointHessian_11(const Vector<Real> &u,
1432 const Vector<Real> &z,
1433 const Vector<Real> &p,
1434 const Vector<Real> &v,
1435 const Vector<Real> &hv,
1436 const bool printToStream = true,
1437 std::ostream & outStream = std::cout,
1438 const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1439 const int order = 1 ) {
1440 std::vector<Real> steps(numSteps);
1441 for(int i=0;i<numSteps;++i) {
1442 steps[i] = pow(10,-i);
1443 }
1444
1445 return checkApplyAdjointHessian_11(u,z,p,v,hv,steps,printToStream,outStream,order);
1446
1447 }
1448
1449 std::vector<std::vector<Real>> checkApplyAdjointHessian_11(const Vector<Real> &u,
1450 const Vector<Real> &z,
1451 const Vector<Real> &p,
1452 const Vector<Real> &v,
1453 const Vector<Real> &hv,
1454 const std::vector<Real> &steps,
1455 const bool printToStream = true,
1456 std::ostream & outStream = std::cout,
1457 const int order = 1 ) {
1458 using Finite_Difference_Arrays::shifts;
1459 using Finite_Difference_Arrays::weights;
1460
1461 Real one(1.0);
1462
1463 Real tol = std::sqrt(ROL_EPSILON<Real>());
1464
1465 int numSteps = steps.size();
1466 int numVals = 4;
1467 std::vector<Real> tmp(numVals);
1468 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1469
1470 // Temporary vectors.
1471 Ptr<Vector<Real>> AJdif = hv.clone();
1472 Ptr<Vector<Real>> AJp = hv.clone();
1473 Ptr<Vector<Real>> AHpv = hv.clone();
1474 Ptr<Vector<Real>> AJnew = hv.clone();
1475 Ptr<Vector<Real>> unew = u.clone();
1476
1477 // Save the format state of the original outStream.
1478 nullstream oldFormatState;
1479 oldFormatState.copyfmt(outStream);
1480
1481 // Apply adjoint Jacobian to p.
1482 this->update(u,z,UpdateType::Temp);
1483 this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1484
1485 // Apply adjoint Hessian at (u,z), in direction v, to p.
1486 this->applyAdjointHessian_11(*AHpv, p, v, u, z, tol);
1487 Real normAHpv = AHpv->norm();
1488
1489 for (int i=0; i<numSteps; i++) {
1490
1491 Real eta = steps[i];
1492
1493 // Apply adjoint Jacobian to p at (u+eta*v,z).
1494 unew->set(u);
1495
1496 AJdif->set(*AJp);
1497 AJdif->scale(weights[order-1][0]);
1498
1499 for(int j=0; j<order; ++j) {
1500
1501 unew->axpy(eta*shifts[order-1][j],v);
1502
1503 if( weights[order-1][j+1] != 0 ) {
1504 this->update(*unew,z,UpdateType::Temp);
1505 this->applyAdjointJacobian_1(*AJnew, p, *unew, z, tol);
1506 AJdif->axpy(weights[order-1][j+1],*AJnew);
1507 }
1508 }
1509
1510 AJdif->scale(one/eta);
1511
1512 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1513 ahpvCheck[i][0] = eta;
1514 ahpvCheck[i][1] = normAHpv;
1515 ahpvCheck[i][2] = AJdif->norm();
1516 AJdif->axpy(-one, *AHpv);
1517 ahpvCheck[i][3] = AJdif->norm();
1518
1519 if (printToStream) {
1520 std::stringstream hist;
1521 if (i==0) {
1522 hist << std::right
1523 << std::setw(20) << "Step size"
1524 << std::setw(20) << "norm(adj(H)(u,v))"
1525 << std::setw(20) << "norm(FD approx)"
1526 << std::setw(20) << "norm(abs error)"
1527 << "\n"
1528 << std::setw(20) << "---------"
1529 << std::setw(20) << "-----------------"
1530 << std::setw(20) << "---------------"
1531 << std::setw(20) << "---------------"
1532 << "\n";
1533 }
1534 hist << std::scientific << std::setprecision(11) << std::right
1535 << std::setw(20) << ahpvCheck[i][0]
1536 << std::setw(20) << ahpvCheck[i][1]
1537 << std::setw(20) << ahpvCheck[i][2]
1538 << std::setw(20) << ahpvCheck[i][3]
1539 << "\n";
1540 outStream << hist.str();
1541 }
1542
1543 }
1544
1545 // Reset format state of outStream.
1546 outStream.copyfmt(oldFormatState);
1547
1548 return ahpvCheck;
1549 } // checkApplyAdjointHessian_11
1550
1554 std::vector<std::vector<Real>> checkApplyAdjointHessian_21(const Vector<Real> &u,
1555 const Vector<Real> &z,
1556 const Vector<Real> &p,
1557 const Vector<Real> &v,
1558 const Vector<Real> &hv,
1559 const bool printToStream = true,
1560 std::ostream & outStream = std::cout,
1561 const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1562 const int order = 1 ) {
1563 std::vector<Real> steps(numSteps);
1564 for(int i=0;i<numSteps;++i) {
1565 steps[i] = pow(10,-i);
1566 }
1567
1568 return checkApplyAdjointHessian_21(u,z,p,v,hv,steps,printToStream,outStream,order);
1569
1570 }
1571
1575 std::vector<std::vector<Real>> checkApplyAdjointHessian_21(const Vector<Real> &u,
1576 const Vector<Real> &z,
1577 const Vector<Real> &p,
1578 const Vector<Real> &v,
1579 const Vector<Real> &hv,
1580 const std::vector<Real> &steps,
1581 const bool printToStream = true,
1582 std::ostream & outStream = std::cout,
1583 const int order = 1 ) {
1584 using Finite_Difference_Arrays::shifts;
1585 using Finite_Difference_Arrays::weights;
1586
1587 Real one(1.0);
1588
1589 Real tol = std::sqrt(ROL_EPSILON<Real>());
1590
1591 int numSteps = steps.size();
1592 int numVals = 4;
1593 std::vector<Real> tmp(numVals);
1594 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1595
1596 // Temporary vectors.
1597 Ptr<Vector<Real>> AJdif = hv.clone();
1598 Ptr<Vector<Real>> AJp = hv.clone();
1599 Ptr<Vector<Real>> AHpv = hv.clone();
1600 Ptr<Vector<Real>> AJnew = hv.clone();
1601 Ptr<Vector<Real>> znew = z.clone();
1602
1603 // Save the format state of the original outStream.
1604 nullstream oldFormatState;
1605 oldFormatState.copyfmt(outStream);
1606
1607 // Apply adjoint Jacobian to p.
1608 this->update(u,z,UpdateType::Temp);
1609 this->applyAdjointJacobian_1(*AJp, p, u, z, tol);
1610
1611 // Apply adjoint Hessian at (u,z), in direction v, to p.
1612 this->applyAdjointHessian_21(*AHpv, p, v, u, z, tol);
1613 Real normAHpv = AHpv->norm();
1614
1615 for (int i=0; i<numSteps; i++) {
1616
1617 Real eta = steps[i];
1618
1619 // Apply adjoint Jacobian to p at (u,z+eta*v).
1620 znew->set(z);
1621
1622 AJdif->set(*AJp);
1623 AJdif->scale(weights[order-1][0]);
1624
1625 for(int j=0; j<order; ++j) {
1626
1627 znew->axpy(eta*shifts[order-1][j],v);
1628
1629 if( weights[order-1][j+1] != 0 ) {
1630 this->update(u,*znew,UpdateType::Temp);
1631 this->applyAdjointJacobian_1(*AJnew, p, u, *znew, tol);
1632 AJdif->axpy(weights[order-1][j+1],*AJnew);
1633 }
1634 }
1635
1636 AJdif->scale(one/eta);
1637
1638 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1639 ahpvCheck[i][0] = eta;
1640 ahpvCheck[i][1] = normAHpv;
1641 ahpvCheck[i][2] = AJdif->norm();
1642 AJdif->axpy(-one, *AHpv);
1643 ahpvCheck[i][3] = AJdif->norm();
1644
1645 if (printToStream) {
1646 std::stringstream hist;
1647 if (i==0) {
1648 hist << std::right
1649 << std::setw(20) << "Step size"
1650 << std::setw(20) << "norm(adj(H)(u,v))"
1651 << std::setw(20) << "norm(FD approx)"
1652 << std::setw(20) << "norm(abs error)"
1653 << "\n"
1654 << std::setw(20) << "---------"
1655 << std::setw(20) << "-----------------"
1656 << std::setw(20) << "---------------"
1657 << std::setw(20) << "---------------"
1658 << "\n";
1659 }
1660 hist << std::scientific << std::setprecision(11) << std::right
1661 << std::setw(20) << ahpvCheck[i][0]
1662 << std::setw(20) << ahpvCheck[i][1]
1663 << std::setw(20) << ahpvCheck[i][2]
1664 << std::setw(20) << ahpvCheck[i][3]
1665 << "\n";
1666 outStream << hist.str();
1667 }
1668
1669 }
1670
1671 // Reset format state of outStream.
1672 outStream.copyfmt(oldFormatState);
1673
1674 return ahpvCheck;
1675 } // checkApplyAdjointHessian_21
1676
1680 std::vector<std::vector<Real>> checkApplyAdjointHessian_12(const Vector<Real> &u,
1681 const Vector<Real> &z,
1682 const Vector<Real> &p,
1683 const Vector<Real> &v,
1684 const Vector<Real> &hv,
1685 const bool printToStream = true,
1686 std::ostream & outStream = std::cout,
1687 const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1688 const int order = 1 ) {
1689 std::vector<Real> steps(numSteps);
1690 for(int i=0;i<numSteps;++i) {
1691 steps[i] = pow(10,-i);
1692 }
1693
1694 return checkApplyAdjointHessian_12(u,z,p,v,hv,steps,printToStream,outStream,order);
1695
1696 }
1697
1698
1699 std::vector<std::vector<Real>> checkApplyAdjointHessian_12(const Vector<Real> &u,
1700 const Vector<Real> &z,
1701 const Vector<Real> &p,
1702 const Vector<Real> &v,
1703 const Vector<Real> &hv,
1704 const std::vector<Real> &steps,
1705 const bool printToStream = true,
1706 std::ostream & outStream = std::cout,
1707 const int order = 1 ) {
1708 using Finite_Difference_Arrays::shifts;
1709 using Finite_Difference_Arrays::weights;
1710
1711 Real one(1.0);
1712
1713 Real tol = std::sqrt(ROL_EPSILON<Real>());
1714
1715 int numSteps = steps.size();
1716 int numVals = 4;
1717 std::vector<Real> tmp(numVals);
1718 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1719
1720 // Temporary vectors.
1721 Ptr<Vector<Real>> AJdif = hv.clone();
1722 Ptr<Vector<Real>> AJp = hv.clone();
1723 Ptr<Vector<Real>> AHpv = hv.clone();
1724 Ptr<Vector<Real>> AJnew = hv.clone();
1725 Ptr<Vector<Real>> unew = u.clone();
1726
1727 // Save the format state of the original outStream.
1728 nullstream oldFormatState;
1729 oldFormatState.copyfmt(outStream);
1730
1731 // Apply adjoint Jacobian to p.
1732 this->update(u,z,UpdateType::Temp);
1733 this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1734
1735 // Apply adjoint Hessian at (u,z), in direction v, to p.
1736 this->applyAdjointHessian_12(*AHpv, p, v, u, z, tol);
1737 Real normAHpv = AHpv->norm();
1738
1739 for (int i=0; i<numSteps; i++) {
1740
1741 Real eta = steps[i];
1742
1743 // Apply adjoint Jacobian to p at (u+eta*v,z).
1744 unew->set(u);
1745
1746 AJdif->set(*AJp);
1747 AJdif->scale(weights[order-1][0]);
1748
1749 for(int j=0; j<order; ++j) {
1750
1751 unew->axpy(eta*shifts[order-1][j],v);
1752
1753 if( weights[order-1][j+1] != 0 ) {
1754 this->update(*unew,z,UpdateType::Temp);
1755 this->applyAdjointJacobian_2(*AJnew, p, *unew, z, tol);
1756 AJdif->axpy(weights[order-1][j+1],*AJnew);
1757 }
1758 }
1759
1760 AJdif->scale(one/eta);
1761
1762 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1763 ahpvCheck[i][0] = eta;
1764 ahpvCheck[i][1] = normAHpv;
1765 ahpvCheck[i][2] = AJdif->norm();
1766 AJdif->axpy(-one, *AHpv);
1767 ahpvCheck[i][3] = AJdif->norm();
1768
1769 if (printToStream) {
1770 std::stringstream hist;
1771 if (i==0) {
1772 hist << std::right
1773 << std::setw(20) << "Step size"
1774 << std::setw(20) << "norm(adj(H)(u,v))"
1775 << std::setw(20) << "norm(FD approx)"
1776 << std::setw(20) << "norm(abs error)"
1777 << "\n"
1778 << std::setw(20) << "---------"
1779 << std::setw(20) << "-----------------"
1780 << std::setw(20) << "---------------"
1781 << std::setw(20) << "---------------"
1782 << "\n";
1783 }
1784 hist << std::scientific << std::setprecision(11) << std::right
1785 << std::setw(20) << ahpvCheck[i][0]
1786 << std::setw(20) << ahpvCheck[i][1]
1787 << std::setw(20) << ahpvCheck[i][2]
1788 << std::setw(20) << ahpvCheck[i][3]
1789 << "\n";
1790 outStream << hist.str();
1791 }
1792
1793 }
1794
1795 // Reset format state of outStream.
1796 outStream.copyfmt(oldFormatState);
1797
1798 return ahpvCheck;
1799 } // checkApplyAdjointHessian_12
1800
1801 std::vector<std::vector<Real>> checkApplyAdjointHessian_22(const Vector<Real> &u,
1802 const Vector<Real> &z,
1803 const Vector<Real> &p,
1804 const Vector<Real> &v,
1805 const Vector<Real> &hv,
1806 const bool printToStream = true,
1807 std::ostream & outStream = std::cout,
1808 const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
1809 const int order = 1 ) {
1810 std::vector<Real> steps(numSteps);
1811 for(int i=0;i<numSteps;++i) {
1812 steps[i] = pow(10,-i);
1813 }
1814
1815 return checkApplyAdjointHessian_22(u,z,p,v,hv,steps,printToStream,outStream,order);
1816
1817 }
1818
1819 std::vector<std::vector<Real>> checkApplyAdjointHessian_22(const Vector<Real> &u,
1820 const Vector<Real> &z,
1821 const Vector<Real> &p,
1822 const Vector<Real> &v,
1823 const Vector<Real> &hv,
1824 const std::vector<Real> &steps,
1825 const bool printToStream = true,
1826 std::ostream & outStream = std::cout,
1827 const int order = 1 ) {
1828 using Finite_Difference_Arrays::shifts;
1829 using Finite_Difference_Arrays::weights;
1830
1831 Real one(1.0);
1832
1833 Real tol = std::sqrt(ROL_EPSILON<Real>());
1834
1835 int numSteps = steps.size();
1836 int numVals = 4;
1837 std::vector<Real> tmp(numVals);
1838 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1839
1840 // Temporary vectors.
1841 Ptr<Vector<Real>> AJdif = hv.clone();
1842 Ptr<Vector<Real>> AJp = hv.clone();
1843 Ptr<Vector<Real>> AHpv = hv.clone();
1844 Ptr<Vector<Real>> AJnew = hv.clone();
1845 Ptr<Vector<Real>> znew = z.clone();
1846
1847 // Save the format state of the original outStream.
1848 nullstream oldFormatState;
1849 oldFormatState.copyfmt(outStream);
1850
1851 // Apply adjoint Jacobian to p.
1852 this->update(u,z,UpdateType::Temp);
1853 this->applyAdjointJacobian_2(*AJp, p, u, z, tol);
1854
1855 // Apply adjoint Hessian at (u,z), in direction v, to p.
1856 this->applyAdjointHessian_22(*AHpv, p, v, u, z, tol);
1857 Real normAHpv = AHpv->norm();
1858
1859 for (int i=0; i<numSteps; i++) {
1860
1861 Real eta = steps[i];
1862
1863 // Apply adjoint Jacobian to p at (u,z+eta*v).
1864 znew->set(z);
1865
1866 AJdif->set(*AJp);
1867 AJdif->scale(weights[order-1][0]);
1868
1869 for(int j=0; j<order; ++j) {
1870
1871 znew->axpy(eta*shifts[order-1][j],v);
1872
1873 if( weights[order-1][j+1] != 0 ) {
1874 this->update(u,*znew,UpdateType::Temp);
1875 this->applyAdjointJacobian_2(*AJnew, p, u, *znew, tol);
1876 AJdif->axpy(weights[order-1][j+1],*AJnew);
1877 }
1878 }
1879
1880 AJdif->scale(one/eta);
1881
1882 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
1883 ahpvCheck[i][0] = eta;
1884 ahpvCheck[i][1] = normAHpv;
1885 ahpvCheck[i][2] = AJdif->norm();
1886 AJdif->axpy(-one, *AHpv);
1887 ahpvCheck[i][3] = AJdif->norm();
1888
1889 if (printToStream) {
1890 std::stringstream hist;
1891 if (i==0) {
1892 hist << std::right
1893 << std::setw(20) << "Step size"
1894 << std::setw(20) << "norm(adj(H)(u,v))"
1895 << std::setw(20) << "norm(FD approx)"
1896 << std::setw(20) << "norm(abs error)"
1897 << "\n"
1898 << std::setw(20) << "---------"
1899 << std::setw(20) << "-----------------"
1900 << std::setw(20) << "---------------"
1901 << std::setw(20) << "---------------"
1902 << "\n";
1903 }
1904 hist << std::scientific << std::setprecision(11) << std::right
1905 << std::setw(20) << ahpvCheck[i][0]
1906 << std::setw(20) << ahpvCheck[i][1]
1907 << std::setw(20) << ahpvCheck[i][2]
1908 << std::setw(20) << ahpvCheck[i][3]
1909 << "\n";
1910 outStream << hist.str();
1911 }
1912
1913 }
1914
1915 // Reset format state of outStream.
1916 outStream.copyfmt(oldFormatState);
1917
1918 return ahpvCheck;
1919 } // checkApplyAdjointHessian_22
1920
1921}; // class Constraint_SimOpt
1922
1923} // namespace ROL
1924
1925#endif
ROL_Constraint.hpp
ROL_NonlinearLeastSquaresObjective.hpp
ROL_Objective_FSsolver.hpp
ROL_SimConstraint.hpp
ROL_TypeE_CompositeStepAlgorithm.hpp
ROL_TypeU_TrustRegionAlgorithm.hpp
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL_NUM_CHECKDERIV_STEPS
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
Definition ROL_Types.hpp:40
ROL_Vector_SimOpt.hpp
ROL::Constraint_SimOpt
Defines the constraint operator interface for simulation-based optimization.
Definition ROL_Constraint_SimOpt.hpp:74
ROL::Constraint_SimOpt::update
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable i...
Definition ROL_Constraint_SimOpt.hpp:129
ROL::Constraint_SimOpt::applyAdjointJacobian_2
virtual void applyAdjointJacobian_2(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
Definition ROL_Constraint_SimOpt.hpp:511
ROL::Constraint_SimOpt::update_2
virtual void update_2(const Vector< Real > &z, UpdateType type, int iter=-1)
Definition ROL_Constraint_SimOpt.hpp:152
ROL::Constraint_SimOpt::checkInverseJacobian_1
virtual Real checkInverseJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Definition ROL_Constraint_SimOpt.hpp:1118
ROL::Constraint_SimOpt::checkAdjointConsistencyJacobian_1
virtual Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
Definition ROL_Constraint_SimOpt.hpp:992
ROL::Constraint_SimOpt::checkApplyJacobian_1
std::vector< std::vector< Real > > checkApplyJacobian_1(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1180
ROL::Constraint_SimOpt::firstSolve_
bool firstSolve_
Definition ROL_Constraint_SimOpt.hpp:105
ROL::Constraint_SimOpt::unew_
Ptr< Vector< Real > > unew_
Definition ROL_Constraint_SimOpt.hpp:77
ROL::Constraint_SimOpt::decr_
Real decr_
Definition ROL_Constraint_SimOpt.hpp:98
ROL::Constraint_SimOpt::rtol_
Real rtol_
Definition ROL_Constraint_SimOpt.hpp:95
ROL::Constraint_SimOpt::applyInverseJacobian_1
virtual void applyInverseJacobian_1(Vector< Real > &ijv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the inverse partial constraint Jacobian at , , to the vector .
Definition ROL_Constraint_SimOpt.hpp:416
ROL::Constraint_SimOpt::checkAdjointConsistencyJacobian_2
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualw, const Vector< Real > &dualv, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the secondary interface,...
Definition ROL_Constraint_SimOpt.hpp:1086
ROL::Constraint_SimOpt::applyInverseAdjointJacobian_1
virtual void applyInverseAdjointJacobian_1(Vector< Real > &iajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector .
Definition ROL_Constraint_SimOpt.hpp:581
ROL::Constraint_SimOpt::DEFAULT_decr_
const Real DEFAULT_decr_
Definition ROL_Constraint_SimOpt.hpp:85
ROL::Constraint_SimOpt::applyJacobian_2
virtual void applyJacobian_2(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
Definition ROL_Constraint_SimOpt.hpp:374
ROL::Constraint_SimOpt::zero_
bool zero_
Definition ROL_Constraint_SimOpt.hpp:101
ROL::Constraint_SimOpt::maxit_
int maxit_
Definition ROL_Constraint_SimOpt.hpp:99
ROL::Constraint_SimOpt::applyAdjointHessian_12
virtual void applyAdjointHessian_12(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
Definition ROL_Constraint_SimOpt.hpp:652
ROL::Constraint_SimOpt::applyAdjointHessian_11
virtual void applyAdjointHessian_11(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
Definition ROL_Constraint_SimOpt.hpp:607
ROL::Constraint_SimOpt::atol_
Real atol_
Definition ROL_Constraint_SimOpt.hpp:94
ROL::Constraint_SimOpt::update_1
virtual void update_1(const Vector< Real > &u, bool flag=true, int iter=-1)
Update constraint functions with respect to Sim variable. x is the optimization variable,...
Definition ROL_Constraint_SimOpt.hpp:143
ROL::Constraint_SimOpt::update_2
virtual void update_2(const Vector< Real > &z, bool flag=true, int iter=-1)
Update constraint functions with respect to Opt variable. x is the optimization variable,...
Definition ROL_Constraint_SimOpt.hpp:151
ROL::Constraint_SimOpt::applyAdjointHessian_22
virtual void applyAdjointHessian_22(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
Definition ROL_Constraint_SimOpt.hpp:741
ROL::Constraint_SimOpt::value
virtual void value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol)=0
Evaluate the constraint operator at .
ROL::Constraint_SimOpt::applyAdjointJacobian_1
virtual void applyAdjointJacobian_1(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualv, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the secondary int...
Definition ROL_Constraint_SimOpt.hpp:466
ROL::Constraint_SimOpt::applyAdjointJacobian_2
virtual void applyAdjointJacobian_2(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualv, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the secondary interfa...
Definition ROL_Constraint_SimOpt.hpp:537
ROL::Constraint_SimOpt::checkApplyJacobian_2
std::vector< std::vector< Real > > checkApplyJacobian_2(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1305
ROL::Constraint_SimOpt::checkApplyAdjointHessian_21
std::vector< std::vector< Real > > checkApplyAdjointHessian_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
, , , ,
Definition ROL_Constraint_SimOpt.hpp:1575
ROL::Constraint_SimOpt::checkAdjointConsistencyJacobian_1
virtual Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualw, const Vector< Real > &dualv, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the secondary interface,...
Definition ROL_Constraint_SimOpt.hpp:1017
ROL::Constraint_SimOpt::update
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update constraint function.
Definition ROL_Constraint_SimOpt.hpp:874
ROL::Constraint_SimOpt::applyAdjointJacobian_1
virtual void applyAdjointJacobian_1(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
Definition ROL_Constraint_SimOpt.hpp:440
ROL::Constraint_SimOpt::solverType_
int solverType_
Definition ROL_Constraint_SimOpt.hpp:102
ROL::Constraint_SimOpt::checkApplyAdjointHessian_22
std::vector< std::vector< Real > > checkApplyAdjointHessian_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1819
ROL::Constraint_SimOpt::DEFAULT_rtol_
const Real DEFAULT_rtol_
Definition ROL_Constraint_SimOpt.hpp:82
ROL::Constraint_SimOpt::applyAdjointJacobian
virtual void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
Definition ROL_Constraint_SimOpt.hpp:905
ROL::Constraint_SimOpt::checkAdjointConsistencyJacobian_2
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
Definition ROL_Constraint_SimOpt.hpp:1062
ROL::Constraint_SimOpt::solveAugmentedSystem
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system
Definition ROL_Constraint_SimOpt.hpp:806
ROL::Constraint_SimOpt::applyAdjointHessian_21
virtual void applyAdjointHessian_21(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
Definition ROL_Constraint_SimOpt.hpp:697
ROL::Constraint_SimOpt::DEFAULT_solverType_
const int DEFAULT_solverType_
Definition ROL_Constraint_SimOpt.hpp:89
ROL::Constraint_SimOpt::checkApplyAdjointHessian_11
std::vector< std::vector< Real > > checkApplyAdjointHessian_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1449
ROL::Constraint_SimOpt::setSolveParameters
virtual void setSolveParameters(ParameterList &parlist)
Set solve parameters.
Definition ROL_Constraint_SimOpt.hpp:303
ROL::Constraint_SimOpt::checkInverseAdjointJacobian_1
virtual Real checkInverseAdjointJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Definition ROL_Constraint_SimOpt.hpp:1148
ROL::Constraint_SimOpt::DEFAULT_factor_
const Real DEFAULT_factor_
Definition ROL_Constraint_SimOpt.hpp:84
ROL::Constraint_SimOpt::checkApplyAdjointHessian_21
std::vector< std::vector< Real > > checkApplyAdjointHessian_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
, , , ,
Definition ROL_Constraint_SimOpt.hpp:1554
ROL::Constraint_SimOpt::checkApplyAdjointHessian_11
std::vector< std::vector< Real > > checkApplyAdjointHessian_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1431
ROL::Constraint_SimOpt::DEFAULT_atol_
const Real DEFAULT_atol_
Definition ROL_Constraint_SimOpt.hpp:81
ROL::Constraint_SimOpt::DEFAULT_zero_
const bool DEFAULT_zero_
Definition ROL_Constraint_SimOpt.hpp:88
ROL::Constraint_SimOpt::print_
bool print_
Definition ROL_Constraint_SimOpt.hpp:100
ROL::Constraint_SimOpt::applyJacobian_1
virtual void applyJacobian_1(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
Definition ROL_Constraint_SimOpt.hpp:331
ROL::Constraint_SimOpt::update_1
virtual void update_1(const Vector< Real > &u, UpdateType type, int iter=-1)
Definition ROL_Constraint_SimOpt.hpp:144
ROL::Constraint_SimOpt::checkApplyJacobian_2
std::vector< std::vector< Real > > checkApplyJacobian_2(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1324
ROL::Constraint_SimOpt::DEFAULT_print_
const bool DEFAULT_print_
Definition ROL_Constraint_SimOpt.hpp:87
ROL::Constraint_SimOpt::update
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable i...
Definition ROL_Constraint_SimOpt.hpp:869
ROL::Constraint_SimOpt::applyAdjointHessian
virtual void applyAdjointHessian(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ,...
Definition ROL_Constraint_SimOpt.hpp:922
ROL::Constraint_SimOpt::checkApplyAdjointHessian_12
std::vector< std::vector< Real > > checkApplyAdjointHessian_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1699
ROL::Constraint_SimOpt::stol_
Real stol_
Definition ROL_Constraint_SimOpt.hpp:96
ROL::Constraint_SimOpt::update
virtual void update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-1)
Definition ROL_Constraint_SimOpt.hpp:133
ROL::Constraint_SimOpt::jv_
Ptr< Vector< Real > > jv_
Definition ROL_Constraint_SimOpt.hpp:78
ROL::Constraint_SimOpt::checkApplyAdjointHessian_22
std::vector< std::vector< Real > > checkApplyAdjointHessian_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1801
ROL::Constraint_SimOpt::solve
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
Definition ROL_Constraint_SimOpt.hpp:191
ROL::Constraint_SimOpt::checkApplyAdjointHessian_12
std::vector< std::vector< Real > > checkApplyAdjointHessian_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
, , , ,
Definition ROL_Constraint_SimOpt.hpp:1680
ROL::Constraint_SimOpt::applyJacobian
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
Definition ROL_Constraint_SimOpt.hpp:889
ROL::Constraint_SimOpt::DEFAULT_maxit_
const int DEFAULT_maxit_
Definition ROL_Constraint_SimOpt.hpp:86
ROL::Constraint_SimOpt::factor_
Real factor_
Definition ROL_Constraint_SimOpt.hpp:97
ROL::Constraint_SimOpt::applyPreconditioner
virtual void applyPreconditioner(Vector< Real > &pv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g, Real &tol)
Apply a constraint preconditioner at , , to vector . In general, this preconditioner satisfies the fo...
Definition ROL_Constraint_SimOpt.hpp:834
ROL::Constraint_SimOpt::checkApplyJacobian_1
std::vector< std::vector< Real > > checkApplyJacobian_1(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Definition ROL_Constraint_SimOpt.hpp:1199
ROL::Constraint_SimOpt::solve_update
virtual void solve_update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-1)
Update SimOpt constraint during solve (disconnected from optimization updates).
Definition ROL_Constraint_SimOpt.hpp:160
ROL::Constraint_SimOpt::DEFAULT_stol_
const Real DEFAULT_stol_
Definition ROL_Constraint_SimOpt.hpp:83
ROL::Constraint_SimOpt::checkSolve
virtual Real checkSolve(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &c, const bool printToStream=true, std::ostream &outStream=std::cout)
Definition ROL_Constraint_SimOpt.hpp:951
ROL::Constraint_SimOpt::Constraint_SimOpt
Constraint_SimOpt()
Definition ROL_Constraint_SimOpt.hpp:108
ROL::Constraint_SimOpt::value
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
Definition ROL_Constraint_SimOpt.hpp:880
ROL::Constraint
Defines the general constraint operator interface.
Definition ROL_Constraint.hpp:52
ROL::Constraint::applyPreconditioner
virtual void applyPreconditioner(Vector< Real > &pv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g, Real &tol)
Apply a constraint preconditioner at , , to vector . Ideally, this preconditioner satisfies the follo...
Definition ROL_Constraint.hpp:249
ROL::Constraint::solveAugmentedSystem
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system
Definition ROL_ConstraintDef.hpp:161
ROL::Vector_SimOpt
Defines the linear algebra or vector space interface for simulation-based optimization.
Definition ROL_Vector_SimOpt.hpp:23
ROL::Vector_SimOpt::get_2
ROL::Ptr< const Vector< Real > > get_2() const
Definition ROL_Vector_SimOpt.hpp:144
ROL::Vector_SimOpt::get_1
ROL::Ptr< const Vector< Real > > get_1() const
Definition ROL_Vector_SimOpt.hpp:140
ROL::Vector_SimOpt::set_1
void set_1(const Vector< Real > &vec)
Definition ROL_Vector_SimOpt.hpp:156
ROL::Vector_SimOpt::set_2
void set_2(const Vector< Real > &vec)
Definition ROL_Vector_SimOpt.hpp:160
ROL::Vector_SimOpt::dual
const Vector< Real > & dual(void) const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
Definition ROL_Vector_SimOpt.hpp:73
ROL::Vector
Defines the linear algebra or vector space interface.
Definition ROL_Vector.hpp:50
ROL::Vector::apply
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition ROL_Vector.hpp:204
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition ROL_Vector.hpp:175
ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Vector::dual
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
Definition ROL_Vector.hpp:192
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition ROL_Vector.hpp:133
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Vector::dimension
virtual int dimension() const
Return dimension of the vector space.
Definition ROL_Vector.hpp:162
ROL::Vector::basis
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition ROL_Vector.hpp:148
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition ROL_Vector.hpp:119
ROL::Finite_Difference_Arrays::shifts
const int shifts[4][4]
Definition ROL_Types.hpp:831
ROL::Finite_Difference_Arrays::weights
const double weights[4][5]
Definition ROL_Types.hpp:838
ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition ROL_Types.hpp:57
Generated by doxygen 1.9.8