45 #ifndef ROL_DAIFLETCHERPROJECTION_DEF_H 46 #define ROL_DAIFLETCHERPROJECTION_DEF_H 50 template<
typename Real>
61 DEFAULT_maxit_ (5000),
62 DEFAULT_verbosity_ (0),
63 atol_ (DEFAULT_atol_),
64 rtol_ (DEFAULT_rtol_),
65 ltol_ (DEFAULT_ltol_),
66 maxit_ (DEFAULT_maxit_),
67 verbosity_ (DEFAULT_verbosity_) {
69 ROL_TEST_FOR_EXCEPTION(
dim_!=1,std::logic_error,
70 ">>> ROL::DaiFletcherProjection : The range of the linear constraint must be one dimensional!");
72 mul1_ =
static_cast<Real
>(0);
73 dlam1_ =
static_cast<Real
>(2);
75 Real tol(std::sqrt(ROL_EPSILON<Real>()));
80 mul_->setScalar(static_cast<Real>(1));
90 Real res0 = std::max(resl,resu);
92 res0 =
static_cast<Real
>(1);
97 template<
typename Real>
106 atol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Absolute Tolerance",
DEFAULT_atol_);
107 rtol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Relative Tolerance",
DEFAULT_rtol_);
108 ltol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Multiplier Tolerance",
DEFAULT_ltol_);
109 maxit_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Iteration Limit",
DEFAULT_maxit_);
113 template<
typename Real>
115 if (con_ == nullPtr) {
119 mul1_ = -residual(x)/cdot_;
121 dlam1_ =
static_cast<Real
>(2);
123 project_df(x, mul1_, dlam1_, stream);
124 mul_->setScalar(mul1_);
128 template<
typename Real>
130 return xprim_->dot(x) + b_;
133 template<
typename Real>
140 template<
typename Real>
142 const Real
zero(0), one(1), two(2), c1(0.1), c2(0.75), c3(0.25);
143 Real lamLower(0), lamUpper(0), lamNew(0), res(0), resLower(0), resUpper(0), s(0);
147 update_primal(*xnew_,x,lam);
148 res = residual(*xnew_);
153 std::ios_base::fmtflags streamFlags(stream.flags());
154 if (verbosity_ > 2) {
155 stream << std::scientific << std::setprecision(6);
157 stream <<
" Polyhedral Projection using the Dai-Fletcher Algorithm" << std::endl;
158 stream <<
" Bracketing Phase" << std::endl;
165 update_primal(*xnew_,x,lam);
166 res = residual(*xnew_);
167 if (verbosity_ > 2) {
169 stream << std::setw(6) << std::left <<
"iter";
170 stream << std::setw(15) << std::left <<
"lam";
171 stream << std::setw(15) << std::left <<
"res";
172 stream << std::setw(15) << std::left <<
"lower lam";
173 stream << std::setw(15) << std::left <<
"lower res";
176 stream << std::setw(6) << std::left << cnt;
177 stream << std::setw(15) << std::left << lam;
178 stream << std::setw(15) << std::left << res;
179 stream << std::setw(15) << std::left << lamLower;
180 stream << std::setw(15) << std::left << resLower;
183 while ( res <
zero && std::abs(res) > rtol && cnt < maxit_ ) {
184 s = std::max(resLower/res-one,c1);
189 update_primal(*xnew_,x,lam);
190 res = residual(*xnew_);
192 if (verbosity_ > 2) {
194 stream << std::setw(6) << std::left << cnt;
195 stream << std::setw(15) << std::left << lam;
196 stream << std::setw(15) << std::left << res;
197 stream << std::setw(15) << std::left << lamLower;
198 stream << std::setw(15) << std::left << resLower;
209 update_primal(*xnew_,x,lam);
210 res = residual(*xnew_);
211 if (verbosity_ > 2) {
213 stream << std::setw(6) << std::left <<
"iter";
214 stream << std::setw(15) << std::left <<
"lam";
215 stream << std::setw(15) << std::left <<
"res";
216 stream << std::setw(15) << std::left <<
"upper lam";
217 stream << std::setw(15) << std::left <<
"upper res";
220 stream << std::setw(6) << std::left << cnt;
221 stream << std::setw(15) << std::left << lam;
222 stream << std::setw(15) << std::left << res;
223 stream << std::setw(15) << std::left << lamUpper;
224 stream << std::setw(15) << std::left << resUpper;
227 while ( res >
zero && std::abs(res) > rtol && cnt < maxit_ ) {
228 s = std::max(resUpper/res-one,c1);
233 update_primal(*xnew_,x,lam);
234 res = residual(*xnew_);
236 if (verbosity_ > 2) {
238 stream << std::setw(6) << std::left << cnt;
239 stream << std::setw(15) << std::left << lam;
240 stream << std::setw(15) << std::left << res;
241 stream << std::setw(15) << std::left << lamUpper;
242 stream << std::setw(15) << std::left << resUpper;
249 if (verbosity_ > 2) {
250 stream <<
" Bracket: ";
251 stream << std::setw(15) << std::left << lamLower;
252 stream << std::setw(15) << std::left << lamUpper;
257 rtol = ctol_*std::max(one,std::min(std::abs(resLower),std::abs(resUpper)));
261 s = (resUpper - resLower) / resUpper;
262 lam = (resUpper * lamLower - resLower * lamUpper) / (resUpper - resLower);
263 dlam = lamUpper - lam;
264 update_primal(*xnew_,x,lam);
265 res = residual(*xnew_);
267 if (verbosity_ > 2) {
269 stream <<
" Secant Phase" << std::endl;
271 stream << std::setw(6) << std::left <<
"iter";
272 stream << std::setw(15) << std::left <<
"lam";
273 stream << std::setw(15) << std::left <<
"res";
274 stream << std::setw(15) << std::left <<
"stepsize";
275 stream << std::setw(15) << std::left <<
"rtol";
276 stream << std::setw(15) << std::left <<
"lbnd";
277 stream << std::setw(15) << std::left <<
"lres";
278 stream << std::setw(15) << std::left <<
"ubnd";
279 stream << std::setw(15) << std::left <<
"ures";
282 stream << std::setw(6) << std::left << cnt;
283 stream << std::setw(15) << std::left << lam;
284 stream << std::setw(15) << std::left << res;
285 stream << std::setw(15) << std::left << dlam;
286 stream << std::setw(15) << std::left << rtol;
287 stream << std::setw(15) << std::left << lamLower;
288 stream << std::setw(15) << std::left << resLower;
289 stream << std::setw(15) << std::left << lamUpper;
290 stream << std::setw(15) << std::left << resUpper;
293 for (cnt = 1; cnt < maxit_; cnt++) {
295 if ( std::abs(res) <= rtol ||
296 std::abs(lamUpper-lamLower) < ltol_*std::max(std::abs(lamUpper),std::abs(lamLower)) ) {
307 s = (resUpper - resLower) / resUpper;
308 lam = (lamLower * resUpper - lamUpper * resLower) / (resUpper - resLower);
309 dlam = lamUpper - lam;
315 if (resUpper <= (c1+one)*res) {
316 dlam = (lamUpper - lam) / c1;
317 lamNew = std::max(lam - dlam, c2*lamLower + c3*lam);
320 lamNew = std::max((lam * resUpper - lamUpper * res) / (resUpper - res),
321 c2*lamLower + c3*lam);
327 s = (lamUpper - lamLower) / (lamUpper - lam);
337 s = (resUpper - resLower) / resUpper;
338 lam = (lamLower * resUpper - lamUpper * resLower) / (resUpper - resLower);
339 dlam = lamUpper - lam;
345 if (resLower >= (c1+one)*res) {
346 dlam = (lam - lamLower) / c1;
347 lamNew = std::max(lam + dlam, c2*lamUpper + c3*lam);
350 lamNew = std::max((lamLower * res - lam * resLower) / (res - resLower),
351 c2*lamUpper + c3*lam);
352 dlam = lamNew - lamLower;
357 s = (lamUpper - lamLower) / (lamUpper - lam);
360 update_primal(*xnew_,x,lam);
361 res = residual(*xnew_);
363 if (verbosity_ > 2) {
365 stream << std::setw(6) << std::left << cnt;
366 stream << std::setw(15) << std::left << lam;
367 stream << std::setw(15) << std::left << res;
368 stream << std::setw(15) << std::left << dlam;
369 stream << std::setw(15) << std::left << rtol;
370 stream << std::setw(15) << std::left << lamLower;
371 stream << std::setw(15) << std::left << resLower;
372 stream << std::setw(15) << std::left << lamUpper;
373 stream << std::setw(15) << std::left << resUpper;
377 if (verbosity_ > 2) {
382 if (std::abs(res) > rtol ) {
384 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
385 stream << std::abs(res) <<
" rtol = " << rtol << std::endl;
387 stream.flags(streamFlags);
Ptr< Vector< Real > > xprim_
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Ptr< Vector< Real > > res_
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Ptr< Vector< Real > > xdual_
const Ptr< BoundConstraint< Real > > bnd_
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Real lam) const
Defines the linear algebra or vector space interface.
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
DaiFletcherProjection(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 project(Vector< Real > &x, std::ostream &stream=std::cout) override
const Ptr< Constraint< Real > > con_
Ptr< Vector< Real > > mul_
virtual int dimension() const
Return dimension of the vector space.
void project_df(Vector< Real > &x, Real &lam, Real &dlam, std::ostream &stream=std::cout) const
Provides the interface to apply upper and lower bound constraints.
Real residual(const Vector< Real > &x) const
virtual void set(const Vector &x)
Set where .
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Ptr< Vector< Real > > xnew_
Defines the general constraint operator interface.