45 #ifndef ROL_EQUALITYCONSTRAINT_DEF_H 46 #define ROL_EQUALITYCONSTRAINT_DEF_H 61 Real h = std::max(1.0,x.
norm()/v.
norm())*tol;
65 Teuchos::RCP<Vector<Real> > c = jv.
clone();
66 this->value(*c,x,ctol);
69 Teuchos::RCP<Vector<Real> > xnew = x.
clone();
76 this->value(jv,*xnew,ctol);
89 applyAdjointJacobian(ajv,v,x,v.
dual(),tol);
110 Teuchos::RCP<Vector<Real> > xnew = x.
clone();
111 Teuchos::RCP<Vector<Real> > ex = x.
clone();
112 Teuchos::RCP<Vector<Real> > eajv = ajv.
clone();
113 Teuchos::RCP<Vector<Real> > cnew = dualv.
clone();
114 Teuchos::RCP<Vector<Real> > c0 = dualv.
clone();
115 this->value(*c0,x,ctol);
118 for (
int i = 0; i < ajv.
dimension(); i++ ) {
121 h = std::max(1.0,x.
norm()/ex->norm())*tol;
125 this->value(*cnew,*xnew,ctol);
126 cnew->axpy(-1.0,*c0);
128 ajv.
axpy(cnew->dot(v.
dual()),*eajv);
165 template <
class Real>
174 template <
class Real>
193 tol = std::sqrt(b1.
dot(b1)+b2.
dot(b2))*tol;
199 Teuchos::RCP<Vector<Real> > r1 = b1.
clone();
200 Teuchos::RCP<Vector<Real> > r2 = b2.
clone();
201 Teuchos::RCP<Vector<Real> > z1 = v1.
clone();
202 Teuchos::RCP<Vector<Real> > z2 = v2.
clone();
203 Teuchos::RCP<Vector<Real> > w1 = b1.
clone();
204 Teuchos::RCP<Vector<Real> > w2 = b2.
clone();
205 std::vector<Teuchos::RCP<Vector<Real> > > V1;
206 std::vector<Teuchos::RCP<Vector<Real> > > V2;
207 Teuchos::RCP<Vector<Real> > V2temp = b2.
clone();
208 std::vector<Teuchos::RCP<Vector<Real> > > Z1;
209 std::vector<Teuchos::RCP<Vector<Real> > > Z2;
210 Teuchos::RCP<Vector<Real> > w1temp = b1.
clone();
211 Teuchos::RCP<Vector<Real> > Z2temp = v2.
clone();
213 std::vector<Real> res(m+1, zero);
214 Teuchos::SerialDenseMatrix<int, Real> H(m+1,m);
215 Teuchos::SerialDenseVector<int, Real> cs(m);
216 Teuchos::SerialDenseVector<int, Real> sn(m);
217 Teuchos::SerialDenseVector<int, Real> s(m+1);
218 Teuchos::SerialDenseVector<int, Real> y(m+1);
219 Teuchos::SerialDenseVector<int, Real> cnorm(m);
220 Teuchos::LAPACK<int, Real> lapack;
223 applyAdjointJacobian(*r1, v2, x, zerotol);
224 r1->scale(-one); r1->axpy(-one, v1.
dual()); r1->plus(b1);
225 applyJacobian(*r2, v1, x, zerotol);
226 r2->scale(-one); r2->plus(b2);
227 res[0] = std::sqrt(r1->dot(*r1) + r2->dot(*r2));
230 if (res[0] == zero) {
235 V1.push_back(b1.
clone()); (V1[0])->
set(*r1); (V1[0])->scale(one/res[0]);
236 V2.push_back(b2.
clone()); (V2[0])->
set(*r2); (V2[0])->scale(one/res[0]);
240 for (i=0; i<m; i++) {
243 V2temp->set(*(V2[i]));
244 applyPreconditioner(*Z2temp, *V2temp, x, b1, zerotol);
245 Z2.push_back(v2.
clone()); (Z2[i])->
set(*Z2temp);
246 Z1.push_back(v1.
clone()); (Z1[i])->
set((V1[i])->dual());
249 applyJacobian(*w2, *(Z1[i]), x, zerotol);
250 applyAdjointJacobian(*w1temp, *Z2temp, x, zerotol);
251 w1->set(*(V1[i])); w1->plus(*w1temp);
254 for (k=0; k<=i; k++) {
255 H(k,i) = w1->dot(*(V1[k])) + w2->dot(*(V2[k]));
256 w1->axpy(-H(k,i), *(V1[k]));
257 w2->axpy(-H(k,i), *(V2[k]));
259 H(i+1,i) = std::sqrt(w1->dot(*w1) + w2->dot(*w2));
261 V1.push_back(b1.
clone()); (V1[i+1])->
set(*w1); (V1[i+1])->scale(one/H(i+1,i));
262 V2.push_back(b2.
clone()); (V2[i+1])->
set(*w2); (V2[i+1])->scale(one/H(i+1,i));
265 for (k=0; k<=i-1; k++) {
266 temp = cs(k)*H(k,i) + sn(k)*H(k+1,i);
267 H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
272 if ( H(i+1,i) == zero ) {
276 else if ( std::abs(H(i+1,i)) > std::abs(H(i,i)) ) {
277 temp = H(i,i) / H(i+1,i);
278 sn(i) = one / std::sqrt( one + temp*temp );
279 cs(i) = temp * sn(i);
282 temp = H(i+1,i) / H(i,i);
283 cs(i) = one / std::sqrt( one + temp*temp );
284 sn(i) = temp * cs(i);
289 s(i+1) = -sn(i)*s(i);
291 H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
293 resnrm = std::abs(s(i+1));
297 const char uplo =
'U';
298 const char trans =
'N';
299 const char diag =
'N';
300 const char normin =
'N';
304 lapack.LATRS(uplo, trans, diag, normin, i+1, H.values(), m+1, y.values(), &scaling, cnorm.values(), &info);
307 for (k=0; k<=i; k++) {
308 z1->axpy(y(k), *(Z1[k]));
309 z2->axpy(y(k), *(Z2[k]));
315 if (res[i+1] <= tol) {
343 template <
class Real>
347 const bool printToStream,
348 std::ostream & outStream,
351 std::vector<Real> steps(numSteps);
352 for(
int i=0;i<numSteps;++i) {
353 steps[i] = pow(10,-i);
356 return checkApplyJacobian(x,v,jv,steps,printToStream,outStream,order);
362 template <
class Real>
366 const std::vector<Real> &steps,
367 const bool printToStream,
368 std::ostream & outStream,
371 TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
372 "Error: finite difference order must be 1,2,3, or 4" );
379 int numSteps = steps.size();
381 std::vector<Real> tmp(numVals);
382 std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
385 Teuchos::oblackholestream oldFormatState;
386 oldFormatState.copyfmt(outStream);
389 Teuchos::RCP<Vector<Real> > c = jv.
clone();
391 this->value(*c, x, tol);
394 Teuchos::RCP<Vector<Real> > Jv = jv.
clone();
395 this->applyJacobian(*Jv, v, x, tol);
396 Real normJv = Jv->norm();
399 Teuchos::RCP<Vector<Real> > cdif = jv.
clone();
400 Teuchos::RCP<Vector<Real> > cnew = jv.
clone();
401 Teuchos::RCP<Vector<Real> > xnew = x.
clone();
403 for (
int i=0; i<numSteps; i++) {
410 cdif->scale(
weights[order-1][0]);
412 for(
int j=0; j<order; ++j) {
414 xnew->axpy(eta*
shifts[order-1][j], v);
416 if(
weights[order-1][j+1] != 0 ) {
418 this->value(*cnew,*xnew,tol);
419 cdif->axpy(
weights[order-1][j+1],*cnew);
424 cdif->scale(1.0/eta);
428 jvCheck[i][1] = normJv;
429 jvCheck[i][2] = cdif->norm();
430 cdif->axpy(-1.0, *Jv);
431 jvCheck[i][3] = cdif->norm();
434 std::stringstream hist;
437 << std::setw(20) <<
"Step size" 438 << std::setw(20) <<
"norm(Jac*vec)" 439 << std::setw(20) <<
"norm(FD approx)" 440 << std::setw(20) <<
"norm(abs error)" 442 << std::setw(20) <<
"---------" 443 << std::setw(20) <<
"-------------" 444 << std::setw(20) <<
"---------------" 445 << std::setw(20) <<
"---------------" 448 hist << std::scientific << std::setprecision(11) << std::right
449 << std::setw(20) << jvCheck[i][0]
450 << std::setw(20) << jvCheck[i][1]
451 << std::setw(20) << jvCheck[i][2]
452 << std::setw(20) << jvCheck[i][3]
454 outStream << hist.str();
460 outStream.copyfmt(oldFormatState);
466 template <
class Real>
471 const bool printToStream,
472 std::ostream & outStream,
473 const int numSteps) {
477 std::vector<Real> tmp(numVals);
478 std::vector<std::vector<Real> > ajvCheck(numSteps, tmp);
479 Real eta_factor = 1e-1;
483 Teuchos::RCP<Vector<Real> > c0 = c.
clone();
484 Teuchos::RCP<Vector<Real> > cnew = c.
clone();
485 Teuchos::RCP<Vector<Real> > xnew = x.
clone();
486 Teuchos::RCP<Vector<Real> > ajv0 = ajv.
clone();
487 Teuchos::RCP<Vector<Real> > ajv1 = ajv.
clone();
488 Teuchos::RCP<Vector<Real> > ex = x.
clone();
489 Teuchos::RCP<Vector<Real> > eajv = ajv.
clone();
492 Teuchos::oblackholestream oldFormatState;
493 oldFormatState.copyfmt(outStream);
497 this->value(*c0, x, tol);
500 this->applyAdjointJacobian(*ajv0, v, x, tol);
501 Real normAJv = ajv0->norm();
503 for (
int i=0; i<numSteps; i++) {
507 for (
int j = 0; j < ajv.
dimension(); j++ ) {
513 this->value(*cnew,*xnew,tol);
514 cnew->axpy(-1.0,*c0);
515 cnew->scale(1.0/eta);
516 ajv1->axpy(cnew->dot(v.
dual()),*eajv);
520 ajvCheck[i][0] = eta;
521 ajvCheck[i][1] = normAJv;
522 ajvCheck[i][2] = ajv1->norm();
523 ajv1->axpy(-1.0, *ajv0);
524 ajvCheck[i][3] = ajv1->norm();
527 std::stringstream hist;
530 << std::setw(20) <<
"Step size" 531 << std::setw(20) <<
"norm(adj(Jac)*vec)" 532 << std::setw(20) <<
"norm(FD approx)" 533 << std::setw(20) <<
"norm(abs error)" 535 << std::setw(20) <<
"---------" 536 << std::setw(20) <<
"------------------" 537 << std::setw(20) <<
"---------------" 538 << std::setw(20) <<
"---------------" 541 hist << std::scientific << std::setprecision(11) << std::right
542 << std::setw(20) << ajvCheck[i][0]
543 << std::setw(20) << ajvCheck[i][1]
544 << std::setw(20) << ajvCheck[i][2]
545 << std::setw(20) << ajvCheck[i][3]
547 outStream << hist.str();
551 eta = eta*eta_factor;
555 outStream.copyfmt(oldFormatState);
560 template <
class Real>
566 const bool printToStream,
567 std::ostream & outStream) {
570 Teuchos::RCP<Vector<Real> > Jv = dualw.
clone();
571 Teuchos::RCP<Vector<Real> > Jw = dualv.
clone();
573 applyJacobian(*Jv,v,x,tol);
574 applyAdjointJacobian(*Jw,w,x,tol);
576 Real vJw = v.
dot(Jw->dual());
577 Real wJv = w.
dot(Jv->dual());
579 Real diff = std::abs(wJv-vJw);
581 if ( printToStream ) {
582 std::stringstream hist;
583 hist << std::scientific << std::setprecision(8);
584 hist <<
"\nTest Consistency of Jacobian and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = " 586 hist <<
" |<w,Jv>| = " << std::abs(wJv) <<
"\n";
587 hist <<
" Relative Error = " << diff / (std::abs(wJv)+
ROL_UNDERFLOW) <<
"\n";
588 outStream << hist.str();
593 template <
class Real>
598 const bool printToStream,
599 std::ostream & outStream,
602 std::vector<Real> steps(numSteps);
603 for(
int i=0;i<numSteps;++i) {
604 steps[i] = pow(10,-i);
607 return checkApplyAdjointHessian(x,u,v,hv,steps,printToStream,outStream,order);
612 template <
class Real>
617 const std::vector<Real> &steps,
618 const bool printToStream,
619 std::ostream & outStream,
626 int numSteps = steps.size();
628 std::vector<Real> tmp(numVals);
629 std::vector<std::vector<Real> > ahuvCheck(numSteps, tmp);
632 Teuchos::RCP<Vector<Real> > AJdif = hv.
clone();
633 Teuchos::RCP<Vector<Real> > AJu = hv.
clone();
634 Teuchos::RCP<Vector<Real> > AHuv = hv.
clone();
635 Teuchos::RCP<Vector<Real> > AJnew = hv.
clone();
636 Teuchos::RCP<Vector<Real> > xnew = x.
clone();
639 Teuchos::oblackholestream oldFormatState;
640 oldFormatState.copyfmt(outStream);
644 this->applyAdjointJacobian(*AJu, u, x, tol);
647 this->applyAdjointHessian(*AHuv, u, v, x, tol);
648 Real normAHuv = AHuv->norm();
650 for (
int i=0; i<numSteps; i++) {
658 AJdif->scale(
weights[order-1][0]);
660 for(
int j=0; j<order; ++j) {
662 xnew->axpy(eta*
shifts[order-1][j],v);
664 if(
weights[order-1][j+1] != 0 ) {
666 this->applyAdjointJacobian(*AJnew, u, *xnew, tol);
667 AJdif->axpy(
weights[order-1][j+1],*AJnew);
671 AJdif->scale(1.0/eta);
674 ahuvCheck[i][0] = eta;
675 ahuvCheck[i][1] = normAHuv;
676 ahuvCheck[i][2] = AJdif->norm();
677 AJdif->axpy(-1.0, *AHuv);
678 ahuvCheck[i][3] = AJdif->norm();
681 std::stringstream hist;
684 << std::setw(20) <<
"Step size" 685 << std::setw(20) <<
"norm(adj(H)(u,v))" 686 << std::setw(20) <<
"norm(FD approx)" 687 << std::setw(20) <<
"norm(abs error)" 689 << std::setw(20) <<
"---------" 690 << std::setw(20) <<
"-----------------" 691 << std::setw(20) <<
"---------------" 692 << std::setw(20) <<
"---------------" 695 hist << std::scientific << std::setprecision(11) << std::right
696 << std::setw(20) << ahuvCheck[i][0]
697 << std::setw(20) << ahuvCheck[i][1]
698 << std::setw(20) << ahuvCheck[i][2]
699 << std::setw(20) << ahuvCheck[i][3]
701 outStream << hist.str();
707 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 scale(const Real alpha)=0
Compute where .
virtual int dimension() const
Return dimension of the vector space.
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 Real checkAdjointConsistencyJacobian(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &x, const bool printToStream=true, std::ostream &outStream=std::cout)
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 > > 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 Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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 .
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 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 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 Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
virtual Real norm() const =0
Returns where .
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. ...
static const double ROL_UNDERFLOW
Platform-dependent minimum double.
static const double ROL_EPSILON
Platform-dependent machine epsilon.