ROL
ROL_ConstraintDef.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_CONSTRAINT_DEF_H
11 #define ROL_CONSTRAINT_DEF_H
12 
13 #include "ROL_LinearAlgebra.hpp"
14 #include "ROL_LAPACK.hpp"
15 
16 namespace ROL {
17 
18 template <class Real>
19 void Constraint<Real>::applyJacobian(Vector<Real> &jv,
20  const Vector<Real> &v,
21  const Vector<Real> &x,
22  Real &tol) {
23  // By default we compute the finite-difference approximation.
24  const Real one(1);
25  Real ctol = std::sqrt(ROL_EPSILON<Real>());
26 
27  // Get step length.
28  Real h = std::max(one,x.norm()/v.norm())*tol;
29  //Real h = 2.0/(v.norm()*v.norm())*tol;
30 
31  // Compute constraint at x.
32  ROL::Ptr<Vector<Real> > c = jv.clone();
33  this->value(*c,x,ctol);
34 
35  // Compute perturbation x + h*v.
36  ROL::Ptr<Vector<Real> > xnew = x.clone();
37  xnew->set(x);
38  xnew->axpy(h,v);
39  this->update(*xnew,UpdateType::Temp);
40 
41  // Compute constraint at x + h*v.
42  jv.zero();
43  this->value(jv,*xnew,ctol);
44 
45  // Compute Newton quotient.
46  jv.axpy(-one,*c);
47  jv.scale(one/h);
48 }
49 
50 
51 template <class Real>
52 void Constraint<Real>::applyAdjointJacobian(Vector<Real> &ajv,
53  const Vector<Real> &v,
54  const Vector<Real> &x,
55  Real &tol) {
56  applyAdjointJacobian(ajv,v,x,v.dual(),tol);
57 }
58 
59 
60 
61 
62 
63 template <class Real>
64 void Constraint<Real>::applyAdjointJacobian(Vector<Real> &ajv,
65  const Vector<Real> &v,
66  const Vector<Real> &x,
67  const Vector<Real> &dualv,
68  Real &tol) {
69  // By default we compute the finite-difference approximation.
70  // This requires the implementation of a vector-space basis for the optimization variables.
71  // The default implementation requires that the constraint space is equal to its dual.
72  const Real one(1);
73  Real h(0), ctol = std::sqrt(ROL_EPSILON<Real>());
74 
75  ROL::Ptr<Vector<Real> > xnew = x.clone();
76  ROL::Ptr<Vector<Real> > ex = x.clone();
77  ROL::Ptr<Vector<Real> > eajv = ajv.clone();
78  ROL::Ptr<Vector<Real> > cnew = dualv.clone(); // in general, should be in the constraint space
79  ROL::Ptr<Vector<Real> > c0 = dualv.clone(); // in general, should be in the constraint space
80  this->value(*c0,x,ctol);
81 
82  ajv.zero();
83  for ( int i = 0; i < ajv.dimension(); i++ ) {
84  ex = x.basis(i);
85  eajv = ajv.basis(i);
86  h = std::max(one,x.norm()/ex->norm())*tol;
87  xnew->set(x);
88  xnew->axpy(h,*ex);
89  this->update(*xnew,UpdateType::Temp);
90  this->value(*cnew,*xnew,ctol);
91  cnew->axpy(-one,*c0);
92  cnew->scale(one/h);
93  //ajv.axpy(cnew->dot(v.dual()),*eajv);
94  ajv.axpy(cnew->apply(v),*eajv);
95  }
96 }
97 
98 
99 /*template <class Real>
100 void Constraint<Real>::applyHessian(Vector<Real> &huv,
101  const Vector<Real> &u,
102  const Vector<Real> &v,
103  const Vector<Real> &x,
104  Real &tol ) {
105  Real jtol = std::sqrt(ROL_EPSILON<Real>());
106 
107  // Get step length.
108  Real h = std::max(1.0,x.norm()/v.norm())*tol;
109  //Real h = 2.0/(v.norm()*v.norm())*tol;
110 
111  // Compute constraint Jacobian at x.
112  ROL::Ptr<Vector<Real> > ju = huv.clone();
113  this->applyJacobian(*ju,u,x,jtol);
114 
115  // Compute new step x + h*v.
116  ROL::Ptr<Vector<Real> > xnew = x.clone();
117  xnew->set(x);
118  xnew->axpy(h,v);
119  this->update(*xnew);
120 
121  // Compute constraint Jacobian at x + h*v.
122  huv.zero();
123  this->applyJacobian(huv,u,*xnew,jtol);
124 
125  // Compute Newton quotient.
126  huv.axpy(-1.0,*ju);
127  huv.scale(1.0/h);
128 }*/
129 
130 
131 template <class Real>
132 void Constraint<Real>::applyAdjointHessian(Vector<Real> &huv,
133  const Vector<Real> &u,
134  const Vector<Real> &v,
135  const Vector<Real> &x,
136  Real &tol ) {
137  // Get step length.
138  Real h = std::max(static_cast<Real>(1),x.norm()/v.norm())*tol;
139 
140  // Compute constraint Jacobian at x.
141  ROL::Ptr<Vector<Real> > aju = huv.clone();
142  applyAdjointJacobian(*aju,u,x,tol);
143 
144  // Compute new step x + h*v.
145  ROL::Ptr<Vector<Real> > xnew = x.clone();
146  xnew->set(x);
147  xnew->axpy(h,v);
148  update(*xnew,UpdateType::Temp);
149 
150  // Compute constraint Jacobian at x + h*v.
151  huv.zero();
152  applyAdjointJacobian(huv,u,*xnew,tol);
153 
154  // Compute Newton quotient.
155  huv.axpy(static_cast<Real>(-1),*aju);
156  huv.scale(static_cast<Real>(1)/h);
157 }
158 
159 
160 template <class Real>
161 std::vector<Real> Constraint<Real>::solveAugmentedSystem(Vector<Real> &v1,
162  Vector<Real> &v2,
163  const Vector<Real> &b1,
164  const Vector<Real> &b2,
165  const Vector<Real> &x,
166  Real &tol) {
167 
168  /*** Initialization. ***/
169  const Real zero(0), one(1);
170  int m = 200; // Krylov space size.
171  Real zerotol = zero;
172  int i = 0;
173  int k = 0;
174  Real temp = zero;
175  Real resnrm = zero;
176 
177  //tol = std::sqrt(b1.dot(b1)+b2.dot(b2))*1e-8;
178  tol = std::sqrt(b1.dot(b1)+b2.dot(b2))*tol;
179 
180  // Set initial guess to zero.
181  v1.zero(); v2.zero();
182 
183  // Allocate static memory.
184  ROL::Ptr<Vector<Real> > r1 = b1.clone();
185  ROL::Ptr<Vector<Real> > r2 = b2.clone();
186  ROL::Ptr<Vector<Real> > z1 = v1.clone();
187  ROL::Ptr<Vector<Real> > z2 = v2.clone();
188  ROL::Ptr<Vector<Real> > w1 = b1.clone();
189  ROL::Ptr<Vector<Real> > w2 = b2.clone();
190  std::vector<ROL::Ptr<Vector<Real> > > V1;
191  std::vector<ROL::Ptr<Vector<Real> > > V2;
192  ROL::Ptr<Vector<Real> > V2temp = b2.clone();
193  std::vector<ROL::Ptr<Vector<Real> > > Z1;
194  std::vector<ROL::Ptr<Vector<Real> > > Z2;
195  ROL::Ptr<Vector<Real> > w1temp = b1.clone();
196  ROL::Ptr<Vector<Real> > Z2temp = v2.clone();
197 
198  std::vector<Real> res(m+1, zero);
199  LA::Matrix<Real> H(m+1,m);
200  LA::Vector<Real> cs(m);
201  LA::Vector<Real> sn(m);
202  LA::Vector<Real> s(m+1);
203  LA::Vector<Real> y(m+1);
204  LA::Vector<Real> cnorm(m);
205  ROL::LAPACK<int, Real> lapack;
206 
207  // Compute initial residual.
208  applyAdjointJacobian(*r1, v2, x, zerotol);
209  r1->scale(-one); r1->axpy(-one, v1.dual()); r1->plus(b1);
210  applyJacobian(*r2, v1, x, zerotol);
211  r2->scale(-one); r2->plus(b2);
212  res[0] = std::sqrt(r1->dot(*r1) + r2->dot(*r2));
213 
214  // Check if residual is identically zero.
215  if (res[0] == zero) {
216  res.resize(0);
217  return res;
218  }
219 
220  V1.push_back(b1.clone()); (V1[0])->set(*r1); (V1[0])->scale(one/res[0]);
221  V2.push_back(b2.clone()); (V2[0])->set(*r2); (V2[0])->scale(one/res[0]);
222 
223  s(0) = res[0];
224 
225  for (i=0; i<m; i++) {
226 
227  // Apply right preconditioner.
228  V2temp->set(*(V2[i]));
229  applyPreconditioner(*Z2temp, *V2temp, x, b1, zerotol);
230  Z2.push_back(v2.clone()); (Z2[i])->set(*Z2temp);
231  Z1.push_back(v1.clone()); (Z1[i])->set((V1[i])->dual());
232 
233  // Apply operator.
234  applyJacobian(*w2, *(Z1[i]), x, zerotol);
235  applyAdjointJacobian(*w1temp, *Z2temp, x, zerotol);
236  w1->set(*(V1[i])); w1->plus(*w1temp);
237 
238  // Evaluate coefficients and orthogonalize using Gram-Schmidt.
239  for (k=0; k<=i; k++) {
240  H(k,i) = w1->dot(*(V1[k])) + w2->dot(*(V2[k]));
241  w1->axpy(-H(k,i), *(V1[k]));
242  w2->axpy(-H(k,i), *(V2[k]));
243  }
244  H(i+1,i) = std::sqrt(w1->dot(*w1) + w2->dot(*w2));
245 
246  V1.push_back(b1.clone()); (V1[i+1])->set(*w1); (V1[i+1])->scale(one/H(i+1,i));
247  V2.push_back(b2.clone()); (V2[i+1])->set(*w2); (V2[i+1])->scale(one/H(i+1,i));
248 
249  // Apply Givens rotations.
250  for (k=0; k<=i-1; k++) {
251  temp = cs(k)*H(k,i) + sn(k)*H(k+1,i);
252  H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
253  H(k,i) = temp;
254  }
255 
256  // Form i-th rotation matrix.
257  if ( H(i+1,i) == zero ) {
258  cs(i) = one;
259  sn(i) = zero;
260  }
261  else if ( std::abs(H(i+1,i)) > std::abs(H(i,i)) ) {
262  temp = H(i,i) / H(i+1,i);
263  sn(i) = one / std::sqrt( one + temp*temp );
264  cs(i) = temp * sn(i);
265  }
266  else {
267  temp = H(i+1,i) / H(i,i);
268  cs(i) = one / std::sqrt( one + temp*temp );
269  sn(i) = temp * cs(i);
270  }
271 
272  // Approximate residual norm.
273  temp = cs(i)*s(i);
274  s(i+1) = -sn(i)*s(i);
275  s(i) = temp;
276  H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
277  H(i+1,i) = zero;
278  resnrm = std::abs(s(i+1));
279  res[i+1] = resnrm;
280 
281  // Update solution approximation.
282  const char uplo = 'U';
283  const char trans = 'N';
284  const char diag = 'N';
285  const char normin = 'N';
286  Real scaling = zero;
287  int info = 0;
288  y = s;
289  lapack.LATRS(uplo, trans, diag, normin, i+1, H.values(), m+1, y.values(), &scaling, cnorm.values(), &info);
290  z1->zero();
291  z2->zero();
292  for (k=0; k<=i; k++) {
293  z1->axpy(y(k), *(Z1[k]));
294  z2->axpy(y(k), *(Z2[k]));
295  }
296 
297  // Evaluate special stopping condition.
298  //tol = ???;
299 
300 // std::cout << " " << i+1 << ": " << res[i+1]/res[0] << std::endl;
301  if (res[i+1] <= tol) {
302 // std::cout << " solved in " << i+1 << " iterations to " << res[i+1] << " (" << res[i+1]/res[0] << ")" << std::endl;
303  // Update solution vector.
304  v1.plus(*z1);
305  v2.plus(*z2);
306  break;
307  }
308 
309  } // for (int i=0; i++; i<m)
310 
311  res.resize(i+2);
312 
313  /*
314  std::stringstream hist;
315  hist << std::scientific << std::setprecision(8);
316  hist << "\n Augmented System Solver:\n";
317  hist << " Iter Residual\n";
318  for (unsigned j=0; j<res.size(); j++) {
319  hist << " " << std::left << std::setw(14) << res[j] << "\n";
320  }
321  hist << "\n";
322  std::cout << hist.str();
323  */
324 
325  return res;
326 }
327 
328 
329 template <class Real>
330 std::vector<std::vector<Real> > Constraint<Real>::checkApplyJacobian(const Vector<Real> &x,
331  const Vector<Real> &v,
332  const Vector<Real> &jv,
333  const bool printToStream,
334  std::ostream & outStream,
335  const int numSteps,
336  const int order) {
337  std::vector<Real> steps(numSteps);
338  for(int i=0;i<numSteps;++i) {
339  steps[i] = pow(10,-i);
340  }
341 
342  return checkApplyJacobian(x,v,jv,steps,printToStream,outStream,order);
343 }
344 
345 
346 
347 
348 template <class Real>
349 std::vector<std::vector<Real> > Constraint<Real>::checkApplyJacobian(const Vector<Real> &x,
350  const Vector<Real> &v,
351  const Vector<Real> &jv,
352  const std::vector<Real> &steps,
353  const bool printToStream,
354  std::ostream & outStream,
355  const int order) {
356  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
357  "Error: finite difference order must be 1,2,3, or 4" );
358 
359  const Real one(1.0);
360 
361  using Finite_Difference_Arrays::shifts;
362  using Finite_Difference_Arrays::weights;
363 
364  Real tol = std::sqrt(ROL_EPSILON<Real>());
365 
366  int numSteps = steps.size();
367  int numVals = 4;
368  std::vector<Real> tmp(numVals);
369  std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
370 
371  // Save the format state of the original outStream.
372  ROL::nullstream oldFormatState;
373  oldFormatState.copyfmt(outStream);
374 
375  // Compute constraint value at x.
376  ROL::Ptr<Vector<Real> > c = jv.clone();
377  this->update(x,UpdateType::Temp);
378  this->value(*c, x, tol);
379 
380  // Compute (Jacobian at x) times (vector v).
381  ROL::Ptr<Vector<Real> > Jv = jv.clone();
382  this->applyJacobian(*Jv, v, x, tol);
383  Real normJv = Jv->norm();
384 
385  // Temporary vectors.
386  ROL::Ptr<Vector<Real> > cdif = jv.clone();
387  ROL::Ptr<Vector<Real> > cnew = jv.clone();
388  ROL::Ptr<Vector<Real> > xnew = x.clone();
389 
390  for (int i=0; i<numSteps; i++) {
391 
392  Real eta = steps[i];
393 
394  xnew->set(x);
395 
396  cdif->set(*c);
397  cdif->scale(weights[order-1][0]);
398 
399  for(int j=0; j<order; ++j) {
400 
401  xnew->axpy(eta*shifts[order-1][j], v);
402 
403  if( weights[order-1][j+1] != 0 ) {
404  this->update(*xnew,UpdateType::Temp);
405  this->value(*cnew,*xnew,tol);
406  cdif->axpy(weights[order-1][j+1],*cnew);
407  }
408 
409  }
410 
411  cdif->scale(one/eta);
412 
413  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
414  jvCheck[i][0] = eta;
415  jvCheck[i][1] = normJv;
416  jvCheck[i][2] = cdif->norm();
417  cdif->axpy(-one, *Jv);
418  jvCheck[i][3] = cdif->norm();
419 
420  if (printToStream) {
421  std::stringstream hist;
422  if (i==0) {
423  hist << std::right
424  << std::setw(20) << "Step size"
425  << std::setw(20) << "norm(Jac*vec)"
426  << std::setw(20) << "norm(FD approx)"
427  << std::setw(20) << "norm(abs error)"
428  << "\n"
429  << std::setw(20) << "---------"
430  << std::setw(20) << "-------------"
431  << std::setw(20) << "---------------"
432  << std::setw(20) << "---------------"
433  << "\n";
434  }
435  hist << std::scientific << std::setprecision(11) << std::right
436  << std::setw(20) << jvCheck[i][0]
437  << std::setw(20) << jvCheck[i][1]
438  << std::setw(20) << jvCheck[i][2]
439  << std::setw(20) << jvCheck[i][3]
440  << "\n";
441  outStream << hist.str();
442  }
443 
444  }
445 
446  // Reset format state of outStream.
447  outStream.copyfmt(oldFormatState);
448 
449  return jvCheck;
450 } // checkApplyJacobian
451 
452 
453 template <class Real>
454 std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointJacobian(const Vector<Real> &x,
455  const Vector<Real> &v,
456  const Vector<Real> &c,
457  const Vector<Real> &ajv,
458  const bool printToStream,
459  std::ostream & outStream,
460  const int numSteps) {
461  Real tol = std::sqrt(ROL_EPSILON<Real>());
462  const Real one(1);
463 
464  int numVals = 4;
465  std::vector<Real> tmp(numVals);
466  std::vector<std::vector<Real> > ajvCheck(numSteps, tmp);
467  Real eta_factor = 1e-1;
468  Real eta = one;
469 
470  // Temporary vectors.
471  ROL::Ptr<Vector<Real> > c0 = c.clone();
472  ROL::Ptr<Vector<Real> > cnew = c.clone();
473  ROL::Ptr<Vector<Real> > xnew = x.clone();
474  ROL::Ptr<Vector<Real> > ajv0 = ajv.clone();
475  ROL::Ptr<Vector<Real> > ajv1 = ajv.clone();
476  ROL::Ptr<Vector<Real> > ex = x.clone();
477  ROL::Ptr<Vector<Real> > eajv = ajv.clone();
478 
479  // Save the format state of the original outStream.
480  ROL::nullstream oldFormatState;
481  oldFormatState.copyfmt(outStream);
482 
483  // Compute constraint value at x.
484  this->update(x,UpdateType::Temp);
485  this->value(*c0, x, tol);
486 
487  // Compute (Jacobian at x) times (vector v).
488  this->applyAdjointJacobian(*ajv0, v, x, tol);
489  Real normAJv = ajv0->norm();
490 
491  for (int i=0; i<numSteps; i++) {
492 
493  ajv1->zero();
494 
495  for ( int j = 0; j < ajv.dimension(); j++ ) {
496  ex = x.basis(j);
497  eajv = ajv.basis(j);
498  xnew->set(x);
499  xnew->axpy(eta,*ex);
500  this->update(*xnew,UpdateType::Temp);
501  this->value(*cnew,*xnew,tol);
502  cnew->axpy(-one,*c0);
503  cnew->scale(one/eta);
504  //ajv1->axpy(cnew->dot(v.dual()),*eajv);
505  ajv1->axpy(cnew->apply(v),*eajv);
506  }
507 
508  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
509  ajvCheck[i][0] = eta;
510  ajvCheck[i][1] = normAJv;
511  ajvCheck[i][2] = ajv1->norm();
512  ajv1->axpy(-one, *ajv0);
513  ajvCheck[i][3] = ajv1->norm();
514 
515  if (printToStream) {
516  std::stringstream hist;
517  if (i==0) {
518  hist << std::right
519  << std::setw(20) << "Step size"
520  << std::setw(20) << "norm(adj(Jac)*vec)"
521  << std::setw(20) << "norm(FD approx)"
522  << std::setw(20) << "norm(abs error)"
523  << "\n"
524  << std::setw(20) << "---------"
525  << std::setw(20) << "------------------"
526  << std::setw(20) << "---------------"
527  << std::setw(20) << "---------------"
528  << "\n";
529  }
530  hist << std::scientific << std::setprecision(11) << std::right
531  << std::setw(20) << ajvCheck[i][0]
532  << std::setw(20) << ajvCheck[i][1]
533  << std::setw(20) << ajvCheck[i][2]
534  << std::setw(20) << ajvCheck[i][3]
535  << "\n";
536  outStream << hist.str();
537  }
538 
539  // Update eta.
540  eta = eta*eta_factor;
541  }
542 
543  // Reset format state of outStream.
544  outStream.copyfmt(oldFormatState);
545 
546  return ajvCheck;
547 } // checkApplyAdjointJacobian
548 
549 template <class Real>
550 Real Constraint<Real>::checkAdjointConsistencyJacobian(const Vector<Real> &w,
551  const Vector<Real> &v,
552  const Vector<Real> &x,
553  const Vector<Real> &dualw,
554  const Vector<Real> &dualv,
555  const bool printToStream,
556  std::ostream & outStream) {
557  Real tol = ROL_EPSILON<Real>();
558 
559  ROL::Ptr<Vector<Real> > Jv = dualw.clone();
560  ROL::Ptr<Vector<Real> > Jw = dualv.clone();
561 
562  this->update(x,UpdateType::Temp);
563  applyJacobian(*Jv,v,x,tol);
564  applyAdjointJacobian(*Jw,w,x,tol);
565 
566  //Real vJw = v.dot(Jw->dual());
567  Real vJw = v.apply(*Jw);
568  //Real wJv = w.dot(Jv->dual());
569  Real wJv = w.apply(*Jv);
570 
571  Real diff = std::abs(wJv-vJw);
572 
573  if ( printToStream ) {
574  std::stringstream hist;
575  hist << std::scientific << std::setprecision(8);
576  hist << "\nTest Consistency of Jacobian and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
577  << diff << "\n";
578  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
579  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
580  outStream << hist.str();
581  }
582  return diff;
583 } // checkAdjointConsistencyJacobian
584 
585 template <class Real>
586 std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointHessian(const Vector<Real> &x,
587  const Vector<Real> &u,
588  const Vector<Real> &v,
589  const Vector<Real> &hv,
590  const bool printToStream,
591  std::ostream & outStream,
592  const int numSteps,
593  const int order) {
594  std::vector<Real> steps(numSteps);
595  for(int i=0;i<numSteps;++i) {
596  steps[i] = pow(10,-i);
597  }
598 
599  return checkApplyAdjointHessian(x,u,v,hv,steps,printToStream,outStream,order);
600 }
601 
602 
603 template <class Real>
604 std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointHessian(const Vector<Real> &x,
605  const Vector<Real> &u,
606  const Vector<Real> &v,
607  const Vector<Real> &hv,
608  const std::vector<Real> &steps,
609  const bool printToStream,
610  std::ostream & outStream,
611  const int order) {
612  using Finite_Difference_Arrays::shifts;
613  using Finite_Difference_Arrays::weights;
614 
615  const Real one(1);
616  Real tol = std::sqrt(ROL_EPSILON<Real>());
617 
618  int numSteps = steps.size();
619  int numVals = 4;
620  std::vector<Real> tmp(numVals);
621  std::vector<std::vector<Real> > ahuvCheck(numSteps, tmp);
622 
623  // Temporary vectors.
624  ROL::Ptr<Vector<Real> > AJdif = hv.clone();
625  ROL::Ptr<Vector<Real> > AJu = hv.clone();
626  ROL::Ptr<Vector<Real> > AHuv = hv.clone();
627  ROL::Ptr<Vector<Real> > AJnew = hv.clone();
628  ROL::Ptr<Vector<Real> > xnew = x.clone();
629 
630  // Save the format state of the original outStream.
631  ROL::nullstream oldFormatState;
632  oldFormatState.copyfmt(outStream);
633 
634  // Apply adjoint Jacobian to u.
635  this->update(x,UpdateType::Temp);
636  this->applyAdjointJacobian(*AJu, u, x, tol);
637 
638  // Apply adjoint Hessian at x, in direction v, to u.
639  this->applyAdjointHessian(*AHuv, u, v, x, tol);
640  Real normAHuv = AHuv->norm();
641 
642  for (int i=0; i<numSteps; i++) {
643 
644  Real eta = steps[i];
645 
646  // Apply adjoint Jacobian to u at x+eta*v.
647  xnew->set(x);
648 
649  AJdif->set(*AJu);
650  AJdif->scale(weights[order-1][0]);
651 
652  for(int j=0; j<order; ++j) {
653 
654  xnew->axpy(eta*shifts[order-1][j],v);
655 
656  if( weights[order-1][j+1] != 0 ) {
657  this->update(*xnew,UpdateType::Temp);
658  this->applyAdjointJacobian(*AJnew, u, *xnew, tol);
659  AJdif->axpy(weights[order-1][j+1],*AJnew);
660  }
661  }
662 
663  AJdif->scale(one/eta);
664 
665  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
666  ahuvCheck[i][0] = eta;
667  ahuvCheck[i][1] = normAHuv;
668  ahuvCheck[i][2] = AJdif->norm();
669  AJdif->axpy(-one, *AHuv);
670  ahuvCheck[i][3] = AJdif->norm();
671 
672  if (printToStream) {
673  std::stringstream hist;
674  if (i==0) {
675  hist << std::right
676  << std::setw(20) << "Step size"
677  << std::setw(20) << "norm(adj(H)(u,v))"
678  << std::setw(20) << "norm(FD approx)"
679  << std::setw(20) << "norm(abs error)"
680  << "\n"
681  << std::setw(20) << "---------"
682  << std::setw(20) << "-----------------"
683  << std::setw(20) << "---------------"
684  << std::setw(20) << "---------------"
685  << "\n";
686  }
687  hist << std::scientific << std::setprecision(11) << std::right
688  << std::setw(20) << ahuvCheck[i][0]
689  << std::setw(20) << ahuvCheck[i][1]
690  << std::setw(20) << ahuvCheck[i][2]
691  << std::setw(20) << ahuvCheck[i][3]
692  << "\n";
693  outStream << hist.str();
694  }
695 
696  }
697 
698  // Reset format state of outStream.
699  outStream.copyfmt(oldFormatState);
700 
701  return ahuvCheck;
702 } // checkApplyAdjointHessian
703 
704 } // namespace ROL
705 
706 #endif
ROL::Vector::dual
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
ROL::Constraint::applyJacobian
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
Definition: ROL_ConstraintDef.hpp:19
ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Constraint::checkAdjointConsistencyJacobian
virtual Real checkAdjointConsistencyJacobian(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &x, const bool printToStream=true, std::ostream &outStream=std::cout)
Definition: ROL_Constraint.hpp:320
ROL::Vector::dimension
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:162
ROL::Vector::apply
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: ROL_Vector.hpp:204
ROL::Vector::basis
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:148
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Finite_Difference_Arrays::weights
const double weights[4][5]
Definition: ROL_Types.hpp:834
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
ROL::update
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1) override
Definition: ROL_Objective_TimeSimOpt.hpp:140
ROL::Constraint::checkApplyAdjointJacobian
virtual std::vector< std::vector< Real > > checkApplyAdjointJacobian(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &c, const Vector< Real > &ajv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS)
Finite-difference check for the application of the adjoint of constraint Jacobian.
Definition: ROL_ConstraintDef.hpp:454
ROL::Finite_Difference_Arrays::shifts
const int shifts[4][4]
Definition: ROL_Types.hpp:827
ROL::value
ROL::Objective_SimOpt value
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL::Vector::dot
virtual Real dot(const Vector &x) const =0
Compute where .
applyJacobian
void applyJacobian(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &sol)
Definition: step/interiorpoint/test_02.cpp:118
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:111
ROL::Objective_SerialSimOpt::zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:46
ROL::Constraint::applyAdjointHessian
virtual void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ...
Definition: ROL_ConstraintDef.hpp:132
ROL::Constraint::checkApplyAdjointHessian
virtual std::vector< std::vector< Real > > checkApplyAdjointHessian(const Vector< Real > &x, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &step, const bool printToScreen=true, std::ostream &outStream=std::cout, const int order=1)
Finite-difference check for the application of the adjoint of constraint Hessian. ...
Definition: ROL_ConstraintDef.hpp:604
ROL::Constraint::checkApplyJacobian
virtual std::vector< std::vector< Real > > checkApplyJacobian(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &jv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Finite-difference check for the constraint Jacobian application.
Definition: ROL_ConstraintDef.hpp:349
ROL::Constraint::applyAdjointJacobian
virtual void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
Definition: ROL_ConstraintDef.hpp:52
ROL::details::nullstream
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
ROL::Constraint::solveAugmentedSystem
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system where , , , , is an identity or Riesz operator...
Definition: ROL_ConstraintDef.hpp:161
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
Generated by   doxygen 1.8.5