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  Ptr<Vector<Real>> r1 = b1.clone();
185  Ptr<Vector<Real>> r2 = b2.clone();
186  Ptr<Vector<Real>> z1 = v1.clone();
187  Ptr<Vector<Real>> z2 = v2.clone();
188  Ptr<Vector<Real>> w1 = b1.clone();
189  Ptr<Vector<Real>> w2 = b2.clone();
190  std::vector<Ptr<Vector<Real>>> V1;
191  std::vector<Ptr<Vector<Real>>> V2;
192  Ptr<Vector<Real>> V2temp = b2.clone();
193  std::vector<Ptr<Vector<Real>>> Z1;
194  std::vector<Ptr<Vector<Real>>> Z2;
195  Ptr<Vector<Real>> w1temp = b1.clone();
196  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  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  r1->set(b1); r2->set(b2);
213  res[0] = std::sqrt(r1->dot(*r1) + r2->dot(*r2));
214 
215  // Check if residual is identically zero.
216  if (res[0] == zero) {
217  res.resize(0);
218  return res;
219  }
220 
221  V1.push_back(b1.clone()); (V1[0])->set(*r1); (V1[0])->scale(one/res[0]);
222  V2.push_back(b2.clone()); (V2[0])->set(*r2); (V2[0])->scale(one/res[0]);
223 
224  s(0) = res[0];
225 
226  for (i=0; i<m; i++) {
227 
228  // Apply right preconditioner.
229  V2temp->set(*(V2[i]));
230  applyPreconditioner(*Z2temp, *V2temp, x, b1, zerotol);
231  Z2.push_back(v2.clone()); (Z2[i])->set(*Z2temp);
232  Z1.push_back(v1.clone()); (Z1[i])->set((V1[i])->dual());
233 
234  // Apply operator.
235  applyJacobian(*w2, *(Z1[i]), x, zerotol);
236  applyAdjointJacobian(*w1temp, *Z2temp, x, zerotol);
237  w1->set(*(V1[i])); w1->plus(*w1temp);
238 
239  // Evaluate coefficients and orthogonalize using Gram-Schmidt.
240  for (k=0; k<=i; k++) {
241  H(k,i) = w1->dot(*(V1[k])) + w2->dot(*(V2[k]));
242  w1->axpy(-H(k,i), *(V1[k]));
243  w2->axpy(-H(k,i), *(V2[k]));
244  }
245  H(i+1,i) = std::sqrt(w1->dot(*w1) + w2->dot(*w2));
246 
247  V1.push_back(b1.clone()); (V1[i+1])->set(*w1); (V1[i+1])->scale(one/H(i+1,i));
248  V2.push_back(b2.clone()); (V2[i+1])->set(*w2); (V2[i+1])->scale(one/H(i+1,i));
249 
250  // Apply Givens rotations.
251  for (k=0; k<=i-1; k++) {
252  temp = cs(k)*H(k,i) + sn(k)*H(k+1,i);
253  H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
254  H(k,i) = temp;
255  }
256 
257  // Form i-th rotation matrix.
258  if ( H(i+1,i) == zero ) {
259  cs(i) = one;
260  sn(i) = zero;
261  }
262  else if ( std::abs(H(i+1,i)) > std::abs(H(i,i)) ) {
263  temp = H(i,i) / H(i+1,i);
264  sn(i) = one / std::sqrt( one + temp*temp );
265  cs(i) = temp * sn(i);
266  }
267  else {
268  temp = H(i+1,i) / H(i,i);
269  cs(i) = one / std::sqrt( one + temp*temp );
270  sn(i) = temp * cs(i);
271  }
272 
273  // Approximate residual norm.
274  temp = cs(i)*s(i);
275  s(i+1) = -sn(i)*s(i);
276  s(i) = temp;
277  H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
278  H(i+1,i) = zero;
279  resnrm = std::abs(s(i+1));
280  res[i+1] = resnrm;
281 
282  // Update solution approximation.
283  const char uplo = 'U';
284  const char trans = 'N';
285  const char diag = 'N';
286  const char normin = 'N';
287  Real scaling = zero;
288  int info = 0;
289  y = s;
290  lapack.LATRS(uplo, trans, diag, normin, i+1, H.values(), m+1, y.values(), &scaling, cnorm.values(), &info);
291  z1->zero();
292  z2->zero();
293  for (k=0; k<=i; k++) {
294  z1->axpy(y(k), *(Z1[k]));
295  z2->axpy(y(k), *(Z2[k]));
296  }
297 
298  // Evaluate special stopping condition.
299  //tol = ???;
300 
301 // std::cout << " " << i+1 << ": " << res[i+1]/res[0] << std::endl;
302  if (res[i+1] <= tol) {
303 // std::cout << " solved in " << i+1 << " iterations to " << res[i+1] << " (" << res[i+1]/res[0] << ")" << std::endl;
304  // Update solution vector.
305  //v1.plus(*z1);
306  //v2.plus(*z2);
307  break;
308  }
309 
310  } // for (int i=0; i++; i<m)
311  v1.set(*z1); v2.set(*z2);
312 
313  res.resize(i+2);
314 
315  /*
316  std::stringstream hist;
317  hist << std::scientific << std::setprecision(8);
318  hist << "\n Augmented System Solver:\n";
319  hist << " Iter Residual\n";
320  for (unsigned j=0; j<res.size(); j++) {
321  hist << " " << std::left << std::setw(14) << res[j] << "\n";
322  }
323  hist << "\n";
324  std::cout << hist.str();
325  */
326 
327  return res;
328 }
329 
330 
331 template <class Real>
332 std::vector<std::vector<Real> > Constraint<Real>::checkApplyJacobian(const Vector<Real> &x,
333  const Vector<Real> &v,
334  const Vector<Real> &jv,
335  const bool printToStream,
336  std::ostream & outStream,
337  const int numSteps,
338  const int order) {
339  std::vector<Real> steps(numSteps);
340  for(int i=0;i<numSteps;++i) {
341  steps[i] = pow(10,-i);
342  }
343 
344  return checkApplyJacobian(x,v,jv,steps,printToStream,outStream,order);
345 }
346 
347 
348 
349 
350 template <class Real>
351 std::vector<std::vector<Real> > Constraint<Real>::checkApplyJacobian(const Vector<Real> &x,
352  const Vector<Real> &v,
353  const Vector<Real> &jv,
354  const std::vector<Real> &steps,
355  const bool printToStream,
356  std::ostream & outStream,
357  const int order) {
358  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
359  "Error: finite difference order must be 1,2,3, or 4" );
360 
361  const Real one(1.0);
362 
363  using Finite_Difference_Arrays::shifts;
364  using Finite_Difference_Arrays::weights;
365 
366  Real tol = std::sqrt(ROL_EPSILON<Real>());
367 
368  int numSteps = steps.size();
369  int numVals = 4;
370  std::vector<Real> tmp(numVals);
371  std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
372 
373  // Save the format state of the original outStream.
374  ROL::nullstream oldFormatState;
375  oldFormatState.copyfmt(outStream);
376 
377  // Compute constraint value at x.
378  ROL::Ptr<Vector<Real> > c = jv.clone();
379  this->update(x,UpdateType::Temp);
380  this->value(*c, x, tol);
381 
382  // Compute (Jacobian at x) times (vector v).
383  ROL::Ptr<Vector<Real> > Jv = jv.clone();
384  this->applyJacobian(*Jv, v, x, tol);
385  Real normJv = Jv->norm();
386 
387  // Temporary vectors.
388  ROL::Ptr<Vector<Real> > cdif = jv.clone();
389  ROL::Ptr<Vector<Real> > cnew = jv.clone();
390  ROL::Ptr<Vector<Real> > xnew = x.clone();
391 
392  for (int i=0; i<numSteps; i++) {
393 
394  Real eta = steps[i];
395 
396  xnew->set(x);
397 
398  cdif->set(*c);
399  cdif->scale(weights[order-1][0]);
400 
401  for(int j=0; j<order; ++j) {
402 
403  xnew->axpy(eta*shifts[order-1][j], v);
404 
405  if( weights[order-1][j+1] != 0 ) {
406  this->update(*xnew,UpdateType::Temp);
407  this->value(*cnew,*xnew,tol);
408  cdif->axpy(weights[order-1][j+1],*cnew);
409  }
410 
411  }
412 
413  cdif->scale(one/eta);
414 
415  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
416  jvCheck[i][0] = eta;
417  jvCheck[i][1] = normJv;
418  jvCheck[i][2] = cdif->norm();
419  cdif->axpy(-one, *Jv);
420  jvCheck[i][3] = cdif->norm();
421 
422  if (printToStream) {
423  std::stringstream hist;
424  if (i==0) {
425  hist << std::right
426  << std::setw(20) << "Step size"
427  << std::setw(20) << "norm(Jac*vec)"
428  << std::setw(20) << "norm(FD approx)"
429  << std::setw(20) << "norm(abs error)"
430  << "\n"
431  << std::setw(20) << "---------"
432  << std::setw(20) << "-------------"
433  << std::setw(20) << "---------------"
434  << std::setw(20) << "---------------"
435  << "\n";
436  }
437  hist << std::scientific << std::setprecision(11) << std::right
438  << std::setw(20) << jvCheck[i][0]
439  << std::setw(20) << jvCheck[i][1]
440  << std::setw(20) << jvCheck[i][2]
441  << std::setw(20) << jvCheck[i][3]
442  << "\n";
443  outStream << hist.str();
444  }
445 
446  }
447 
448  // Reset format state of outStream.
449  outStream.copyfmt(oldFormatState);
450 
451  return jvCheck;
452 } // checkApplyJacobian
453 
454 
455 template <class Real>
456 std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointJacobian(const Vector<Real> &x,
457  const Vector<Real> &v,
458  const Vector<Real> &c,
459  const Vector<Real> &ajv,
460  const bool printToStream,
461  std::ostream & outStream,
462  const int numSteps) {
463  Real tol = std::sqrt(ROL_EPSILON<Real>());
464  const Real one(1);
465 
466  int numVals = 4;
467  std::vector<Real> tmp(numVals);
468  std::vector<std::vector<Real> > ajvCheck(numSteps, tmp);
469  Real eta_factor = 1e-1;
470  Real eta = one;
471 
472  // Temporary vectors.
473  ROL::Ptr<Vector<Real> > c0 = c.clone();
474  ROL::Ptr<Vector<Real> > cnew = c.clone();
475  ROL::Ptr<Vector<Real> > xnew = x.clone();
476  ROL::Ptr<Vector<Real> > ajv0 = ajv.clone();
477  ROL::Ptr<Vector<Real> > ajv1 = ajv.clone();
478  ROL::Ptr<Vector<Real> > ex = x.clone();
479  ROL::Ptr<Vector<Real> > eajv = ajv.clone();
480 
481  // Save the format state of the original outStream.
482  ROL::nullstream oldFormatState;
483  oldFormatState.copyfmt(outStream);
484 
485  // Compute constraint value at x.
486  this->update(x,UpdateType::Temp);
487  this->value(*c0, x, tol);
488 
489  // Compute (Jacobian at x) times (vector v).
490  this->applyAdjointJacobian(*ajv0, v, x, tol);
491  Real normAJv = ajv0->norm();
492 
493  for (int i=0; i<numSteps; i++) {
494 
495  ajv1->zero();
496 
497  for ( int j = 0; j < ajv.dimension(); j++ ) {
498  ex = x.basis(j);
499  eajv = ajv.basis(j);
500  xnew->set(x);
501  xnew->axpy(eta,*ex);
502  this->update(*xnew,UpdateType::Temp);
503  this->value(*cnew,*xnew,tol);
504  cnew->axpy(-one,*c0);
505  cnew->scale(one/eta);
506  //ajv1->axpy(cnew->dot(v.dual()),*eajv);
507  ajv1->axpy(cnew->apply(v),*eajv);
508  }
509 
510  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
511  ajvCheck[i][0] = eta;
512  ajvCheck[i][1] = normAJv;
513  ajvCheck[i][2] = ajv1->norm();
514  ajv1->axpy(-one, *ajv0);
515  ajvCheck[i][3] = ajv1->norm();
516 
517  if (printToStream) {
518  std::stringstream hist;
519  if (i==0) {
520  hist << std::right
521  << std::setw(20) << "Step size"
522  << std::setw(20) << "norm(adj(Jac)*vec)"
523  << std::setw(20) << "norm(FD approx)"
524  << std::setw(20) << "norm(abs error)"
525  << "\n"
526  << std::setw(20) << "---------"
527  << std::setw(20) << "------------------"
528  << std::setw(20) << "---------------"
529  << std::setw(20) << "---------------"
530  << "\n";
531  }
532  hist << std::scientific << std::setprecision(11) << std::right
533  << std::setw(20) << ajvCheck[i][0]
534  << std::setw(20) << ajvCheck[i][1]
535  << std::setw(20) << ajvCheck[i][2]
536  << std::setw(20) << ajvCheck[i][3]
537  << "\n";
538  outStream << hist.str();
539  }
540 
541  // Update eta.
542  eta = eta*eta_factor;
543  }
544 
545  // Reset format state of outStream.
546  outStream.copyfmt(oldFormatState);
547 
548  return ajvCheck;
549 } // checkApplyAdjointJacobian
550 
551 template <class Real>
552 Real Constraint<Real>::checkAdjointConsistencyJacobian(const Vector<Real> &w,
553  const Vector<Real> &v,
554  const Vector<Real> &x,
555  const Vector<Real> &dualw,
556  const Vector<Real> &dualv,
557  const bool printToStream,
558  std::ostream & outStream) {
559  Real tol = ROL_EPSILON<Real>();
560 
561  ROL::Ptr<Vector<Real> > Jv = dualw.clone();
562  ROL::Ptr<Vector<Real> > Jw = dualv.clone();
563 
564  this->update(x,UpdateType::Temp);
565  applyJacobian(*Jv,v,x,tol);
566  applyAdjointJacobian(*Jw,w,x,tol);
567 
568  //Real vJw = v.dot(Jw->dual());
569  Real vJw = v.apply(*Jw);
570  //Real wJv = w.dot(Jv->dual());
571  Real wJv = w.apply(*Jv);
572 
573  Real diff = std::abs(wJv-vJw);
574 
575  if ( printToStream ) {
576  std::stringstream hist;
577  hist << std::scientific << std::setprecision(8);
578  hist << "\nTest Consistency of Jacobian and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
579  << diff << "\n";
580  hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
581  hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
582  outStream << hist.str();
583  }
584  return diff;
585 } // checkAdjointConsistencyJacobian
586 
587 template <class Real>
588 std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointHessian(const Vector<Real> &x,
589  const Vector<Real> &u,
590  const Vector<Real> &v,
591  const Vector<Real> &hv,
592  const bool printToStream,
593  std::ostream & outStream,
594  const int numSteps,
595  const int order) {
596  std::vector<Real> steps(numSteps);
597  for(int i=0;i<numSteps;++i) {
598  steps[i] = pow(10,-i);
599  }
600 
601  return checkApplyAdjointHessian(x,u,v,hv,steps,printToStream,outStream,order);
602 }
603 
604 
605 template <class Real>
606 std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointHessian(const Vector<Real> &x,
607  const Vector<Real> &u,
608  const Vector<Real> &v,
609  const Vector<Real> &hv,
610  const std::vector<Real> &steps,
611  const bool printToStream,
612  std::ostream & outStream,
613  const int order) {
614  using Finite_Difference_Arrays::shifts;
615  using Finite_Difference_Arrays::weights;
616 
617  const Real one(1);
618  Real tol = std::sqrt(ROL_EPSILON<Real>());
619 
620  int numSteps = steps.size();
621  int numVals = 4;
622  std::vector<Real> tmp(numVals);
623  std::vector<std::vector<Real> > ahuvCheck(numSteps, tmp);
624 
625  // Temporary vectors.
626  ROL::Ptr<Vector<Real> > AJdif = hv.clone();
627  ROL::Ptr<Vector<Real> > AJu = hv.clone();
628  ROL::Ptr<Vector<Real> > AHuv = hv.clone();
629  ROL::Ptr<Vector<Real> > AJnew = hv.clone();
630  ROL::Ptr<Vector<Real> > xnew = x.clone();
631 
632  // Save the format state of the original outStream.
633  ROL::nullstream oldFormatState;
634  oldFormatState.copyfmt(outStream);
635 
636  // Apply adjoint Jacobian to u.
637  this->update(x,UpdateType::Temp);
638  this->applyAdjointJacobian(*AJu, u, x, tol);
639 
640  // Apply adjoint Hessian at x, in direction v, to u.
641  this->applyAdjointHessian(*AHuv, u, v, x, tol);
642  Real normAHuv = AHuv->norm();
643 
644  for (int i=0; i<numSteps; i++) {
645 
646  Real eta = steps[i];
647 
648  // Apply adjoint Jacobian to u at x+eta*v.
649  xnew->set(x);
650 
651  AJdif->set(*AJu);
652  AJdif->scale(weights[order-1][0]);
653 
654  for(int j=0; j<order; ++j) {
655 
656  xnew->axpy(eta*shifts[order-1][j],v);
657 
658  if( weights[order-1][j+1] != 0 ) {
659  this->update(*xnew,UpdateType::Temp);
660  this->applyAdjointJacobian(*AJnew, u, *xnew, tol);
661  AJdif->axpy(weights[order-1][j+1],*AJnew);
662  }
663  }
664 
665  AJdif->scale(one/eta);
666 
667  // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
668  ahuvCheck[i][0] = eta;
669  ahuvCheck[i][1] = normAHuv;
670  ahuvCheck[i][2] = AJdif->norm();
671  AJdif->axpy(-one, *AHuv);
672  ahuvCheck[i][3] = AJdif->norm();
673 
674  if (printToStream) {
675  std::stringstream hist;
676  if (i==0) {
677  hist << std::right
678  << std::setw(20) << "Step size"
679  << std::setw(20) << "norm(adj(H)(u,v))"
680  << std::setw(20) << "norm(FD approx)"
681  << std::setw(20) << "norm(abs error)"
682  << "\n"
683  << std::setw(20) << "---------"
684  << std::setw(20) << "-----------------"
685  << std::setw(20) << "---------------"
686  << std::setw(20) << "---------------"
687  << "\n";
688  }
689  hist << std::scientific << std::setprecision(11) << std::right
690  << std::setw(20) << ahuvCheck[i][0]
691  << std::setw(20) << ahuvCheck[i][1]
692  << std::setw(20) << ahuvCheck[i][2]
693  << std::setw(20) << ahuvCheck[i][3]
694  << "\n";
695  outStream << hist.str();
696  }
697 
698  }
699 
700  // Reset format state of outStream.
701  outStream.copyfmt(oldFormatState);
702 
703  return ahuvCheck;
704 } // checkApplyAdjointHessian
705 
706 } // namespace ROL
707 
708 #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::Finite_Difference_Arrays::weights
const double weights[4][5]
Definition: ROL_Types.hpp:838
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:456
ROL::Finite_Difference_Arrays::shifts
const int shifts[4][4]
Definition: ROL_Types.hpp:831
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::details::nullstream
basic_nullstream< char, std::char_traits< char >> nullstream
Definition: ROL_Stream.hpp:36
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:606
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:351
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::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::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
Generated by   doxygen 1.8.5