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
51namespace ROL {
52
53template <class Real>
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
86template <class Real>
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
98template <class Real>
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>
135void 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
166template <class Real>
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);
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
195template <class Real>
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
364template <class Real>
365std::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
383template <class Real>
384std::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
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
488template <class Real>
489std::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
584template <class Real>
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
620template <class Real>
621std::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
638template <class Real>
639std::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) {
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
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
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 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 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 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 Real checkAdjointConsistencyJacobian(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &x, const bool printToStream=true, std::ostream &outStream=std::cout)
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.
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
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: ROL_Vector.hpp:238
virtual Real norm() const =0
Returns where .
virtual void scale(const Real alpha)=0
Compute where .
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
virtual void plus(const Vector &x)=0
Compute , where .
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:196
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:182
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
virtual Real dot(const Vector &x) const =0
Compute where .
const double weights[4][5]
Definition: ROL_Types.hpp:868
ROL::Objective_SerialSimOpt Objective_SimOpt value(const V &u, const V &z, Real &tol) override
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1) override
void applyJacobian(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &sol)