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",
"Conjugate Gradients");
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 template<
typename Real>
108 DEFAULT_decr_ (1e-4),
109 DEFAULT_factor_ (0.5),
110 DEFAULT_regscale_ (1e-4),
111 DEFAULT_errscale_ (1e-2),
112 DEFAULT_maxit_ (5000),
114 DEFAULT_verbosity_ (0),
115 DEFAULT_useproj_ (false),
116 atol_ (DEFAULT_atol_),
117 rtol_ (DEFAULT_rtol_),
118 stol_ (DEFAULT_stol_),
119 decr_ (DEFAULT_decr_),
120 factor_ (DEFAULT_factor_),
121 regscale_ (DEFAULT_regscale_),
122 errscale_ (DEFAULT_errscale_),
123 maxit_ (DEFAULT_maxit_),
124 lstype_ (DEFAULT_lstype_),
125 verbosity_ (DEFAULT_verbosity_),
126 useproj_ (DEFAULT_useproj_) {
131 ParameterList &ppl = list.sublist(
"General").sublist(
"Polyhedral Projection");
135 decr_ = ppl.sublist(
"Semismooth Newton").get(
"Sufficient Decrease Tolerance",
DEFAULT_decr_);
145 klist.sublist(
"General").sublist(
"Krylov") = ppl.sublist(
"Semismooth Newton").sublist(
"Krylov");
146 klist.sublist(
"General").set(
"Inexact Hessian-Times-A-Vector",
false);
147 krylov_ = KrylovFactory<Real>(klist);
152 template<
typename Real>
154 if (con_ == nullPtr) {
158 project_ssn(x, *mul_, *dlam_, stream);
162 template<
typename Real>
164 Real tol(std::sqrt(ROL_EPSILON<Real>()));
166 con_->value(r,y,tol);
170 template<
typename Real>
178 Ptr<Precond> M = makePtr<Precond>(mu);
179 Ptr<Jacobian> J = makePtr<Jacobian>(con_,bnd_,makePtrFromRef(y),xdual_,xprim_,mu);
180 krylov_->run(s,*J,r,*M,iter,flag);
183 template<
typename Real>
187 Real tol(std::sqrt(ROL_EPSILON<Real>()));
190 con_->applyAdjointJacobian(*xdual_,lam,x,tol);
191 y.
plus(xdual_->dual());
195 template<
typename Real>
199 std::ostream &stream)
const {
200 const Real
zero(0), half(0.5), one(1);
202 update_primal(*xnew_,x,lam);
203 Real rnorm = residual(*res_,*xnew_);
208 Real alpha(1), tmp(0), mu(0), rho(1), dd(0);
209 int iter(0), flag(0);
210 std::ios_base::fmtflags streamFlags(stream.flags());
211 if (verbosity_ > 2) {
213 stream << std::scientific << std::setprecision(6);
214 stream <<
" Polyhedral Projection using Dual Semismooth Newton" << std::endl;
216 stream << std::setw(6) << std::left <<
"iter";
217 stream << std::setw(15) << std::left <<
"rnorm";
218 stream << std::setw(15) << std::left <<
"alpha";
219 stream << std::setw(15) << std::left <<
"mu";
220 stream << std::setw(15) << std::left <<
"rho";
221 stream << std::setw(15) << std::left <<
"rtol";
222 stream << std::setw(8) << std::left <<
"kiter";
223 stream << std::setw(8) << std::left <<
"kflag";
226 for (
int cnt = 0; cnt < maxit_; ++cnt) {
228 mu = regscale_*std::max(rnorm,std::sqrt(rnorm));
229 rho = std::min(half,errscale_*std::min(std::sqrt(rnorm),rnorm));
230 solve_newton_system(dlam,*res_,*xnew_,mu,rho,iter,flag);
231 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
232 update_primal(*xnew_,x,*lnew_);
235 tmp = residual(*res_,*xnew_);
236 while ( tmp > (one-decr_*alpha)*rnorm && alpha > stol_ ) {
238 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
239 update_primal(*xnew_,x,*lnew_);
240 tmp = residual(*res_,*xnew_);
245 rnorm = residual(*res_,*xnew_);
247 tmp = dlam.
apply(*res_);
249 while ( tmp < decr_*(one-rho)*mu*dd && alpha > stol_ ) {
251 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
252 update_primal(*xnew_,x,*lnew_);
253 rnorm = residual(*res_,*xnew_);
255 tmp = dlam.
apply(*res_);
262 lam.
axpy(-alpha*tmp/(rnorm*rnorm),res_->dual());
263 update_primal(*xnew_,x,lam);
264 rnorm = residual(*res_,*xnew_);
266 if (verbosity_ > 2) {
268 stream << std::setw(6) << std::left << cnt;
269 stream << std::setw(15) << std::left << rnorm;
270 stream << std::setw(15) << std::left << alpha;
271 stream << std::setw(15) << std::left << mu;
272 stream << std::setw(15) << std::left << rho;
273 stream << std::setw(15) << std::left << ctol_;
274 stream << std::setw(8) << std::left << iter;
275 stream << std::setw(8) << std::left << flag;
278 if (rnorm <= ctol_)
break;
281 if (verbosity_ > 2) {
286 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
287 stream << rnorm <<
" rtol = " << ctol_ << std::endl;
290 stream.flags(streamFlags);
293 template<
typename Real>
296 Real resl = ROL_INF<Real>(), resu = ROL_INF<Real>();
297 if (bnd_->isLowerActivated()) resl = residual(*res_,*bnd_->getLowerBound());
298 if (bnd_->isUpperActivated()) resu = residual(*res_,*bnd_->getUpperBound());
299 Real res0 = std::max(resl,resu);
300 if (res0 < atol_) res0 =
static_cast<Real
>(1);
301 return std::min(atol_,rtol_*res0);
Ptr< Vector< Real > > lnew_
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
Ptr< Krylov< Real > > krylov_
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
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)
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
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
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()
Ptr< Vector< Real > > dlam_
Provides the interface to apply upper and lower bound constraints.
Real residual(Vector< Real > &r, const Vector< Real > &y) const
Real compute_tolerance() const
Ptr< Vector< Real > > xnew_
void project_ssn(Vector< Real > &x, Vector< Real > &lam, Vector< Real > &dlam, std::ostream &stream=std::cout) const
virtual void set(const Vector &x)
Set where .
virtual Real norm() const =0
Returns where .
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Vector< Real > &lam) const
Defines the general constraint operator interface.