44 #ifndef ROL_EQUALITY_CONSTRAINT_SIMOPT_H 45 #define ROL_EQUALITY_CONSTRAINT_SIMOPT_H 98 Teuchos::RCP<Vector<Real> >
unew_;
99 Teuchos::RCP<Vector<Real> >
jv_;
125 unew_(Teuchos::null),
jv_(Teuchos::null),
202 Real cnorm = c.
norm(), alpha = 1.0, tmp = 0.0;
203 const Real ctol = std::min(
atol_,
rtol_*cnorm);
206 std::cout <<
"\n Default EqualityConstraint_SimOpt::solve\n";
208 std::cout << std::setw(6) << std::left <<
"iter";
209 std::cout << std::setw(15) << std::left <<
"rnorm";
210 std::cout << std::setw(15) << std::left <<
"alpha";
213 for (cnt = 0; cnt <
maxit_; ++cnt) {
222 while ( tmp > (1.0-
decr_*alpha)*cnorm &&
233 std::cout << std::setw(6) << std::left << cnt;
234 std::cout << std::scientific << std::setprecision(6);
235 std::cout << std::setw(15) << std::left << tmp;
236 std::cout << std::scientific << std::setprecision(6);
237 std::cout << std::setw(15) << std::left << alpha;
270 Teuchos::ParameterList & list = parlist.sublist(
"SimOpt").sublist(
"Solve");
300 Real ctol = std::sqrt(ROL_EPSILON<Real>());
303 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
304 h = std::max(1.0,u.
norm()/v.
norm())*tol;
307 Teuchos::RCP<Vector<Real> > unew = u.
clone();
312 value(jv,*unew,z,ctol);
314 Teuchos::RCP<Vector<Real> > cold = jv.
clone();
316 value(*cold,u,z,ctol);
343 Real ctol = std::sqrt(ROL_EPSILON<Real>());
346 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
347 h = std::max(1.0,u.
norm()/v.
norm())*tol;
350 Teuchos::RCP<Vector<Real> > znew = z.
clone();
355 value(jv,u,*znew,ctol);
357 Teuchos::RCP<Vector<Real> > cold = jv.
clone();
359 value(*cold,u,z,ctol);
385 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
386 "The method applyInverseJacobian_1 is used but not implemented!\n");
436 Real ctol = std::sqrt(ROL_EPSILON<Real>());
438 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
439 h = std::max(1.0,u.
norm()/v.
norm())*tol;
441 Teuchos::RCP<Vector<Real> > cold = dualv.
clone();
442 Teuchos::RCP<Vector<Real> > cnew = dualv.
clone();
444 value(*cold,u,z,ctol);
445 Teuchos::RCP<Vector<Real> > unew = u.
clone();
447 for (
int i = 0; i < u.
dimension(); i++) {
449 unew->axpy(h,*(u.
basis(i)));
451 value(*cnew,*unew,z,ctol);
452 cnew->axpy(-1.0,*cold);
454 ajv.
axpy(cnew->dot(v),*((u.
dual()).basis(i)));
507 Real ctol = std::sqrt(ROL_EPSILON<Real>());
509 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
510 h = std::max(1.0,u.
norm()/v.
norm())*tol;
512 Teuchos::RCP<Vector<Real> > cold = dualv.
clone();
513 Teuchos::RCP<Vector<Real> > cnew = dualv.
clone();
515 value(*cold,u,z,ctol);
516 Teuchos::RCP<Vector<Real> > znew = z.
clone();
518 for (
int i = 0; i < z.
dimension(); i++) {
520 znew->axpy(h,*(z.
basis(i)));
522 value(*cnew,u,*znew,ctol);
523 cnew->axpy(-1.0,*cold);
525 ajv.
axpy(cnew->dot(v),*((z.
dual()).basis(i)));
550 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
551 "The method applyInverseAdjointJacobian_1 is used but not implemented!\n");
577 Real jtol = std::sqrt(ROL_EPSILON<Real>());
580 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
581 h = std::max(1.0,u.
norm()/v.
norm())*tol;
584 Teuchos::RCP<Vector<Real> > unew = u.
clone();
590 Teuchos::RCP<Vector<Real> > jv = ahwv.
clone();
622 Real jtol = std::sqrt(ROL_EPSILON<Real>());
625 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
626 h = std::max(1.0,u.
norm()/v.
norm())*tol;
629 Teuchos::RCP<Vector<Real> > unew = u.
clone();
635 Teuchos::RCP<Vector<Real> > jv = ahwv.
clone();
667 Real jtol = std::sqrt(ROL_EPSILON<Real>());
670 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
671 h = std::max(1.0,u.
norm()/v.
norm())*tol;
674 Teuchos::RCP<Vector<Real> > znew = z.
clone();
680 Teuchos::RCP<Vector<Real> > jv = ahwv.
clone();
711 Real jtol = std::sqrt(ROL_EPSILON<Real>());
714 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
715 h = std::max(1.0,u.
norm()/v.
norm())*tol;
718 Teuchos::RCP<Vector<Real> > znew = z.
clone();
724 Teuchos::RCP<Vector<Real> > jv = ahwv.
clone();
804 Teuchos::RCP<ROL::Vector<Real> > ijv = (xs.
get_1())->clone();
809 catch (
const std::logic_error &e) {
815 Teuchos::RCP<ROL::Vector<Real> > ijv_dual = (gs.
get_1())->clone();
816 ijv_dual->
set(ijv->dual());
821 catch (
const std::logic_error &e) {
864 Teuchos::RCP<Vector<Real> > jv2 = jv.
clone();
878 Teuchos::RCP<Vector<Real> > ajv1 = (ajvs.
get_1())->clone();
881 Teuchos::RCP<Vector<Real> > ajv2 = (ajvs.
get_2())->clone();
899 Teuchos::RCP<Vector<Real> > C11 = (ahwvs.
get_1())->clone();
900 Teuchos::RCP<Vector<Real> > C21 = (ahwvs.
get_1())->clone();
906 Teuchos::RCP<Vector<Real> > C12 = (ahwvs.
get_2())->clone();
907 Teuchos::RCP<Vector<Real> > C22 = (ahwvs.
get_2())->clone();
919 const bool printToStream =
true,
920 std::ostream & outStream = std::cout) {
922 Real tol = ROL_EPSILON<Real>();
923 Teuchos::RCP<ROL::Vector<Real> > r = c.
clone();
924 Teuchos::RCP<ROL::Vector<Real> > s = u.
clone();
927 Teuchos::RCP<ROL::Vector<Real> > cs = c.
clone();
931 Real rnorm = r->norm();
932 Real cnorm = cs->norm();
933 if ( printToStream ) {
934 std::stringstream hist;
935 hist << std::scientific << std::setprecision(8);
936 hist <<
"\nTest SimOpt solve at feasible (u,z):\n";
937 hist <<
" Solver Residual = " << rnorm <<
"\n";
938 hist <<
" ||c(u,z)|| = " << cnorm <<
"\n";
939 outStream << hist.str();
961 const bool printToStream =
true,
962 std::ostream & outStream = std::cout) {
988 const bool printToStream =
true,
989 std::ostream & outStream = std::cout) {
990 Real tol = ROL_EPSILON<Real>();
991 Teuchos::RCP<Vector<Real> > Jv = dualw.
clone();
994 Real wJv = w.
dot(Jv->dual());
995 Teuchos::RCP<Vector<Real> > Jw = dualv.
clone();
998 Real vJw = v.
dot(Jw->dual());
999 Real diff = std::abs(wJv-vJw);
1000 if ( printToStream ) {
1001 std::stringstream hist;
1002 hist << std::scientific << std::setprecision(8);
1003 hist <<
"\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = " 1005 hist <<
" |<w,Jv>| = " << std::abs(wJv) <<
"\n";
1006 hist <<
" Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) <<
"\n";
1007 outStream << hist.str();
1029 const bool printToStream =
true,
1030 std::ostream & outStream = std::cout) {
1055 const bool printToStream =
true,
1056 std::ostream & outStream = std::cout) {
1057 Real tol = ROL_EPSILON<Real>();
1058 Teuchos::RCP<Vector<Real> > Jv = dualw.
clone();
1061 Real wJv = w.
dot(Jv->dual());
1062 Teuchos::RCP<Vector<Real> > Jw = dualv.
clone();
1065 Real vJw = v.
dot(Jw->dual());
1066 Real diff = std::abs(wJv-vJw);
1067 if ( printToStream ) {
1068 std::stringstream hist;
1069 hist << std::scientific << std::setprecision(8);
1070 hist <<
"\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = " 1072 hist <<
" |<w,Jv>| = " << std::abs(wJv) <<
"\n";
1073 hist <<
" Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) <<
"\n";
1074 outStream << hist.str();
1083 const bool printToStream =
true,
1084 std::ostream & outStream = std::cout) {
1085 Real tol = ROL_EPSILON<Real>();
1086 Teuchos::RCP<Vector<Real> > Jv = jv.
clone();
1089 Teuchos::RCP<Vector<Real> > iJJv = u.
clone();
1092 Teuchos::RCP<Vector<Real> > diff = v.
clone();
1094 diff->axpy(-1.0,*iJJv);
1095 Real dnorm = diff->norm();
1096 Real vnorm = v.
norm();
1097 if ( printToStream ) {
1098 std::stringstream hist;
1099 hist << std::scientific << std::setprecision(8);
1100 hist <<
"\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = " 1102 hist <<
" ||v|| = " << vnorm <<
"\n";
1103 hist <<
" Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) <<
"\n";
1104 outStream << hist.str();
1113 const bool printToStream =
true,
1114 std::ostream & outStream = std::cout) {
1115 Real tol = ROL_EPSILON<Real>();
1116 Teuchos::RCP<Vector<Real> > Jv = jv.
clone();
1119 Teuchos::RCP<Vector<Real> > iJJv = v.
clone();
1122 Teuchos::RCP<Vector<Real> > diff = v.
clone();
1124 diff->axpy(-1.0,*iJJv);
1125 Real dnorm = diff->norm();
1126 Real vnorm = v.
norm();
1127 if ( printToStream ) {
1128 std::stringstream hist;
1129 hist << std::scientific << std::setprecision(8);
1130 hist <<
"\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = " 1132 hist <<
" ||v|| = " << vnorm <<
"\n";
1133 hist <<
" Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) <<
"\n";
1134 outStream << hist.str();
1145 const bool printToStream =
true,
1146 std::ostream & outStream = std::cout,
1148 const int order = 1) {
1149 std::vector<Real> steps(numSteps);
1150 for(
int i=0;i<numSteps;++i) {
1151 steps[i] = pow(10,-i);
1164 const std::vector<Real> &steps,
1165 const bool printToStream =
true,
1166 std::ostream & outStream = std::cout,
1167 const int order = 1) {
1169 TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1170 "Error: finite difference order must be 1,2,3, or 4" );
1177 Real tol = std::sqrt(ROL_EPSILON<Real>());
1179 int numSteps = steps.size();
1181 std::vector<Real> tmp(numVals);
1182 std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
1185 Teuchos::oblackholestream oldFormatState;
1186 oldFormatState.copyfmt(outStream);
1189 Teuchos::RCP<Vector<Real> > c = jv.
clone();
1191 this->
value(*c, u, z, tol);
1194 Teuchos::RCP<Vector<Real> > Jv = jv.
clone();
1196 Real normJv = Jv->norm();
1199 Teuchos::RCP<Vector<Real> > cdif = jv.
clone();
1200 Teuchos::RCP<Vector<Real> > cnew = jv.
clone();
1201 Teuchos::RCP<Vector<Real> > unew = u.
clone();
1203 for (
int i=0; i<numSteps; i++) {
1205 Real eta = steps[i];
1210 cdif->scale(
weights[order-1][0]);
1212 for(
int j=0; j<order; ++j) {
1214 unew->axpy(eta*
shifts[order-1][j], v);
1216 if(
weights[order-1][j+1] != 0 ) {
1218 this->
value(*cnew,*unew,z,tol);
1219 cdif->axpy(
weights[order-1][j+1],*cnew);
1224 cdif->scale(one/eta);
1227 jvCheck[i][0] = eta;
1228 jvCheck[i][1] = normJv;
1229 jvCheck[i][2] = cdif->norm();
1230 cdif->axpy(-one, *Jv);
1231 jvCheck[i][3] = cdif->norm();
1233 if (printToStream) {
1234 std::stringstream hist;
1237 << std::setw(20) <<
"Step size" 1238 << std::setw(20) <<
"norm(Jac*vec)" 1239 << std::setw(20) <<
"norm(FD approx)" 1240 << std::setw(20) <<
"norm(abs error)" 1242 << std::setw(20) <<
"---------" 1243 << std::setw(20) <<
"-------------" 1244 << std::setw(20) <<
"---------------" 1245 << std::setw(20) <<
"---------------" 1248 hist << std::scientific << std::setprecision(11) << std::right
1249 << std::setw(20) << jvCheck[i][0]
1250 << std::setw(20) << jvCheck[i][1]
1251 << std::setw(20) << jvCheck[i][2]
1252 << std::setw(20) << jvCheck[i][3]
1254 outStream << hist.str();
1260 outStream.copyfmt(oldFormatState);
1270 const bool printToStream =
true,
1271 std::ostream & outStream = std::cout,
1273 const int order = 1) {
1274 std::vector<Real> steps(numSteps);
1275 for(
int i=0;i<numSteps;++i) {
1276 steps[i] = pow(10,-i);
1289 const std::vector<Real> &steps,
1290 const bool printToStream =
true,
1291 std::ostream & outStream = std::cout,
1292 const int order = 1) {
1294 TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1295 "Error: finite difference order must be 1,2,3, or 4" );
1302 Real tol = std::sqrt(ROL_EPSILON<Real>());
1304 int numSteps = steps.size();
1306 std::vector<Real> tmp(numVals);
1307 std::vector<std::vector<Real> > jvCheck(numSteps, tmp);
1310 Teuchos::oblackholestream oldFormatState;
1311 oldFormatState.copyfmt(outStream);
1314 Teuchos::RCP<Vector<Real> > c = jv.
clone();
1316 this->
value(*c, u, z, tol);
1319 Teuchos::RCP<Vector<Real> > Jv = jv.
clone();
1321 Real normJv = Jv->norm();
1324 Teuchos::RCP<Vector<Real> > cdif = jv.
clone();
1325 Teuchos::RCP<Vector<Real> > cnew = jv.
clone();
1326 Teuchos::RCP<Vector<Real> > znew = z.
clone();
1328 for (
int i=0; i<numSteps; i++) {
1330 Real eta = steps[i];
1335 cdif->scale(
weights[order-1][0]);
1337 for(
int j=0; j<order; ++j) {
1339 znew->axpy(eta*
shifts[order-1][j], v);
1341 if(
weights[order-1][j+1] != 0 ) {
1343 this->
value(*cnew,u,*znew,tol);
1344 cdif->axpy(
weights[order-1][j+1],*cnew);
1349 cdif->scale(one/eta);
1352 jvCheck[i][0] = eta;
1353 jvCheck[i][1] = normJv;
1354 jvCheck[i][2] = cdif->norm();
1355 cdif->axpy(-one, *Jv);
1356 jvCheck[i][3] = cdif->norm();
1358 if (printToStream) {
1359 std::stringstream hist;
1362 << std::setw(20) <<
"Step size" 1363 << std::setw(20) <<
"norm(Jac*vec)" 1364 << std::setw(20) <<
"norm(FD approx)" 1365 << std::setw(20) <<
"norm(abs error)" 1367 << std::setw(20) <<
"---------" 1368 << std::setw(20) <<
"-------------" 1369 << std::setw(20) <<
"---------------" 1370 << std::setw(20) <<
"---------------" 1373 hist << std::scientific << std::setprecision(11) << std::right
1374 << std::setw(20) << jvCheck[i][0]
1375 << std::setw(20) << jvCheck[i][1]
1376 << std::setw(20) << jvCheck[i][2]
1377 << std::setw(20) << jvCheck[i][3]
1379 outStream << hist.str();
1385 outStream.copyfmt(oldFormatState);
1397 const bool printToStream =
true,
1398 std::ostream & outStream = std::cout,
1400 const int order = 1 ) {
1401 std::vector<Real> steps(numSteps);
1402 for(
int i=0;i<numSteps;++i) {
1403 steps[i] = pow(10,-i);
1415 const std::vector<Real> &steps,
1416 const bool printToStream =
true,
1417 std::ostream & outStream = std::cout,
1418 const int order = 1 ) {
1424 Real tol = std::sqrt(ROL_EPSILON<Real>());
1426 int numSteps = steps.size();
1428 std::vector<Real> tmp(numVals);
1429 std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1432 Teuchos::RCP<Vector<Real> > AJdif = hv.
clone();
1433 Teuchos::RCP<Vector<Real> > AJp = hv.
clone();
1434 Teuchos::RCP<Vector<Real> > AHpv = hv.
clone();
1435 Teuchos::RCP<Vector<Real> > AJnew = hv.
clone();
1436 Teuchos::RCP<Vector<Real> > unew = u.
clone();
1439 Teuchos::oblackholestream oldFormatState;
1440 oldFormatState.copyfmt(outStream);
1448 Real normAHpv = AHpv->norm();
1450 for (
int i=0; i<numSteps; i++) {
1452 Real eta = steps[i];
1458 AJdif->scale(
weights[order-1][0]);
1460 for(
int j=0; j<order; ++j) {
1462 unew->axpy(eta*
shifts[order-1][j],v);
1464 if(
weights[order-1][j+1] != 0 ) {
1467 AJdif->axpy(
weights[order-1][j+1],*AJnew);
1471 AJdif->scale(one/eta);
1474 ahpvCheck[i][0] = eta;
1475 ahpvCheck[i][1] = normAHpv;
1476 ahpvCheck[i][2] = AJdif->norm();
1477 AJdif->axpy(-one, *AHpv);
1478 ahpvCheck[i][3] = AJdif->norm();
1480 if (printToStream) {
1481 std::stringstream hist;
1484 << std::setw(20) <<
"Step size" 1485 << std::setw(20) <<
"norm(adj(H)(u,v))" 1486 << std::setw(20) <<
"norm(FD approx)" 1487 << std::setw(20) <<
"norm(abs error)" 1489 << std::setw(20) <<
"---------" 1490 << std::setw(20) <<
"-----------------" 1491 << std::setw(20) <<
"---------------" 1492 << std::setw(20) <<
"---------------" 1495 hist << std::scientific << std::setprecision(11) << std::right
1496 << std::setw(20) << ahpvCheck[i][0]
1497 << std::setw(20) << ahpvCheck[i][1]
1498 << std::setw(20) << ahpvCheck[i][2]
1499 << std::setw(20) << ahpvCheck[i][3]
1501 outStream << hist.str();
1507 outStream.copyfmt(oldFormatState);
1517 const bool printToStream =
true,
1518 std::ostream & outStream = std::cout,
1520 const int order = 1 ) {
1521 std::vector<Real> steps(numSteps);
1522 for(
int i=0;i<numSteps;++i) {
1523 steps[i] = pow(10,-i);
1535 const std::vector<Real> &steps,
1536 const bool printToStream =
true,
1537 std::ostream & outStream = std::cout,
1538 const int order = 1 ) {
1544 Real tol = std::sqrt(ROL_EPSILON<Real>());
1546 int numSteps = steps.size();
1548 std::vector<Real> tmp(numVals);
1549 std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1552 Teuchos::RCP<Vector<Real> > AJdif = hv.
clone();
1553 Teuchos::RCP<Vector<Real> > AJp = hv.
clone();
1554 Teuchos::RCP<Vector<Real> > AHpv = hv.
clone();
1555 Teuchos::RCP<Vector<Real> > AJnew = hv.
clone();
1556 Teuchos::RCP<Vector<Real> > znew = z.
clone();
1559 Teuchos::oblackholestream oldFormatState;
1560 oldFormatState.copyfmt(outStream);
1568 Real normAHpv = AHpv->norm();
1570 for (
int i=0; i<numSteps; i++) {
1572 Real eta = steps[i];
1578 AJdif->scale(
weights[order-1][0]);
1580 for(
int j=0; j<order; ++j) {
1582 znew->axpy(eta*
shifts[order-1][j],v);
1584 if(
weights[order-1][j+1] != 0 ) {
1587 AJdif->axpy(
weights[order-1][j+1],*AJnew);
1591 AJdif->scale(one/eta);
1594 ahpvCheck[i][0] = eta;
1595 ahpvCheck[i][1] = normAHpv;
1596 ahpvCheck[i][2] = AJdif->norm();
1597 AJdif->axpy(-one, *AHpv);
1598 ahpvCheck[i][3] = AJdif->norm();
1600 if (printToStream) {
1601 std::stringstream hist;
1604 << std::setw(20) <<
"Step size" 1605 << std::setw(20) <<
"norm(adj(H)(u,v))" 1606 << std::setw(20) <<
"norm(FD approx)" 1607 << std::setw(20) <<
"norm(abs error)" 1609 << std::setw(20) <<
"---------" 1610 << std::setw(20) <<
"-----------------" 1611 << std::setw(20) <<
"---------------" 1612 << std::setw(20) <<
"---------------" 1615 hist << std::scientific << std::setprecision(11) << std::right
1616 << std::setw(20) << ahpvCheck[i][0]
1617 << std::setw(20) << ahpvCheck[i][1]
1618 << std::setw(20) << ahpvCheck[i][2]
1619 << std::setw(20) << ahpvCheck[i][3]
1621 outStream << hist.str();
1627 outStream.copyfmt(oldFormatState);
1637 const bool printToStream =
true,
1638 std::ostream & outStream = std::cout,
1640 const int order = 1 ) {
1641 std::vector<Real> steps(numSteps);
1642 for(
int i=0;i<numSteps;++i) {
1643 steps[i] = pow(10,-i);
1655 const std::vector<Real> &steps,
1656 const bool printToStream =
true,
1657 std::ostream & outStream = std::cout,
1658 const int order = 1 ) {
1664 Real tol = std::sqrt(ROL_EPSILON<Real>());
1666 int numSteps = steps.size();
1668 std::vector<Real> tmp(numVals);
1669 std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1672 Teuchos::RCP<Vector<Real> > AJdif = hv.
clone();
1673 Teuchos::RCP<Vector<Real> > AJp = hv.
clone();
1674 Teuchos::RCP<Vector<Real> > AHpv = hv.
clone();
1675 Teuchos::RCP<Vector<Real> > AJnew = hv.
clone();
1676 Teuchos::RCP<Vector<Real> > unew = u.
clone();
1679 Teuchos::oblackholestream oldFormatState;
1680 oldFormatState.copyfmt(outStream);
1688 Real normAHpv = AHpv->norm();
1690 for (
int i=0; i<numSteps; i++) {
1692 Real eta = steps[i];
1698 AJdif->scale(
weights[order-1][0]);
1700 for(
int j=0; j<order; ++j) {
1702 unew->axpy(eta*
shifts[order-1][j],v);
1704 if(
weights[order-1][j+1] != 0 ) {
1707 AJdif->axpy(
weights[order-1][j+1],*AJnew);
1711 AJdif->scale(one/eta);
1714 ahpvCheck[i][0] = eta;
1715 ahpvCheck[i][1] = normAHpv;
1716 ahpvCheck[i][2] = AJdif->norm();
1717 AJdif->axpy(-one, *AHpv);
1718 ahpvCheck[i][3] = AJdif->norm();
1720 if (printToStream) {
1721 std::stringstream hist;
1724 << std::setw(20) <<
"Step size" 1725 << std::setw(20) <<
"norm(adj(H)(u,v))" 1726 << std::setw(20) <<
"norm(FD approx)" 1727 << std::setw(20) <<
"norm(abs error)" 1729 << std::setw(20) <<
"---------" 1730 << std::setw(20) <<
"-----------------" 1731 << std::setw(20) <<
"---------------" 1732 << std::setw(20) <<
"---------------" 1735 hist << std::scientific << std::setprecision(11) << std::right
1736 << std::setw(20) << ahpvCheck[i][0]
1737 << std::setw(20) << ahpvCheck[i][1]
1738 << std::setw(20) << ahpvCheck[i][2]
1739 << std::setw(20) << ahpvCheck[i][3]
1741 outStream << hist.str();
1747 outStream.copyfmt(oldFormatState);
1757 const bool printToStream =
true,
1758 std::ostream & outStream = std::cout,
1760 const int order = 1 ) {
1761 std::vector<Real> steps(numSteps);
1762 for(
int i=0;i<numSteps;++i) {
1763 steps[i] = pow(10,-i);
1775 const std::vector<Real> &steps,
1776 const bool printToStream =
true,
1777 std::ostream & outStream = std::cout,
1778 const int order = 1 ) {
1784 Real tol = std::sqrt(ROL_EPSILON<Real>());
1786 int numSteps = steps.size();
1788 std::vector<Real> tmp(numVals);
1789 std::vector<std::vector<Real> > ahpvCheck(numSteps, tmp);
1792 Teuchos::RCP<Vector<Real> > AJdif = hv.
clone();
1793 Teuchos::RCP<Vector<Real> > AJp = hv.
clone();
1794 Teuchos::RCP<Vector<Real> > AHpv = hv.
clone();
1795 Teuchos::RCP<Vector<Real> > AJnew = hv.
clone();
1796 Teuchos::RCP<Vector<Real> > znew = z.
clone();
1799 Teuchos::oblackholestream oldFormatState;
1800 oldFormatState.copyfmt(outStream);
1808 Real normAHpv = AHpv->norm();
1810 for (
int i=0; i<numSteps; i++) {
1812 Real eta = steps[i];
1818 AJdif->scale(
weights[order-1][0]);
1820 for(
int j=0; j<order; ++j) {
1822 znew->axpy(eta*
shifts[order-1][j],v);
1824 if(
weights[order-1][j+1] != 0 ) {
1827 AJdif->axpy(
weights[order-1][j+1],*AJnew);
1831 AJdif->scale(one/eta);
1834 ahpvCheck[i][0] = eta;
1835 ahpvCheck[i][1] = normAHpv;
1836 ahpvCheck[i][2] = AJdif->norm();
1837 AJdif->axpy(-one, *AHpv);
1838 ahpvCheck[i][3] = AJdif->norm();
1840 if (printToStream) {
1841 std::stringstream hist;
1844 << std::setw(20) <<
"Step size" 1845 << std::setw(20) <<
"norm(adj(H)(u,v))" 1846 << std::setw(20) <<
"norm(FD approx)" 1847 << std::setw(20) <<
"norm(abs error)" 1849 << std::setw(20) <<
"---------" 1850 << std::setw(20) <<
"-----------------" 1851 << std::setw(20) <<
"---------------" 1852 << std::setw(20) <<
"---------------" 1855 hist << std::scientific << std::setprecision(11) << std::right
1856 << std::setw(20) << ahpvCheck[i][0]
1857 << std::setw(20) << ahpvCheck[i][1]
1858 << std::setw(20) << ahpvCheck[i][2]
1859 << std::setw(20) << ahpvCheck[i][3]
1861 outStream << hist.str();
1867 outStream.copyfmt(oldFormatState);
std::vector< std::vector< Real > > checkApplyAdjointHessian_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
virtual Real checkInverseJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
virtual void applyAdjointJacobian_2(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualv, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the secondary interfa...
virtual void applyAdjointJacobian_1(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
virtual void applyAdjointJacobian_2(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
virtual void scale(const Real alpha)=0
Compute where .
virtual void applyJacobian_2(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
virtual void setSolveParameters(Teuchos::ParameterList &parlist)
Set solve parameters.
std::vector< std::vector< Real > > checkApplyAdjointHessian_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
virtual void plus(const Vector &x)=0
Compute , where .
const double weights[4][5]
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Defines the linear algebra or vector space interface for simulation-based optimization.
virtual void value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol)=0
Evaluate the constraint operator at .
std::vector< std::vector< Real > > checkApplyJacobian_2(const Vector< Real > &u, const Vector< Real > &z, 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)
virtual void applyInverseAdjointJacobian_1(Vector< Real > &iajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector ...
Teuchos::RCP< const Vector< Real > > get_2() const
virtual void applyAdjointHessian_11(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
virtual void applyAdjointHessian_22(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
Contains definitions of custom data types in ROL.
std::vector< std::vector< Real > > checkApplyAdjointHessian_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
const Real DEFAULT_factor_
void set_1(const Vector< Real > &vec)
virtual void applyAdjointHessian_21(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void applyAdjointHessian(Vector< Real > &ahwv, const Vector< Real > &w, 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 void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
Teuchos::RCP< Vector< Real > > jv_
const bool DEFAULT_print_
Defines the equality constraint operator interface for simulation-based optimization.
virtual Real dot(const Vector &x) const =0
Compute where .
virtual Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
Teuchos::RCP< Vector< Real > > unew_
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
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 . In general, this preconditioner satisfies the fo...
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable i...
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 const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Defines the equality constraint operator interface.
virtual Real checkSolve(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, const ROL::Vector< Real > &c, const bool printToStream=true, std::ostream &outStream=std::cout)
virtual void applyInverseJacobian_1(Vector< Real > &ijv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the inverse partial constraint Jacobian at , , to the vector .
virtual void applyAdjointHessian_12(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
EqualityConstraint_SimOpt()
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
virtual int dimension() const
Return dimension of the vector space.
virtual void applyJacobian_1(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
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 .
Teuchos::RCP< const Vector< Real > > get_1() const
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualw, const Vector< Real > &dualv, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the secondary interface, for use with dual spaces where the user does not define the dual() operation.
std::vector< std::vector< Real > > checkApplyJacobian_2(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
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...
virtual void update_1(const Vector< Real > &u, bool flag=true, int iter=-1)
Update constraint functions with respect to Sim variable. x is the optimization variable...
std::vector< std::vector< Real > > checkApplyAdjointHessian_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
std::vector< std::vector< Real > > checkApplyAdjointHessian_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
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 update_2(const Vector< Real > &z, bool flag=true, int iter=-1)
Update constraint functions with respect to Opt variable. x is the optimization variable, flag = true if optimization variable is changed, iter is the outer algorithm iterations count.
virtual Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualw, const Vector< Real > &dualv, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the secondary interface, for use with dual spaces where the user does not define the dual() operation.
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable i...
virtual Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
virtual void set(const Vector &x)
Set where .
virtual void applyAdjointJacobian_1(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualv, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the secondary int...
virtual Real checkInverseAdjointJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
virtual Real norm() const =0
Returns where .
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 operator, and is a zero operator.
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
std::vector< std::vector< Real > > checkApplyJacobian_1(const Vector< Real > &u, const Vector< Real > &z, 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)
std::vector< std::vector< Real > > checkApplyAdjointHessian_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
void set_2(const Vector< Real > &vec)
std::vector< std::vector< Real > > checkApplyAdjointHessian_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
std::vector< std::vector< Real > > checkApplyJacobian_1(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &jv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
std::vector< std::vector< Real > > checkApplyAdjointHessian_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &hv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)