45 #ifndef ROL_SEMISMOOTHNEWTONPROJECTION_DEF_H 46 #define ROL_SEMISMOOTHNEWTONPROJECTION_DEF_H 50 template<
typename Real>
62 DEFAULT_factor_ (0.5),
63 DEFAULT_regscale_ (1e-4),
64 DEFAULT_errscale_ (1e-2),
65 DEFAULT_maxit_ (5000),
67 DEFAULT_verbosity_ (0),
68 DEFAULT_useproj_ (false),
69 atol_ (DEFAULT_atol_),
70 rtol_ (DEFAULT_rtol_),
71 stol_ (DEFAULT_stol_),
72 decr_ (DEFAULT_decr_),
73 factor_ (DEFAULT_factor_),
74 regscale_ (DEFAULT_regscale_),
75 errscale_ (DEFAULT_errscale_),
76 maxit_ (DEFAULT_maxit_),
77 lstype_ (DEFAULT_lstype_),
78 verbosity_ (DEFAULT_verbosity_),
79 useproj_ (DEFAULT_useproj_) {
86 list.sublist(
"General").sublist(
"Krylov").set(
"Type",
"CG");
87 list.sublist(
"General").sublist(
"Krylov").set(
"Absolute Tolerance", 1e-6);
88 list.sublist(
"General").sublist(
"Krylov").set(
"Relative Tolerance", 1e-4);
89 list.sublist(
"General").sublist(
"Krylov").set(
"Iteration Limit",
dim_);
90 list.sublist(
"General").set(
"Inexact Hessian-Times-A-Vector",
false);
91 krylov_ = KrylovFactory<Real>(list);
96 Real res0 = std::max(resl,resu);
98 res0 =
static_cast<Real
>(1);
103 template<
typename Real>
112 ParameterList &ppl = list.sublist(
"General").sublist(
"Polyhedral Projection");
116 decr_ = ppl.sublist(
"Semismooth Newton").get(
"Sufficient Decrease Tolerance",
DEFAULT_decr_);
126 klist.sublist(
"General").sublist(
"Krylov") = ppl.sublist(
"Semismooth Newton").sublist(
"Krylov");
127 klist.sublist(
"General").set(
"Inexact Hessian-Times-A-Vector",
false);
128 krylov_ = KrylovFactory<Real>(klist);
131 template<
typename Real>
133 if (con_ == nullPtr) {
137 project_ssn(x, *mul_, *dlam_, stream);
141 template<
typename Real>
143 Real tol(std::sqrt(ROL_EPSILON<Real>()));
145 con_->value(r,y,tol);
149 template<
typename Real>
157 Ptr<Precond> M = makePtr<Precond>();
158 Ptr<Jacobian> J = makePtr<Jacobian>(con_,bnd_,makePtrFromRef(y),xdual_,xprim_,mu);
159 krylov_->run(s,*J,r,*M,iter,flag);
162 template<
typename Real>
166 Real tol(std::sqrt(ROL_EPSILON<Real>()));
169 con_->applyAdjointJacobian(*xdual_,lam,x,tol);
170 y.
plus(xdual_->dual());
174 template<
typename Real>
178 std::ostream &stream)
const {
179 const Real
zero(0), half(0.5), one(1);
181 update_primal(*xnew_,x,lam);
182 Real rnorm = residual(*res_,*xnew_);
187 Real alpha(1), tmp(0), mu(0), rho(1), dd(0);
188 int iter(0), flag(0);
189 std::ios_base::fmtflags streamFlags(stream.flags());
190 if (verbosity_ > 2) {
192 stream << std::scientific << std::setprecision(6);
193 stream <<
" Polyhedral Projection using Dual Semismooth Newton" << std::endl;
195 stream << std::setw(6) << std::left <<
"iter";
196 stream << std::setw(15) << std::left <<
"rnorm";
197 stream << std::setw(15) << std::left <<
"alpha";
198 stream << std::setw(15) << std::left <<
"mu";
199 stream << std::setw(15) << std::left <<
"rho";
200 stream << std::setw(15) << std::left <<
"rtol";
201 stream << std::setw(8) << std::left <<
"kiter";
202 stream << std::setw(8) << std::left <<
"kflag";
205 for (
int cnt = 0; cnt < maxit_; ++cnt) {
207 mu = regscale_*std::max(rnorm,std::sqrt(rnorm));
208 rho = std::min(half,errscale_*std::min(std::sqrt(rnorm),rnorm));
209 solve_newton_system(dlam,*res_,*xnew_,mu,rho,iter,flag);
210 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
211 update_primal(*xnew_,x,*lnew_);
214 tmp = residual(*res_,*xnew_);
215 while ( tmp > (one-decr_*alpha)*rnorm && alpha > stol_ ) {
217 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
218 update_primal(*xnew_,x,*lnew_);
219 tmp = residual(*res_,*xnew_);
224 rnorm = residual(*res_,*xnew_);
226 tmp = dlam.
apply(*res_);
228 while ( tmp < decr_*(one-rho)*mu*dd && alpha > stol_ ) {
230 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
231 update_primal(*xnew_,x,*lnew_);
232 rnorm = residual(*res_,*xnew_);
234 tmp = dlam.
apply(*res_);
241 lam.
axpy(-alpha*tmp/(rnorm*rnorm),res_->dual());
242 update_primal(*xnew_,x,lam);
243 rnorm = residual(*res_,*xnew_);
245 if (verbosity_ > 2) {
247 stream << std::setw(6) << std::left << cnt;
248 stream << std::setw(15) << std::left << rnorm;
249 stream << std::setw(15) << std::left << alpha;
250 stream << std::setw(15) << std::left << mu;
251 stream << std::setw(15) << std::left << rho;
252 stream << std::setw(15) << std::left << ctol_;
253 stream << std::setw(8) << std::left << iter;
254 stream << std::setw(8) << std::left << flag;
257 if (rnorm <= ctol_)
break;
260 if (verbosity_ > 2) {
265 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
266 stream << rnorm <<
" rtol = " << ctol_ << std::endl;
269 stream.flags(streamFlags);
Ptr< Vector< Real > > lnew_
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Ptr< Krylov< Real > > krylov_
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
Ptr< Vector< Real > > res_
SemismoothNewtonProjection(const Vector< Real > &xprim, const Vector< Real > &xdual, const Ptr< BoundConstraint< Real >> &bnd, const Ptr< Constraint< Real >> &con, const Vector< Real > &mul, const Vector< Real > &res)
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Vector< Real > &lam) const
const Ptr< BoundConstraint< Real > > bnd_
Real residual(Vector< Real > &r, const Vector< Real > &y) const
void project_ssn(Vector< Real > &x, Vector< Real > &lam, Vector< Real > &dlam, std::ostream &stream=std::cout) const
Defines the linear algebra or vector space interface.
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
void solve_newton_system(Vector< Real > &s, const Vector< Real > &r, const Vector< Real > &y, const Real mu, const Real rho, int &iter, int &flag) const
Ptr< Vector< Real > > dlam_
virtual int dimension() const
Return dimension of the vector space.
Provides the interface to apply upper and lower bound constraints.
Ptr< Vector< Real > > xnew_
virtual void set(const Vector &x)
Set where .
virtual Real norm() const =0
Returns where .
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Defines the general constraint operator interface.