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