ROL
Loading...
Searching...
No Matches
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 const bool breakdown = (H(i+1,i) == zero);
248 if (!breakdown) {
249 V1.push_back(b1.clone()); (V1[i+1])->set(*w1); (V1[i+1])->scale(one/H(i+1,i));
250 V2.push_back(b2.clone()); (V2[i+1])->set(*w2); (V2[i+1])->scale(one/H(i+1,i));
251 }
252
253 // Apply Givens rotations.
254 for (k=0; k<=i-1; k++) {
255 temp = cs(k)*H(k,i) + sn(k)*H(k+1,i);
256 H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
257 H(k,i) = temp;
258 }
259
260 // Form i-th rotation matrix.
261 if ( H(i+1,i) == zero ) {
262 cs(i) = one;
263 sn(i) = zero;
264 }
265 else if ( std::abs(H(i+1,i)) > std::abs(H(i,i)) ) {
266 temp = H(i,i) / H(i+1,i);
267 sn(i) = one / std::sqrt( one + temp*temp );
268 cs(i) = temp * sn(i);
269 }
270 else {
271 temp = H(i+1,i) / H(i,i);
272 cs(i) = one / std::sqrt( one + temp*temp );
273 sn(i) = temp * cs(i);
274 }
275
276 // Approximate residual norm.
277 temp = cs(i)*s(i);
278 s(i+1) = -sn(i)*s(i);
279 s(i) = temp;
280 H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
281 H(i+1,i) = zero;
282 resnrm = std::abs(s(i+1));
283 res[i+1] = resnrm;
284
285 // Update solution approximation.
286 const char uplo = 'U';
287 const char trans = 'N';
288 const char diag = 'N';
289 const char normin = 'N';
290 Real scaling = zero;
291 int info = 0;
292 y = s;
293 lapack.LATRS(uplo, trans, diag, normin, i+1, H.values(), m+1, y.values(), &scaling, cnorm.values(), &info);
294 z1->zero();
295 z2->zero();
296 for (k=0; k<=i; k++) {
297 z1->axpy(y(k), *(Z1[k]));
298 z2->axpy(y(k), *(Z2[k]));
299 }
300
301 if (res[i+1] <= tol || breakdown) {
302 break;
303 }
304
305 } // for (int i=0; i++; i<m)
306 v1.set(*z1); v2.set(*z2);
307
308 res.resize(i+2);
309
310 /*
311 std::stringstream hist;
312 hist << std::scientific << std::setprecision(8);
313 hist << "\n Augmented System Solver:\n";
314 hist << " Iter Residual\n";
315 for (unsigned j=0; j<res.size(); j++) {
316 hist << " " << std::left << std::setw(14) << res[j] << "\n";
317 }
318 hist << "\n";
319 std::cout << hist.str();
320 */
321
322 return res;
323}
324
325
326template <class Real>
327std::vector<std::vector<Real> > Constraint<Real>::checkApplyJacobian(const Vector<Real> &x,
328 const Vector<Real> &v,
329 const Vector<Real> &jv,
330 const bool printToStream,
331 std::ostream & outStream,
332 const int numSteps,
333 const int order) {
334 std::vector<Real> steps(numSteps);
335 for(int i=0;i<numSteps;++i) {
336 steps[i] = pow(10,-i);
337 }
338
339 return checkApplyJacobian(x,v,jv,steps,printToStream,outStream,order);
340}
341
342
343
344
345template <class Real>
346std::vector<std::vector<Real> > Constraint<Real>::checkApplyJacobian(const Vector<Real> &x,
347 const Vector<Real> &v,
348 const Vector<Real> &jv,
349 const std::vector<Real> &steps,
350 const bool printToStream,
351 std::ostream & outStream,
352 const int order) {
353 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
354 "Error: finite difference order must be 1,2,3, or 4" );
355
356 const Real one(1.0);
357
358 using Finite_Difference_Arrays::shifts;
359 using Finite_Difference_Arrays::weights;
360
361 Real tol = std::sqrt(ROL_EPSILON<Real>());
362
363 int numSteps = steps.size();
364 int numVals = 4;
365 std::vector<Real> tmp(numVals);
366 std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
367
368 // Save the format state of the original outStream.
369 ROL::nullstream oldFormatState;
370 oldFormatState.copyfmt(outStream);
371
372 // Compute constraint value at x.
373 ROL::Ptr<Vector<Real> > c = jv.clone();
374 this->update(x,UpdateType::Temp);
375 this->value(*c, x, tol);
376
377 // Compute (Jacobian at x) times (vector v).
378 ROL::Ptr<Vector<Real> > Jv = jv.clone();
379 this->applyJacobian(*Jv, v, x, tol);
380 Real normJv = Jv->norm();
381
382 // Temporary vectors.
383 ROL::Ptr<Vector<Real> > cdif = jv.clone();
384 ROL::Ptr<Vector<Real> > cnew = jv.clone();
385 ROL::Ptr<Vector<Real> > xnew = x.clone();
386
387 for (int i=0; i<numSteps; i++) {
388
389 Real eta = steps[i];
390
391 xnew->set(x);
392
393 cdif->set(*c);
394 cdif->scale(weights[order-1][0]);
395
396 for(int j=0; j<order; ++j) {
397
398 xnew->axpy(eta*shifts[order-1][j], v);
399
400 if( weights[order-1][j+1] != 0 ) {
401 this->update(*xnew,UpdateType::Temp);
402 this->value(*cnew,*xnew,tol);
403 cdif->axpy(weights[order-1][j+1],*cnew);
404 }
405
406 }
407
408 cdif->scale(one/eta);
409
410 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
411 jvCheck[i][0] = eta;
412 jvCheck[i][1] = normJv;
413 jvCheck[i][2] = cdif->norm();
414 cdif->axpy(-one, *Jv);
415 jvCheck[i][3] = cdif->norm();
416
417 if (printToStream) {
418 std::stringstream hist;
419 if (i==0) {
420 hist << std::right
421 << std::setw(20) << "Step size"
422 << std::setw(20) << "norm(Jac*vec)"
423 << std::setw(20) << "norm(FD approx)"
424 << std::setw(20) << "norm(abs error)"
425 << "\n"
426 << std::setw(20) << "---------"
427 << std::setw(20) << "-------------"
428 << std::setw(20) << "---------------"
429 << std::setw(20) << "---------------"
430 << "\n";
431 }
432 hist << std::scientific << std::setprecision(11) << std::right
433 << std::setw(20) << jvCheck[i][0]
434 << std::setw(20) << jvCheck[i][1]
435 << std::setw(20) << jvCheck[i][2]
436 << std::setw(20) << jvCheck[i][3]
437 << "\n";
438 outStream << hist.str();
439 }
440
441 }
442
443 // Reset format state of outStream.
444 outStream.copyfmt(oldFormatState);
445
446 return jvCheck;
447} // checkApplyJacobian
448
449
450template <class Real>
451std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointJacobian(const Vector<Real> &x,
452 const Vector<Real> &v,
453 const Vector<Real> &c,
454 const Vector<Real> &ajv,
455 const bool printToStream,
456 std::ostream & outStream,
457 const int numSteps) {
458 Real tol = std::sqrt(ROL_EPSILON<Real>());
459 const Real one(1);
460
461 int numVals = 4;
462 std::vector<Real> tmp(numVals);
463 std::vector<std::vector<Real> > ajvCheck(numSteps, tmp);
464 Real eta_factor = 1e-1;
465 Real eta = one;
466
467 // Temporary vectors.
468 ROL::Ptr<Vector<Real> > c0 = c.clone();
469 ROL::Ptr<Vector<Real> > cnew = c.clone();
470 ROL::Ptr<Vector<Real> > xnew = x.clone();
471 ROL::Ptr<Vector<Real> > ajv0 = ajv.clone();
472 ROL::Ptr<Vector<Real> > ajv1 = ajv.clone();
473 ROL::Ptr<Vector<Real> > ex = x.clone();
474 ROL::Ptr<Vector<Real> > eajv = ajv.clone();
475
476 // Save the format state of the original outStream.
477 ROL::nullstream oldFormatState;
478 oldFormatState.copyfmt(outStream);
479
480 // Compute constraint value at x.
481 this->update(x,UpdateType::Temp);
482 this->value(*c0, x, tol);
483
484 // Compute (Jacobian at x) times (vector v).
485 this->applyAdjointJacobian(*ajv0, v, x, tol);
486 Real normAJv = ajv0->norm();
487
488 for (int i=0; i<numSteps; i++) {
489
490 ajv1->zero();
491
492 for ( int j = 0; j < ajv.dimension(); j++ ) {
493 ex = x.basis(j);
494 eajv = ajv.basis(j);
495 xnew->set(x);
496 xnew->axpy(eta,*ex);
497 this->update(*xnew,UpdateType::Temp);
498 this->value(*cnew,*xnew,tol);
499 cnew->axpy(-one,*c0);
500 cnew->scale(one/eta);
501 //ajv1->axpy(cnew->dot(v.dual()),*eajv);
502 ajv1->axpy(cnew->apply(v),*eajv);
503 }
504
505 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
506 ajvCheck[i][0] = eta;
507 ajvCheck[i][1] = normAJv;
508 ajvCheck[i][2] = ajv1->norm();
509 ajv1->axpy(-one, *ajv0);
510 ajvCheck[i][3] = ajv1->norm();
511
512 if (printToStream) {
513 std::stringstream hist;
514 if (i==0) {
515 hist << std::right
516 << std::setw(20) << "Step size"
517 << std::setw(20) << "norm(adj(Jac)*vec)"
518 << std::setw(20) << "norm(FD approx)"
519 << std::setw(20) << "norm(abs error)"
520 << "\n"
521 << std::setw(20) << "---------"
522 << std::setw(20) << "------------------"
523 << std::setw(20) << "---------------"
524 << std::setw(20) << "---------------"
525 << "\n";
526 }
527 hist << std::scientific << std::setprecision(11) << std::right
528 << std::setw(20) << ajvCheck[i][0]
529 << std::setw(20) << ajvCheck[i][1]
530 << std::setw(20) << ajvCheck[i][2]
531 << std::setw(20) << ajvCheck[i][3]
532 << "\n";
533 outStream << hist.str();
534 }
535
536 // Update eta.
537 eta = eta*eta_factor;
538 }
539
540 // Reset format state of outStream.
541 outStream.copyfmt(oldFormatState);
542
543 return ajvCheck;
544} // checkApplyAdjointJacobian
545
546template <class Real>
547Real Constraint<Real>::checkAdjointConsistencyJacobian(const Vector<Real> &w,
548 const Vector<Real> &v,
549 const Vector<Real> &x,
550 const Vector<Real> &dualw,
551 const Vector<Real> &dualv,
552 const bool printToStream,
553 std::ostream & outStream) {
554 Real tol = ROL_EPSILON<Real>();
555
556 ROL::Ptr<Vector<Real> > Jv = dualw.clone();
557 ROL::Ptr<Vector<Real> > Jw = dualv.clone();
558
559 this->update(x,UpdateType::Temp);
560 applyJacobian(*Jv,v,x,tol);
561 applyAdjointJacobian(*Jw,w,x,tol);
562
563 //Real vJw = v.dot(Jw->dual());
564 Real vJw = v.apply(*Jw);
565 //Real wJv = w.dot(Jv->dual());
566 Real wJv = w.apply(*Jv);
567
568 Real diff = std::abs(wJv-vJw);
569
570 if ( printToStream ) {
571 std::stringstream hist;
572 hist << std::scientific << std::setprecision(8);
573 hist << "\nTest Consistency of Jacobian and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
574 << diff << "\n";
575 hist << " |<w,Jv>| = " << std::abs(wJv) << "\n";
576 hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n";
577 outStream << hist.str();
578 }
579 return diff;
580} // checkAdjointConsistencyJacobian
581
582template <class Real>
583std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointHessian(const Vector<Real> &x,
584 const Vector<Real> &u,
585 const Vector<Real> &v,
586 const Vector<Real> &hv,
587 const bool printToStream,
588 std::ostream & outStream,
589 const int numSteps,
590 const int order) {
591 std::vector<Real> steps(numSteps);
592 for(int i=0;i<numSteps;++i) {
593 steps[i] = pow(10,-i);
594 }
595
596 return checkApplyAdjointHessian(x,u,v,hv,steps,printToStream,outStream,order);
597}
598
599
600template <class Real>
601std::vector<std::vector<Real> > Constraint<Real>::checkApplyAdjointHessian(const Vector<Real> &x,
602 const Vector<Real> &u,
603 const Vector<Real> &v,
604 const Vector<Real> &hv,
605 const std::vector<Real> &steps,
606 const bool printToStream,
607 std::ostream & outStream,
608 const int order) {
609 using Finite_Difference_Arrays::shifts;
610 using Finite_Difference_Arrays::weights;
611
612 const Real one(1);
613 Real tol = std::sqrt(ROL_EPSILON<Real>());
614
615 int numSteps = steps.size();
616 int numVals = 4;
617 std::vector<Real> tmp(numVals);
618 std::vector<std::vector<Real> > ahuvCheck(numSteps, tmp);
619
620 // Temporary vectors.
621 ROL::Ptr<Vector<Real> > AJdif = hv.clone();
622 ROL::Ptr<Vector<Real> > AJu = hv.clone();
623 ROL::Ptr<Vector<Real> > AHuv = hv.clone();
624 ROL::Ptr<Vector<Real> > AJnew = hv.clone();
625 ROL::Ptr<Vector<Real> > xnew = x.clone();
626
627 // Save the format state of the original outStream.
628 ROL::nullstream oldFormatState;
629 oldFormatState.copyfmt(outStream);
630
631 // Apply adjoint Jacobian to u.
632 this->update(x,UpdateType::Temp);
633 this->applyAdjointJacobian(*AJu, u, x, tol);
634
635 // Apply adjoint Hessian at x, in direction v, to u.
636 this->applyAdjointHessian(*AHuv, u, v, x, tol);
637 Real normAHuv = AHuv->norm();
638
639 for (int i=0; i<numSteps; i++) {
640
641 Real eta = steps[i];
642
643 // Apply adjoint Jacobian to u at x+eta*v.
644 xnew->set(x);
645
646 AJdif->set(*AJu);
647 AJdif->scale(weights[order-1][0]);
648
649 for(int j=0; j<order; ++j) {
650
651 xnew->axpy(eta*shifts[order-1][j],v);
652
653 if( weights[order-1][j+1] != 0 ) {
654 this->update(*xnew,UpdateType::Temp);
655 this->applyAdjointJacobian(*AJnew, u, *xnew, tol);
656 AJdif->axpy(weights[order-1][j+1],*AJnew);
657 }
658 }
659
660 AJdif->scale(one/eta);
661
662 // Compute norms of Jacobian-vector products, finite-difference approximations, and error.
663 ahuvCheck[i][0] = eta;
664 ahuvCheck[i][1] = normAHuv;
665 ahuvCheck[i][2] = AJdif->norm();
666 AJdif->axpy(-one, *AHuv);
667 ahuvCheck[i][3] = AJdif->norm();
668
669 if (printToStream) {
670 std::stringstream hist;
671 if (i==0) {
672 hist << std::right
673 << std::setw(20) << "Step size"
674 << std::setw(20) << "norm(adj(H)(u,v))"
675 << std::setw(20) << "norm(FD approx)"
676 << std::setw(20) << "norm(abs error)"
677 << "\n"
678 << std::setw(20) << "---------"
679 << std::setw(20) << "-----------------"
680 << std::setw(20) << "---------------"
681 << std::setw(20) << "---------------"
682 << "\n";
683 }
684 hist << std::scientific << std::setprecision(11) << std::right
685 << std::setw(20) << ahuvCheck[i][0]
686 << std::setw(20) << ahuvCheck[i][1]
687 << std::setw(20) << ahuvCheck[i][2]
688 << std::setw(20) << ahuvCheck[i][3]
689 << "\n";
690 outStream << hist.str();
691 }
692
693 }
694
695 // Reset format state of outStream.
696 outStream.copyfmt(oldFormatState);
697
698 return ahuvCheck;
699} // checkApplyAdjointHessian
700
701} // namespace ROL
702
703#endif
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Definition ROL_Objective_SerialSimOpt.hpp:77
update
virtual void update(const Vector< Real > &u_old, const Vector< Real > &u_new, const Vector< Real > &z, bool flag=true, int iter=-1)
Definition ROL_Objective_TimeSimOpt.hpp:79
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.
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.
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 ,...
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 .
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::applyJacobian
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
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 .
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.
ROL::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::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::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:346
ROL::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:601
ROL::ROL::Constraint::applyPreconditioner
virtual void applyPreconditioner(Vector< Real > &pv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g, Real &tol)
Apply a constraint preconditioner at , , to vector . Ideally, this preconditioner satisfies the follo...
Definition ROL_Constraint_SerialSimOpt.hpp:249
ROL::ROL::Vector::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition ROL_Constraint_SerialSimOpt.hpp:176
ROL::ROL::Vector::basis
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition ROL_Constraint_SerialSimOpt.hpp:149
ROL::ROL::Vector::dot
virtual Real dot(const Vector &x) const =0
Compute where .
ROL::ROL::Vector::dimension
virtual int dimension() const
Return dimension of the vector space.
Definition ROL_Constraint_SerialSimOpt.hpp:163
ROL::ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition ROL_Constraint_SerialSimOpt.hpp:120
ROL::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_Constraint_SerialSimOpt.hpp:193
ROL::ROL::Vector::norm
virtual Real norm() const =0
Returns where .
ROL::ROL::Vector::apply
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition ROL_Constraint_SerialSimOpt.hpp:205
ROL::ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition ROL_Constraint_SerialSimOpt.hpp:134
ROL::ROL::Finite_Difference_Arrays::shifts
const int shifts[4][4]
Definition ROL_Constraint_SerialSimOpt.hpp:831
ROL::ROL::Finite_Difference_Arrays::weights
const double weights[4][5]
Definition ROL_Constraint_SerialSimOpt.hpp:838
ROL::ROL::details::nullstream
basic_nullstream< char, std::char_traits< char > > nullstream
Definition ROL_Constraint_SerialSimOpt.hpp:37
ROL::ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition ROL_Constraint_SerialSimOpt.hpp:57
ROL::value
ROL::Objective_SerialSimOpt Objective_SimOpt value(const V &u, const V &z, Real &tol) override
Definition ROL_Objective_SerialSimOpt.hpp:50
ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition ROL_Types.hpp:57
ROL::ROL_UNDERFLOW
Real ROL_UNDERFLOW(void)
Platform-dependent minimum double.
Definition ROL_Types.hpp:79
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.15.0