44 #ifndef ROL_EQUALITY_CONSTRAINT_SIMOPT_H 45 #define ROL_EQUALITY_CONSTRAINT_SIMOPT_H 98 Teuchos::RCP<Vector<Real> >
unew_;
99 Teuchos::RCP<Vector<Real> >
jv_;
123 unew_(Teuchos::null), jv_(Teuchos::null),
126 DEFAULT_factor_(0.5),
127 DEFAULT_decr_(1.e-4),
129 DEFAULT_print_(false),
130 rtol_(DEFAULT_rtol_), stol_(DEFAULT_stol_), factor_(DEFAULT_factor_),
131 decr_(DEFAULT_decr_), maxit_(DEFAULT_maxit_), print_(DEFAULT_print_),
182 Real cnorm = c.
norm(), alpha = 1.0, tmp = 0.0;
185 std::cout <<
"\n Default EqualityConstraint_SimOpt::solve\n";
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";
192 while ( cnorm > rtol_ && cnt < maxit_) {
196 unew_->axpy(-alpha, *jv_);
198 value(c,*unew_,z,tol);
201 while ( tmp > (1.0-decr_*alpha)*cnorm &&
205 unew_->axpy(-alpha,*jv_);
207 value(c,*unew_,z,tol);
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;
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_);
279 h = std::max(1.0,u.
norm()/v.
norm())*tol;
282 Teuchos::RCP<Vector<Real> > unew = u.
clone();
287 value(jv,*unew,z,ctol);
289 Teuchos::RCP<Vector<Real> > cold = jv.
clone();
291 value(*cold,u,z,ctol);
322 h = std::max(1.0,u.
norm()/v.
norm())*tol;
325 Teuchos::RCP<Vector<Real> > znew = z.
clone();
330 value(jv,u,*znew,ctol);
332 Teuchos::RCP<Vector<Real> > cold = jv.
clone();
334 value(*cold,u,z,ctol);
360 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
361 "The method applyInverseJacobian_1 is used but not implemented!\n");
414 h = std::max(1.0,u.
norm()/v.
norm())*tol;
416 Teuchos::RCP<Vector<Real> > cold = dualv.
clone();
417 Teuchos::RCP<Vector<Real> > cnew = dualv.
clone();
419 value(*cold,u,z,ctol);
420 Teuchos::RCP<Vector<Real> > unew = u.
clone();
422 for (
int i = 0; i < u.
dimension(); i++) {
424 unew->axpy(h,*(u.
basis(i)));
426 value(*cnew,*unew,z,ctol);
427 cnew->axpy(-1.0,*cold);
429 ajv.
axpy(cnew->dot(v),*((u.
dual()).basis(i)));
485 h = std::max(1.0,u.
norm()/v.
norm())*tol;
487 Teuchos::RCP<Vector<Real> > cold = dualv.
clone();
488 Teuchos::RCP<Vector<Real> > cnew = dualv.
clone();
490 value(*cold,u,z,ctol);
491 Teuchos::RCP<Vector<Real> > znew = z.
clone();
493 for (
int i = 0; i < z.
dimension(); i++) {
495 znew->axpy(h,*(z.
basis(i)));
497 value(*cnew,u,*znew,ctol);
498 cnew->axpy(-1.0,*cold);
500 ajv.
axpy(cnew->dot(v),*((z.
dual()).basis(i)));
525 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
526 "The method applyInverseAdjointJacobian_1 is used but not implemented!\n");
555 h = std::max(1.0,u.
norm()/v.
norm())*tol;
558 Teuchos::RCP<Vector<Real> > unew = u.
clone();
564 Teuchos::RCP<Vector<Real> > jv = ahwv.
clone();
599 h = std::max(1.0,u.
norm()/v.
norm())*tol;
602 Teuchos::RCP<Vector<Real> > unew = u.
clone();
608 Teuchos::RCP<Vector<Real> > jv = ahwv.
clone();
643 h = std::max(1.0,u.
norm()/v.
norm())*tol;
646 Teuchos::RCP<Vector<Real> > znew = z.
clone();
652 Teuchos::RCP<Vector<Real> > jv = ahwv.
clone();
686 h = std::max(1.0,u.
norm()/v.
norm())*tol;
689 Teuchos::RCP<Vector<Real> > znew = z.
clone();
695 Teuchos::RCP<Vector<Real> > jv = ahwv.
clone();
775 Teuchos::RCP<ROL::Vector<Real> > ijv = (xs.
get_1())->clone();
780 catch (
const std::logic_error &e) {
786 Teuchos::RCP<ROL::Vector<Real> > ijv_dual = (gs.
get_1())->clone();
787 ijv_dual->
set(ijv->dual());
792 catch (
const std::logic_error &e) {
832 const Vector_SimOpt<Real> &vs = Teuchos::dyn_cast<
const Vector_SimOpt<Real> >(
835 Teuchos::RCP<Vector<Real> > jv2 = jv.
clone();
847 const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<
const Vector_SimOpt<Real> >(
849 Teuchos::RCP<Vector<Real> > ajv1 = (ajvs.
get_1())->clone();
852 Teuchos::RCP<Vector<Real> > ajv2 = (ajvs.
get_2())->clone();
865 const Vector_SimOpt<Real> &xs = Teuchos::dyn_cast<
const Vector_SimOpt<Real> >(
867 const Vector_SimOpt<Real> &vs = Teuchos::dyn_cast<
const Vector_SimOpt<Real> >(
870 Teuchos::RCP<Vector<Real> > C11 = (ahwvs.
get_1())->clone();
871 Teuchos::RCP<Vector<Real> > C21 = (ahwvs.
get_1())->clone();
877 Teuchos::RCP<Vector<Real> > C12 = (ahwvs.
get_2())->clone();
878 Teuchos::RCP<Vector<Real> > C22 = (ahwvs.
get_2())->clone();
890 const bool printToStream =
true,
891 std::ostream & outStream = std::cout) {
894 Teuchos::RCP<ROL::Vector<Real> > r = c.
clone();
895 Teuchos::RCP<ROL::Vector<Real> > s = u.
clone();
898 Teuchos::RCP<ROL::Vector<Real> > cs = c.
clone();
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();
932 const bool printToStream =
true,
933 std::ostream & outStream = std::cout) {
959 const bool printToStream =
true,
960 std::ostream & outStream = std::cout) {
962 Teuchos::RCP<Vector<Real> > Jv = dualw.
clone();
965 Real wJv = w.
dot(Jv->dual());
966 Teuchos::RCP<Vector<Real> > Jw = dualv.
clone();
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>| = " 976 hist <<
" |<w,Jv>| = " << std::abs(wJv) <<
"\n";
977 hist <<
" Relative Error = " << diff / (std::abs(wJv)+
ROL_UNDERFLOW) <<
"\n";
978 outStream << hist.str();
1000 const bool printToStream =
true,
1001 std::ostream & outStream = std::cout) {
1026 const bool printToStream =
true,
1027 std::ostream & outStream = std::cout) {
1029 Teuchos::RCP<Vector<Real> > Jv = dualw.
clone();
1032 Real wJv = w.
dot(Jv->dual());
1033 Teuchos::RCP<Vector<Real> > Jw = dualv.
clone();
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>| = " 1043 hist <<
" |<w,Jv>| = " << std::abs(wJv) <<
"\n";
1044 hist <<
" Relative Error = " << diff / (std::abs(wJv)+
ROL_UNDERFLOW) <<
"\n";
1045 outStream << hist.str();
1054 const bool printToStream =
true,
1055 std::ostream & outStream = std::cout) {
1057 Teuchos::RCP<Vector<Real> > Jv = jv.
clone();
1060 Teuchos::RCP<Vector<Real> > iJJv = u.
clone();
1063 Teuchos::RCP<Vector<Real> > diff = v.
clone();
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|| = " 1073 hist <<
" ||v|| = " << vnorm <<
"\n";
1074 hist <<
" Relative Error = " << dnorm / (vnorm+
ROL_UNDERFLOW) <<
"\n";
1075 outStream << hist.str();
1084 const bool printToStream =
true,
1085 std::ostream & outStream = std::cout) {
1087 Teuchos::RCP<Vector<Real> > Jv = jv.
clone();
1090 Teuchos::RCP<Vector<Real> > iJJv = v.
clone();
1093 Teuchos::RCP<Vector<Real> > diff = v.
clone();
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|| = " 1103 hist <<
" ||v|| = " << vnorm <<
"\n";
1104 hist <<
" Relative Error = " << dnorm / (vnorm+
ROL_UNDERFLOW) <<
"\n";
1105 outStream << hist.str();
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...
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.
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 .
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.
const Real DEFAULT_factor_
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.
Defines the linear algebra or vector space interface.
Teuchos::RCP< Vector< Real > > jv_
const bool DEFAULT_print_
Defines the equality constraint operator interface for simulation-based optimization.
virtual Real dot(const Vector &x) const =0
Compute where .
Teuchos::RCP< Vector< Real > > unew_
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 ...
EqualityConstraint_SimOpt()
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 .
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.
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.
static const double ROL_EPSILON
Platform-dependent machine epsilon.