10 #ifndef ROL_RIDDERSPROJECTION_DEF_H
11 #define ROL_RIDDERSPROJECTION_DEF_H
15 template<
typename Real>
26 DEFAULT_maxit_ (5000),
27 DEFAULT_verbosity_ (0),
28 atol_ (DEFAULT_atol_),
29 rtol_ (DEFAULT_rtol_),
30 ltol_ (DEFAULT_ltol_),
31 maxit_ (DEFAULT_maxit_),
32 verbosity_ (DEFAULT_verbosity_) {
36 template<
typename Real>
48 DEFAULT_maxit_ (5000),
49 DEFAULT_verbosity_ (0),
50 atol_ (DEFAULT_atol_),
51 rtol_ (DEFAULT_rtol_),
52 ltol_ (DEFAULT_ltol_),
53 maxit_ (DEFAULT_maxit_),
54 verbosity_ (DEFAULT_verbosity_) {
55 atol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Absolute Tolerance",
DEFAULT_atol_);
56 rtol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Relative Tolerance",
DEFAULT_rtol_);
57 ltol_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Multiplier Tolerance",
DEFAULT_ltol_);
58 maxit_ = list.sublist(
"General").sublist(
"Polyhedral Projection").get(
"Iteration Limit",
DEFAULT_maxit_);
63 template<
typename Real>
71 ROL_TEST_FOR_EXCEPTION(dim_!=1,std::logic_error,
72 ">>> ROL::RiddersProjection : The range of the linear constraint must be one dimensional!");
73 xnew_ = xprim.
clone();
75 mul1_ =
static_cast<Real
>(0);
76 dlam1_ =
static_cast<Real
>(2);
78 Real tol(std::sqrt(ROL_EPSILON<Real>()));
81 con_->value(*res_,*xprim_,tol);
82 b_ = res_->dot(*res_->basis(0));
83 mul_->setScalar(static_cast<Real>(1));
84 con_->applyAdjointJacobian(*xdual_,*mul_,xprim,tol);
85 xprim_->set(xdual_->dual());
86 cdot_ = xprim_->dot(*xprim_);
91 Real resl = ROL_INF<Real>(), resu = ROL_INF<Real>();
92 if (bnd_->isLowerActivated()) resl = residual(*bnd_->getLowerBound());
93 if (bnd_->isUpperActivated()) resu = residual(*bnd_->getUpperBound());
94 Real res0 = std::max(resl,resu);
95 if (res0 < atol_) res0 =
static_cast<Real
>(1);
96 ctol_ = std::min(atol_,rtol_*res0);
99 template<
typename Real>
101 if (con_ == nullPtr) {
105 mul1_ = -residual(x)/cdot_;
107 dlam1_ =
static_cast<Real
>(2);
109 project_df(x, mul1_, dlam1_, stream);
110 mul_->setScalar(mul1_);
114 template<
typename Real>
116 return xprim_->dot(x) + b_;
119 template<
typename Real>
126 template<
typename Real>
128 const Real
zero(0), one(1), c1(0.1);
129 Real lamLower(0), lamUpper(0), res(0), resLower(0), resUpper(0), s(0);
133 update_primal(*xnew_,x,lam);
134 res = residual(*xnew_);
139 std::ios_base::fmtflags streamFlags(stream.flags());
140 if (verbosity_ > 2) {
141 stream << std::scientific << std::setprecision(6);
143 stream <<
" Polyhedral Projection using Ridders' Algorithm" << std::endl;
144 stream <<
" Bracketing Phase" << std::endl;
151 update_primal(*xnew_,x,lam);
152 res = residual(*xnew_);
153 if (verbosity_ > 2) {
155 stream << std::setw(6) << std::left <<
"iter";
156 stream << std::setw(15) << std::left <<
"lam";
157 stream << std::setw(15) << std::left <<
"res";
158 stream << std::setw(15) << std::left <<
"lower lam";
159 stream << std::setw(15) << std::left <<
"lower res";
162 stream << std::setw(6) << std::left << cnt;
163 stream << std::setw(15) << std::left << lam;
164 stream << std::setw(15) << std::left << res;
165 stream << std::setw(15) << std::left << lamLower;
166 stream << std::setw(15) << std::left << resLower;
169 while ( res <
zero && std::abs(res) > rtol && cnt < maxit_ ) {
170 s = std::max(resLower/res-one,c1);
175 update_primal(*xnew_,x,lam);
176 res = residual(*xnew_);
178 if (verbosity_ > 2) {
180 stream << std::setw(6) << std::left << cnt;
181 stream << std::setw(15) << std::left << lam;
182 stream << std::setw(15) << std::left << res;
183 stream << std::setw(15) << std::left << lamLower;
184 stream << std::setw(15) << std::left << resLower;
195 update_primal(*xnew_,x,lam);
196 res = residual(*xnew_);
197 if (verbosity_ > 2) {
199 stream << std::setw(6) << std::left <<
"iter";
200 stream << std::setw(15) << std::left <<
"lam";
201 stream << std::setw(15) << std::left <<
"res";
202 stream << std::setw(15) << std::left <<
"upper lam";
203 stream << std::setw(15) << std::left <<
"upper res";
206 stream << std::setw(6) << std::left << cnt;
207 stream << std::setw(15) << std::left << lam;
208 stream << std::setw(15) << std::left << res;
209 stream << std::setw(15) << std::left << lamUpper;
210 stream << std::setw(15) << std::left << resUpper;
213 while ( res >
zero && std::abs(res) > rtol && cnt < maxit_ ) {
214 s = std::max(resUpper/res-one,c1);
219 update_primal(*xnew_,x,lam);
220 res = residual(*xnew_);
222 if (verbosity_ > 2) {
224 stream << std::setw(6) << std::left << cnt;
225 stream << std::setw(15) << std::left << lam;
226 stream << std::setw(15) << std::left << res;
227 stream << std::setw(15) << std::left << lamUpper;
228 stream << std::setw(15) << std::left << resUpper;
235 if (verbosity_ > 2) {
236 stream <<
" Bracket: ";
237 stream << std::setw(15) << std::left << lamLower;
238 stream << std::setw(15) << std::left << lamUpper;
245 if (verbosity_ > 2) {
247 stream <<
" Ridders' Phase" << std::endl;
249 stream << std::setw(6) << std::left <<
"iter";
250 stream << std::setw(15) << std::left <<
"rtol";
251 stream << std::setw(15) << std::left <<
"lam";
252 stream << std::setw(15) << std::left <<
"res";
253 stream << std::setw(15) << std::left <<
"lam mid";
254 stream << std::setw(15) << std::left <<
"res mid";
255 stream << std::setw(15) << std::left <<
"lam low";
256 stream << std::setw(15) << std::left <<
"res low";
257 stream << std::setw(15) << std::left <<
"lam up";
258 stream << std::setw(15) << std::left <<
"res up";
261 const Real half(0.5);
262 Real lamMid(0), resMid(0);
263 for (cnt = 0; cnt < maxit_; cnt++) {
266 if (std::abs(lamUpper-lamLower) < ltol_)
break;
268 lamMid = half*(lamUpper+lamLower);
269 update_primal(*xnew_,x,lamMid);
270 resMid = residual(*xnew_);
271 if (std::abs(resMid) <= rtol) {
277 lam = lamMid-(lamMid-lamLower)*resMid/std::sqrt(resMid*resMid-resLower*resUpper);
278 update_primal(*xnew_,x,lam);
279 res = residual(*xnew_);
280 if (std::abs(res) <= rtol)
break;
283 if (resMid < -rtol) {
284 resLower = (lam < lamMid ? resMid : res);
285 lamLower = (lam < lamMid ? lamMid : lam);
295 if (resMid < -rtol) {
302 resUpper = (lam < lamMid ? res : resMid);
303 lamUpper = (lam < lamMid ? lam : lamMid);
308 if (verbosity_ > 2) {
310 stream << std::setw(6) << std::left << cnt;
311 stream << std::setw(15) << std::left << rtol;
312 stream << std::setw(15) << std::left << lam;
313 stream << std::setw(15) << std::left << res;
314 stream << std::setw(15) << std::left << lamMid;
315 stream << std::setw(15) << std::left << resMid;
316 stream << std::setw(15) << std::left << lamLower;
317 stream << std::setw(15) << std::left << resLower;
318 stream << std::setw(15) << std::left << lamUpper;
319 stream << std::setw(15) << std::left << resUpper;
323 if (verbosity_ > 2) {
326 stream << std::setw(6) << std::left << cnt;
327 stream << std::setw(15) << std::left << rtol;
328 stream << std::setw(15) << std::left << lam;
329 stream << std::setw(15) << std::left << res;
330 stream << std::setw(15) << std::left << lamMid;
331 stream << std::setw(15) << std::left << resMid;
332 stream << std::setw(15) << std::left << lamLower;
333 stream << std::setw(15) << std::left << resLower;
334 stream << std::setw(15) << std::left << lamUpper;
335 stream << std::setw(15) << std::left << resUpper;
342 if (std::abs(res) > rtol ) {
344 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
345 stream << std::abs(res) <<
" rtol = " << rtol << std::endl;
347 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.