ROL
ROL_CompositeStep.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_COMPOSITESTEP_H
11 #define ROL_COMPOSITESTEP_H
12 
13 #include "ROL_Types.hpp"
14 #include "ROL_Step.hpp"
15 #include "ROL_LAPACK.hpp"
16 #include "ROL_LinearAlgebra.hpp"
17 #include <sstream>
18 #include <iomanip>
19 
26 namespace ROL {
27 
28 template <class Real>
29 class CompositeStep : public Step<Real> {
30 private:
31 
32  // Vectors used for cloning.
33  ROL::Ptr<Vector<Real> > xvec_;
34  ROL::Ptr<Vector<Real> > gvec_;
35  ROL::Ptr<Vector<Real> > cvec_;
36  ROL::Ptr<Vector<Real> > lvec_;
37 
38  // Diagnostic return flags for subalgorithms.
39  int flagCG_;
40  int flagAC_;
41  int iterCG_;
42 
43  // Stopping conditions.
44  int maxiterCG_;
45  int maxiterOSS_;
46  Real tolCG_;
47  Real tolOSS_;
48  bool tolOSSfixed_;
49 
50  // Tolerances and stopping conditions for subalgorithms.
51  Real lmhtol_;
52  Real qntol_;
53  Real pgtol_;
54  Real projtol_;
55  Real tangtol_;
56  Real tntmax_;
57 
58  // Trust-region parameters.
59  Real zeta_;
60  Real Delta_;
61  Real penalty_;
62  Real eta_;
63  bool useConHess_;
64 
65  Real ared_;
66  Real pred_;
67  Real snorm_;
68  Real nnorm_;
69  Real tnorm_;
70 
71  // Output flags.
72  bool infoQN_;
73  bool infoLM_;
74  bool infoTS_;
75  bool infoAC_;
76  bool infoLS_;
77  bool infoALL_;
78 
79  // Performance summary.
80  int totalIterCG_;
81  int totalProj_;
82  int totalNegCurv_;
83  int totalRef_;
84  int totalCallLS_;
85  int totalIterLS_;
86 
87  template <typename T> int sgn(T val) {
88  return (T(0) < val) - (val < T(0));
89  }
90 
91  void printInfoLS(const std::vector<Real> &res) const {
92  if (infoLS_) {
93  std::stringstream hist;
94  hist << std::scientific << std::setprecision(8);
95  hist << "\n Augmented System Solver:\n";
96  hist << " True Residual\n";
97  for (unsigned j=0; j<res.size(); j++) {
98  hist << " " << std::left << std::setw(14) << res[j] << "\n";
99  }
100  hist << "\n";
101  std::cout << hist.str();
102  }
103  }
104 
105  Real setTolOSS(const Real intol) const {
106  return tolOSSfixed_ ? tolOSS_ : intol;
107  }
108 
109 public:
110  using Step<Real>::initialize;
111  using Step<Real>::compute;
112  using Step<Real>::update;
113 
114  virtual ~CompositeStep() {}
115 
116  CompositeStep( ROL::ParameterList & parlist ) : Step<Real>() {
117  //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
118  flagCG_ = 0;
119  flagAC_ = 0;
120  iterCG_ = 0;
121 
122  ROL::ParameterList& steplist = parlist.sublist("Step").sublist("Composite Step");
123 
124  //maxiterOSS_ = steplist.sublist("Optimality System Solver").get("Iteration Limit", 50);
125  tolOSS_ = steplist.sublist("Optimality System Solver").get("Nominal Relative Tolerance", 1e-8);
126  tolOSSfixed_ = steplist.sublist("Optimality System Solver").get("Fix Tolerance", true);
127 
128  maxiterCG_ = steplist.sublist("Tangential Subproblem Solver").get("Iteration Limit", 20);
129  tolCG_ = steplist.sublist("Tangential Subproblem Solver").get("Relative Tolerance", 1e-2);
130  Delta_ = steplist.get("Initial Radius", 1e2);
131  useConHess_ = steplist.get("Use Constraint Hessian", true);
132 
133  int outLvl = steplist.get("Output Level", 0);
134 
135  lmhtol_ = tolOSS_;
136  qntol_ = tolOSS_;
137  pgtol_ = tolOSS_;
138  projtol_ = tolOSS_;
139  tangtol_ = tolOSS_;
140  tntmax_ = 2.0;
141 
142  zeta_ = 0.8;
143  penalty_ = 1.0;
144  eta_ = 1e-8;
145 
146  snorm_ = 0.0;
147  nnorm_ = 0.0;
148  tnorm_ = 0.0;
149 
150  infoALL_ = false;
151  if (outLvl > 0) {
152  infoALL_ = true;
153  }
154  infoQN_ = false;
155  infoLM_ = false;
156  infoTS_ = false;
157  infoAC_ = false;
158  infoLS_ = false;
159  infoQN_ = infoQN_ || infoALL_;
160  infoLM_ = infoLM_ || infoALL_;
161  infoTS_ = infoTS_ || infoALL_;
162  infoAC_ = infoAC_ || infoALL_;
163  infoLS_ = infoLS_ || infoALL_;
164 
165  totalIterCG_ = 0;
166  totalProj_ = 0;
167  totalNegCurv_ = 0;
168  totalRef_ = 0;
169  totalCallLS_ = 0;
170  totalIterLS_ = 0;
171  }
172 
175  void initialize( Vector<Real> &x, const Vector<Real> &g, Vector<Real> &l, const Vector<Real> &c,
176  Objective<Real> &obj, Constraint<Real> &con,
177  AlgorithmState<Real> &algo_state ) {
178  //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
179  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
180  state->descentVec = x.clone();
181  state->gradientVec = g.clone();
182  state->constraintVec = c.clone();
183 
184  xvec_ = x.clone();
185  gvec_ = g.clone();
186  lvec_ = l.clone();
187  cvec_ = c.clone();
188 
189  ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
190  ROL::Ptr<Vector<Real> > gl = gvec_->clone();
191 
192  algo_state.nfval = 0;
193  algo_state.ncval = 0;
194  algo_state.ngrad = 0;
195 
196  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
197 
198  // Update objective and constraint.
199  obj.update(x,true,algo_state.iter);
200  algo_state.value = obj.value(x, zerotol);
201  algo_state.nfval++;
202  con.update(x,true,algo_state.iter);
203  con.value(*cvec_, x, zerotol);
204  algo_state.cnorm = cvec_->norm();
205  algo_state.ncval++;
206  obj.gradient(*gvec_, x, zerotol);
207 
208  // Compute gradient of Lagrangian at new multiplier guess.
209  computeLagrangeMultiplier(l, x, *gvec_, con);
210  con.applyAdjointJacobian(*ajl, l, x, zerotol);
211  gl->set(*gvec_); gl->plus(*ajl);
212  algo_state.ngrad++;
213  algo_state.gnorm = gl->norm();
214  }
215 
218  void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
219  Objective<Real> &obj, Constraint<Real> &con,
220  AlgorithmState<Real> &algo_state ) {
221  //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
222  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
223  Real f = 0.0;
224  ROL::Ptr<Vector<Real> > n = xvec_->clone();
225  ROL::Ptr<Vector<Real> > c = cvec_->clone();
226  ROL::Ptr<Vector<Real> > t = xvec_->clone();
227  ROL::Ptr<Vector<Real> > tCP = xvec_->clone();
228  ROL::Ptr<Vector<Real> > g = gvec_->clone();
229  ROL::Ptr<Vector<Real> > gf = gvec_->clone();
230  ROL::Ptr<Vector<Real> > Wg = xvec_->clone();
231  ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
232 
233  Real f_new = 0.0;
234  ROL::Ptr<Vector<Real> > l_new = lvec_->clone();
235  ROL::Ptr<Vector<Real> > c_new = cvec_->clone();
236  ROL::Ptr<Vector<Real> > g_new = gvec_->clone();
237  ROL::Ptr<Vector<Real> > gf_new = gvec_->clone();
238 
239  // Evaluate objective ... should have been stored.
240  f = obj.value(x, zerotol);
241  algo_state.nfval++;
242  // Compute gradient of objective ... should have been stored.
243  obj.gradient(*gf, x, zerotol);
244  // Evaluate constraint ... should have been stored.
245  con.value(*c, x, zerotol);
246 
247  // Compute quasi-normal step.
248  computeQuasinormalStep(*n, *c, x, zeta_*Delta_, con);
249 
250  // Compute gradient of Lagrangian ... should have been stored.
251  con.applyAdjointJacobian(*ajl, l, x, zerotol);
252  g->set(*gf);
253  g->plus(*ajl);
254  algo_state.ngrad++;
255 
256  // Solve tangential subproblem.
257  solveTangentialSubproblem(*t, *tCP, *Wg, x, *g, *n, l, Delta_, obj, con);
258  totalIterCG_ += iterCG_;
259 
260  // Check acceptance of subproblem solutions, adjust merit function penalty parameter, ensure global convergence.
261  accept(s, *n, *t, f_new, *c_new, *gf_new, *l_new, *g_new, x, l, f, *gf, *c, *g, *tCP, *Wg, obj, con, algo_state);
262  }
263 
266  void update( Vector<Real> &x, Vector<Real> &l, const Vector<Real> &s,
267  Objective<Real> &obj, Constraint<Real> &con,
268  AlgorithmState<Real> &algo_state ) {
269  //ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
270 
271  Real zero(0);
272  Real one(1);
273  Real two(2);
274  Real seven(7);
275  Real half(0.5);
276  Real zp9(0.9);
277  Real zp8(0.8);
278  Real em12(1e-12);
279  Real zerotol = std::sqrt(ROL_EPSILON<Real>());//zero;
280  Real ratio(zero);
281 
282  ROL::Ptr<Vector<Real> > g = gvec_->clone();
283  ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
284  ROL::Ptr<Vector<Real> > gl = gvec_->clone();
285  ROL::Ptr<Vector<Real> > c = cvec_->clone();
286 
287  // Determine if the step gives sufficient reduction in the merit function,
288  // update the trust-region radius.
289  ratio = ared_/pred_;
290  if ((std::abs(ared_) < em12) && std::abs(pred_) < em12) {
291  ratio = one;
292  }
293  if (ratio >= eta_) {
294  x.plus(s);
295  if (ratio >= zp9) {
296  Delta_ = std::max(seven*snorm_, Delta_);
297  }
298  else if (ratio >= zp8) {
299  Delta_ = std::max(two*snorm_, Delta_);
300  }
301  obj.update(x,true,algo_state.iter);
302  con.update(x,true,algo_state.iter);
303  flagAC_ = 1;
304  }
305  else {
306  Delta_ = half*std::max(nnorm_, tnorm_);
307  obj.update(x,false,algo_state.iter);
308  con.update(x,false,algo_state.iter);
309  flagAC_ = 0;
310  } // if (ratio >= eta)
311 
312  Real val = obj.value(x, zerotol);
313  algo_state.nfval++;
314  obj.gradient(*g, x, zerotol);
315  computeLagrangeMultiplier(l, x, *g, con);
316  con.applyAdjointJacobian(*ajl, l, x, zerotol);
317  gl->set(*g); gl->plus(*ajl);
318  algo_state.ngrad++;
319  con.value(*c, x, zerotol);
320 
321  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
322  state->gradientVec->set(*gl);
323  state->constraintVec->set(*c);
324 
325  algo_state.value = val;
326  algo_state.gnorm = gl->norm();
327  algo_state.cnorm = c->norm();
328  algo_state.iter++;
329  algo_state.snorm = snorm_;
330 
331  // Update algorithm state
332  //(algo_state.iterateVec)->set(x);
333  }
334 
338  void compute( Vector<Real> &s, const Vector<Real> &x, Objective<Real> &obj,
339  BoundConstraint<Real> &con,
340  AlgorithmState<Real> &algo_state ) {}
341 
345  void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj,
346  BoundConstraint<Real> &con,
347  AlgorithmState<Real> &algo_state ) {}
348 
351  std::string printHeader( void ) const {
352  std::stringstream hist;
353  hist << " ";
354  hist << std::setw(6) << std::left << "iter";
355  hist << std::setw(15) << std::left << "fval";
356  hist << std::setw(15) << std::left << "cnorm";
357  hist << std::setw(15) << std::left << "gLnorm";
358  hist << std::setw(15) << std::left << "snorm";
359  hist << std::setw(10) << std::left << "delta";
360  hist << std::setw(10) << std::left << "nnorm";
361  hist << std::setw(10) << std::left << "tnorm";
362  hist << std::setw(8) << std::left << "#fval";
363  hist << std::setw(8) << std::left << "#grad";
364  hist << std::setw(8) << std::left << "iterCG";
365  hist << std::setw(8) << std::left << "flagCG";
366  hist << std::setw(8) << std::left << "accept";
367  hist << std::setw(8) << std::left << "linsys";
368  hist << "\n";
369  return hist.str();
370  }
371 
372  std::string printName( void ) const {
373  std::stringstream hist;
374  hist << "\n" << " Composite-step trust-region solver";
375  hist << "\n";
376  return hist.str();
377  }
378 
381  std::string print( AlgorithmState<Real> & algo_state, bool pHeader = false ) const {
382  //const ROL::Ptr<const StepState<Real> >& step_state = Step<Real>::getStepState();
383 
384  std::stringstream hist;
385  hist << std::scientific << std::setprecision(6);
386  if ( algo_state.iter == 0 ) {
387  hist << printName();
388  }
389  if ( pHeader ) {
390  hist << printHeader();
391  }
392  if ( algo_state.iter == 0 ) {
393  hist << " ";
394  hist << std::setw(6) << std::left << algo_state.iter;
395  hist << std::setw(15) << std::left << algo_state.value;
396  hist << std::setw(15) << std::left << algo_state.cnorm;
397  hist << std::setw(15) << std::left << algo_state.gnorm;
398  hist << "\n";
399  }
400  else {
401  hist << " ";
402  hist << std::setw(6) << std::left << algo_state.iter;
403  hist << std::setw(15) << std::left << algo_state.value;
404  hist << std::setw(15) << std::left << algo_state.cnorm;
405  hist << std::setw(15) << std::left << algo_state.gnorm;
406  hist << std::setw(15) << std::left << algo_state.snorm;
407  hist << std::scientific << std::setprecision(2);
408  hist << std::setw(10) << std::left << Delta_;
409  hist << std::setw(10) << std::left << nnorm_;
410  hist << std::setw(10) << std::left << tnorm_;
411  hist << std::scientific << std::setprecision(6);
412  hist << std::setw(8) << std::left << algo_state.nfval;
413  hist << std::setw(8) << std::left << algo_state.ngrad;
414  hist << std::setw(8) << std::left << iterCG_;
415  hist << std::setw(8) << std::left << flagCG_;
416  hist << std::setw(8) << std::left << flagAC_;
417  hist << std::left << totalCallLS_ << "/" << totalIterLS_;
418  hist << "\n";
419  }
420  return hist.str();
421  }
422 
434  void computeLagrangeMultiplier(Vector<Real> &l, const Vector<Real> &x, const Vector<Real> &gf, Constraint<Real> &con) {
435 
436  Real one(1);
437  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
438  std::vector<Real> augiters;
439 
440  if (infoLM_) {
441  std::stringstream hist;
442  hist << "\n Lagrange multiplier step\n";
443  std::cout << hist.str();
444  }
445 
446  /* Apply adjoint of constraint Jacobian to current multiplier. */
447  ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
448  con.applyAdjointJacobian(*ajl, l, x, zerotol);
449 
450  /* Form right-hand side of the augmented system. */
451  ROL::Ptr<Vector<Real> > b1 = gvec_->clone();
452  ROL::Ptr<Vector<Real> > b2 = cvec_->clone();
453  // b1 is the negative gradient of the Lagrangian
454  b1->set(gf); b1->plus(*ajl); b1->scale(-one);
455  // b2 is zero
456  b2->zero();
457 
458  /* Declare left-hand side of augmented system. */
459  ROL::Ptr<Vector<Real> > v1 = xvec_->clone();
460  ROL::Ptr<Vector<Real> > v2 = lvec_->clone();
461 
462  /* Compute linear solver tolerance. */
463  Real b1norm = b1->norm();
464  Real tol = setTolOSS(lmhtol_*b1norm);
465 
466  /* Solve augmented system. */
467  augiters = con.solveAugmentedSystem(*v1, *v2, *b1, *b2, x, tol);
468  totalCallLS_++;
469  totalIterLS_ = totalIterLS_ + augiters.size();
470  printInfoLS(augiters);
471 
472  /* Return updated Lagrange multiplier. */
473  // v2 is the multiplier update
474  l.plus(*v2);
475 
476  } // computeLagrangeMultiplier
477 
478 
501  void computeQuasinormalStep(Vector<Real> &n, const Vector<Real> &c, const Vector<Real> &x, Real delta, Constraint<Real> &con) {
502 
503  if (infoQN_) {
504  std::stringstream hist;
505  hist << "\n Quasi-normal step\n";
506  std::cout << hist.str();
507  }
508 
509  Real zero(0);
510  Real one(1);
511  Real zerotol = std::sqrt(ROL_EPSILON<Real>()); //zero;
512  std::vector<Real> augiters;
513 
514  /* Compute Cauchy step nCP. */
515  ROL::Ptr<Vector<Real> > nCP = xvec_->clone();
516  ROL::Ptr<Vector<Real> > nCPdual = gvec_->clone();
517  ROL::Ptr<Vector<Real> > nN = xvec_->clone();
518  ROL::Ptr<Vector<Real> > ctemp = cvec_->clone();
519  ROL::Ptr<Vector<Real> > dualc0 = lvec_->clone();
520  dualc0->set(c.dual());
521  con.applyAdjointJacobian(*nCPdual, *dualc0, x, zerotol);
522  nCP->set(nCPdual->dual());
523  con.applyJacobian(*ctemp, *nCP, x, zerotol);
524 
525  Real normsquare_ctemp = ctemp->dot(*ctemp);
526  if (normsquare_ctemp != zero) {
527  nCP->scale( -(nCP->dot(*nCP))/normsquare_ctemp );
528  }
529 
530  /* If the Cauchy step nCP is outside the trust region,
531  return the scaled Cauchy step. */
532  Real norm_nCP = nCP->norm();
533  if (norm_nCP >= delta) {
534  n.set(*nCP);
535  n.scale( delta/norm_nCP );
536  if (infoQN_) {
537  std::stringstream hist;
538  hist << " taking partial Cauchy step\n";
539  std::cout << hist.str();
540  }
541  return;
542  }
543 
544  /* Compute 'Newton' step, for example, by solving a problem
545  related to finding the minimum norm solution of min || c(x_k)*s + c ||^2. */
546  // Compute tolerance for linear solver.
547  con.applyJacobian(*ctemp, *nCP, x, zerotol);
548  ctemp->plus(c);
549  Real tol = setTolOSS(qntol_*ctemp->norm());
550  // Form right-hand side.
551  ctemp->scale(-one);
552  nCPdual->set(nCP->dual());
553  nCPdual->scale(-one);
554  // Declare left-hand side of augmented system.
555  ROL::Ptr<Vector<Real> > dn = xvec_->clone();
556  ROL::Ptr<Vector<Real> > y = lvec_->clone();
557  // Solve augmented system.
558  augiters = con.solveAugmentedSystem(*dn, *y, *nCPdual, *ctemp, x, tol);
559  totalCallLS_++;
560  totalIterLS_ = totalIterLS_ + augiters.size();
561  printInfoLS(augiters);
562 
563  nN->set(*dn);
564  nN->plus(*nCP);
565 
566  /* Either take full or partial Newton step, depending on
567  the trust-region constraint. */
568  Real norm_nN = nN->norm();
569  if (norm_nN <= delta) {
570  // Take full feasibility step.
571  n.set(*nN);
572  if (infoQN_) {
573  std::stringstream hist;
574  hist << " taking full Newton step\n";
575  std::cout << hist.str();
576  }
577  }
578  else {
579  // Take convex combination n = nCP+tau*(nN-nCP),
580  // so that ||n|| = delta. In other words, solve
581  // scalar quadratic equation: ||nCP+tau*(nN-nCP)||^2 = delta^2.
582  Real aa = dn->dot(*dn);
583  Real bb = dn->dot(*nCP);
584  Real cc = norm_nCP*norm_nCP - delta*delta;
585  Real tau = (-bb+sqrt(bb*bb-aa*cc))/aa;
586  n.set(*nCP);
587  n.axpy(tau, *dn);
588  if (infoQN_) {
589  std::stringstream hist;
590  hist << " taking dogleg step\n";
591  std::cout << hist.str();
592  }
593  }
594 
595  } // computeQuasinormalStep
596 
597 
611  void solveTangentialSubproblem(Vector<Real> &t, Vector<Real> &tCP, Vector<Real> &Wg,
612  const Vector<Real> &x, const Vector<Real> &g, const Vector<Real> &n, const Vector<Real> &l,
613  Real delta, Objective<Real> &obj, Constraint<Real> &con) {
614 
615  /* Initialization of the CG step. */
616  bool orthocheck = true; // set to true if want to check orthogonality
617  // of Wr and r, otherwise set to false
618  Real tol_ortho = 0.5; // orthogonality measure; represets a bound on norm( \hat{S}, 2), where
619  // \hat{S} is defined in Heinkenschloss/Ridzal., "A Matrix-Free Trust-Region SQP Method"
620  Real S_max = 1.0; // another orthogonality measure; norm(S) needs to be bounded by
621  // a modest constant; norm(S) is small if the approximation of
622  // the null space projector is good
623  Real zero(0);
624  Real one(1);
625  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
626  std::vector<Real> augiters;
627  iterCG_ = 1;
628  flagCG_ = 0;
629  t.zero();
630  tCP.zero();
631  ROL::Ptr<Vector<Real> > r = gvec_->clone();
632  ROL::Ptr<Vector<Real> > pdesc = xvec_->clone();
633  ROL::Ptr<Vector<Real> > tprev = xvec_->clone();
634  ROL::Ptr<Vector<Real> > Wr = xvec_->clone();
635  ROL::Ptr<Vector<Real> > Hp = gvec_->clone();
636  ROL::Ptr<Vector<Real> > xtemp = xvec_->clone();
637  ROL::Ptr<Vector<Real> > gtemp = gvec_->clone();
638  ROL::Ptr<Vector<Real> > ltemp = lvec_->clone();
639  ROL::Ptr<Vector<Real> > czero = cvec_->clone();
640  czero->zero();
641  r->set(g);
642  obj.hessVec(*gtemp, n, x, zerotol);
643  r->plus(*gtemp);
644  if (useConHess_) {
645  con.applyAdjointHessian(*gtemp, l, n, x, zerotol);
646  r->plus(*gtemp);
647  }
648  Real normg = r->norm();
649  Real normWg = zero;
650  Real pHp = zero;
651  Real rp = zero;
652  Real alpha = zero;
653  Real normp = zero;
654  Real normr = zero;
655  Real normt = zero;
656  std::vector<Real> normWr(maxiterCG_+1, zero);
657 
658  std::vector<ROL::Ptr<Vector<Real > > > p; // stores search directions
659  std::vector<ROL::Ptr<Vector<Real > > > Hps; // stores duals of hessvec's applied to p's
660  std::vector<ROL::Ptr<Vector<Real > > > rs; // stores duals of residuals
661  std::vector<ROL::Ptr<Vector<Real > > > Wrs; // stores duals of projected residuals
662 
663  Real rptol(1e-12);
664 
665  if (infoTS_) {
666  std::stringstream hist;
667  hist << "\n Tangential subproblem\n";
668  hist << std::setw(6) << std::right << "iter" << std::setw(18) << "||Wr||/||Wr0||" << std::setw(15) << "||s||";
669  hist << std::setw(15) << "delta" << std::setw(15) << "||c'(x)s||" << "\n";
670  std::cout << hist.str();
671  }
672 
673  if (normg == 0) {
674  if (infoTS_) {
675  std::stringstream hist;
676  hist << " >>> Tangential subproblem: Initial gradient is zero! \n";
677  std::cout << hist.str();
678  }
679  iterCG_ = 0; Wg.zero(); flagCG_ = 0;
680  return;
681  }
682 
683  /* Start CG loop. */
684  while (iterCG_ < maxiterCG_) {
685 
686  // Store tangential Cauchy point (which is the current iterate in the second iteration).
687  if (iterCG_ == 2) {
688  tCP.set(t);
689  }
690 
691  // Compute (inexact) projection W*r.
692  if (iterCG_ == 1) {
693  // Solve augmented system.
694  Real tol = setTolOSS(pgtol_);
695  augiters = con.solveAugmentedSystem(*Wr, *ltemp, *r, *czero, x, tol);
696  totalCallLS_++;
697  totalIterLS_ = totalIterLS_ + augiters.size();
698  printInfoLS(augiters);
699 
700  Wg.set(*Wr);
701  normWg = Wg.norm();
702  if (orthocheck) {
703  Wrs.push_back(xvec_->clone());
704  (Wrs[iterCG_-1])->set(*Wr);
705  }
706  // Check if done (small initial projected residual).
707  if (normWg == zero) {
708  flagCG_ = 0;
709  iterCG_--;
710  if (infoTS_) {
711  std::stringstream hist;
712  hist << " Initial projected residual is close to zero! \n";
713  std::cout << hist.str();
714  }
715  return;
716  }
717  // Set first residual to projected gradient.
718  // change r->set(Wg);
719  r->set(Wg.dual());
720  if (orthocheck) {
721  rs.push_back(xvec_->clone());
722  // change (rs[0])->set(*r);
723  (rs[0])->set(r->dual());
724  }
725  }
726  else {
727  // Solve augmented system.
728  Real tol = setTolOSS(projtol_);
729  augiters = con.solveAugmentedSystem(*Wr, *ltemp, *r, *czero, x, tol);
730  totalCallLS_++;
731  totalIterLS_ = totalIterLS_ + augiters.size();
732  printInfoLS(augiters);
733 
734  if (orthocheck) {
735  Wrs.push_back(xvec_->clone());
736  (Wrs[iterCG_-1])->set(*Wr);
737  }
738  }
739 
740  normWr[iterCG_-1] = Wr->norm();
741 
742  if (infoTS_) {
743  ROL::Ptr<Vector<Real> > ct = cvec_->clone();
744  con.applyJacobian(*ct, t, x, zerotol);
745  Real linc = ct->norm();
746  std::stringstream hist;
747  hist << std::scientific << std::setprecision(6);
748  hist << std::setw(6) << std::right << iterCG_-1 << std::setw(18) << normWr[iterCG_-1]/normWg << std::setw(15) << t.norm();
749  hist << std::setw(15) << delta << std::setw(15) << linc << "\n";
750  std::cout << hist.str();
751  }
752 
753  // Check if done (small relative residual).
754  if (normWr[iterCG_-1]/normWg < tolCG_) {
755  flagCG_ = 0;
756  iterCG_ = iterCG_-1;
757  if (infoTS_) {
758  std::stringstream hist;
759  hist << " || W(g + H*(n+s)) || <= cgtol*|| W(g + H*n)|| \n";
760  std::cout << hist.str();
761  }
762  return;
763  }
764 
765  // Check nonorthogonality, one-norm of (WR*R/diag^2 - I)
766  if (orthocheck) {
767  ROL::LA::Matrix<Real> Wrr(iterCG_,iterCG_); // holds matrix Wrs'*rs
768  ROL::LA::Matrix<Real> T(iterCG_,iterCG_); // holds matrix T=(1/diag)*Wrs'*rs*(1/diag)
769  ROL::LA::Matrix<Real> Tm1(iterCG_,iterCG_); // holds matrix Tm1=T-I
770  for (int i=0; i<iterCG_; i++) {
771  for (int j=0; j<iterCG_; j++) {
772  Wrr(i,j) = (Wrs[i])->dot(*rs[j]);
773  T(i,j) = Wrr(i,j)/(normWr[i]*normWr[j]);
774  Tm1(i,j) = T(i,j);
775  if (i==j) {
776  Tm1(i,j) = Tm1(i,j) - one;
777  }
778  }
779  }
780  if (Tm1.normOne() >= tol_ortho) {
781  ROL::LAPACK<int,Real> lapack;
782  std::vector<int> ipiv(iterCG_);
783  int info;
784  std::vector<Real> work(3*iterCG_);
785  // compute inverse of T
786  lapack.GETRF(iterCG_, iterCG_, T.values(), T.stride(), &ipiv[0], &info);
787  lapack.GETRI(iterCG_, T.values(), T.stride(), &ipiv[0], &work[0], 3*iterCG_, &info);
788  Tm1 = T;
789  for (int i=0; i<iterCG_; i++) {
790  Tm1(i,i) = Tm1(i,i) - one;
791  }
792  if (Tm1.normOne() > S_max) {
793  flagCG_ = 4;
794  if (infoTS_) {
795  std::stringstream hist;
796  hist << " large nonorthogonality in W(R)'*R detected \n";
797  std::cout << hist.str();
798  }
799  return;
800  }
801  }
802  }
803 
804  // Full orthogonalization.
805  p.push_back(xvec_->clone());
806  (p[iterCG_-1])->set(*Wr);
807  (p[iterCG_-1])->scale(-one);
808  for (int j=1; j<iterCG_; j++) {
809  Real scal = (p[iterCG_-1])->dot(*(Hps[j-1])) / (p[j-1])->dot(*(Hps[j-1]));
810  ROL::Ptr<Vector<Real> > pj = xvec_->clone();
811  pj->set(*p[j-1]);
812  pj->scale(-scal);
813  (p[iterCG_-1])->plus(*pj);
814  }
815 
816  // change Hps.push_back(gvec_->clone());
817  Hps.push_back(xvec_->clone());
818  // change obj.hessVec(*(Hps[iterCG_-1]), *(p[iterCG_-1]), x, zerotol);
819  obj.hessVec(*Hp, *(p[iterCG_-1]), x, zerotol);
820  if (useConHess_) {
821  con.applyAdjointHessian(*gtemp, l, *(p[iterCG_-1]), x, zerotol);
822  // change (Hps[iterCG_-1])->plus(*gtemp);
823  Hp->plus(*gtemp);
824  }
825  // "Preconditioning" step.
826  (Hps[iterCG_-1])->set(Hp->dual());
827 
828  pHp = (p[iterCG_-1])->dot(*(Hps[iterCG_-1]));
829  // change rp = (p[iterCG_-1])->dot(*r);
830  rp = (p[iterCG_-1])->dot(*(rs[iterCG_-1]));
831 
832  normp = (p[iterCG_-1])->norm();
833  normr = r->norm();
834 
835  // Negative curvature stopping condition.
836  if (pHp <= 0) {
837  pdesc->set(*(p[iterCG_-1])); // p is the descent direction
838  if ((std::abs(rp) >= rptol*normp*normr) && (sgn(rp) == 1)) {
839  pdesc->scale(-one); // -p is the descent direction
840  }
841  flagCG_ = 2;
842  Real a = pdesc->dot(*pdesc);
843  Real b = pdesc->dot(t);
844  Real c = t.dot(t) - delta*delta;
845  // Positive root of a*theta^2 + 2*b*theta + c = 0.
846  Real theta = (-b + std::sqrt(b*b - a*c)) / a;
847  xtemp->set(*(p[iterCG_-1]));
848  xtemp->scale(theta);
849  t.plus(*xtemp);
850  // Store as tangential Cauchy point if terminating in first iteration.
851  if (iterCG_ == 1) {
852  tCP.set(t);
853  }
854  if (infoTS_) {
855  std::stringstream hist;
856  hist << " negative curvature detected \n";
857  std::cout << hist.str();
858  }
859  return;
860  }
861 
862  // Want to enforce nonzero alpha's.
863  if (std::abs(rp) < rptol*normp*normr) {
864  flagCG_ = 5;
865  if (infoTS_) {
866  std::stringstream hist;
867  hist << " Zero alpha due to inexactness. \n";
868  std::cout << hist.str();
869  }
870  return;
871  }
872 
873  alpha = - rp/pHp;
874 
875  // Iterate update.
876  tprev->set(t);
877  xtemp->set(*(p[iterCG_-1]));
878  xtemp->scale(alpha);
879  t.plus(*xtemp);
880 
881  // Trust-region stopping condition.
882  normt = t.norm();
883  if (normt >= delta) {
884  pdesc->set(*(p[iterCG_-1])); // p is the descent direction
885  if (sgn(rp) == 1) {
886  pdesc->scale(-one); // -p is the descent direction
887  }
888  Real a = pdesc->dot(*pdesc);
889  Real b = pdesc->dot(*tprev);
890  Real c = tprev->dot(*tprev) - delta*delta;
891  // Positive root of a*theta^2 + 2*b*theta + c = 0.
892  Real theta = (-b + std::sqrt(b*b - a*c)) / a;
893  xtemp->set(*(p[iterCG_-1]));
894  xtemp->scale(theta);
895  t.set(*tprev);
896  t.plus(*xtemp);
897  // Store as tangential Cauchy point if terminating in first iteration.
898  if (iterCG_ == 1) {
899  tCP.set(t);
900  }
901  flagCG_ = 3;
902  if (infoTS_) {
903  std::stringstream hist;
904  hist << " trust-region condition active \n";
905  std::cout << hist.str();
906  }
907  return;
908  }
909 
910  // Residual update.
911  xtemp->set(*(Hps[iterCG_-1]));
912  xtemp->scale(alpha);
913  // change r->plus(*gtemp);
914  r->plus(xtemp->dual());
915  if (orthocheck) {
916  // change rs.push_back(gvec_->clone());
917  rs.push_back(xvec_->clone());
918  // change (rs[iterCG_])->set(*r);
919  (rs[iterCG_])->set(r->dual());
920  }
921 
922  iterCG_++;
923 
924  } // while (iterCG_ < maxiterCG_)
925 
926  flagCG_ = 1;
927  if (infoTS_) {
928  std::stringstream hist;
929  hist << " maximum number of iterations reached \n";
930  std::cout << hist.str();
931  }
932 
933  } // solveTangentialSubproblem
934 
935 
938  void accept(Vector<Real> &s, Vector<Real> &n, Vector<Real> &t, Real f_new, Vector<Real> &c_new,
939  Vector<Real> &gf_new, Vector<Real> &l_new, Vector<Real> &g_new,
940  const Vector<Real> &x, const Vector<Real> &l, Real f, const Vector<Real> &gf, const Vector<Real> &c,
941  const Vector<Real> &g, Vector<Real> &tCP, Vector<Real> &Wg,
942  Objective<Real> &obj, Constraint<Real> &con, AlgorithmState<Real> &algo_state) {
943 
944  Real beta = 1e-8; // predicted reduction parameter
945  Real tol_red_tang = 1e-3; // internal reduction factor for tangtol
946  Real tol_red_all = 1e-1; // internal reduction factor for qntol, lmhtol, pgtol, projtol, tangtol
947  //bool glob_refine = true; // true - if subsolver tolerances are adjusted in this routine, keep adjusted values globally
948  // false - if subsolver tolerances are adjusted in this routine, discard adjusted values
949  Real tol_fdiff = 1e-12; // relative objective function difference for ared computation
950  int ct_max = 10; // maximum number of globalization tries
951  Real mintol = 1e-16; // smallest projection tolerance value
952 
953  // Determines max value of |rpred|/pred.
954  Real rpred_over_pred = 0.5*(1-eta_);
955 
956  if (infoAC_) {
957  std::stringstream hist;
958  hist << "\n Composite step acceptance\n";
959  std::cout << hist.str();
960  }
961 
962  Real zero = 0.0;
963  Real one = 1.0;
964  Real two = 2.0;
965  Real half = one/two;
966  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
967  std::vector<Real> augiters;
968 
969  Real pred = zero;
970  Real ared = zero;
971  Real rpred = zero;
972  Real part_pred = zero;
973  Real linc_preproj = zero;
974  Real linc_postproj = zero;
975  Real tangtol_start = zero;
976  Real tangtol = tangtol_;
977  //Real projtol = projtol_;
978  bool flag = false;
979  int num_proj = 0;
980  bool try_tCP = false;
981  Real fdiff = zero;
982 
983  ROL::Ptr<Vector<Real> > xtrial = xvec_->clone();
984  ROL::Ptr<Vector<Real> > Jl = gvec_->clone();
985  ROL::Ptr<Vector<Real> > gfJl = gvec_->clone();
986  ROL::Ptr<Vector<Real> > Jnc = cvec_->clone();
987  ROL::Ptr<Vector<Real> > t_orig = xvec_->clone();
988  ROL::Ptr<Vector<Real> > t_dual = gvec_->clone();
989  ROL::Ptr<Vector<Real> > Jt_orig = cvec_->clone();
990  ROL::Ptr<Vector<Real> > t_m_tCP = xvec_->clone();
991  ROL::Ptr<Vector<Real> > ltemp = lvec_->clone();
992  ROL::Ptr<Vector<Real> > xtemp = xvec_->clone();
993  ROL::Ptr<Vector<Real> > rt = cvec_->clone();
994  ROL::Ptr<Vector<Real> > Hn = gvec_->clone();
995  ROL::Ptr<Vector<Real> > Hto = gvec_->clone();
996  ROL::Ptr<Vector<Real> > cxxvec = gvec_->clone();
997  ROL::Ptr<Vector<Real> > czero = cvec_->clone();
998  czero->zero();
999  Real Jnc_normsquared = zero;
1000  Real c_normsquared = zero;
1001 
1002  // Compute and store some quantities for later use. Necessary
1003  // because of the function and constraint updates below.
1004  con.applyAdjointJacobian(*Jl, l, x, zerotol);
1005  con.applyJacobian(*Jnc, n, x, zerotol);
1006  Jnc->plus(c);
1007  Jnc_normsquared = Jnc->dot(*Jnc);
1008  c_normsquared = c.dot(c);
1009 
1010  for (int ct=0; ct<ct_max; ct++) {
1011 
1012  try_tCP = true;
1013  t_m_tCP->set(t);
1014  t_m_tCP->scale(-one);
1015  t_m_tCP->plus(tCP);
1016  if (t_m_tCP->norm() == zero) {
1017  try_tCP = false;
1018  }
1019 
1020  t_orig->set(t);
1021  con.applyJacobian(*Jt_orig, *t_orig, x, zerotol);
1022  linc_preproj = Jt_orig->norm();
1023  pred = one;
1024  rpred = two*rpred_over_pred*pred;
1025  flag = false;
1026  num_proj = 1;
1027  tangtol_start = tangtol;
1028 
1029  while (std::abs(rpred)/pred > rpred_over_pred) {
1030  // Compute projected tangential step.
1031  if (flag) {
1032  tangtol = tol_red_tang*tangtol;
1033  num_proj++;
1034  if (tangtol < mintol) {
1035  if (infoAC_) {
1036  std::stringstream hist;
1037  hist << "\n The projection of the tangential step cannot be done with sufficient precision.\n";
1038  hist << " Is the quasi-normal step very small? Continuing with no global convergence guarantees.\n";
1039  std::cout << hist.str();
1040  }
1041  break;
1042  }
1043  }
1044  // Solve augmented system.
1045  Real tol = setTolOSS(tangtol);
1046  // change augiters = con.solveAugmentedSystem(t, *ltemp, *t_orig, *czero, x, tol);
1047  t_dual->set(t_orig->dual());
1048  augiters = con.solveAugmentedSystem(t, *ltemp, *t_dual, *czero, x, tol);
1049  totalCallLS_++;
1050  totalIterLS_ = totalIterLS_ + augiters.size();
1051  printInfoLS(augiters);
1052  totalProj_++;
1053  con.applyJacobian(*rt, t, x, zerotol);
1054  linc_postproj = rt->norm();
1055 
1056  // Compute composite step.
1057  s.set(t);
1058  s.plus(n);
1059 
1060  // Compute some quantities before updating the objective and the constraint.
1061  obj.hessVec(*Hn, n, x, zerotol);
1062  if (useConHess_) {
1063  con.applyAdjointHessian(*cxxvec, l, n, x, zerotol);
1064  Hn->plus(*cxxvec);
1065  }
1066  obj.hessVec(*Hto, *t_orig, x, zerotol);
1067  if (useConHess_) {
1068  con.applyAdjointHessian(*cxxvec, l, *t_orig, x, zerotol);
1069  Hto->plus(*cxxvec);
1070  }
1071 
1072  // Compute objective, constraint, etc. values at the trial point.
1073  xtrial->set(x);
1074  xtrial->plus(s);
1075  obj.update(*xtrial,false,algo_state.iter);
1076  con.update(*xtrial,false,algo_state.iter);
1077  f_new = obj.value(*xtrial, zerotol);
1078  obj.gradient(gf_new, *xtrial, zerotol);
1079  con.value(c_new, *xtrial, zerotol);
1080  l_new.set(l);
1081  computeLagrangeMultiplier(l_new, *xtrial, gf_new, con);
1082 
1083  // Penalty parameter update.
1084  part_pred = - Wg.dot(*t_orig);
1085  gfJl->set(gf);
1086  gfJl->plus(*Jl);
1087  // change part_pred -= gfJl->dot(n);
1088  part_pred -= n.dot(gfJl->dual());
1089  // change part_pred -= half*Hn->dot(n);
1090  part_pred -= half*n.dot(Hn->dual());
1091  // change part_pred -= half*Hto->dot(*t_orig);
1092  part_pred -= half*t_orig->dot(Hto->dual());
1093  ltemp->set(l_new);
1094  ltemp->axpy(-one, l);
1095  // change part_pred -= Jnc->dot(*ltemp);
1096  part_pred -= Jnc->dot(ltemp->dual());
1097 
1098  if ( part_pred < -half*penalty_*(c_normsquared-Jnc_normsquared) ) {
1099  penalty_ = ( -two * part_pred / (c_normsquared-Jnc_normsquared) ) + beta;
1100  }
1101 
1102  pred = part_pred + penalty_*(c_normsquared-Jnc_normsquared);
1103 
1104  // Computation of rpred.
1105  // change rpred = - ltemp->dot(*rt) - penalty_ * rt->dot(*rt) - two * penalty_ * rt->dot(*Jnc);
1106  rpred = - rt->dot(ltemp->dual()) - penalty_ * rt->dot(*rt) - two * penalty_ * rt->dot(*Jnc);
1107  // change ROL::Ptr<Vector<Real> > lrt = lvec_->clone();
1108  //lrt->set(*rt);
1109  //rpred = - ltemp->dot(*rt) - penalty_ * std::pow(rt->norm(), 2) - two * penalty_ * lrt->dot(*Jnc);
1110  flag = 1;
1111 
1112  } // while (std::abs(rpred)/pred > rpred_over_pred)
1113 
1114  tangtol = tangtol_start;
1115 
1116  // Check if the solution of the tangential subproblem is
1117  // disproportionally large compared to total trial step.
1118  xtemp->set(n);
1119  xtemp->plus(t);
1120  if ( t_orig->norm()/xtemp->norm() < tntmax_ ) {
1121  break;
1122  }
1123  else {
1124  t_m_tCP->set(*t_orig);
1125  t_m_tCP->scale(-one);
1126  t_m_tCP->plus(tCP);
1127  if ((t_m_tCP->norm() > 0) && try_tCP) {
1128  if (infoAC_) {
1129  std::stringstream hist;
1130  hist << " ---> now trying tangential Cauchy point\n";
1131  std::cout << hist.str();
1132  }
1133  t.set(tCP);
1134  }
1135  else {
1136  if (infoAC_) {
1137  std::stringstream hist;
1138  hist << " ---> recomputing quasi-normal step and re-solving tangential subproblem\n";
1139  std::cout << hist.str();
1140  }
1141  totalRef_++;
1142  // Reset global quantities.
1143  obj.update(x);
1144  con.update(x);
1145  /*lmhtol = tol_red_all*lmhtol;
1146  qntol = tol_red_all*qntol;
1147  pgtol = tol_red_all*pgtol;
1148  projtol = tol_red_all*projtol;
1149  tangtol = tol_red_all*tangtol;
1150  if (glob_refine) {
1151  lmhtol_ = lmhtol;
1152  qntol_ = qntol;
1153  pgtol_ = pgtol;
1154  projtol_ = projtol;
1155  tangtol_ = tangtol;
1156  }*/
1157  if (!tolOSSfixed_) {
1158  lmhtol_ *= tol_red_all;
1159  qntol_ *= tol_red_all;
1160  pgtol_ *= tol_red_all;
1161  projtol_ *= tol_red_all;
1162  tangtol_ *= tol_red_all;
1163  }
1164  // Recompute the quasi-normal step.
1165  computeQuasinormalStep(n, c, x, zeta_*Delta_, con);
1166  // Solve tangential subproblem.
1167  solveTangentialSubproblem(t, tCP, Wg, x, g, n, l, Delta_, obj, con);
1168  totalIterCG_ += iterCG_;
1169  if (flagCG_ == 1) {
1170  totalNegCurv_++;
1171  }
1172  }
1173  } // else w.r.t. ( t_orig->norm()/xtemp->norm() < tntmax )
1174 
1175  } // for (int ct=0; ct<ct_max; ct++)
1176 
1177  // Compute actual reduction;
1178  fdiff = f - f_new;
1179  // Heuristic 1: If fdiff is very small compared to f, set it to 0,
1180  // in order to prevent machine precision issues.
1181  Real em24(1e-24);
1182  Real em14(1e-14);
1183  if (std::abs(fdiff / (f+em24)) < tol_fdiff) {
1184  fdiff = em14;
1185  }
1186  // change ared = fdiff + (l.dot(c) - l_new.dot(c_new)) + penalty_*(c.dot(c) - c_new.dot(c_new));
1187  // change ared = fdiff + (l.dot(c) - l_new.dot(c_new)) + penalty_*(std::pow(c.norm(),2) - std::pow(c_new.norm(),2));
1188  ared = fdiff + (c.dot(l.dual()) - c_new.dot(l_new.dual())) + penalty_*(c.dot(c) - c_new.dot(c_new));
1189 
1190  // Store actual and predicted reduction.
1191  ared_ = ared;
1192  pred_ = pred;
1193 
1194  // Store step and vector norms.
1195  snorm_ = s.norm();
1196  nnorm_ = n.norm();
1197  tnorm_ = t.norm();
1198 
1199  // Print diagnostics.
1200  if (infoAC_) {
1201  std::stringstream hist;
1202  hist << "\n Trial step info ...\n";
1203  hist << " n_norm = " << nnorm_ << "\n";
1204  hist << " t_norm = " << tnorm_ << "\n";
1205  hist << " s_norm = " << snorm_ << "\n";
1206  hist << " xtrial_norm = " << xtrial->norm() << "\n";
1207  hist << " f_old = " << f << "\n";
1208  hist << " f_trial = " << f_new << "\n";
1209  hist << " f_old-f_trial = " << f-f_new << "\n";
1210  hist << " ||c_old|| = " << c.norm() << "\n";
1211  hist << " ||c_trial|| = " << c_new.norm() << "\n";
1212  hist << " ||Jac*t_preproj|| = " << linc_preproj << "\n";
1213  hist << " ||Jac*t_postproj|| = " << linc_postproj << "\n";
1214  hist << " ||t_tilde||/||t|| = " << t_orig->norm() / t.norm() << "\n";
1215  hist << " ||t_tilde||/||n+t|| = " << t_orig->norm() / snorm_ << "\n";
1216  hist << " # projections = " << num_proj << "\n";
1217  hist << " penalty param = " << penalty_ << "\n";
1218  hist << " ared = " << ared_ << "\n";
1219  hist << " pred = " << pred_ << "\n";
1220  hist << " ared/pred = " << ared_/pred_ << "\n";
1221  std::cout << hist.str();
1222  }
1223 
1224  } // accept
1225 
1226 }; // class CompositeStep
1227 
1228 } // namespace ROL
1229 
1230 #endif
ROL::Objective
Provides the interface to evaluate objective functions.
Definition: ROL_Objective.hpp:44
ROL::CompositeStep::cvec_
ROL::Ptr< Vector< Real > > cvec_
Definition: ROL_CompositeStep.hpp:35
ROL::CompositeStep::compute
void compute(Vector< Real > &s, const Vector< Real > &x, const Vector< Real > &l, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step.
Definition: ROL_CompositeStep.hpp:218
ROL::Vector::dual
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
ROL::Constraint::applyJacobian
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
Definition: ROL_ConstraintDef.hpp:19
ROL::CompositeStep::printHeader
std::string printHeader(void) const
Print iterate header.
Definition: ROL_CompositeStep.hpp:351
ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::CompositeStep::setTolOSS
Real setTolOSS(const Real intol) const
Definition: ROL_CompositeStep.hpp:105
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Constraint::update
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update constraint function.
Definition: ROL_Constraint.hpp:68
ROL::CompositeStep::accept
void accept(Vector< Real > &s, Vector< Real > &n, Vector< Real > &t, Real f_new, Vector< Real > &c_new, Vector< Real > &gf_new, Vector< Real > &l_new, Vector< Real > &g_new, const Vector< Real > &x, const Vector< Real > &l, Real f, const Vector< Real > &gf, const Vector< Real > &c, const Vector< Real > &g, Vector< Real > &tCP, Vector< Real > &Wg, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Check acceptance of subproblem solutions, adjust merit function penalty parameter, ensure global convergence.
Definition: ROL_CompositeStep.hpp:938
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
ROL::CompositeStep::compute
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step for bound constraints; here only to satisfy the interface requirements, does nothing, needs refactoring.
Definition: ROL_CompositeStep.hpp:338
ROL::CompositeStep::xvec_
ROL::Ptr< Vector< Real > > xvec_
Definition: ROL_CompositeStep.hpp:33
ROL::CompositeStep::gvec_
ROL::Ptr< Vector< Real > > gvec_
Definition: ROL_CompositeStep.hpp:34
ROL::Objective::value
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
ROL::Step
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:34
ROL::Objective::hessVec
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Definition: ROL_ObjectiveDef.hpp:61
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL::CompositeStep::initialize
void initialize(Vector< Real > &x, const Vector< Real > &g, Vector< Real > &l, const Vector< Real > &c, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Initialize step.
Definition: ROL_CompositeStep.hpp:175
ROL::CompositeStep::CompositeStep
CompositeStep(ROL::ParameterList &parlist)
Definition: ROL_CompositeStep.hpp:116
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL::Constraint::value
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)=0
Evaluate the constraint operator at .
ROL::Objective::update
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Definition: ROL_Objective.hpp:62
ROL::Vector::dot
virtual Real dot(const Vector &x) const =0
Compute where .
ROL::CompositeStep::update
void update(Vector< Real > &x, Vector< Real > &l, const Vector< Real > &s, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, if successful.
Definition: ROL_CompositeStep.hpp:266
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:111
ROL::CompositeStep::print
std::string print(AlgorithmState< Real > &algo_state, bool pHeader=false) const
Print iterate status.
Definition: ROL_CompositeStep.hpp:381
ROL::AlgorithmState
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:109
ROL::Objective::gradient
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ObjectiveDef.hpp:40
ROL::CompositeStep::computeQuasinormalStep
void computeQuasinormalStep(Vector< Real > &n, const Vector< Real > &c, const Vector< Real > &x, Real delta, Constraint< Real > &con)
Compute quasi-normal step by minimizing the norm of the linearized constraint.
Definition: ROL_CompositeStep.hpp:501
ROL::Constraint::applyAdjointHessian
virtual void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ...
Definition: ROL_ConstraintDef.hpp:132
ROL::Step::getState
ROL::Ptr< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:39
ROL::CompositeStep::update
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, for bound constraints; here only to satisfy the interface requirements, does nothing, needs refactoring.
Definition: ROL_CompositeStep.hpp:345
ROL::CompositeStep::lvec_
ROL::Ptr< Vector< Real > > lvec_
Definition: ROL_CompositeStep.hpp:36
ROL::CompositeStep
Implements the computation of optimization steps with composite-step trust-region methods...
Definition: ROL_CompositeStep.hpp:29
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:39
ROL::Constraint::applyAdjointJacobian
virtual void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
Definition: ROL_ConstraintDef.hpp:52
ROL::Constraint::solveAugmentedSystem
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system where , , , , is an identity or Riesz operator...
Definition: ROL_ConstraintDef.hpp:161
ROL::CompositeStep::solveTangentialSubproblem
void solveTangentialSubproblem(Vector< Real > &t, Vector< Real > &tCP, Vector< Real > &Wg, const Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &n, const Vector< Real > &l, Real delta, Objective< Real > &obj, Constraint< Real > &con)
Solve tangential subproblem.
Definition: ROL_CompositeStep.hpp:611
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
ROL::CompositeStep::printInfoLS
void printInfoLS(const std::vector< Real > &res) const
Definition: ROL_CompositeStep.hpp:91
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
ROL::CompositeStep::computeLagrangeMultiplier
void computeLagrangeMultiplier(Vector< Real > &l, const Vector< Real > &x, const Vector< Real > &gf, Constraint< Real > &con)
Compute Lagrange multipliers by solving the least-squares problem minimizing the gradient of the Lagr...
Definition: ROL_CompositeStep.hpp:434
ROL_Step.hpp
ROL::Constraint
Defines the general constraint operator interface.
Definition: ROL_Constraint.hpp:52
ROL::CompositeStep::printName
std::string printName(void) const
Print step name.
Definition: ROL_CompositeStep.hpp:372
Generated by   doxygen 1.8.5