ROL
ROL_BoundFletcher.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
45#ifndef ROL_BOUNDFLETCHER_H
46#define ROL_BOUNDFLETCHER_H
47
48#include "ROL_FletcherBase.hpp"
49#include "ROL_Objective.hpp"
50#include "ROL_Constraint.hpp"
51#include "ROL_Vector.hpp"
52#include "ROL_Types.hpp"
53#include "ROL_Ptr.hpp"
54#include "ROL_Krylov.hpp"
56#include <iostream>
57
58namespace ROL {
59
60template <class Real>
61class BoundFletcher : public FletcherBase<Real> {
62private:
63 // Required for Fletcher penalty function definition
64 using FletcherBase<Real>::obj_;
65 using FletcherBase<Real>::con_;
66 Ptr<const Vector<Real> > low_;
67 Ptr<const Vector<Real> > upp_;
68
71
72 // Evaluation counters
73 using FletcherBase<Real>::nfval_;
74 using FletcherBase<Real>::ngval_;
75 using FletcherBase<Real>::ncval_;
76
77 using FletcherBase<Real>::fPhi_; // value of penalty function
78 using FletcherBase<Real>::gPhi_; // gradient of penalty function
79
80 using FletcherBase<Real>::y_; // multiplier estimate
81
82 using FletcherBase<Real>::fval_; // value of objective function
83 using FletcherBase<Real>::g_; // gradient of objective value
84 using FletcherBase<Real>::c_; // constraint value
85 using FletcherBase<Real>::scaledc_; // penaltyParameter_ * c_
86 using FletcherBase<Real>::gL_; // gradient of Lagrangian (g - A*y)
87
88 using FletcherBase<Real>::cnorm_; // norm of constraint violation
89
96
97 using FletcherBase<Real>::delta_; // regularization parameter
98
99 Ptr<Vector<Real> > Q_;
100 Ptr<Vector<Real> > umx_;
101 Ptr<Vector<Real> > DQ_;
102 Ptr<Vector<Real> > Qsqrt_;
103
105 // Flag to determine type of augmented system to solve
106 // AugSolve_ = 0 : Symmetric system
107 // [ I Q^{1/2} A][w] = [b1]
108 // [A'Q^{1/2} ][v] [b2]
109 // AugSolve_ = 1 : Nonsymmetric system
110 // [ I A][w] = [b1]
111 // [A'Q ][v] [b2]
113
114 Ptr<Vector<Real> > QsgL_; // scaled gradient of Lagrangian Q^{1/2}*(g-A*y)
115 Ptr<Vector<Real> > QgL_; // scaled gradient of Lagrangian Q*(g-A*y)
116 Ptr<Vector<Real> > Qsg_; // Scaled gradient of objective Q^{1/2}*g
117 Ptr<Vector<Real> > DQgL_; // gradient of Lagrangian scaled by DQ, DQ*(g-A*y)
118 Ptr<Vector<Real> > Qv_; // used in augmented system solve
119
120 // Temporaries
121 Ptr<Vector<Real> > Tv_; // Temporary for matvecs
122 Ptr<Vector<Real> > w_; // first component of augmented system solve solution
123 Ptr<Vector<Real> > v_; // second component of augmented system solve solution
124 Ptr<Vector<Real> > htmp1_; // Temporary for rhs
125 Ptr<Vector<Real> > htmp2_; // Temporary for rhs
126
127 Ptr<Vector<Real> > xzeros_; // zero vector
128 Ptr<Vector<Real> > czeros_; // zero vector
129
130 using FletcherBase<Real>::useInexact_;
131
134
135 using FletcherBase<Real>::multSolverError_; // Error from augmented system solve in value()
136 using FletcherBase<Real>::gradSolveError_; // Error from augmented system solve in gradient()
137
138 // For Augmented system solves
139 using FletcherBase<Real>::krylov_;
140 using FletcherBase<Real>::iterKrylov_;
141 using FletcherBase<Real>::flagKrylov_;
142 using FletcherBase<Real>::v1_;
143 using FletcherBase<Real>::v2_;
144 using FletcherBase<Real>::vv_;
145 using FletcherBase<Real>::w1_;
146 using FletcherBase<Real>::w2_;
147 using FletcherBase<Real>::ww_;
148 using FletcherBase<Real>::b1_;
149 using FletcherBase<Real>::b2_;
150 using FletcherBase<Real>::bb_;
151
152 class DiffLower : public Elementwise::BinaryFunction<Real> {
153 public:
154 DiffLower(void) {}
155 Real apply(const Real& x, const Real& y) const {
156 const Real NINF(ROL_NINF<Real>());
157 return (y <= NINF ? static_cast<Real>(-1.) : x - y);
158 }
159 };
160
161 class DiffUpper : public Elementwise::BinaryFunction<Real> {
162 public:
163 DiffUpper(void) {}
164 Real apply(const Real& x, const Real& y) const {
165 const Real INF(ROL_INF<Real>());
166 return (y >= INF ? static_cast<Real>(-1.) : y - x);
167 }
168 };
169
170 class FormQ : public Elementwise::BinaryFunction<Real> {
171 public:
172 FormQ(void) {}
173 Real apply(const Real& x, const Real& y) const {
174 Real zero(0.);
175 if( x < zero && y < zero) {
176 return static_cast<Real>(1);
177 }
178 if( x < zero && y >= zero ) {
179 return y;
180 }
181 if( x >= zero && y < zero ) {
182 return x;
183 }
184 return std::min(x, y);
185 }
186 };
187
188 class FormDQ : public Elementwise::BinaryFunction<Real> {
189 public:
190 FormDQ(void) {}
191 Real apply(const Real& x, const Real& y) const {
192 Real zero(0.), one(1.), mone(-1.);
193 if( x < zero && y < zero) {
194 return zero;
195 }
196 if( x < zero && y >= zero ) {
197 return mone;
198 }
199 if( x >= zero && y < zero ) {
200 return one;
201 }
202 if( x < y ) {
203 return one;
204 } else if( y < x) {
205 return mone;
206 } else {
207 return zero;
208 }
209 }
210 };
211
212 class AugSystemSym : public LinearOperator<Real> {
213 private:
214 const Ptr<Constraint<Real> > con_;
215 const Ptr<const Vector<Real> > x_;
216 const Ptr<Vector<Real> > Qsqrt_;
217 const Ptr<Vector<Real> > Qv_;
218 const Real delta_;
219 public:
221 const Ptr<const Vector<Real> > &x,
222 const Ptr<Vector<Real> > &Qsqrt,
223 const Ptr<Vector<Real> > &Qv,
224 const Real delta) : con_(con), x_(x), Qsqrt_(Qsqrt), Qv_(Qv), delta_(delta) {}
225
226 void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
227 PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
228 const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);
229
230 con_->applyAdjointJacobian(*(Hvp.get(0)), *(vp.get(1)), *x_, tol);
231 Hvp.get(0)->applyBinary(Elementwise::Multiply<Real>(), *Qsqrt_);
232 Hvp.get(0)->plus(*(vp.get(0)));
233
234 Qv_->set(*(vp.get(0)));
235 Qv_->applyBinary(Elementwise::Multiply<Real>(), *Qsqrt_);
236 con_->applyJacobian(*(Hvp.get(1)), *(Qv_), *x_, tol);
237 Hvp.get(1)->axpy(-delta_*delta_, *(vp.get(1)));
238 }
239 };
240
241 class AugSystemNonSym : public LinearOperator<Real> {
242 private:
243 const Ptr<Constraint<Real> > con_;
244 const Ptr<const Vector<Real> > x_;
245 const Ptr<Vector<Real> > Q_;
246 const Ptr<Vector<Real> > Qv_;
247 const Real delta_;
248 public:
250 const Ptr<const Vector<Real> > &x,
251 const Ptr<Vector<Real> > &Q,
252 const Ptr<Vector<Real> > &Qv,
253 const Real delta) : con_(con), x_(x), Q_(Q), Qv_(Qv), delta_(delta) {}
254
255 void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
256 PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
257 const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);
258
259 con_->applyAdjointJacobian(*(Hvp.get(0)), *(vp.get(1)), *x_, tol);
260 Hvp.get(0)->plus(*(vp.get(0)));
261
262 Qv_->set(*(vp.get(0)));
263 Qv_->applyBinary(Elementwise::Multiply<Real>(), *Q_);
264 con_->applyJacobian(*(Hvp.get(1)), *(Qv_), *x_, tol);
265 Hvp.get(1)->axpy(-delta_*delta_, *(vp.get(1)));
266 }
267 };
268
269 class AugSystemPrecond : public LinearOperator<Real> {
270 private:
271 const Ptr<Constraint<Real> > con_;
272 const Ptr<const Vector<Real> > x_;
273 public:
275 const Ptr<const Vector<Real> > x) : con_(con), x_(x) {}
276
277 void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
278 Hv.set(v.dual());
279 }
280 void applyInverse(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
281 Real zero(0);
282 PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
283 const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);
284
285 Hvp.set(0, *(vp.get(0)));
286 // Second x should be dual, but unused?
287 con_->applyPreconditioner(*(Hvp.get(1)),*(vp.get(1)),*x_,*x_, zero);
288 }
289 };
290
291public:
292 BoundFletcher(const ROL::Ptr<Objective<Real> > &obj,
293 const ROL::Ptr<Constraint<Real> > &con,
294 const ROL::Ptr<BoundConstraint<Real> > &bnd,
295 const Vector<Real> &optVec,
296 const Vector<Real> &conVec,
297 ROL::ParameterList &parlist)
298 : FletcherBase<Real>(obj, con), isQComputed_(false), isDQComputed_(false) {
299
300 low_ = bnd->getLowerBound();
301 upp_ = bnd->getUpperBound();
302
303 gPhi_ = optVec.dual().clone();
304 y_ = conVec.dual().clone();
305 g_ = optVec.dual().clone();
306 gL_ = optVec.dual().clone();
307 c_ = conVec.clone();
308 scaledc_ = conVec.clone();
309
310 Q_ = optVec.clone();
311 DQ_ = optVec.clone();
312 umx_ = optVec.clone();
313 Qsqrt_ = optVec.clone();
314 Qv_ = optVec.dual().clone();
315 QsgL_ = optVec.dual().clone();
316 QgL_ = optVec.dual().clone();
317 Qsg_ = optVec.dual().clone();
318 DQgL_ = optVec.dual().clone();
319
320 Tv_ = optVec.dual().clone();
321 w_ = optVec.dual().clone();
322 v_ = conVec.dual().clone();
323 htmp1_ = optVec.dual().clone();
324 htmp2_ = conVec.clone();
325
326 xzeros_ = optVec.dual().clone();
327 xzeros_->zero();
328 czeros_ = conVec.clone();
329 czeros_->zero();
330
331 v1_ = optVec.dual().clone();
332 v2_ = conVec.dual().clone();
333 vv_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({v1_, v2_}));
334
335 w1_ = optVec.dual().clone();
336 w2_ = conVec.dual().clone();
337 ww_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({w1_, w2_}));
338
339 b1_ = optVec.dual().clone();
340 b2_ = conVec.clone();
341 bb_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({b1_, b2_}));
342
343 ROL::ParameterList& sublist = parlist.sublist("Step").sublist("Fletcher");
344 HessianApprox_ = sublist.get("Level of Hessian Approximation", 0);
345
346 AugSolve_ = sublist.get("Type of Augmented System Solve", 0);
347 AugSolve_ = (0 < AugSolve_ && AugSolve_ < 2) ? AugSolve_ : 0;
348
349 penaltyParameter_ = sublist.get("Penalty Parameter", 1.0);
350 quadPenaltyParameter_ = sublist.get("Quadratic Penalty Parameter", 0.0);
351
352 delta_ = sublist.get("Regularization Parameter", 0.0);
353
354 useInexact_ = sublist.get("Inexact Solves", false);
355
356 ROL::ParameterList krylovList;
357 Real atol = static_cast<Real>(1e-12);
358 Real rtol = static_cast<Real>(1e-2);
359 krylovList.sublist("General").sublist("Krylov").set("Type", "GMRES");
360 krylovList.sublist("General").sublist("Krylov").set("Absolute Tolerance", atol);
361 krylovList.sublist("General").sublist("Krylov").set("Relative Tolerance", rtol);
362 krylovList.sublist("General").sublist("Krylov").set("Iteration Limit", 200);
363 krylov_ = KrylovFactory<Real>(krylovList);
364 }
365
366 void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
367 obj_->update(x,flag,iter);
368 con_->update(x,flag,iter);
369 isValueComputed_ = (flag ? false : isValueComputed_);
370 isGradientComputed_ = (flag ? false : isGradientComputed_);
372 isObjValueComputed_ = (flag ? false : isObjValueComputed_);
373 isObjGradComputed_ = (flag ? false : isObjGradComputed_);
374 isConValueComputed_ = (flag ? false : isConValueComputed_);
375 isQComputed_ = (flag ? false : isQComputed_);
376 isDQComputed_ = (flag ? false : isDQComputed_);
377 multSolverError_ = (flag ? ROL_INF<Real>() : multSolverError_);
378 gradSolveError_ = (flag ? ROL_INF<Real>() : gradSolveError_);
379 }
380
381 Real value( const Vector<Real> &x, Real &tol ) {
384 return fPhi_;
385 }
386
387 Real zero(0);
388
389 // Reset tolerances
390 Real origTol = tol;
391 Real tol2 = origTol;
392
393 FletcherBase<Real>::objValue(x, tol2); tol2 = origTol;
394 multSolverError_ = origTol / (static_cast<Real>(2) * std::max(static_cast<Real>(1), cnorm_));
396 tol = multSolverError_;
397
398 fPhi_ = fval_ - c_->dot(y_->dual());
399
401 fPhi_ = fPhi_ + Real(0.5)*quadPenaltyParameter_*(c_->dot(c_->dual()));
402 }
403
404 isValueComputed_ = true;
405
406 return fPhi_;
407 }
408
409 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
410 if( isGradientComputed_ && gradSolveError_ <= tol) {
411 tol = gradSolveError_;
412 g.set(*gPhi_);
413 return;
414 }
415
416 Real zero(0);
417
418 // Reset tolerances
419 Real origTol = tol;
420 Real tol2 = origTol;
421
422 gradSolveError_ = origTol / static_cast<Real>(2);
424
425 bool refine = isGradientComputed_;
426
427 switch( AugSolve_ ) {
428 case 0: {
429 solveAugmentedSystem( *w_, *v_, *xzeros_, *c_, x, gradSolveError_, refine );
431 tol = gradSolveError_;
432
433 w_->applyBinary(Elementwise::Multiply<Real>(), *Qsqrt_);
434 con_->applyAdjointHessian( *gPhi_, *y_, *w_, x, tol2 ); tol2 = origTol;
435 obj_->hessVec( *Tv_, *w_, x, tol2 ); tol2 = origTol;
436 gPhi_->axpy( static_cast<Real>(-1), *Tv_ );
437
438 con_->applyAdjointJacobian( *Tv_, *v_, x, tol2); tol2 = origTol;
439 gPhi_->axpy( -penaltyParameter_, *Tv_);
440
441 Tv_->applyBinary(Elementwise::Multiply<Real>(), *DQgL_);
442 gPhi_->plus( *Tv_ );
443
444 con_->applyAdjointHessian( *Tv_, *v_, *QgL_, x, tol2 ); tol2 = origTol;
445 gPhi_->plus( *Tv_ );
446
447 gPhi_->plus( *gL_ );
448 break;
449 }
450 case 1: {
451 solveAugmentedSystem( *w_, *v_, *xzeros_, *c_, x, gradSolveError_, refine );
453 tol = gradSolveError_;
454
455 gPhi_->set( *w_ );
456 gPhi_->scale( penaltyParameter_ );
457 Tv_->set( *w_ );
458 Tv_->applyBinary( Elementwise::Multiply<Real>(), *DQgL_ );
459 gPhi_->axpy(static_cast<Real>(-1), *Tv_);
460
461 w_->applyBinary( Elementwise::Multiply<Real>(), *Q_ );
462 obj_->hessVec( *Tv_, *w_, x, tol2); tol2 = origTol;
463 gPhi_->axpy( static_cast<Real>(-1), *Tv_ );
464 con_->applyAdjointHessian( *Tv_, *y_, *w_, x, tol2 ); tol2 = origTol;
465 gPhi_->plus( *Tv_ );
466
467 con_->applyAdjointHessian( *Tv_, *v_, *QgL_, x, tol2 ); tol2 = origTol;
468 gPhi_->plus( *Tv_ );
469
470 gPhi_->plus( *gL_ );
471
472 break;
473 }
474 }
475
477 con_->applyAdjointJacobian( *Tv_, *c_, x, tol2 ); tol2 = origTol;
479 }
480
481 g.set(*gPhi_);
482 isGradientComputed_ = true;
483 }
484
485 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
486 Real zero(0);
487
488 // Reset tolerances
489 Real origTol = tol;
490 Real tol2 = origTol;
491
492 // Make sure everything is already computed
493 value(x, tol2); tol2 = origTol;
494 computeMultipliers(x, tol2); tol2 = origTol;
495 gradient(*Tv_, x, tol2); tol2 = origTol;
496
497 switch( AugSolve_ ) {
498 case 0: {
499 // hv <- HL*v
500 obj_->hessVec( hv, v, x, tol2 ); tol2 = origTol;
501 con_->applyAdjointHessian( *Tv_, *y_, v, x, tol2 ); tol2 = origTol;
502 hv.axpy(static_cast<Real>(-1), *Tv_ );
503
504 // htmp1_ <- Q^{1/2}*hv
505 htmp1_->set(hv);
506 htmp1_->applyBinary(Elementwise::Multiply<Real>(), *Qsqrt_);
507 htmp1_->scale(static_cast<Real>(-1));
508 // htmp2_ <- A'*(R(gL) - sigma*I)*v
509 Tv_->set( *DQgL_ );
510 Tv_->applyBinary( Elementwise::Multiply<Real>(), v );
511 Tv_->axpy(-penaltyParameter_, v);
512 con_->applyJacobian( *htmp2_, *Tv_, x, tol2); tol2 = origTol;
513 // v_ <- - (A'QA)^-1 [A' (R(gL)-sigma I) u + A' Q HL u]
514 solveAugmentedSystem( *w_, *v_, *htmp1_, *htmp2_, x, tol2 ); tol2 = origTol;
515 con_->applyAdjointJacobian( *Tv_, *v_, x, tol2 ); tol2 = origTol;
516 hv.plus(*Tv_);
517
518 con_->applyJacobian( *htmp2_, v, x, tol2 ); tol2 = origTol;
519 solveAugmentedSystem( *w_, *v_, *xzeros_, *htmp2_, x, tol2 ); tol2 = origTol;
520 con_->applyAdjointJacobian( *Tv_, *v_, x, tol2 ); tol2 = origTol;
521 hv.axpy( -penaltyParameter_, *Tv_);
522 Tv_->applyBinary( Elementwise::Multiply<Real>(), *DQgL_ );
523 hv.plus( *Tv_ );
524
525 w_->applyBinary( Elementwise::Multiply<Real>(), *Qsqrt_ );
526 obj_->hessVec( *Tv_, *w_, x, tol2 ); tol2 = origTol;
527 hv.axpy( static_cast<Real>(-1), *Tv_);
528 con_->applyAdjointHessian( *Tv_, *y_, *w_, x, tol2 ); tol2 = origTol;
529 hv.plus( *Tv_ );
530 break;
531 }
532 case 1: {
533 // hv <- HL*v
534 obj_->hessVec( hv, v, x, tol2 ); tol2 = origTol;
535 con_->applyAdjointHessian( *Tv_, *y_, v, x, tol2 ); tol2 = origTol;
536 hv.axpy(static_cast<Real>(-1), *Tv_ );
537
538 htmp1_->set(hv);
539 Tv_->set(v);
540 Tv_->applyBinary( Elementwise::Multiply<Real>(), *DQgL_ );
541 Tv_->axpy( -penaltyParameter_, v );
542 Tv_->scale( static_cast<Real>(-1) );
543 con_->applyJacobian( *htmp2_, *Tv_, x, tol2 ); tol2 = origTol;
544 solveAugmentedSystem( *w_, *v_, *htmp1_, *htmp2_, x, tol2 ); tol2 = origTol;
545 hv.set( *w_ );
546
547 con_->applyJacobian( *htmp2_, v, x, tol2 ); tol2 = origTol;
548 solveAugmentedSystem( *w_, *v_, *xzeros_, *htmp2_, x, tol2 ); tol2 = origTol;
549 hv.axpy( penaltyParameter_, *w_ );
550 Tv_->set( *w_ );
551 Tv_->applyBinary( Elementwise::Multiply<Real>(), *DQgL_ );
552 hv.axpy( static_cast<Real>(-1), *Tv_ );
553
554 w_->applyBinary( Elementwise::Multiply<Real>(), *Q_ );
555 obj_->hessVec( *Tv_, *w_, x, tol2 ); tol2 = tol;
556 hv.axpy( static_cast<Real>(-1), *Tv_);
557 con_->applyAdjointHessian( *Tv_, *y_, *w_, x, tol2 ); tol2 = origTol;
558 hv.plus( *Tv_ );
559 break;
560 }
561 }
562
564 con_->applyJacobian( *b2_, v, x, tol2 ); tol2 = origTol;
565 con_->applyAdjointJacobian( *Tv_, *b2_, x, tol2 ); tol2 = origTol;
567 con_->applyAdjointHessian( *Tv_, *c_, v, x, tol2); tol2 = origTol;
569 }
570
571 }
572
574 Vector<Real> &v2,
575 const Vector<Real> &b1,
576 const Vector<Real> &b2,
577 const Vector<Real> &x,
578 Real &tol,
579 bool refine = false) {
580 // Ignore tol for now
581 ROL::Ptr<LinearOperator<Real> > K;
582 switch( AugSolve_ ) {
583 case 0: {
584 K = ROL::makePtr<AugSystemSym>(con_, makePtrFromRef(x), Qsqrt_, Qv_, delta_);
585 break;
586 }
587 case 1: {
588 K = ROL::makePtr<AugSystemNonSym>(con_, makePtrFromRef(x), Q_, Qv_, delta_);
589 break;
590 }
591 }
592 ROL::Ptr<LinearOperator<Real> > P
593 = ROL::makePtr<AugSystemPrecond>(con_, makePtrFromRef(x));
594
595 b1_->set(b1);
596 b2_->set(b2);
597
598 if( refine ) {
599 // TODO: Make sure this tol is actually ok...
600 Real origTol = tol;
601 w1_->set(v1);
602 w2_->set(v2);
603 K->apply(*vv_, *ww_, tol); tol = origTol;
604
605 b1_->axpy( static_cast<Real>(-1), *v1_ );
606 b2_->axpy( static_cast<Real>(-1), *v2_ );
607 }
608
609 v1_->zero();
610 v2_->zero();
611
612 // If inexact, change tolerance
613 if( useInexact_ ) {
614 krylov_->resetAbsoluteTolerance(tol);
615 }
616
617 flagKrylov_ = 0;
618 tol = krylov_->run(*vv_,*K,*bb_,*P,iterKrylov_,flagKrylov_);
619
620 if( refine ) {
621 v1.plus(*v1_);
622 v2.plus(*v2_);
623 } else {
624 v1.set(*v1_);
625 v2.set(*v2_);
626 }
627 }
628
629 void computeMultipliers(const Vector<Real>& x, const Real tol) {
631 return;
632 }
633
634 if( !isMultiplierComputed_ ) {
635 Real tol2 = tol;
636 FletcherBase<Real>::objGrad(x, tol2); tol2 = tol;
637 FletcherBase<Real>::conValue(x, tol2); tol2 = tol;
638 cnorm_ = c_->norm();
639 computeQ(x);
640 computeDQ(x);
641 }
642
643 bool refine = isMultiplierComputed_;
644
645 switch( AugSolve_ ) {
646 case 0: {
647 Qsg_->set(*g_);
648 Qsg_->applyBinary(Elementwise::Multiply<Real>(), *Qsqrt_);
649
650 multSolverError_ = tol;
652
653 gL_->set(*QsgL_);
654 gL_->applyBinary(Elementwise::Divide<Real>(), *Qsqrt_);
655 QgL_->set(*QsgL_);
656 QgL_->applyBinary(Elementwise::Multiply<Real>(), *Qsqrt_);
657 break;
658 }
659 case 1: {
660 multSolverError_ = tol;
662 QgL_->set(*gL_);
663 QgL_->applyBinary(Elementwise::Multiply<Real>(), *Q_);
664 break;
665 }
666 }
667
668 DQgL_->set(*gL_);
669 DQgL_->applyBinary(Elementwise::Multiply<Real>(), *DQ_);
670
672 }
673
674 void computeQ(const Vector<Real>& x) {
675 if( isQComputed_ ) {
676 return;
677 }
678
679 Q_->set(x);
680 Q_->applyBinary(DiffLower(), *low_);
681 umx_->set(x);
682 umx_->applyBinary(DiffUpper(), *upp_);
683 Q_->applyBinary(FormQ(), *umx_);
684 Qsqrt_->set(*Q_);
685 Qsqrt_->applyUnary(Elementwise::SquareRoot<Real>());
686
687 isQComputed_ = true;
688 }
689
690 void computeDQ(const Vector<Real> &x) {
691 if( isDQComputed_ ) {
692 return;
693 }
694
695 DQ_->set(x);
696 DQ_->applyBinary(DiffLower(), *low_);
697 umx_->set(x);
698 umx_->applyBinary(DiffUpper(), *upp_);
699 DQ_->applyBinary(FormDQ(), *umx_);
700
701 isDQComputed_ = true;
702 }
703
704}; // class Fletcher
705
706} // namespace ROL
707
708#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Contains definitions of custom data types in ROL.
Provides the interface to apply upper and lower bound constraints.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
const Ptr< const Vector< Real > > x_
AugSystemNonSym(const Ptr< Constraint< Real > > &con, const Ptr< const Vector< Real > > &x, const Ptr< Vector< Real > > &Q, const Ptr< Vector< Real > > &Qv, const Real delta)
const Ptr< Constraint< Real > > con_
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
const Ptr< const Vector< Real > > x_
const Ptr< Constraint< Real > > con_
AugSystemPrecond(const Ptr< Constraint< Real > > con, const Ptr< const Vector< Real > > x)
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
const Ptr< const Vector< Real > > x_
const Ptr< Vector< Real > > Qsqrt_
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
const Ptr< Vector< Real > > Qv_
AugSystemSym(const Ptr< Constraint< Real > > &con, const Ptr< const Vector< Real > > &x, const Ptr< Vector< Real > > &Qsqrt, const Ptr< Vector< Real > > &Qv, const Real delta)
const Ptr< Constraint< Real > > con_
Real apply(const Real &x, const Real &y) const
Real apply(const Real &x, const Real &y) const
Real apply(const Real &x, const Real &y) const
Real apply(const Real &x, const Real &y) const
Ptr< Vector< Real > > DQ_
void computeQ(const Vector< Real > &x)
Ptr< Vector< Real > > QsgL_
Ptr< Vector< Real > > Q_
Ptr< Vector< Real > > v_
BoundFletcher(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< Constraint< Real > > &con, const ROL::Ptr< BoundConstraint< Real > > &bnd, const Vector< Real > &optVec, const Vector< Real > &conVec, ROL::ParameterList &parlist)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Ptr< Vector< Real > > Qv_
Real value(const Vector< Real > &x, Real &tol)
Compute value.
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
void computeDQ(const Vector< Real > &x)
Ptr< Vector< Real > > QgL_
Ptr< Vector< Real > > xzeros_
Ptr< Vector< Real > > czeros_
Ptr< Vector< Real > > Qsqrt_
Ptr< Vector< Real > > DQgL_
Ptr< const Vector< Real > > low_
Ptr< Vector< Real > > htmp2_
Ptr< Vector< Real > > w_
Ptr< const Vector< Real > > upp_
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Ptr< Vector< Real > > Tv_
Ptr< Vector< Real > > htmp1_
Ptr< Vector< Real > > umx_
void computeMultipliers(const Vector< Real > &x, const Real tol)
Ptr< Vector< Real > > Qsg_
void solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol, bool refine=false)
Defines the general constraint operator interface.
Ptr< Vector< Real > > v1_
Ptr< Vector< Real > > b2_
Ptr< Vector< Real > > c_
Ptr< Vector< Real > > w1_
const Ptr< Constraint< Real > > con_
void conValue(const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > scaledc_
Ptr< Vector< Real > > v2_
const Ptr< Objective< Real > > obj_
Ptr< Vector< Real > > gPhi_
Ptr< Vector< Real > > w2_
void objGrad(const Vector< Real > &x, Real &tol)
void objValue(const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > g_
Ptr< Vector< Real > > y_
Ptr< PartitionedVector< Real > > vv_
Ptr< Vector< Real > > gL_
Ptr< PartitionedVector< Real > > ww_
Ptr< PartitionedVector< Real > > bb_
Ptr< Krylov< Real > > krylov_
Ptr< Vector< Real > > b1_
Provides the interface to apply a linear operator.
Provides the interface to evaluate objective functions.
Defines the linear algebra of vector space on a generic partitioned vector.
ROL::Ptr< const Vector< Real > > get(size_type i) const
void set(const V &x)
Set where .
Defines the linear algebra or vector space interface.
virtual void set(const Vector &x)
Set where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual void plus(const Vector &x)=0
Compute , where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .