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"
55 #include "ROL_PartitionedVector.hpp"
56 #include <iostream>
57 
58 namespace ROL {
59 
60 template <class Real>
61 class BoundFletcher : public FletcherBase<Real> {
62 private:
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 
69  using FletcherBase<Real>::penaltyParameter_;
70  using FletcherBase<Real>::quadPenaltyParameter_;
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 
90  using FletcherBase<Real>::isValueComputed_;
91  using FletcherBase<Real>::isGradientComputed_;
92  using FletcherBase<Real>::isMultiplierComputed_;
93  using FletcherBase<Real>::isObjValueComputed_;
94  using FletcherBase<Real>::isObjGradComputed_;
95  using FletcherBase<Real>::isConValueComputed_;
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 
104  int HessianApprox_;
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]
112  int AugSolve_;
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 
132  bool isQComputed_;
133  bool isDQComputed_;
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:
220  AugSystemSym(const Ptr<Constraint<Real> > &con,
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:
249  AugSystemNonSym(const Ptr<Constraint<Real> > &con,
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:
274  AugSystemPrecond(const Ptr<Constraint<Real> > con,
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 
291 public:
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_);
371  isMultiplierComputed_ = (flag ? false : isMultiplierComputed_);
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 ) {
382  if( isValueComputed_ && multSolverError_*cnorm_ <= tol) {
383  tol = multSolverError_*cnorm_;
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_));
395  computeMultipliers(x, multSolverError_);
396  tol = multSolverError_;
397 
398  fPhi_ = fval_ - c_->dot(y_->dual());
399 
400  if( quadPenaltyParameter_ > zero ) {
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);
423  computeMultipliers(x, gradSolveError_);
424 
425  bool refine = isGradientComputed_;
426 
427  switch( AugSolve_ ) {
428  case 0: {
429  solveAugmentedSystem( *w_, *v_, *xzeros_, *c_, x, gradSolveError_, refine );
430  gradSolveError_ += multSolverError_;
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 );
452  gradSolveError_ += multSolverError_;
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 
476  if( quadPenaltyParameter_ > zero ) {
477  con_->applyAdjointJacobian( *Tv_, *c_, x, tol2 ); tol2 = origTol;
478  gPhi_->axpy( quadPenaltyParameter_, *Tv_ );
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 
563  if( quadPenaltyParameter_ > zero ) {
564  con_->applyJacobian( *b2_, v, x, tol2 ); tol2 = origTol;
565  con_->applyAdjointJacobian( *Tv_, *b2_, x, tol2 ); tol2 = origTol;
566  hv.axpy( quadPenaltyParameter_, *Tv_ );
567  con_->applyAdjointHessian( *Tv_, *c_, v, x, tol2); tol2 = origTol;
568  hv.axpy( -quadPenaltyParameter_, *Tv_ );
569  }
570 
571  }
572 
573  void solveAugmentedSystem(Vector<Real> &v1,
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) {
630  if( isMultiplierComputed_ && multSolverError_ <= 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;
651  solveAugmentedSystem(*QsgL_, *y_, *Qsg_, *scaledc_, x, multSolverError_, refine);
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;
661  solveAugmentedSystem(*gL_, *y_, *g_, *scaledc_, x, multSolverError_, refine);
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 
671  isMultiplierComputed_ = true;
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
ROL::FletcherBase::gPhi_
Ptr< Vector< Real > > gPhi_
Definition: ROL_FletcherBase.hpp:75
ROL::FletcherBase::fPhi_
Real fPhi_
Definition: ROL_FletcherBase.hpp:74
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::BoundFletcher::Q_
Ptr< Vector< Real > > Q_
Definition: ROL_BoundFletcher.hpp:99
ROL::BoundFletcher::htmp1_
Ptr< Vector< Real > > htmp1_
Definition: ROL_BoundFletcher.hpp:124
ROL::BoundFletcher::computeDQ
void computeDQ(const Vector< Real > &x)
Definition: ROL_BoundFletcher.hpp:690
ROL::BoundFletcher::w_
Ptr< Vector< Real > > w_
Definition: ROL_BoundFletcher.hpp:122
ROL::BoundFletcher::AugSystemPrecond::con_
const Ptr< Constraint< Real > > con_
Definition: ROL_BoundFletcher.hpp:271
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
ROL::FletcherBase::gL_
Ptr< Vector< Real > > gL_
Definition: ROL_FletcherBase.hpp:83
ROL::FletcherBase::isValueComputed_
bool isValueComputed_
Definition: ROL_FletcherBase.hpp:87
ROL::FletcherBase::c_
Ptr< Vector< Real > > c_
Definition: ROL_FletcherBase.hpp:81
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:112
ROL::BoundFletcher::AugSystemPrecond::AugSystemPrecond
AugSystemPrecond(const Ptr< Constraint< Real > > con, const Ptr< const Vector< Real > > x)
Definition: ROL_BoundFletcher.hpp:274
ROL::BoundFletcher::DQ_
Ptr< Vector< Real > > DQ_
Definition: ROL_BoundFletcher.hpp:101
ROL::FletcherBase::isGradientComputed_
bool isGradientComputed_
Definition: ROL_FletcherBase.hpp:88
ROL::FletcherBase::flagKrylov_
int flagKrylov_
Definition: ROL_FletcherBase.hpp:104
ROL::BoundFletcher::Tv_
Ptr< Vector< Real > > Tv_
Definition: ROL_BoundFletcher.hpp:121
ROL::PartitionedVector::set
void set(const V &x)
Set where .
Definition: ROL_PartitionedVector.hpp:80
ROL::BoundFletcher::AugSystemSym::Qv_
const Ptr< Vector< Real > > Qv_
Definition: ROL_BoundFletcher.hpp:217
ROL::FletcherBase::con_
const Ptr< Constraint< Real > > con_
Definition: ROL_FletcherBase.hpp:64
ROL::BoundFletcher::AugSystemSym::apply
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Definition: ROL_BoundFletcher.hpp:226
ROL::FletcherBase::isMultiplierComputed_
bool isMultiplierComputed_
Definition: ROL_FletcherBase.hpp:89
ROL::BoundFletcher::AugSystemNonSym::delta_
const Real delta_
Definition: ROL_BoundFletcher.hpp:247
ROL::BoundFletcher::low_
Ptr< const Vector< Real > > low_
Definition: ROL_BoundFletcher.hpp:66
ROL::BoundFletcher::AugSystemSym::Qsqrt_
const Ptr< Vector< Real > > Qsqrt_
Definition: ROL_BoundFletcher.hpp:216
ROL::BoundFletcher::AugSystemNonSym::apply
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Definition: ROL_BoundFletcher.hpp:255
ROL::BoundFletcher::isQComputed_
bool isQComputed_
Definition: ROL_BoundFletcher.hpp:132
ROL::BoundFletcher::AugSystemPrecond
Definition: ROL_BoundFletcher.hpp:269
ROL::BoundFletcher::FormQ::apply
Real apply(const Real &x, const Real &y) const
Definition: ROL_BoundFletcher.hpp:173
ROL::BoundFletcher::AugSystemSym::delta_
const Real delta_
Definition: ROL_BoundFletcher.hpp:218
ROL::FletcherBase::conValue
void conValue(const Vector< Real > &x, Real &tol)
Definition: ROL_FletcherBase.hpp:129
ROL::PartitionedVector
Defines the linear algebra of vector space on a generic partitioned vector.
Definition: ROL_PartitionedVector.hpp:60
ROL::BoundFletcher::BoundFletcher
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)
Definition: ROL_BoundFletcher.hpp:292
ROL::BoundFletcher::AugSystemNonSym::con_
const Ptr< Constraint< Real > > con_
Definition: ROL_BoundFletcher.hpp:243
ROL::BoundFletcher::Qv_
Ptr< Vector< Real > > Qv_
Definition: ROL_BoundFletcher.hpp:118
ROL::FletcherBase::multSolverError_
Real multSolverError_
Definition: ROL_FletcherBase.hpp:94
ROL::BoundFletcher::DiffLower::apply
Real apply(const Real &x, const Real &y) const
Definition: ROL_BoundFletcher.hpp:155
ROL::BoundFletcher::QgL_
Ptr< Vector< Real > > QgL_
Definition: ROL_BoundFletcher.hpp:115
ROL::BoundFletcher::Qsg_
Ptr< Vector< Real > > Qsg_
Definition: ROL_BoundFletcher.hpp:116
ROL::FletcherBase::gradSolveError_
Real gradSolveError_
Definition: ROL_FletcherBase.hpp:95
ROL::BoundFletcher::computeMultipliers
void computeMultipliers(const Vector< Real > &x, const Real tol)
Definition: ROL_BoundFletcher.hpp:629
ROL::FletcherBase::b1_
Ptr< Vector< Real > > b1_
Definition: ROL_FletcherBase.hpp:108
ROL::BoundFletcher::DQgL_
Ptr< Vector< Real > > DQgL_
Definition: ROL_BoundFletcher.hpp:117
ROL::FletcherBase::quadPenaltyParameter_
Real quadPenaltyParameter_
Definition: ROL_FletcherBase.hpp:67
ROL::BoundFletcher::xzeros_
Ptr< Vector< Real > > xzeros_
Definition: ROL_BoundFletcher.hpp:127
ROL::FletcherBase::g_
Ptr< Vector< Real > > g_
Definition: ROL_FletcherBase.hpp:80
ROL::FletcherBase::isObjValueComputed_
bool isObjValueComputed_
Definition: ROL_FletcherBase.hpp:90
ROL::BoundFletcher::FormDQ::apply
Real apply(const Real &x, const Real &y) const
Definition: ROL_BoundFletcher.hpp:191
ROL_Constraint.hpp
ROL::BoundFletcher::FormDQ
Definition: ROL_BoundFletcher.hpp:188
ROL::BoundFletcher::FormQ::FormQ
FormQ(void)
Definition: ROL_BoundFletcher.hpp:172
ROL::BoundFletcher::AugSystemSym::con_
const Ptr< Constraint< Real > > con_
Definition: ROL_BoundFletcher.hpp:214
ROL::BoundFletcher::DiffUpper::DiffUpper
DiffUpper(void)
Definition: ROL_BoundFletcher.hpp:163
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::BoundFletcher::update
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Definition: ROL_BoundFletcher.hpp:366
ROL::BoundFletcher::AugSystemNonSym::AugSystemNonSym
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)
Definition: ROL_BoundFletcher.hpp:249
ROL::BoundFletcher::AugSystemSym::x_
const Ptr< const Vector< Real > > x_
Definition: ROL_BoundFletcher.hpp:215
ROL::FletcherBase::v1_
Ptr< Vector< Real > > v1_
Definition: ROL_FletcherBase.hpp:105
ROL::BoundFletcher::czeros_
Ptr< Vector< Real > > czeros_
Definition: ROL_BoundFletcher.hpp:128
ROL::FletcherBase::useInexact_
bool useInexact_
Definition: ROL_FletcherBase.hpp:99
ROL::FletcherBase::objValue
void objValue(const Vector< Real > &x, Real &tol)
Definition: ROL_FletcherBase.hpp:115
ROL::FletcherBase::bb_
Ptr< PartitionedVector< Real > > bb_
Definition: ROL_FletcherBase.hpp:110
ROL::BoundFletcher::DiffUpper
Definition: ROL_BoundFletcher.hpp:161
ROL::BoundFletcher::v_
Ptr< Vector< Real > > v_
Definition: ROL_BoundFletcher.hpp:123
ROL::BoundFletcher::AugSystemNonSym::x_
const Ptr< const Vector< Real > > x_
Definition: ROL_BoundFletcher.hpp:244
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
ROL::FletcherBase::cnorm_
Real cnorm_
Definition: ROL_FletcherBase.hpp:85
ROL::FletcherBase::fval_
Real fval_
Definition: ROL_FletcherBase.hpp:79
ROL::Constraint
Defines the general constraint operator interface.
Definition: ROL_Constraint.hpp:85
ROL::FletcherBase::isConValueComputed_
bool isConValueComputed_
Definition: ROL_FletcherBase.hpp:92
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL::BoundFletcher::AugSystemSym::AugSystemSym
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)
Definition: ROL_BoundFletcher.hpp:220
ROL::BoundFletcher::gradient
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_BoundFletcher.hpp:409
ROL::BoundFletcher::isDQComputed_
bool isDQComputed_
Definition: ROL_BoundFletcher.hpp:133
ROL::BoundFletcher::DiffLower::DiffLower
DiffLower(void)
Definition: ROL_BoundFletcher.hpp:154
ROL::FletcherBase::iterKrylov_
int iterKrylov_
Definition: ROL_FletcherBase.hpp:103
ROL::BoundFletcher::FormQ
Definition: ROL_BoundFletcher.hpp:170
ROL::BoundFletcher::htmp2_
Ptr< Vector< Real > > htmp2_
Definition: ROL_BoundFletcher.hpp:125
ROL::BoundFletcher::AugSystemNonSym
Definition: ROL_BoundFletcher.hpp:241
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
ROL::FletcherBase
Definition: ROL_FletcherBase.hpp:60
ROL::FletcherBase::v2_
Ptr< Vector< Real > > v2_
Definition: ROL_FletcherBase.hpp:106
ROL::BoundFletcher
Definition: ROL_BoundFletcher.hpp:61
ROL
Definition: ROL_ElementwiseVector.hpp:61
ROL::FletcherBase::b2_
Ptr< Vector< Real > > b2_
Definition: ROL_FletcherBase.hpp:109
ROL::BoundFletcher::AugSystemNonSym::Q_
const Ptr< Vector< Real > > Q_
Definition: ROL_BoundFletcher.hpp:245
ROL::FletcherBase::w2_
Ptr< Vector< Real > > w2_
Definition: ROL_FletcherBase.hpp:112
ROL::BoundFletcher::computeQ
void computeQ(const Vector< Real > &x)
Definition: ROL_BoundFletcher.hpp:674
ROL::FletcherBase::y_
Ptr< Vector< Real > > y_
Definition: ROL_FletcherBase.hpp:77
ROL::Vector::dual
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
Definition: ROL_Vector.hpp:226
ROL_FletcherBase.hpp
ROL::BoundFletcher::AugSystemNonSym::Qv_
const Ptr< Vector< Real > > Qv_
Definition: ROL_BoundFletcher.hpp:246
ROL_Krylov.hpp
ROL::BoundFletcher::AugSystemPrecond::applyInverse
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
Definition: ROL_BoundFletcher.hpp:280
ROL::BoundFletcher::value
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Definition: ROL_BoundFletcher.hpp:381
ROL::FletcherBase::isObjGradComputed_
bool isObjGradComputed_
Definition: ROL_FletcherBase.hpp:91
ROL::PartitionedVector::get
ROL::Ptr< const Vector< Real > > get(size_type i) const
Definition: ROL_PartitionedVector.hpp:243
ROL::LinearOperator
Provides the interface to apply a linear operator.
Definition: ROL_LinearOperator.hpp:71
ROL::BoundFletcher::AugSystemPrecond::apply
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Definition: ROL_BoundFletcher.hpp:277
ROL::BoundFletcher::upp_
Ptr< const Vector< Real > > upp_
Definition: ROL_BoundFletcher.hpp:67
ROL::FletcherBase::ww_
Ptr< PartitionedVector< Real > > ww_
Definition: ROL_FletcherBase.hpp:113
ROL::FletcherBase::scaledc_
Ptr< Vector< Real > > scaledc_
Definition: ROL_FletcherBase.hpp:82
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:72
ROL::BoundFletcher::Qsqrt_
Ptr< Vector< Real > > Qsqrt_
Definition: ROL_BoundFletcher.hpp:102
ROL::BoundFletcher::AugSystemPrecond::x_
const Ptr< const Vector< Real > > x_
Definition: ROL_BoundFletcher.hpp:272
ROL::BoundFletcher::AugSolve_
int AugSolve_
Definition: ROL_BoundFletcher.hpp:112
ROL::FletcherBase::obj_
const Ptr< Objective< Real > > obj_
Definition: ROL_FletcherBase.hpp:63
ROL_PartitionedVector.hpp
ROL::BoundFletcher::hessVec
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Definition: ROL_BoundFletcher.hpp:485
ROL::BoundFletcher::QsgL_
Ptr< Vector< Real > > QsgL_
Definition: ROL_BoundFletcher.hpp:114
ROL::FletcherBase::delta_
Real delta_
Definition: ROL_FletcherBase.hpp:97
ROL::BoundFletcher::FormDQ::FormDQ
FormDQ(void)
Definition: ROL_BoundFletcher.hpp:190
ROL_Vector.hpp
ROL::BoundFletcher::AugSystemSym
Definition: ROL_BoundFletcher.hpp:212
ROL_Objective.hpp
ROL::Objective
Provides the interface to evaluate objective functions.
Definition: ROL_Objective.hpp:77
ROL::BoundFletcher::solveAugmentedSystem
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)
Definition: ROL_BoundFletcher.hpp:573
ROL::BoundFletcher::umx_
Ptr< Vector< Real > > umx_
Definition: ROL_BoundFletcher.hpp:100
ROL::BoundFletcher::DiffUpper::apply
Real apply(const Real &x, const Real &y) const
Definition: ROL_BoundFletcher.hpp:164
ROL::BoundFletcher::DiffLower
Definition: ROL_BoundFletcher.hpp:152
ROL::FletcherBase::w1_
Ptr< Vector< Real > > w1_
Definition: ROL_FletcherBase.hpp:111
ROL::BoundFletcher::HessianApprox_
int HessianApprox_
Definition: ROL_BoundFletcher.hpp:104
ROL::FletcherBase::vv_
Ptr< PartitionedVector< Real > > vv_
Definition: ROL_FletcherBase.hpp:107
ROL::FletcherBase::penaltyParameter_
Real penaltyParameter_
Definition: ROL_FletcherBase.hpp:66
ROL::FletcherBase::krylov_
Ptr< Krylov< Real > > krylov_
Definition: ROL_FletcherBase.hpp:102
ROL::FletcherBase::objGrad
void objGrad(const Vector< Real > &x, Real &tol)
Definition: ROL_FletcherBase.hpp:122
Generated by doxygen 1.8.20