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