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