45 #ifndef ROL_RIDDERSPROJECTION_DEF_H
46 #define ROL_RIDDERSPROJECTION_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_) {
71 template<
typename Real>
83 DEFAULT_maxit_ (5000),
84 DEFAULT_verbosity_ (0),
85 atol_ (DEFAULT_atol_),
86 rtol_ (DEFAULT_rtol_),
87 ltol_ (DEFAULT_ltol_),
88 maxit_ (DEFAULT_maxit_),
89 verbosity_ (DEFAULT_verbosity_) {
90 atol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Absolute Tolerance",
DEFAULT_atol_);
91 rtol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Relative Tolerance",
DEFAULT_rtol_);
92 ltol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Multiplier Tolerance",
DEFAULT_ltol_);
93 maxit_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Iteration Limit",
DEFAULT_maxit_);
98 template<
typename Real>
106 ROL_TEST_FOR_EXCEPTION(dim_!=1,std::logic_error,
107 ">>> ROL::RiddersProjection : The range of the linear constraint must be one dimensional!");
108 xnew_ = xprim.
clone();
110 mul1_ =
static_cast<Real
>(0);
111 dlam1_ =
static_cast<Real
>(2);
113 Real tol(std::sqrt(ROL_EPSILON<Real>()));
116 con_->value(*res_,*xprim_,tol);
117 b_ = res_->dot(*res_->basis(0));
118 mul_->setScalar(static_cast<Real>(1));
119 con_->applyAdjointJacobian(*xdual_,*mul_,xprim,tol);
120 xprim_->set(xdual_->dual());
121 cdot_ = xprim_->dot(*xprim_);
126 Real resl = ROL_INF<Real>(), resu = ROL_INF<Real>();
127 if (bnd_->isLowerActivated()) resl = residual(*bnd_->getLowerBound());
128 if (bnd_->isUpperActivated()) resu = residual(*bnd_->getUpperBound());
129 Real res0 = std::max(resl,resu);
130 if (res0 < atol_) res0 =
static_cast<Real
>(1);
131 ctol_ = std::min(atol_,rtol_*res0);
134 template<
typename Real>
136 if (con_ == nullPtr) {
140 mul1_ = -residual(x)/cdot_;
142 dlam1_ =
static_cast<Real
>(2);
144 project_df(x, mul1_, dlam1_, stream);
145 mul_->setScalar(mul1_);
149 template<
typename Real>
151 return xprim_->dot(x) + b_;
154 template<
typename Real>
161 template<
typename Real>
163 const Real
zero(0), one(1), c1(0.1);
164 Real lamLower(0), lamUpper(0), res(0), resLower(0), resUpper(0), s(0);
168 update_primal(*xnew_,x,lam);
169 res = residual(*xnew_);
174 std::ios_base::fmtflags streamFlags(stream.flags());
175 if (verbosity_ > 2) {
176 stream << std::scientific << std::setprecision(6);
178 stream <<
" Polyhedral Projection using Ridders' Algorithm" << std::endl;
179 stream <<
" Bracketing Phase" << std::endl;
186 update_primal(*xnew_,x,lam);
187 res = residual(*xnew_);
188 if (verbosity_ > 2) {
190 stream << std::setw(6) << std::left <<
"iter";
191 stream << std::setw(15) << std::left <<
"lam";
192 stream << std::setw(15) << std::left <<
"res";
193 stream << std::setw(15) << std::left <<
"lower lam";
194 stream << std::setw(15) << std::left <<
"lower res";
197 stream << std::setw(6) << std::left << cnt;
198 stream << std::setw(15) << std::left << lam;
199 stream << std::setw(15) << std::left << res;
200 stream << std::setw(15) << std::left << lamLower;
201 stream << std::setw(15) << std::left << resLower;
204 while ( res <
zero && std::abs(res) > rtol && cnt < maxit_ ) {
205 s = std::max(resLower/res-one,c1);
210 update_primal(*xnew_,x,lam);
211 res = residual(*xnew_);
213 if (verbosity_ > 2) {
215 stream << std::setw(6) << std::left << cnt;
216 stream << std::setw(15) << std::left << lam;
217 stream << std::setw(15) << std::left << res;
218 stream << std::setw(15) << std::left << lamLower;
219 stream << std::setw(15) << std::left << resLower;
230 update_primal(*xnew_,x,lam);
231 res = residual(*xnew_);
232 if (verbosity_ > 2) {
234 stream << std::setw(6) << std::left <<
"iter";
235 stream << std::setw(15) << std::left <<
"lam";
236 stream << std::setw(15) << std::left <<
"res";
237 stream << std::setw(15) << std::left <<
"upper lam";
238 stream << std::setw(15) << std::left <<
"upper res";
241 stream << std::setw(6) << std::left << cnt;
242 stream << std::setw(15) << std::left << lam;
243 stream << std::setw(15) << std::left << res;
244 stream << std::setw(15) << std::left << lamUpper;
245 stream << std::setw(15) << std::left << resUpper;
248 while ( res >
zero && std::abs(res) > rtol && cnt < maxit_ ) {
249 s = std::max(resUpper/res-one,c1);
254 update_primal(*xnew_,x,lam);
255 res = residual(*xnew_);
257 if (verbosity_ > 2) {
259 stream << std::setw(6) << std::left << cnt;
260 stream << std::setw(15) << std::left << lam;
261 stream << std::setw(15) << std::left << res;
262 stream << std::setw(15) << std::left << lamUpper;
263 stream << std::setw(15) << std::left << resUpper;
270 if (verbosity_ > 2) {
271 stream <<
" Bracket: ";
272 stream << std::setw(15) << std::left << lamLower;
273 stream << std::setw(15) << std::left << lamUpper;
280 if (verbosity_ > 2) {
282 stream <<
" Ridders' Phase" << std::endl;
284 stream << std::setw(6) << std::left <<
"iter";
285 stream << std::setw(15) << std::left <<
"rtol";
286 stream << std::setw(15) << std::left <<
"lam";
287 stream << std::setw(15) << std::left <<
"res";
288 stream << std::setw(15) << std::left <<
"lam mid";
289 stream << std::setw(15) << std::left <<
"res mid";
290 stream << std::setw(15) << std::left <<
"lam low";
291 stream << std::setw(15) << std::left <<
"res low";
292 stream << std::setw(15) << std::left <<
"lam up";
293 stream << std::setw(15) << std::left <<
"res up";
296 const Real half(0.5);
297 Real lamMid(0), resMid(0);
298 for (cnt = 0; cnt < maxit_; cnt++) {
301 if (std::abs(lamUpper-lamLower) < ltol_)
break;
303 lamMid = half*(lamUpper+lamLower);
304 update_primal(*xnew_,x,lamMid);
305 resMid = residual(*xnew_);
306 if (std::abs(resMid) <= rtol) {
312 lam = lamMid-(lamMid-lamLower)*resMid/std::sqrt(resMid*resMid-resLower*resUpper);
313 update_primal(*xnew_,x,lam);
314 res = residual(*xnew_);
315 if (std::abs(res) <= rtol)
break;
318 if (resMid < -rtol) {
319 resLower = (lam < lamMid ? resMid : res);
320 lamLower = (lam < lamMid ? lamMid : lam);
330 if (resMid < -rtol) {
337 resUpper = (lam < lamMid ? res : resMid);
338 lamUpper = (lam < lamMid ? lam : lamMid);
343 if (verbosity_ > 2) {
345 stream << std::setw(6) << std::left << cnt;
346 stream << std::setw(15) << std::left << rtol;
347 stream << std::setw(15) << std::left << lam;
348 stream << std::setw(15) << std::left << res;
349 stream << std::setw(15) << std::left << lamMid;
350 stream << std::setw(15) << std::left << resMid;
351 stream << std::setw(15) << std::left << lamLower;
352 stream << std::setw(15) << std::left << resLower;
353 stream << std::setw(15) << std::left << lamUpper;
354 stream << std::setw(15) << std::left << resUpper;
358 if (verbosity_ > 2) {
361 stream << std::setw(6) << std::left << cnt;
362 stream << std::setw(15) << std::left << rtol;
363 stream << std::setw(15) << std::left << lam;
364 stream << std::setw(15) << std::left << res;
365 stream << std::setw(15) << std::left << lamMid;
366 stream << std::setw(15) << std::left << resMid;
367 stream << std::setw(15) << std::left << lamLower;
368 stream << std::setw(15) << std::left << resLower;
369 stream << std::setw(15) << std::left << lamUpper;
370 stream << std::setw(15) << std::left << resUpper;
377 if (std::abs(res) > rtol ) {
379 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
380 stream << std::abs(res) <<
" rtol = " << rtol << std::endl;
382 stream.flags(streamFlags);
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Real residual(const Vector< Real > &x) const
Defines the linear algebra or vector space interface.
void project_df(Vector< Real > &x, Real &lam, Real &dlam, std::ostream &stream=std::cout) const
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
RiddersProjection(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)
Provides the interface to apply upper and lower bound constraints.
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Real lam) const
virtual void set(const Vector &x)
Set where .
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
void initialize(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)
Defines the general constraint operator interface.