45 #ifndef ROL_CONSTRAINT_DEF_H
46 #define ROL_CONSTRAINT_DEF_H
48 #include "ROL_LinearAlgebra.hpp"
49 #include "ROL_LAPACK.hpp"
60 Real ctol = std::sqrt(ROL_EPSILON<Real>());
63 Real h = std::max(one,x.
norm()/v.
norm())*tol;
67 ROL::Ptr<Vector<Real> > c = jv.
clone();
68 this->
value(*c,x,ctol);
71 ROL::Ptr<Vector<Real> > xnew = x.
clone();
78 this->
value(jv,*xnew,ctol);
91 applyAdjointJacobian(ajv,v,x,v.
dual(),tol);
108 Real h(0), ctol = std::sqrt(ROL_EPSILON<Real>());
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();
114 ROL::Ptr<Vector<Real> > c0 = dualv.
clone();
115 this->
value(*c0,x,ctol);
118 for (
int i = 0; i < ajv.
dimension(); i++ ) {
121 h = std::max(one,x.
norm()/ex->norm())*tol;
125 this->
value(*cnew,*xnew,ctol);
126 cnew->axpy(-one,*c0);
128 ajv.
axpy(cnew->dot(v.
dual()),*eajv);
165 template <
class Real>
172 Real h = std::max(static_cast<Real>(1),x.
norm()/v.
norm())*tol;
175 ROL::Ptr<Vector<Real> > aju = huv.
clone();
176 applyAdjointJacobian(*aju,u,x,tol);
179 ROL::Ptr<Vector<Real> > xnew = x.
clone();
186 applyAdjointJacobian(huv,u,*xnew,tol);
189 huv.
axpy(static_cast<Real>(-1),*aju);
190 huv.
scale(static_cast<Real>(1)/h);
194 template <
class Real>
203 const Real
zero(0), one(1);
212 tol = std::sqrt(b1.
dot(b1)+b2.
dot(b2))*tol;
218 ROL::Ptr<Vector<Real> > r1 = b1.
clone();
219 ROL::Ptr<Vector<Real> > r2 = b2.
clone();
220 ROL::Ptr<Vector<Real> > z1 = v1.
clone();
221 ROL::Ptr<Vector<Real> > z2 = v2.
clone();
222 ROL::Ptr<Vector<Real> > w1 = b1.
clone();
223 ROL::Ptr<Vector<Real> > w2 = b2.
clone();
224 std::vector<ROL::Ptr<Vector<Real> > > V1;
225 std::vector<ROL::Ptr<Vector<Real> > > V2;
226 ROL::Ptr<Vector<Real> > V2temp = b2.
clone();
227 std::vector<ROL::Ptr<Vector<Real> > > Z1;
228 std::vector<ROL::Ptr<Vector<Real> > > Z2;
229 ROL::Ptr<Vector<Real> > w1temp = b1.
clone();
230 ROL::Ptr<Vector<Real> > Z2temp = v2.
clone();
232 std::vector<Real> res(m+1,
zero);
233 LA::Matrix<Real> H(m+1,m);
234 LA::Vector<Real> cs(m);
235 LA::Vector<Real> sn(m);
236 LA::Vector<Real> s(m+1);
237 LA::Vector<Real> y(m+1);
238 LA::Vector<Real> cnorm(m);
239 ROL::LAPACK<int, Real> lapack;
242 applyAdjointJacobian(*r1, v2, x, zerotol);
243 r1->scale(-one); r1->axpy(-one, v1.
dual()); r1->plus(b1);
245 r2->scale(-one); r2->plus(b2);
246 res[0] = std::sqrt(r1->dot(*r1) + r2->dot(*r2));
249 if (res[0] ==
zero) {
254 V1.push_back(b1.
clone()); (V1[0])->set(*r1); (V1[0])->scale(one/res[0]);
255 V2.push_back(b2.
clone()); (V2[0])->set(*r2); (V2[0])->scale(one/res[0]);
259 for (i=0; i<m; i++) {
262 V2temp->set(*(V2[i]));
263 applyPreconditioner(*Z2temp, *V2temp, x, b1, zerotol);
264 Z2.push_back(v2.
clone()); (Z2[i])->set(*Z2temp);
265 Z1.push_back(v1.
clone()); (Z1[i])->set((V1[i])->dual());
269 applyAdjointJacobian(*w1temp, *Z2temp, x, zerotol);
270 w1->set(*(V1[i])); w1->plus(*w1temp);
273 for (k=0; k<=i; k++) {
274 H(k,i) = w1->dot(*(V1[k])) + w2->dot(*(V2[k]));
275 w1->axpy(-H(k,i), *(V1[k]));
276 w2->axpy(-H(k,i), *(V2[k]));
278 H(i+1,i) = std::sqrt(w1->dot(*w1) + w2->dot(*w2));
280 V1.push_back(b1.
clone()); (V1[i+1])->set(*w1); (V1[i+1])->scale(one/H(i+1,i));
281 V2.push_back(b2.
clone()); (V2[i+1])->set(*w2); (V2[i+1])->scale(one/H(i+1,i));
284 for (k=0; k<=i-1; k++) {
285 temp = cs(k)*H(k,i) + sn(k)*H(k+1,i);
286 H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
291 if ( H(i+1,i) ==
zero ) {
295 else if ( std::abs(H(i+1,i)) > std::abs(H(i,i)) ) {
296 temp = H(i,i) / H(i+1,i);
297 sn(i) = one / std::sqrt( one + temp*temp );
298 cs(i) = temp * sn(i);
301 temp = H(i+1,i) / H(i,i);
302 cs(i) = one / std::sqrt( one + temp*temp );
303 sn(i) = temp * cs(i);
308 s(i+1) = -sn(i)*s(i);
310 H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
312 resnrm = std::abs(s(i+1));
316 const char uplo =
'U';
317 const char trans =
'N';
318 const char diag =
'N';
319 const char normin =
'N';
323 lapack.LATRS(uplo, trans, diag, normin, i+1, H.values(), m+1, y.values(), &scaling, cnorm.values(), &info);
326 for (k=0; k<=i; k++) {
327 z1->axpy(y(k), *(Z1[k]));
328 z2->axpy(y(k), *(Z2[k]));
334 std::cout <<
" " << i+1 <<
": " << res[i+1]/res[0] << std::endl;
335 if (res[i+1] <= tol) {
336 std::cout <<
" solved in " << i+1 <<
" iterations to " << res[i+1] <<
" (" << res[i+1]/res[0] <<
")" << std::endl;
363 template <
class Real>
367 const bool printToStream,
368 std::ostream & outStream,
371 std::vector<Real> steps(numSteps);
372 for(
int i=0;i<numSteps;++i) {
373 steps[i] = pow(10,-i);
376 return checkApplyJacobian(x,v,jv,steps,printToStream,outStream,order);
382 template <
class Real>
386 const std::vector<Real> &steps,
387 const bool printToStream,
388 std::ostream & outStream,
390 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
391 "Error: finite difference order must be 1,2,3, or 4" );
398 Real tol = std::sqrt(ROL_EPSILON<Real>());
400 int numSteps = steps.size();
402 std::vector<Real> tmp(numVals);
403 std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
407 oldFormatState.copyfmt(outStream);
410 ROL::Ptr<Vector<Real> > c = jv.
clone();
412 this->
value(*c, x, tol);
415 ROL::Ptr<Vector<Real> > Jv = jv.
clone();
417 Real normJv = Jv->norm();
420 ROL::Ptr<Vector<Real> > cdif = jv.
clone();
421 ROL::Ptr<Vector<Real> > cnew = jv.
clone();
422 ROL::Ptr<Vector<Real> > xnew = x.
clone();
424 for (
int i=0; i<numSteps; i++) {
431 cdif->scale(
weights[order-1][0]);
433 for(
int j=0; j<order; ++j) {
435 xnew->axpy(eta*
shifts[order-1][j], v);
437 if(
weights[order-1][j+1] != 0 ) {
439 this->
value(*cnew,*xnew,tol);
440 cdif->axpy(
weights[order-1][j+1],*cnew);
445 cdif->scale(one/eta);
449 jvCheck[i][1] = normJv;
450 jvCheck[i][2] = cdif->norm();
451 cdif->axpy(-one, *Jv);
452 jvCheck[i][3] = cdif->norm();
455 std::stringstream hist;
458 << std::setw(20) <<
"Step size"
459 << std::setw(20) <<
"norm(Jac*vec)"
460 << std::setw(20) <<
"norm(FD approx)"
461 << std::setw(20) <<
"norm(abs error)"
463 << std::setw(20) <<
"---------"
464 << std::setw(20) <<
"-------------"
465 << std::setw(20) <<
"---------------"
466 << std::setw(20) <<
"---------------"
469 hist << std::scientific << std::setprecision(11) << std::right
470 << std::setw(20) << jvCheck[i][0]
471 << std::setw(20) << jvCheck[i][1]
472 << std::setw(20) << jvCheck[i][2]
473 << std::setw(20) << jvCheck[i][3]
475 outStream << hist.str();
481 outStream.copyfmt(oldFormatState);
487 template <
class Real>
492 const bool printToStream,
493 std::ostream & outStream,
494 const int numSteps) {
495 Real tol = std::sqrt(ROL_EPSILON<Real>());
499 std::vector<Real> tmp(numVals);
500 std::vector<std::vector<Real> > ajvCheck(numSteps, tmp);
501 Real eta_factor = 1e-1;
505 ROL::Ptr<Vector<Real> > c0 = c.
clone();
506 ROL::Ptr<Vector<Real> > cnew = c.
clone();
507 ROL::Ptr<Vector<Real> > xnew = x.
clone();
508 ROL::Ptr<Vector<Real> > ajv0 = ajv.
clone();
509 ROL::Ptr<Vector<Real> > ajv1 = ajv.
clone();
510 ROL::Ptr<Vector<Real> > ex = x.
clone();
511 ROL::Ptr<Vector<Real> > eajv = ajv.
clone();
515 oldFormatState.copyfmt(outStream);
519 this->
value(*c0, x, tol);
522 this->applyAdjointJacobian(*ajv0, v, x, tol);
523 Real normAJv = ajv0->norm();
525 for (
int i=0; i<numSteps; i++) {
529 for (
int j = 0; j < ajv.
dimension(); j++ ) {
535 this->
value(*cnew,*xnew,tol);
536 cnew->axpy(-one,*c0);
537 cnew->scale(one/eta);
538 ajv1->axpy(cnew->dot(v.
dual()),*eajv);
542 ajvCheck[i][0] = eta;
543 ajvCheck[i][1] = normAJv;
544 ajvCheck[i][2] = ajv1->norm();
545 ajv1->axpy(-one, *ajv0);
546 ajvCheck[i][3] = ajv1->norm();
549 std::stringstream hist;
552 << std::setw(20) <<
"Step size"
553 << std::setw(20) <<
"norm(adj(Jac)*vec)"
554 << std::setw(20) <<
"norm(FD approx)"
555 << std::setw(20) <<
"norm(abs error)"
557 << std::setw(20) <<
"---------"
558 << std::setw(20) <<
"------------------"
559 << std::setw(20) <<
"---------------"
560 << std::setw(20) <<
"---------------"
563 hist << std::scientific << std::setprecision(11) << std::right
564 << std::setw(20) << ajvCheck[i][0]
565 << std::setw(20) << ajvCheck[i][1]
566 << std::setw(20) << ajvCheck[i][2]
567 << std::setw(20) << ajvCheck[i][3]
569 outStream << hist.str();
573 eta = eta*eta_factor;
577 outStream.copyfmt(oldFormatState);
582 template <
class Real>
588 const bool printToStream,
589 std::ostream & outStream) {
590 Real tol = ROL_EPSILON<Real>();
592 ROL::Ptr<Vector<Real> > Jv = dualw.
clone();
593 ROL::Ptr<Vector<Real> > Jw = dualv.
clone();
596 applyAdjointJacobian(*Jw,w,x,tol);
598 Real vJw = v.
dot(Jw->dual());
599 Real wJv = w.
dot(Jv->dual());
601 Real diff = std::abs(wJv-vJw);
603 if ( printToStream ) {
604 std::stringstream hist;
605 hist << std::scientific << std::setprecision(8);
606 hist <<
"\nTest Consistency of Jacobian and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
608 hist <<
" |<w,Jv>| = " << std::abs(wJv) <<
"\n";
609 hist <<
" Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) <<
"\n";
610 outStream << hist.str();
615 template <
class Real>
620 const bool printToStream,
621 std::ostream & outStream,
624 std::vector<Real> steps(numSteps);
625 for(
int i=0;i<numSteps;++i) {
626 steps[i] = pow(10,-i);
629 return checkApplyAdjointHessian(x,u,v,hv,steps,printToStream,outStream,order);
633 template <
class Real>
638 const std::vector<Real> &steps,
639 const bool printToStream,
640 std::ostream & outStream,
646 Real tol = std::sqrt(ROL_EPSILON<Real>());
648 int numSteps = steps.size();
650 std::vector<Real> tmp(numVals);
651 std::vector<std::vector<Real> > ahuvCheck(numSteps, tmp);
654 ROL::Ptr<Vector<Real> > AJdif = hv.
clone();
655 ROL::Ptr<Vector<Real> > AJu = hv.
clone();
656 ROL::Ptr<Vector<Real> > AHuv = hv.
clone();
657 ROL::Ptr<Vector<Real> > AJnew = hv.
clone();
658 ROL::Ptr<Vector<Real> > xnew = x.
clone();
662 oldFormatState.copyfmt(outStream);
666 this->applyAdjointJacobian(*AJu, u, x, tol);
669 this->applyAdjointHessian(*AHuv, u, v, x, tol);
670 Real normAHuv = AHuv->norm();
672 for (
int i=0; i<numSteps; i++) {
680 AJdif->scale(
weights[order-1][0]);
682 for(
int j=0; j<order; ++j) {
684 xnew->axpy(eta*
shifts[order-1][j],v);
686 if(
weights[order-1][j+1] != 0 ) {
688 this->applyAdjointJacobian(*AJnew, u, *xnew, tol);
689 AJdif->axpy(
weights[order-1][j+1],*AJnew);
693 AJdif->scale(one/eta);
696 ahuvCheck[i][0] = eta;
697 ahuvCheck[i][1] = normAHuv;
698 ahuvCheck[i][2] = AJdif->norm();
699 AJdif->axpy(-one, *AHuv);
700 ahuvCheck[i][3] = AJdif->norm();
703 std::stringstream hist;
706 << std::setw(20) <<
"Step size"
707 << std::setw(20) <<
"norm(adj(H)(u,v))"
708 << std::setw(20) <<
"norm(FD approx)"
709 << std::setw(20) <<
"norm(abs error)"
711 << std::setw(20) <<
"---------"
712 << std::setw(20) <<
"-----------------"
713 << std::setw(20) <<
"---------------"
714 << std::setw(20) <<
"---------------"
717 hist << std::scientific << std::setprecision(11) << std::right
718 << std::setw(20) << ahuvCheck[i][0]
719 << std::setw(20) << ahuvCheck[i][1]
720 << std::setw(20) << ahuvCheck[i][2]
721 << std::setw(20) << ahuvCheck[i][3]
723 outStream << hist.str();
729 outStream.copyfmt(oldFormatState);
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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 int dimension() const
Return dimension of the vector space.
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
virtual void plus(const Vector &x)=0
Compute , where .
const double weights[4][5]
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1) override
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::Objective_SimOpt value
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
virtual Real dot(const Vector &x) const =0
Compute where .
void applyJacobian(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &sol)
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
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 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< 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 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 .
basic_nullstream< char, char_traits< char >> nullstream
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...
virtual Real norm() const =0
Returns where .