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 Brents' 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 <<
" Brents' 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 low";
289 stream << std::setw(15) << std::left <<
"res low";
290 stream << std::setw(15) << std::left <<
"lam up";
291 stream << std::setw(15) << std::left <<
"res up";
294 const Real half(0.5), two(2), three(3);
295 const Real eps(ROL_EPSILON<Real>()), tol0(rtol);
296 Real d1(1), d2(1), tol(1);
297 Real p(0), q(0), r(0), m(0);
298 lam = lamUpper; res = resUpper;
299 update_primal(*xnew_,x,lamUpper);
300 for (cnt = 0; cnt < maxit_; cnt++) {
302 lam = lamLower; res = resLower;
303 d1 = lamUpper-lamLower; d2 = d1;
305 if (std::abs(res) < std::abs(resUpper)) {
306 lamLower = lamUpper; lamUpper = lam; lam = lamLower;
307 resLower = resUpper; resUpper = res; res = resLower;
309 tol = two*eps*std::abs(lamUpper) + half*tol0;
310 m = half*(lam - lamUpper);
311 if (std::abs(m) <= tol || std::abs(resUpper) <= rtol)
break;
312 if (std::abs(d2) < tol || std::abs(resLower) <= std::abs(resUpper)) {
316 s = resUpper/resLower;
317 if (lamLower == lam) {
324 p = s*(two*m*q*(q-r)-(lamUpper-lamLower)*(r-one));
325 q = (q-one)*(r-one)*(s-one);
327 if (p >
zero) q = -q;
329 if (two*p < three*m*q-std::abs(tol*q) && p < std::abs(half*d2*q)) {
336 lamLower = lamUpper; resLower = resUpper;
337 if (std::abs(d1) > tol) lamUpper += d1;
338 else if (m >
zero) lamUpper += tol;
339 else lamUpper -= tol;
340 update_primal(*xnew_,x,lamUpper);
341 resUpper = residual(*xnew_);
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 << lamLower;
350 stream << std::setw(15) << std::left << resLower;
351 stream << std::setw(15) << std::left << lamUpper;
352 stream << std::setw(15) << std::left << resUpper;
356 if (verbosity_ > 2) {
359 stream << std::setw(6) << std::left << cnt;
360 stream << std::setw(15) << std::left << rtol;
361 stream << std::setw(15) << std::left << lam;
362 stream << std::setw(15) << std::left << res;
363 stream << std::setw(15) << std::left << lamLower;
364 stream << std::setw(15) << std::left << resLower;
365 stream << std::setw(15) << std::left << lamUpper;
366 stream << std::setw(15) << std::left << resUpper;
374 if (std::abs(res) > rtol ) {
376 stream <<
">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
377 stream << std::abs(res) <<
" rtol = " << rtol << std::endl;
379 stream.flags(streamFlags);