45 #ifndef ROL_PRIMALDUALINTERIORPOINTSTEP_H
46 #define ROL_PRIMALDUALINTERIORPOINTSTEP_H
48 #include "ROL_VectorsNorms.hpp"
113 namespace InteriorPoint {
114 template <
class Real>
115 class PrimalDualInteriorPointStep :
public Step<Real> {
117 typedef Vector<Real>
V;
118 typedef PartitionedVector<Real>
PV;
119 typedef Objective<Real> OBJ;
120 typedef BoundConstraint<Real> BND;
121 typedef Krylov<Real> KRYLOV;
122 typedef LinearOperator<Real> LINOP;
123 typedef LinearOperatorFromConstraint<Real> LOPEC;
124 typedef Constraint<Real> EQCON;
125 typedef StepState<Real> STATE;
126 typedef InteriorPointPenalty<Real> PENALTY;
127 typedef PrimalDualInteriorPointResidual<Real> RESIDUAL;
131 ROL::Ptr<KRYLOV> krylov_;
132 ROL::Ptr<LINOP> precond_;
146 ROL::Ptr<V> xscratch_;
147 ROL::Ptr<V> lscratch_;
149 ROL::Ptr<V> zlscratch_;
150 ROL::Ptr<V> zuscratch_;
160 Elementwise::Multiply<Real> mult_;
188 void updateState(
const V& x,
const V &l, OBJ &obj,
189 EQCON &con, BND &bnd, ALGO &algo_state ) {
191 Real tol = std::sqrt(ROL_EPSILON<Real>());
194 obj.update(x,
true,algo_state.iter);
195 con.update(x,
true,algo_state.iter);
197 algo_state.value = obj.value(tol);
198 con.value(*(state->constraintVec),x,tol);
200 obj.gradient(*(state->gradientVec),x,tol);
201 con.applyAdjointJacobian(*g_,l,x,tol);
203 state->gradientVec->plus(*g_);
205 state->gradientVec->axpy(-1.0,*zl_);
206 state->gradientVec->axpy(-1.0,*zu_);
209 algo_state.gnorm =
normLinf(*(state->gradientVec))/sd_;
212 algo_state.cnorm =
normLinf(*(state->constraintVec));
214 Elementwise::Multiply<Real> mult;
221 xscratch_->applyBinary(mult,*zl_);
223 exactLowerViolation =
normLinf(*xscratch_)/scl_;
226 xscratch_->applyBunary(mult,*zu_);
228 exactUpperBound =
normLinf(*xscratch_)/scu_;
231 algo_state.aggregateModelError = std::max(exactLowerViolation,
232 exactUpperViolation);
242 void updateScalingParameters(
void) {
248 sd_ = (nl+nzl+nzu)/(diml_+2*dimx);
249 sd_ = std::max(smax_,sd_)
254 scl_ = std::max(smax_,scl_);
258 scu_ = std::max(smax_,scu_);
264 using Step<Real>::initialize;
265 using Step<Real>::compute;
268 PrimalDualInteriorPointStep( ROL::ParameterList &parlist,
269 const ROL::Ptr<Krylov<Real> > &krylov = ROL::nullPtr,
270 const ROL::Ptr<LinearOperator<Real> > &precond = ROL::nullPtr ) :
271 Step<Real>(), krylov_(krylov), precond_(precond), iterKrylov_(0), flagKrylov_(0) {
273 typedef ROL::ParameterList PL;
275 PL &iplist = parlist.sublist(
"Step").sublist(
"Primal Dual Interior Point");
277 kappa1_ = iplist.get(
"Bound Perturbation Coefficient 1", 1.e-2);
278 kappa2_ = iplist.get(
"Bound Perturbation Coefficient 2", 1.e-2);
279 lambda_max_ = iplist.get(
" Multiplier Maximum Value", 1.e3 );
280 smax_ = iplist.get(
"Maximum Scaling Parameter", 1.e2 );
281 tau_min_ = iplist.get(
"Minimum Fraction-to-Boundary Parameter", 0.99 );
282 kappa_mu_ = iplist.get(
"Multiplicative Penalty Reduction Factor", 0.2 );
283 theta_mu_ = iplist.get(
"Penalty Update Power", 1.5 );
284 eps_tol_ = iplist.get(
"Error Tolerance", 1.e-8);
285 symmetrize_ = iplist.get(
"Symmetrize Primal Dual System",
true );
287 PL &filter = iplist.sublist(
"Filter Parameters");
289 if(krylov_ == ROL::nullPtr) {
290 krylov_ = KrylovFactory<Real>(parlist);
293 if( precond_ == ROL::nullPtr) {
296 apply(
V& Hv,
const V &v, Real tol )
const {
301 precond_ = ROL::makePtr<IdentityOperator<Real>>();
308 ~PrimalDualInteriorPointStep() {
311 void initialize( Vector<Real> &x,
const Vector<Real> &g, Vector<Real> &l,
const Vector<Real> &c,
312 Objective<Real> &obj, Constraint<Real> &con, BoundConstraint<Real> &bnd,
313 AlgorithmState<Real> &algo_state ) {
316 using Elementwise::ValueSet;
318 ROL::Ptr<PENALTY> &ipPen =
dynamic_cast<PENALTY&
>(obj);
322 state->descentVec = x.clone();
323 state->gradientVec = g.clone();
324 state->constaintVec = c.clone();
326 diml_ = l.dimension();
327 dimx_ = x.dimension();
331 Real tol = std::sqrt(ROL_EPSILON<Real>());
338 xscratch_ = x.clone();
339 lscratch_ = l.clone();
341 zlscratch_ = x.clone();
342 zuscratch_ = x.clone();
352 xl_ = bnd.getLowerBound();
353 xu_ = bnd.getUpperBound();
355 maskl_ = ipPen.getLowerMask();
356 masku_ = ipPen.getUpperMask();
366 ROL::Ptr<V> xdiff = xu_->clone();
368 xdiff->axpy(-1.0,*xl_);
369 xdiff->scale(kappa2_);
371 class Max1X :
public Elementwise::UnaryFunction<Real> {
373 Real
apply(
const Real &x )
const {
374 return std::max(1.0,x);
379 Elementwise::AbsoluteValue<Real> absval;
380 Elementwise::Min min;
383 ROL::Ptr<V> pl = xl_->clone();
384 pl->applyUnary(absval);
385 pl->applyUnary(max1x);
387 pl->applyBinary(min,xdiff);
390 ROL::Ptr<V> pu = xu_->clone();
391 pu->applyUnary(absval);
392 pu->applyUnary(max1x);
394 pu->applyBinary(min,xdiff);
401 pbnd_ = ROL::makePtr<BoundConstraint<Real>>(pl,pu)
418 con.solveAugmentedSystem(*xscratch_,*l_,*g_,*lscratch_,x,tol);
423 if( normInf(l_) > lambda_max_ ) {
428 algo_state.nfval = 0;
429 algo_state.ncval = 0;
430 algo_state.ngrad = 0;
432 updateState(x,l,obj,con,bnd,algo_state);
437 void compute( Vector<Real> &s,
const Vector<Real> &x,
const Vector<Real> &l,
438 Objective<Real> &obj, Constraint<Real> &con,
439 BoundConstraint<Real> &bnd, AlgorithmState<Real> &algo_state ) {
444 Elementwise::Fill<Real> minus_mu(-mu_);
445 Elementwise::Divide<Real> div;
446 Elementwise::Multiply<Real> mult;
450 ROL::Ptr<OBJ> obj_ptr = ROL::makePtrFromRef(obj);
451 ROL::Ptr<EQCON> con_ptr = ROL::makePtrFromRef(con);
452 ROL::Ptr<BND> bnd_ptr = ROL::makePtrFromRef(bnd);
461 state->constraintVec,
468 ROL::Ptr<RESIDUAL> residual = ROL::makePtr<RESIDUAL>(obj,con,bnd,*sol,maskL_,maskU_,w_,mu_,symmetrize_);
470 residual->value(*rhs,*sysvec,tol);
474 LOPEC jacobian( sysvec, residual );
477 krylov_->run(*sol,jacobian,*residual,*precond_,iterKrylov_,flagKrylov_);
489 void update( Vector<Real> &x, Vector<Real> &l,
const Vector<Real> &s,
490 Objective<Real> &obj, Constraint<Real> &con,
491 BoundConstraint<Real> &bnd, AlgorithmState<Real> &algo_state ) {
494 Elementwise::Shift<Real> minus_mu(-mu_);
497 xscratch_->applyBinary(mult,*zl_);
498 xscratch_->applyUnary(minus_mu);
500 lowerViolation =
normLinf(*xscratch_)/scl_;
503 xscratch_->applyBinary(mult,*zu_);
504 xscratch_->applyUnary(minus_mu);
506 upperBound =
normLinf(*xscratch_)/scu_;
509 Real Emu = algo_state.gnorm;
510 Emu = std::max(Emu,algo_state.cnorm);
511 Emu = std::max(Emu,upperBound);
512 Emu = std::max(Emu,lowerBound);
515 if(Emu < (kappa_epsilon_*mu_) ) {
521 mu_ = std::min(kappa_mu_*mu_old,std::pow(mu_old,theta_mu_));
522 mu_ = std::max(eps_tol_/10.0,mu_);
525 tau_ = std::max(tau_min_,1.0-mu_);
537 std::string printHeader(
void )
const {
538 std::string head(
"");
543 std::string printName(
void )
const {
544 std::string name(
"");
555 #endif // ROL_PRIMALDUALINTERIORPOINTSTEP_H
PartitionedVector< Real > PV
Real normLinf(const Vector< Real > &x)
ROL::Ptr< Vector< Real > > CreatePartitionedVector(const ROL::Ptr< Vector< Real >> &a)
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1) override
Defines the linear algebra or vector space interface.
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Real normL1(const Vector< Real > &x)
ROL::Ptr< StepState< Real > > getState(void)
ROL::DiagonalOperator apply
virtual void set(const Vector &x)
Set where .