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