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