204 const Real
zero(0), one(1);
213 tol = std::sqrt(b1.
dot(b1)+b2.
dot(b2))*tol;
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();
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;
243 applyAdjointJacobian(*r1, v2, x, zerotol);
244 r1->scale(-one); r1->axpy(-one, v1.
dual()); r1->plus(b1);
246 r2->scale(-one); r2->plus(b2);
247 res[0] = std::sqrt(r1->dot(*r1) + r2->dot(*r2));
250 if (res[0] ==
zero) {
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]);
260 for (i=0; i<m; i++) {
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());
270 applyAdjointJacobian(*w1temp, *Z2temp, x, zerotol);
271 w1->set(*(V1[i])); w1->plus(*w1temp);
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]));
279 H(i+1,i) = std::sqrt(w1->dot(*w1) + w2->dot(*w2));
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));
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);
292 if ( H(i+1,i) ==
zero ) {
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);
302 temp = H(i+1,i) / H(i,i);
303 cs(i) = one / std::sqrt( one + temp*temp );
304 sn(i) = temp * cs(i);
309 s(i+1) = -sn(i)*s(i);
311 H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
313 resnrm = std::abs(s(i+1));
317 const char uplo =
'U';
318 const char trans =
'N';
319 const char diag =
'N';
320 const char normin =
'N';
324 lapack.LATRS(uplo, trans, diag, normin, i+1, H.values(), m+1, y.values(), &scaling, cnorm.values(), &info);
327 for (k=0; k<=i; k++) {
328 z1->axpy(y(k), *(Z1[k]));
329 z2->axpy(y(k), *(Z2[k]));
336 if (res[i+1] <= tol) {
387 const std::vector<Real> &steps,
388 const bool printToStream,
389 std::ostream & outStream,
391 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
392 "Error: finite difference order must be 1,2,3, or 4" );
399 Real tol = std::sqrt(ROL_EPSILON<Real>());
401 int numSteps = steps.size();
403 std::vector<Real> tmp(numVals);
404 std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
407 ROL::nullstream oldFormatState;
408 oldFormatState.copyfmt(outStream);
411 ROL::Ptr<Vector<Real> > c = jv.
clone();
413 this->
value(*c, x, tol);
416 ROL::Ptr<Vector<Real> > Jv = jv.
clone();
418 Real normJv = Jv->norm();
421 ROL::Ptr<Vector<Real> > cdif = jv.
clone();
422 ROL::Ptr<Vector<Real> > cnew = jv.
clone();
423 ROL::Ptr<Vector<Real> > xnew = x.
clone();
425 for (
int i=0; i<numSteps; i++) {
432 cdif->scale(weights[order-1][0]);
434 for(
int j=0; j<order; ++j) {
436 xnew->axpy(eta*shifts[order-1][j], v);
438 if( weights[order-1][j+1] != 0 ) {
440 this->
value(*cnew,*xnew,tol);
441 cdif->axpy(weights[order-1][j+1],*cnew);
446 cdif->scale(one/eta);
450 jvCheck[i][1] = normJv;
451 jvCheck[i][2] = cdif->norm();
452 cdif->axpy(-one, *Jv);
453 jvCheck[i][3] = cdif->norm();
456 std::stringstream hist;
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)"
464 << std::setw(20) <<
"---------"
465 << std::setw(20) <<
"-------------"
466 << std::setw(20) <<
"---------------"
467 << std::setw(20) <<
"---------------"
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]
476 outStream << hist.str();
482 outStream.copyfmt(oldFormatState);
493 const bool printToStream,
494 std::ostream & outStream,
495 const int numSteps) {
496 Real tol = std::sqrt(ROL_EPSILON<Real>());
500 std::vector<Real> tmp(numVals);
501 std::vector<std::vector<Real> > ajvCheck(numSteps, tmp);
502 Real eta_factor = 1e-1;
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();
515 ROL::nullstream oldFormatState;
516 oldFormatState.copyfmt(outStream);
520 this->
value(*c0, x, tol);
523 this->applyAdjointJacobian(*ajv0, v, x, tol);
524 Real normAJv = ajv0->norm();
526 for (
int i=0; i<numSteps; i++) {
530 for (
int j = 0; j < ajv.
dimension(); j++ ) {
536 this->
value(*cnew,*xnew,tol);
537 cnew->axpy(-one,*c0);
538 cnew->scale(one/eta);
540 ajv1->axpy(cnew->apply(v),*eajv);
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();
551 std::stringstream hist;
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)"
559 << std::setw(20) <<
"---------"
560 << std::setw(20) <<
"------------------"
561 << std::setw(20) <<
"---------------"
562 << std::setw(20) <<
"---------------"
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]
571 outStream << hist.str();
575 eta = eta*eta_factor;
579 outStream.copyfmt(oldFormatState);
643 const std::vector<Real> &steps,
644 const bool printToStream,
645 std::ostream & outStream,
651 Real tol = std::sqrt(ROL_EPSILON<Real>());
653 int numSteps = steps.size();
655 std::vector<Real> tmp(numVals);
656 std::vector<std::vector<Real> > ahuvCheck(numSteps, tmp);
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();
666 ROL::nullstream oldFormatState;
667 oldFormatState.copyfmt(outStream);
671 this->applyAdjointJacobian(*AJu, u, x, tol);
674 this->applyAdjointHessian(*AHuv, u, v, x, tol);
675 Real normAHuv = AHuv->norm();
677 for (
int i=0; i<numSteps; i++) {
685 AJdif->scale(weights[order-1][0]);
687 for(
int j=0; j<order; ++j) {
689 xnew->axpy(eta*shifts[order-1][j],v);
691 if( weights[order-1][j+1] != 0 ) {
693 this->applyAdjointJacobian(*AJnew, u, *xnew, tol);
694 AJdif->axpy(weights[order-1][j+1],*AJnew);
698 AJdif->scale(one/eta);
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();
708 std::stringstream hist;
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)"
716 << std::setw(20) <<
"---------"
717 << std::setw(20) <<
"-----------------"
718 << std::setw(20) <<
"---------------"
719 << std::setw(20) <<
"---------------"
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]
728 outStream << hist.str();
734 outStream.copyfmt(oldFormatState);
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.