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