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