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