ROL
ROL_LinMore.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_LINMORE_H
45#define ROL_LINMORE_H
46
51#include "ROL_TrustRegion.hpp"
52#include "ROL_LinMoreModel.hpp"
53#include "ROL_Elementwise_Function.hpp"
54#include "ROL_Elementwise_Reduce.hpp"
55
56namespace ROL {
57
58template<class Real>
59class LinMore : public TrustRegion<Real> {
60private:
61
62 Ptr<Vector<Real>> x_, s_, g_;
63 Ptr<Vector<Real>> pwa1_, pwa2_, dwa1_, dwa2_;
64
66 int maxit_;
67
68 unsigned verbosity_;
69
70 class LowerBreakPoint : public Elementwise::BinaryFunction<Real> {
71 public:
72 Real apply( const Real &x, const Real &y) const {
73 const Real zero(0), one(1);
74 return (x > zero && y < zero) ? -x/y : -one;
75 }
77
78 class UpperBreakPoint : public Elementwise::BinaryFunction<Real> {
79 public:
80 Real apply( const Real &x, const Real &y) const {
81 const Real zero(0), one(1);
82 return (x > zero && y > zero) ? x/y : -one;
83 }
85
86 class PositiveMin : public Elementwise::ReductionOp<Real> {
87 public:
88 void reduce(const Real &input, Real &output) const {
89 const Real zero(0);
90 output = (input<output && input>zero) ? input : output;
91 }
92 void reduce( const volatile Real &input, Real volatile &output ) const {
93 const Real zero(0);
94 output = (input<output && input>zero) ? input : output;
95 }
96 Real initialValue() const {
97 return ROL_INF<Real>();
98 }
99 Elementwise::EReductionType reductionType() const {
100 return Elementwise::REDUCE_MIN;
101 }
103
104 class PositiveMax : public Elementwise::ReductionOp<Real> {
105 public:
106 void reduce(const Real &input, Real &output) const {
107 const Real zero(0);
108 output = (input>output && input>zero) ? input : output;
109 }
110 void reduce( const volatile Real &input, Real volatile &output ) const {
111 const Real zero(0);
112 output = (input>output && input>zero) ? input : output;
113 }
114 Real initialValue() const {
115 return static_cast<Real>(0);
116 }
117 Elementwise::EReductionType reductionType() const {
118 return Elementwise::REDUCE_MAX;
119 }
121
122public:
123
124 // Constructor
125 LinMore( ROL::ParameterList &parlist ) : TrustRegion<Real>(parlist), alpha_(1) {
126 // Unravel Parameter List
127 Real em4(1e-4), em2(1e-2);
128 maxit_ = parlist.sublist("General").sublist("Krylov").get("Iteration Limit",20);
129 tol1_ = parlist.sublist("General").sublist("Krylov").get("Absolute Tolerance",em4);
130 tol2_ = parlist.sublist("General").sublist("Krylov").get("Relative Tolerance",em2);
131 // Get verbosity level
132 verbosity_ = parlist.sublist("General").get("Print Verbosity", 0);
133 }
134
135 void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
137 x_ = x.clone(); s_ = x.clone(); g_ = g.clone();
138 pwa1_ = x.clone(); pwa2_ = x.clone();
139 dwa1_ = g.clone(); dwa2_ = g.clone();
140 }
141
143 Real &snorm,
144 int &iflag,
145 int &iter,
146 const Real del,
147 TrustRegionModel<Real> &model ) {
148 const Real zero(0), half(0.5), one(1);
149 Real tol0 = std::sqrt(ROL_EPSILON<Real>());
150 Real gfnorm(0), gfnormf(0), tol(0), stol(0);
151 int dim = s.dimension();
152 // Compute Cauchy point (TRON notation: x_ = x[1])
153 snorm = dcauchy(*s_,alpha_,*model.getIterate(),model.getGradient()->dual(),
154 del,model,*pwa1_,*pwa2_,*dwa1_); // Solve 1D optimization problem for alpha_
155 x_->set(*model.getIterate()); // TRON notation: model.getIterate() = x[0]
156 x_->plus(*s_); // Set x_ = x[0] + alpha_*g
157 model.getBoundConstraint()->project(*x_); // Project x_ onto bounds
158
159 // Model gradient at s = x[1] - x[0]
160 s.set(*x_); s.axpy(-one,*model.getIterate()); // s_ = x[i+1]-x[0]
161 dynamic_cast<LinMoreModel<Real>&>(model).applyFullHessian(*g_,s,tol0);
162 g_->plus(*model.getGradient());
163 model.getBoundConstraint()->pruneActive(*g_,*x_,zero);
164 gfnorm = g_->norm();
165 if (verbosity_ > 0) {
166 std::cout << std::endl;
167 std::cout << " Computation of Cauchy point" << std::endl;
168 std::cout << " Step length (alpha): " << alpha_ << std::endl;
169 std::cout << " Step length (alpha*g): " << snorm << std::endl;
170 std::cout << " Norm of Cauchy point: " << x_->norm() << std::endl;
171 std::cout << " Norm of free gradient components: " << gfnorm << std::endl;
172 }
173
174 // Main step computation loop
175 // There are at most dim iterations since at least one
176 // face becomes active at each iteration
177 iter = 0;
178 for (int i = 0; i < dim; ++i) {
179 // Run Truncated CG
180 int flagCG = 0, iterCG = 0;
181 tol = tol2_*gfnorm;
182 stol = zero;
183 snorm = dtrpcg(*s_,flagCG,iterCG,*g_,*x_,del,model,
184 tol,stol,maxit_,
185 *pwa1_,*dwa1_,*pwa2_,*dwa2_);
186 iter += iterCG;
187 if (verbosity_ > 0) {
188 std::cout << std::endl;
189 std::cout << " Computation of CG step" << std::endl;
190 std::cout << " Number of faces: " << dim << std::endl;
191 std::cout << " Current face (i): " << i << std::endl;
192 std::cout << " CG step length: " << snorm << std::endl;
193 std::cout << " Number of CG iterations: " << iterCG << std::endl;
194 std::cout << " CG flag: " << flagCG << std::endl;
195 std::cout << " Total number of iterations: " << iter << std::endl;
196 }
197
198 // Projected search
199 snorm = dprsrch(*x_,*s_,g_->dual(),model,*pwa1_,*dwa1_);
200 if (verbosity_ > 0) {
201 std::cout << " Step length (beta*s): " << snorm << std::endl;
202 std::cout << " Iterate length: " << x_->norm() << std::endl;
203 }
204
205 // Model gradient at s = x[i+1] - x[0]
206 s.set(*x_); s.axpy(-one,*model.getIterate()); // s_ = x[i+1]-x[0]
207 dynamic_cast<LinMoreModel<Real>&>(model).applyFullHessian(*g_,s,tol0);
208 g_->plus(*model.getGradient());
209 model.getBoundConstraint()->pruneActive(*g_,*x_,zero);
210 gfnormf = g_->norm();
211 if (verbosity_ > 0) {
212 std::cout << std::endl;
213 std::cout << " Update model gradient" << std::endl;
214 std::cout << " Step length: " << s.norm() << std::endl;
215 std::cout << " Norm of free gradient components: " << gfnormf << std::endl;
216 std::cout << std::endl;
217 }
218
219 // Termination check
220 if (gfnormf <= tol2_*gfnorm) {
221 iflag = 0;
222 break;
223 }
224 else if (iter >= maxit_) {
225 iflag = 1;
226 break;
227 }
228 else if (flagCG == 2) {
229 iflag = 2;
230 break;
231 }
232 else if (flagCG == 3) {
233 iflag = 3;
234 break;
235 }
236
237 // Update free gradient norm
238 gfnorm = gfnormf;
239 }
240 // Update norm of step and update model predicted reduction
241 snorm = s.norm();
242 Real gs(0), q(0);
243 dynamic_cast<LinMoreModel<Real>&>(model).applyFullHessian(*dwa1_,s,tol0);
244 gs = s.dot(model.getGradient()->dual());
245 q = half * s.dot(dwa1_->dual()) + gs;
247 }
248
249private:
250
251 // Compute the projected step s = P(x + alpha*w) - x
252 // Returns the norm of the projected step s
253 // s -- The projected step upon return
254 // w -- The direction vector w (unchanged)
255 // x -- The anchor vector x (unchanged)
256 // alpha -- The step size (unchanged)
257 // model -- Contains the bound constraint information
259 const Vector<Real> &x, const Real alpha,
260 TrustRegionModel<Real> &model) const {
261 s.set(x); s.axpy(alpha,w);
262 model.getBoundConstraint()->project(s);
263 s.axpy(static_cast<Real>(-1),x);
264 return s.norm();
265 }
266
267 // Compute minimal and maximal break points of x+alpha*s
268 // with in the interval [xl,xu] specified by the bound constraint
269 // x -- The anchor vector x (unchanged)
270 // s -- The descent vector s (unchanged)
271 // model -- Contains the bound constraint information
272 // bpmin -- The minimum break point
273 // bpmax -- The maximum break point
274 // pwa -- A primal working vector
275 void dbreakpt(const Vector<Real> &x, const Vector<Real> &s,
277 Real &bpmin, Real &bpmax,
278 Vector<Real> &pwa) {
279 const Real zero(0), one(1);
280 bpmin = one; bpmax = zero;
281 Real lbpmin = one, lbpmax = zero, ubpmin = one, ubpmax = zero;
282 // Compute lower break points
283 if (model.getBoundConstraint()->isLowerActivated()) {
284 pwa.set(x);
285 pwa.axpy(-one,*model.getBoundConstraint()->getLowerBound());
286 pwa.applyBinary(lbp_,s);
287 if (pwa.norm() != zero) {
288 lbpmin = pwa.reduce(pmin_);
289 lbpmax = pwa.reduce(pmax_);
290 }
291 }
292 // Compute upper break points
293 if (model.getBoundConstraint()->isUpperActivated()) {
294 pwa.set(*model.getBoundConstraint()->getUpperBound());
295 pwa.axpy(-one,x);
296 pwa.applyBinary(ubp_,s);
297 if (pwa.norm() != zero) {
298 ubpmin = pwa.reduce(pmin_);
299 ubpmax = pwa.reduce(pmax_);
300 }
301 }
302 bpmin = std::min(lbpmin,ubpmin);
303 bpmax = std::max(lbpmax,ubpmax);
304 if (bpmin > bpmax) {
305 bpmin = zero;
306 bpmax = zero;
307 }
308 if (verbosity_ > 0) {
309 std::cout << std::endl;
310 std::cout << " Computation of break points" << std::endl;
311 std::cout << " Minimum break point: " << bpmin << std::endl;
312 std::cout << " Maximum break point: " << bpmax << std::endl;
313 }
314 }
315
316 // Compute Cauchy point, i.e., the minimizer of q(P(x - alpha*g)-x)
317 // subject to the trust region constraint ||P(x - alpha*g)-x|| <= del
318 // s -- The Cauchy point upon return
319 // alpha -- The step length for the Cauchy point upon return
320 // x -- The anchor vector x (unchanged)
321 // g -- The (dual) gradient vector g (unchanged)
322 // del -- The trust region radius (unchanged)
323 // model -- Contains the objective and bound constraint information
324 // pwa1 -- Primal working array
325 // pwa2 -- Primal working array
326 // dwa -- Dual working array
327 Real dcauchy(Vector<Real> &s, Real &alpha,
328 const Vector<Real> &x, const Vector<Real> &g,
329 const Real del, TrustRegionModel<Real> &model,
330 Vector<Real> &pwa1, Vector<Real> &pwa2, Vector<Real> &dwa) {
331 const Real half(0.5), one(1), mu0(0.01), interpf(0.1), extrapf(10);
332 // const Real zero(0); // Unused
333 Real tol = std::sqrt(ROL_EPSILON<Real>());
334 bool interp = false;
335 Real q(0), gs(0), bpmin(0), bpmax(0), snorm(0);
336 // Compute minimal and maximal break points of x[0] - alpha g[0]
337 pwa1.set(g); pwa1.scale(-one);
338 dbreakpt(x,pwa1,model,bpmin,bpmax,pwa2);
339 // Compute s = P(x[0] - alpha g[0].dual())
340 snorm = dgpstep(s,g,x,-alpha,model);
341 if (snorm > del) {
342 interp = true;
343 }
344 else {
345 dynamic_cast<LinMoreModel<Real>&>(model).applyFullHessian(dwa,s,tol);
346 gs = s.dot(g);
347 q = half * s.dot(dwa.dual()) + gs;
348 interp = (q > mu0*gs);
349 }
350 // Either increase or decrease alpha to find approximate Cauchy point
351 if (interp) {
352 bool search = true;
353 while (search) {
354 alpha *= interpf;
355 snorm = dgpstep(s,g,x,-alpha,model);
356 if (snorm <= del) {
357 dynamic_cast<LinMoreModel<Real>&>(model).applyFullHessian(dwa,s,tol);
358 gs = s.dot(g);
359 q = half * s.dot(dwa.dual()) + gs;
360 search = (q > mu0*gs);
361 }
362 }
363 }
364 else {
365 bool search = true;
366 Real alphas = alpha;
367 while (search && alpha <= bpmax) {
368 alpha *= extrapf;
369 snorm = dgpstep(s,g,x,-alpha,model);
370 if (snorm <= del) {
371 dynamic_cast<LinMoreModel<Real>&>(model).applyFullHessian(dwa,s,tol);
372 gs = s.dot(g);
373 q = half * s.dot(dwa.dual()) + gs;
374 if (q < mu0*gs) {
375 search = true;
376 alphas = alpha;
377 }
378 }
379 else {
380 search = false;
381 }
382 }
383 alpha = alphas;
384 snorm = dgpstep(s,g,x,-alpha,model);
385 }
386 return snorm;
387 }
388
389 // Perform projected search to determine beta such that
390 // q(P(x + beta*s)-x) <= mu0*g'(P(x + beta*s)-x) for mu0 in (0,1)
391 // x -- The anchor vector x, upon return x = P(x + beta*s)
392 // s -- The direction vector s, upon return s = P(x + beta*s) - x
393 // g -- The free components of the gradient vector g (unchanged)
394 // model -- Contains objective and bound constraint information
395 // pwa -- Primal working array
396 // dwa -- Dual working array
398 const Vector<Real> &g, TrustRegionModel<Real> &model,
399 Vector<Real> &pwa, Vector<Real> &dwa) {
400 const Real half(0.5), one(1), mu0(0.01), interpf(0.5);
401 Real tol = std::sqrt(ROL_EPSILON<Real>());
402 Real beta(1), snorm(0), q(0), gs(0), bpmin(0), bpmax(0);
403 int nsteps = 0;
404 // Compute break points of x+beta*s;
405 dbreakpt(x,s,model,bpmin,bpmax,pwa);
406 // Reduce beta until sufficient decrease is satisfied
407 bool search = true;
408 while (search && beta > bpmin) {
409 nsteps++;
410 snorm = dgpstep(pwa,s,x,beta,model);
411 dynamic_cast<LinMoreModel<Real>&>(model).applyFreeHessian(dwa,pwa,x,tol);
412 gs = pwa.dot(g);
413 q = half * s.dot(dwa.dual()) + gs;
414 if (q <= mu0*gs) {
415 search = false;
416 }
417 else {
418 beta *= interpf;
419 }
420 }
421 if (beta < one && beta < bpmin) {
422 beta = bpmin;
423 }
424 snorm = dgpstep(pwa,s,x,beta,model);
425 s.set(pwa);
426 x.plus(s);
427 if (verbosity_ > 0) {
428 std::cout << std::endl;
429 std::cout << " Projected search" << std::endl;
430 std::cout << " Step length (beta): " << beta << std::endl;
431 }
432 return snorm;
433 }
434
435 // Compute sigma such that ||x+sigma*p||_inv(M) = del. This is called
436 // if dtrpcg detects negative curvature or if the step violates
437 // the trust region bound
438 // xtx -- The dot product <x, inv(M)x> (unchanged)
439 // ptp -- The dot product <p, inv(M)p> (unchanged)
440 // ptx -- The dot product <p, inv(M)x> (unchanged)
441 // del -- Trust region radius (unchanged)
442 Real dtrqsol(const Real xtx, const Real ptp, const Real ptx, const Real del) {
443 const Real zero(0);
444 Real dsq = del*del;
445 Real rad = ptx*ptx + ptp*(dsq-xtx);
446 rad = std::sqrt(std::max(rad,zero));
447 Real sigma(0);
448 if (ptx > zero) {
449 sigma = (dsq-xtx)/(ptx+rad);
450 }
451 else if (rad > zero) {
452 sigma = (rad-ptx)/ptp;
453 }
454 else {
455 sigma = zero;
456 }
457 return sigma;
458 }
459
460 // Solve the trust region subproblem: minimize q(w) subject to the
461 // trust region constraint ||w||_inv(M) <= del using the Steihaug-Toint
462 // Conjugate Gradients algorithm
463 // w -- The step vector to be computed
464 // iflag -- Termination flag
465 // iter -- Number of CG iterations
466 // del -- Trust region radius (unchanged)
467 // model -- Contains the objective and bound constraint information
468 // tol -- Residual stopping tolerance (unchanged)
469 // stol -- Preconditioned residual stopping tolerance (unchanged)
470 // itermax -- Maximum number of iterations
471 // p -- Primal working array that stores the CG step
472 // q -- Dual working array that stores the Hessian applied to p
473 // r -- Primal working array that stores the preconditioned residual
474 // t -- Dual working array that stores the residual
475 Real dtrpcg(Vector<Real> &w, int &iflag, int &iter,
476 const Vector<Real> &g, const Vector<Real> &x,
477 const Real del, TrustRegionModel<Real> &model,
478 const Real tol, const Real stol, const int itermax,
480 Vector<Real> &t) {
481 Real tol0 = std::sqrt(ROL_EPSILON<Real>());
482 const Real zero(0), one(1), two(2);
483 Real rho(0), tnorm(0), rnorm(0), rnorm0(0), kappa(0), beta(0), sigma(0), alpha(0), rtr(0);
484 Real sMs(0), pMp(0), sMp(0);
485 iter = 0; iflag = 0;
486 // Initialize step
487 w.zero();
488 // Compute residual
489 t.set(g); t.scale(-one);
490 // Preconditioned residual
491 dynamic_cast<LinMoreModel<Real>&>(model).applyFreePrecond(r,t,x,tol0);
492 rho = r.dot(t.dual());
493 rnorm0 = std::sqrt(rho);
494 if ( rnorm0 == zero ) {
495 return zero;
496 }
497 // Initialize direction
498 p.set(r);
499 pMp = rho;
500 // Iterate CG
501 for (iter = 0; iter < itermax; ++iter) {
502 // Apply Hessian to direction dir
503 dynamic_cast<LinMoreModel<Real>&>(model).applyFreeHessian(q,p,x,tol0);
504 // Compute sigma such that ||s+sigma*dir|| = del
505 kappa = p.dot(q.dual());
506 alpha = (kappa>zero) ? rho/kappa : zero;
507 sigma = dtrqsol(sMs,pMp,sMp,del);
508 // Check for negative curvature or if iterate exceeds trust region
509 if (kappa <= zero || alpha >= sigma) {
510 w.axpy(sigma,p);
511 iflag = (kappa<=zero) ? 2 : 3;
512 break;
513 }
514 // Update iterate and residuals
515 w.axpy(alpha,p);
516 t.axpy(-alpha,q);
517 dynamic_cast<LinMoreModel<Real>&>(model).applyFreePrecond(r,t,x,tol0);
518 // Exit if residual tolerance is met
519 rtr = r.dot(t.dual());
520 rnorm = std::sqrt(rtr);
521 tnorm = t.norm();
522 if (rnorm <= stol || tnorm <= tol) {
523 iflag = 0;
524 break;
525 }
526 // Compute p = r + beta * p
527 beta = rtr/rho;
528 p.scale(beta); p.plus(r);
529 rho = rtr;
530 // Update dot products
531 // sMs = <s, inv(M)s>
532 // sMp = <s, inv(M)p>
533 // pMp = <p, inv(M)p>
534 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
535 sMp = beta*(sMp + alpha*pMp);
536 pMp = rho + beta*beta*pMp;
537 }
538 // Check iteration count
539 if (iter == itermax) {
540 iflag = 1;
541 }
542 if (iflag != 1) {
543 iter++;
544 }
545 return w.norm();
546 }
547
548};
549
550}
551
552#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the interface to evaluate projected trust-region model functions from the Kelley-Sachs bound...
Real apply(const Real &x, const Real &y) const
void reduce(const volatile Real &input, Real volatile &output) const
Elementwise::EReductionType reductionType() const
void reduce(const Real &input, Real &output) const
void reduce(const Real &input, Real &output) const
void reduce(const volatile Real &input, Real volatile &output) const
Elementwise::EReductionType reductionType() const
Real apply(const Real &x, const Real &y) const
Provides interface for truncated CG trust-region subproblem solver.
Ptr< Vector< Real > > dwa1_
Ptr< Vector< Real > > pwa2_
Real dtrpcg(Vector< Real > &w, int &iflag, int &iter, const Vector< Real > &g, const Vector< Real > &x, const Real del, TrustRegionModel< Real > &model, const Real tol, const Real stol, const int itermax, Vector< Real > &p, Vector< Real > &q, Vector< Real > &r, Vector< Real > &t)
Real dcauchy(Vector< Real > &s, Real &alpha, const Vector< Real > &x, const Vector< Real > &g, const Real del, TrustRegionModel< Real > &model, Vector< Real > &pwa1, Vector< Real > &pwa2, Vector< Real > &dwa)
Real dtrqsol(const Real xtx, const Real ptp, const Real ptx, const Real del)
ROL::LinMore::UpperBreakPoint ubp_
unsigned verbosity_
void dbreakpt(const Vector< Real > &x, const Vector< Real > &s, TrustRegionModel< Real > &model, Real &bpmin, Real &bpmax, Vector< Real > &pwa)
Ptr< Vector< Real > > g_
LinMore(ROL::ParameterList &parlist)
void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, TrustRegionModel< Real > &model) const
Real dprsrch(Vector< Real > &x, Vector< Real > &s, const Vector< Real > &g, TrustRegionModel< Real > &model, Vector< Real > &pwa, Vector< Real > &dwa)
Ptr< Vector< Real > > x_
Ptr< Vector< Real > > s_
Ptr< Vector< Real > > dwa2_
void run(Vector< Real > &s, Real &snorm, int &iflag, int &iter, const Real del, TrustRegionModel< Real > &model)
Ptr< Vector< Real > > pwa1_
ROL::LinMore::PositiveMin pmin_
ROL::LinMore::LowerBreakPoint lbp_
ROL::LinMore::PositiveMax pmax_
Provides the interface to evaluate trust-region model functions.
virtual const Ptr< const Vector< Real > > getGradient(void) const
virtual const Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
virtual const Ptr< const Vector< Real > > getIterate(void) const
Provides interface for and implements trust-region subproblem solvers.
virtual void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
void setPredictedReduction(const Real pRed)
Defines the linear algebra or vector space interface.
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
virtual void scale(const Real alpha)=0
Compute where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual void plus(const Vector &x)=0
Compute , where .
virtual void zero()
Set to zero vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual Real reduce(const Elementwise::ReductionOp< Real > &r) const
virtual int dimension() const
Return dimension of the vector space.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
virtual Real dot(const Vector &x) const =0
Compute where .
constexpr auto dim