10 #ifndef ROL_SEMISMOOTHNEWTONPROJECTION_DEF_H
11 #define ROL_SEMISMOOTHNEWTONPROJECTION_DEF_H
15 template<
typename Real>
27 DEFAULT_factor_ (0.5),
28 DEFAULT_regscale_ (1e-4),
29 DEFAULT_errscale_ (1e-2),
30 DEFAULT_maxit_ (5000),
32 DEFAULT_verbosity_ (0),
33 DEFAULT_useproj_ (false),
34 atol_ (DEFAULT_atol_),
35 rtol_ (DEFAULT_rtol_),
36 stol_ (DEFAULT_stol_),
37 decr_ (DEFAULT_decr_),
38 factor_ (DEFAULT_factor_),
39 regscale_ (DEFAULT_regscale_),
40 errscale_ (DEFAULT_errscale_),
41 maxit_ (DEFAULT_maxit_),
42 lstype_ (DEFAULT_lstype_),
43 verbosity_ (DEFAULT_verbosity_),
44 useproj_ (DEFAULT_useproj_) {
51 list.sublist(
"General").sublist(
"Krylov").set(
"Type",
"Conjugate Gradients");
52 list.sublist(
"General").sublist(
"Krylov").set(
"Absolute Tolerance", 1e-6);
53 list.sublist(
"General").sublist(
"Krylov").set(
"Relative Tolerance", 1e-4);
54 list.sublist(
"General").sublist(
"Krylov").set(
"Iteration Limit",
dim_);
55 list.sublist(
"General").set(
"Inexact Hessian-Times-A-Vector",
false);
56 krylov_ = KrylovFactory<Real>(list);
61 template<
typename Real>
74 DEFAULT_factor_ (0.5),
75 DEFAULT_regscale_ (1e-4),
76 DEFAULT_errscale_ (1e-2),
77 DEFAULT_maxit_ (5000),
79 DEFAULT_verbosity_ (0),
80 DEFAULT_useproj_ (false),
81 atol_ (DEFAULT_atol_),
82 rtol_ (DEFAULT_rtol_),
83 stol_ (DEFAULT_stol_),
84 decr_ (DEFAULT_decr_),
85 factor_ (DEFAULT_factor_),
86 regscale_ (DEFAULT_regscale_),
87 errscale_ (DEFAULT_errscale_),
88 maxit_ (DEFAULT_maxit_),
89 lstype_ (DEFAULT_lstype_),
90 verbosity_ (DEFAULT_verbosity_),
91 useproj_ (DEFAULT_useproj_) {
96 ParameterList &ppl = list.sublist(
"General").sublist(
"Polyhedral Projection");
100 decr_ = ppl.sublist(
"Semismooth Newton").get(
"Sufficient Decrease Tolerance",
DEFAULT_decr_);
110 klist.sublist(
"General").sublist(
"Krylov") = ppl.sublist(
"Semismooth Newton").sublist(
"Krylov");
111 klist.sublist(
"General").set(
"Inexact Hessian-Times-A-Vector",
false);
112 krylov_ = KrylovFactory<Real>(klist);
117 template<
typename Real>
119 if (con_ == nullPtr) {
123 project_ssn(x, *mul_, *dlam_, stream);
127 template<
typename Real>
129 Real tol(std::sqrt(ROL_EPSILON<Real>()));
131 con_->value(r,y,tol);
135 template<
typename Real>
143 Ptr<Precond> M = makePtr<Precond>(mu);
144 Ptr<Jacobian> J = makePtr<Jacobian>(con_,bnd_,makePtrFromRef(y),xdual_,xprim_,mu);
145 krylov_->run(s,*J,r,*M,iter,flag);
148 template<
typename Real>
152 Real tol(std::sqrt(ROL_EPSILON<Real>()));
155 con_->applyAdjointJacobian(*xdual_,lam,x,tol);
156 y.
plus(xdual_->dual());
160 template<
typename Real>
164 std::ostream &stream)
const {
165 const Real
zero(0), half(0.5), one(1);
167 update_primal(*xnew_,x,lam);
168 Real rnorm = residual(*res_,*xnew_);
173 Real alpha(1), tmp(0), mu(0), rho(1), dd(0);
174 int iter(0), flag(0);
175 std::ios_base::fmtflags streamFlags(stream.flags());
176 if (verbosity_ > 2) {
178 stream << std::scientific << std::setprecision(6);
179 stream <<
" Polyhedral Projection using Dual Semismooth Newton" << std::endl;
181 stream << std::setw(6) << std::left <<
"iter";
182 stream << std::setw(15) << std::left <<
"rnorm";
183 stream << std::setw(15) << std::left <<
"alpha";
184 stream << std::setw(15) << std::left <<
"mu";
185 stream << std::setw(15) << std::left <<
"rho";
186 stream << std::setw(15) << std::left <<
"rtol";
187 stream << std::setw(8) << std::left <<
"kiter";
188 stream << std::setw(8) << std::left <<
"kflag";
191 for (
int cnt = 0; cnt < maxit_; ++cnt) {
193 mu = regscale_*std::max(rnorm,std::sqrt(rnorm));
194 rho = std::min(half,errscale_*std::min(std::sqrt(rnorm),rnorm));
195 solve_newton_system(dlam,*res_,*xnew_,mu,rho,iter,flag);
196 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
197 update_primal(*xnew_,x,*lnew_);
200 tmp = residual(*res_,*xnew_);
201 while ( tmp > (one-decr_*alpha)*rnorm && alpha > stol_ ) {
203 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
204 update_primal(*xnew_,x,*lnew_);
205 tmp = residual(*res_,*xnew_);
210 rnorm = residual(*res_,*xnew_);
212 tmp = dlam.
apply(*res_);
214 while ( tmp < decr_*(one-rho)*mu*dd && alpha > stol_ ) {
216 lnew_->set(lam); lnew_->axpy(-alpha, dlam);
217 update_primal(*xnew_,x,*lnew_);
218 rnorm = residual(*res_,*xnew_);
220 tmp = dlam.
apply(*res_);
227 lam.
axpy(-alpha*tmp/(rnorm*rnorm),res_->dual());
228 update_primal(*xnew_,x,lam);
229 rnorm = residual(*res_,*xnew_);
231 if (verbosity_ > 2) {
233 stream << std::setw(6) << std::left << cnt;
234 stream << std::setw(15) << std::left << rnorm;
235 stream << std::setw(15) << std::left << alpha;
236 stream << std::setw(15) << std::left << mu;
237 stream << std::setw(15) << std::left << rho;
238 stream << std::setw(15) << std::left << ctol_;
239 stream << std::setw(8) << std::left << iter;
240 stream << std::setw(8) << std::left << flag;
243 if (rnorm <= ctol_)
break;
246 if (verbosity_ > 2) {
251 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
252 stream << rnorm <<
" rtol = " << ctol_ << std::endl;
255 stream.flags(streamFlags);
258 template<
typename Real>
261 Real resl = ROL_INF<Real>(), resu = ROL_INF<Real>();
262 if (bnd_->isLowerActivated()) resl = residual(*res_,*bnd_->getLowerBound());
263 if (bnd_->isUpperActivated()) resu = residual(*res_,*bnd_->getUpperBound());
264 Real res0 = std::max(resl,resu);
265 if (res0 < atol_) res0 =
static_cast<Real
>(1);
266 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.