10 #ifndef ROL_BOUNDFLETCHER_H
11 #define ROL_BOUNDFLETCHER_H
18 #include "ROL_Ptr.hpp"
31 Ptr<const Vector<Real> >
low_;
32 Ptr<const Vector<Real> >
upp_;
64 Ptr<Vector<Real> >
Q_;
66 Ptr<Vector<Real> >
DQ_;
83 Ptr<Vector<Real> >
Qv_;
86 Ptr<Vector<Real> >
Tv_;
87 Ptr<Vector<Real> >
w_;
88 Ptr<Vector<Real> >
v_;
117 class DiffLower :
public Elementwise::BinaryFunction<Real> {
120 Real
apply(
const Real& x,
const Real& y)
const {
121 const Real NINF(ROL_NINF<Real>());
122 return (y <= NINF ? static_cast<Real>(-1.) : x - y);
126 class DiffUpper :
public Elementwise::BinaryFunction<Real> {
129 Real
apply(
const Real& x,
const Real& y)
const {
130 const Real INF(ROL_INF<Real>());
131 return (y >= INF ? static_cast<Real>(-1.) : y - x);
135 class FormQ :
public Elementwise::BinaryFunction<Real> {
138 Real
apply(
const Real& x,
const Real& y)
const {
140 if( x < zero && y < zero) {
141 return static_cast<Real
>(1);
143 if( x < zero && y >= zero ) {
146 if( x >= zero && y < zero ) {
149 return std::min(x, y);
153 class FormDQ :
public Elementwise::BinaryFunction<Real> {
156 Real
apply(
const Real& x,
const Real& y)
const {
157 Real
zero(0.), one(1.), mone(-1.);
161 if( x < zero && y >=
zero ) {
179 const Ptr<Constraint<Real> >
con_;
180 const Ptr<const Vector<Real> >
x_;
182 const Ptr<Vector<Real> >
Qv_;
195 con_->applyAdjointJacobian(*(Hvp.
get(0)), *(vp.
get(1)), *
x_, tol);
196 Hvp.
get(0)->applyBinary(Elementwise::Multiply<Real>(), *
Qsqrt_);
197 Hvp.
get(0)->plus(*(vp.
get(0)));
200 Qv_->applyBinary(Elementwise::Multiply<Real>(), *
Qsqrt_);
208 const Ptr<Constraint<Real> >
con_;
209 const Ptr<const Vector<Real> >
x_;
210 const Ptr<Vector<Real> >
Q_;
211 const Ptr<Vector<Real> >
Qv_;
224 con_->applyAdjointJacobian(*(Hvp.
get(0)), *(vp.
get(1)), *
x_, tol);
225 Hvp.
get(0)->plus(*(vp.
get(0)));
228 Qv_->applyBinary(Elementwise::Multiply<Real>(), *
Q_);
236 const Ptr<Constraint<Real> >
con_;
237 const Ptr<const Vector<Real> >
x_;
252 con_->applyPreconditioner(*(Hvp.
get(1)),*(vp.
get(1)),*
x_,*
x_, zero);
262 ROL::ParameterList &parlist)
265 low_ = bnd->getLowerBound();
266 upp_ = bnd->getUpperBound();
269 y_ = conVec.
dual().clone();
270 g_ = optVec.
dual().clone();
286 w_ = optVec.
dual().clone();
287 v_ = conVec.
dual().clone();
298 vv_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({
v1_,
v2_}));
302 ww_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({
w1_,
w2_}));
306 bb_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({
b1_,
b2_}));
308 ROL::ParameterList& sublist = parlist.sublist(
"Step").sublist(
"Fletcher");
309 HessianApprox_ = sublist.get(
"Level of Hessian Approximation", 0);
311 AugSolve_ = sublist.get(
"Type of Augmented System Solve", 0);
317 delta_ = sublist.get(
"Regularization Parameter", 0.0);
319 useInexact_ = sublist.get(
"Inexact Solves",
false);
321 ROL::ParameterList krylovList;
322 Real atol =
static_cast<Real
>(1e-12);
323 Real rtol =
static_cast<Real
>(1e-2);
324 krylovList.sublist(
"General").sublist(
"Krylov").set(
"Type",
"GMRES");
325 krylovList.sublist(
"General").sublist(
"Krylov").set(
"Absolute Tolerance", atol);
326 krylovList.sublist(
"General").sublist(
"Krylov").set(
"Relative Tolerance", rtol);
327 krylovList.sublist(
"General").sublist(
"Krylov").set(
"Iteration Limit", 200);
328 krylov_ = KrylovFactory<Real>(krylovList);
332 obj_->update(x,flag,iter);
333 con_->update(x,flag,iter);
398 w_->applyBinary(Elementwise::Multiply<Real>(), *
Qsqrt_);
399 con_->applyAdjointHessian( *
gPhi_, *
y_, *
w_, x, tol2 ); tol2 = origTol;
400 obj_->hessVec( *
Tv_, *
w_, x, tol2 ); tol2 = origTol;
401 gPhi_->axpy( static_cast<Real>(-1), *
Tv_ );
403 con_->applyAdjointJacobian( *
Tv_, *
v_, x, tol2); tol2 = origTol;
406 Tv_->applyBinary(Elementwise::Multiply<Real>(), *
DQgL_);
409 con_->applyAdjointHessian( *
Tv_, *
v_, *
QgL_, x, tol2 ); tol2 = origTol;
423 Tv_->applyBinary( Elementwise::Multiply<Real>(), *
DQgL_ );
424 gPhi_->axpy(static_cast<Real>(-1), *
Tv_);
426 w_->applyBinary( Elementwise::Multiply<Real>(), *
Q_ );
427 obj_->hessVec( *
Tv_, *
w_, x, tol2); tol2 = origTol;
428 gPhi_->axpy( static_cast<Real>(-1), *
Tv_ );
429 con_->applyAdjointHessian( *
Tv_, *
y_, *
w_, x, tol2 ); tol2 = origTol;
432 con_->applyAdjointHessian( *
Tv_, *
v_, *
QgL_, x, tol2 ); tol2 = origTol;
442 con_->applyAdjointJacobian( *
Tv_, *
c_, x, tol2 ); tol2 = origTol;
458 value(x, tol2); tol2 = origTol;
465 obj_->hessVec( hv, v, x, tol2 ); tol2 = origTol;
466 con_->applyAdjointHessian( *
Tv_, *
y_, v, x, tol2 ); tol2 = origTol;
467 hv.
axpy(static_cast<Real>(-1), *
Tv_ );
471 htmp1_->applyBinary(Elementwise::Multiply<Real>(), *
Qsqrt_);
472 htmp1_->scale(static_cast<Real>(-1));
475 Tv_->applyBinary( Elementwise::Multiply<Real>(), v );
477 con_->applyJacobian( *
htmp2_, *
Tv_, x, tol2); tol2 = origTol;
480 con_->applyAdjointJacobian( *
Tv_, *
v_, x, tol2 ); tol2 = origTol;
483 con_->applyJacobian( *
htmp2_, v, x, tol2 ); tol2 = origTol;
485 con_->applyAdjointJacobian( *
Tv_, *
v_, x, tol2 ); tol2 = origTol;
487 Tv_->applyBinary( Elementwise::Multiply<Real>(), *
DQgL_ );
490 w_->applyBinary( Elementwise::Multiply<Real>(), *
Qsqrt_ );
491 obj_->hessVec( *
Tv_, *
w_, x, tol2 ); tol2 = origTol;
492 hv.
axpy( static_cast<Real>(-1), *
Tv_);
493 con_->applyAdjointHessian( *
Tv_, *
y_, *
w_, x, tol2 ); tol2 = origTol;
499 obj_->hessVec( hv, v, x, tol2 ); tol2 = origTol;
500 con_->applyAdjointHessian( *
Tv_, *
y_, v, x, tol2 ); tol2 = origTol;
501 hv.
axpy(static_cast<Real>(-1), *
Tv_ );
505 Tv_->applyBinary( Elementwise::Multiply<Real>(), *
DQgL_ );
507 Tv_->scale( static_cast<Real>(-1) );
508 con_->applyJacobian( *
htmp2_, *
Tv_, x, tol2 ); tol2 = origTol;
512 con_->applyJacobian( *
htmp2_, v, x, tol2 ); tol2 = origTol;
516 Tv_->applyBinary( Elementwise::Multiply<Real>(), *
DQgL_ );
517 hv.
axpy( static_cast<Real>(-1), *
Tv_ );
519 w_->applyBinary( Elementwise::Multiply<Real>(), *
Q_ );
520 obj_->hessVec( *
Tv_, *
w_, x, tol2 ); tol2 = tol;
521 hv.
axpy( static_cast<Real>(-1), *
Tv_);
522 con_->applyAdjointHessian( *
Tv_, *
y_, *
w_, x, tol2 ); tol2 = origTol;
529 con_->applyJacobian( *
b2_, v, x, tol2 ); tol2 = origTol;
530 con_->applyAdjointJacobian( *
Tv_, *
b2_, x, tol2 ); tol2 = origTol;
532 con_->applyAdjointHessian( *
Tv_, *
c_, v, x, tol2); tol2 = origTol;
544 bool refine =
false) {
546 ROL::Ptr<LinearOperator<Real> > K;
553 K = ROL::makePtr<AugSystemNonSym>(
con_, makePtrFromRef(x),
Q_,
Qv_,
delta_);
557 ROL::Ptr<LinearOperator<Real> > P
558 = ROL::makePtr<AugSystemPrecond>(
con_, makePtrFromRef(x));
568 K->apply(*
vv_, *
ww_, tol); tol = origTol;
570 b1_->axpy( static_cast<Real>(-1), *
v1_ );
571 b2_->axpy( static_cast<Real>(-1), *
v2_ );
579 krylov_->resetAbsoluteTolerance(tol);
613 Qsg_->applyBinary(Elementwise::Multiply<Real>(), *
Qsqrt_);
619 gL_->applyBinary(Elementwise::Divide<Real>(), *
Qsqrt_);
621 QgL_->applyBinary(Elementwise::Multiply<Real>(), *
Qsqrt_);
628 QgL_->applyBinary(Elementwise::Multiply<Real>(), *
Q_);
634 DQgL_->applyBinary(Elementwise::Multiply<Real>(), *
DQ_);
650 Qsqrt_->applyUnary(Elementwise::SquareRoot<Real>());
Provides the interface to evaluate objective functions.
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Ptr< Vector< Real > > w2_
Ptr< Vector< Real > > gPhi_
Ptr< Vector< Real > > Qsg_
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
const Ptr< Constraint< Real > > con_
Real quadPenaltyParameter_
ROL::Ptr< const Vector< Real > > get(size_type i) const
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Defines the linear algebra of vector space on a generic partitioned vector.
AugSystemSym(const Ptr< Constraint< Real > > &con, const Ptr< const Vector< Real > > &x, const Ptr< Vector< Real > > &Qsqrt, const Ptr< Vector< Real > > &Qv, const Real delta)
void objValue(const Vector< Real > &x, Real &tol)
const Ptr< Objective< Real > > obj_
Contains definitions of custom data types in ROL.
const Ptr< Constraint< Real > > con_
Ptr< Vector< Real > > czeros_
Real apply(const Real &x, const Real &y) const
Ptr< Vector< Real > > w1_
Real apply(const Real &x, const Real &y) const
Defines the linear algebra or vector space interface.
Ptr< Krylov< Real > > krylov_
const Ptr< Vector< Real > > Qv_
const Ptr< Vector< Real > > Q_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
const Ptr< const Vector< Real > > x_
BoundFletcher(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< Constraint< Real > > &con, const ROL::Ptr< BoundConstraint< Real > > &bnd, const Vector< Real > &optVec, const Vector< Real > &conVec, ROL::ParameterList &parlist)
Ptr< Vector< Real > > QgL_
Ptr< Vector< Real > > gL_
Ptr< Vector< Real > > DQgL_
Ptr< const Vector< Real > > low_
Ptr< Vector< Real > > Tv_
Ptr< Vector< Real > > v1_
Ptr< Vector< Real > > b1_
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Ptr< Vector< Real > > Qv_
Real value(const Vector< Real > &x, Real &tol)
Compute value.
bool isMultiplierComputed_
void computeMultipliers(const Vector< Real > &x, const Real tol)
const Ptr< const Vector< Real > > x_
Ptr< Vector< Real > > scaledc_
const Ptr< Constraint< Real > > con_
void conValue(const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > Qsqrt_
const Ptr< Constraint< Real > > con_
void computeDQ(const Vector< Real > &x)
void set(const V &x)
Set where .
const Ptr< Vector< Real > > Qsqrt_
Ptr< PartitionedVector< Real > > vv_
AugSystemPrecond(const Ptr< Constraint< Real > > con, const Ptr< const Vector< Real > > x)
Provides the interface to apply a linear operator.
Provides the interface to apply upper and lower bound constraints.
Ptr< const Vector< Real > > upp_
AugSystemNonSym(const Ptr< Constraint< Real > > &con, const Ptr< const Vector< Real > > &x, const Ptr< Vector< Real > > &Q, const Ptr< Vector< Real > > &Qv, const Real delta)
Ptr< Vector< Real > > umx_
Ptr< Vector< Real > > DQ_
const Ptr< Vector< Real > > Qv_
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
virtual void set(const Vector &x)
Set where .
Ptr< Vector< Real > > v2_
Ptr< Vector< Real > > QsgL_
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Ptr< Vector< Real > > b2_
Ptr< PartitionedVector< Real > > bb_
Ptr< PartitionedVector< Real > > ww_
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Ptr< Vector< Real > > htmp1_
const Ptr< const Vector< Real > > x_
void computeQ(const Vector< Real > &x)
Ptr< Vector< Real > > htmp2_
void objGrad(const Vector< Real > &x, Real &tol)
Defines the general constraint operator interface.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
void solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol, bool refine=false)
Ptr< Vector< Real > > xzeros_