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