44#ifndef ROL_CONSTRAINT_SIMOPT_H
45#define ROL_CONSTRAINT_SIMOPT_H
55class Constraint_SimOpt;
230 Ptr<std::ostream> stream = makeStreamPtr(std::cout,
print_);
233 Real cnorm = c.
norm();
234 const Real ctol = std::min(
atol_,
rtol_*cnorm);
242 Real alpha(1), tmp(0);
244 *stream << std::endl;
245 *stream <<
" Default Constraint_SimOpt::solve" << std::endl;
247 *stream << std::setw(6) << std::left <<
"iter";
248 *stream << std::setw(15) << std::left <<
"rnorm";
249 *stream << std::setw(15) << std::left <<
"alpha";
250 *stream << std::endl;
251 for (cnt = 0; cnt <
maxit_; ++cnt) {
260 while ( tmp > (one-
decr_*alpha)*cnorm &&
270 *stream << std::setw(6) << std::left << cnt;
271 *stream << std::scientific << std::setprecision(6);
272 *stream << std::setw(15) << std::left << tmp;
273 *stream << std::scientific << std::setprecision(6);
274 *stream << std::setw(15) << std::left << alpha;
275 *stream << std::endl;
280 if (cnorm <= ctol)
break;
285 Ptr<Constraint<Real>> con = makePtr<SimConstraint<Real>>(makePtrFromRef(*
this),makePtrFromRef(z),
true);
286 Ptr<Objective<Real>> obj = makePtr<NonlinearLeastSquaresObjective<Real>>(con,u,c,
true);
287 ParameterList parlist;
288 parlist.sublist(
"General").set(
"Output Level",1);
289 parlist.sublist(
"General").sublist(
"Krylov").set(
"Iteration Limit",100);
290 parlist.sublist(
"Step").sublist(
"Trust Region").set(
"Subproblem Solver",
"Truncated CG");
291 parlist.sublist(
"Status Test").set(
"Gradient Tolerance",ctol);
292 parlist.sublist(
"Status Test").set(
"Step Tolerance",
stol_);
293 parlist.sublist(
"Status Test").set(
"Iteration Limit",
maxit_);
294 Ptr<TypeU::Algorithm<Real>> algo = makePtr<TypeU::TrustRegionAlgorithm<Real>>(parlist);
295 algo->run(u,*obj,*stream);
299 Ptr<Constraint<Real>> con = makePtr<SimConstraint<Real>>(makePtrFromRef(*
this),makePtrFromRef(z),
true);
300 Ptr<Objective<Real>> obj = makePtr<Objective_FSsolver<Real>>();
302 ParameterList parlist;
303 parlist.sublist(
"General").set(
"Output Level",1);
304 parlist.sublist(
"Status Test").set(
"Constraint Tolerance",ctol);
305 parlist.sublist(
"Status Test").set(
"Step Tolerance",
stol_);
306 parlist.sublist(
"Status Test").set(
"Iteration Limit",
maxit_);
307 Ptr<TypeE::Algorithm<Real>> algo = makePtr<TypeE::AugmentedLagrangianAlgorithm<Real>>(parlist);
308 algo->run(u,*obj,*con,*l,*stream);
312 ROL_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
313 ">>> ERROR (ROL:Constraint_SimOpt:solve): Invalid solver type!");
338 ParameterList & list = parlist.sublist(
"SimOpt").sublist(
"Solve");
370 Real ctol = std::sqrt(ROL_EPSILON<Real>());
373 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
374 h = std::max(1.0,u.
norm()/v.
norm())*tol;
377 Ptr<Vector<Real>> unew = u.
clone();
382 value(jv,*unew,z,ctol);
384 Ptr<Vector<Real>> cold = jv.
clone();
386 value(*cold,u,z,ctol);
413 Real ctol = std::sqrt(ROL_EPSILON<Real>());
416 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
417 h = std::max(1.0,u.
norm()/v.
norm())*tol;
420 Ptr<Vector<Real>> znew = z.
clone();
425 value(jv,u,*znew,ctol);
427 Ptr<Vector<Real>> cold = jv.
clone();
429 value(*cold,u,z,ctol);
455 ROL_TEST_FOR_EXCEPTION(
true, std::logic_error,
456 "The method applyInverseJacobian_1 is used but not implemented!\n");
506 Real ctol = std::sqrt(ROL_EPSILON<Real>());
508 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
509 h = std::max(1.0,u.
norm()/v.
norm())*tol;
511 Ptr<Vector<Real>> cold = dualv.
clone();
512 Ptr<Vector<Real>> cnew = dualv.
clone();
514 value(*cold,u,z,ctol);
515 Ptr<Vector<Real>> unew = u.
clone();
517 for (
int i = 0; i < u.
dimension(); i++) {
519 unew->axpy(h,*(u.
basis(i)));
521 value(*cnew,*unew,z,ctol);
522 cnew->axpy(-1.0,*cold);
524 ajv.
axpy(cnew->dot(v),*((u.
dual()).basis(i)));
577 Real ctol = std::sqrt(ROL_EPSILON<Real>());
579 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
580 h = std::max(1.0,u.
norm()/v.
norm())*tol;
582 Ptr<Vector<Real>> cold = dualv.
clone();
583 Ptr<Vector<Real>> cnew = dualv.
clone();
585 value(*cold,u,z,ctol);
586 Ptr<Vector<Real>> znew = z.
clone();
588 for (
int i = 0; i < z.
dimension(); i++) {
590 znew->axpy(h,*(z.
basis(i)));
592 value(*cnew,u,*znew,ctol);
593 cnew->axpy(-1.0,*cold);
595 ajv.
axpy(cnew->dot(v),*((z.
dual()).basis(i)));
620 ROL_TEST_FOR_EXCEPTION(
true, std::logic_error,
621 "The method applyInverseAdjointJacobian_1 is used but not implemented!\n");
647 Real jtol = std::sqrt(ROL_EPSILON<Real>());
650 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
651 h = std::max(1.0,u.
norm()/v.
norm())*tol;
654 Ptr<Vector<Real>> unew = u.
clone();
660 Ptr<Vector<Real>> jv = ahwv.
clone();
692 Real jtol = std::sqrt(ROL_EPSILON<Real>());
695 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
696 h = std::max(1.0,u.
norm()/v.
norm())*tol;
699 Ptr<Vector<Real>> unew = u.
clone();
705 Ptr<Vector<Real>> jv = ahwv.
clone();
737 Real jtol = std::sqrt(ROL_EPSILON<Real>());
740 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
741 h = std::max(1.0,u.
norm()/v.
norm())*tol;
744 Ptr<Vector<Real>> znew = z.
clone();
750 Ptr<Vector<Real>> jv = ahwv.
clone();
781 Real jtol = std::sqrt(ROL_EPSILON<Real>());
784 if (v.
norm() > std::sqrt(ROL_EPSILON<Real>())) {
785 h = std::max(1.0,u.
norm()/v.
norm())*tol;
788 Ptr<Vector<Real>> znew = z.
clone();
794 Ptr<Vector<Real>> jv = ahwv.
clone();
874 Ptr<Vector<Real>> ijv = (xs.
get_1())->clone();
879 catch (
const std::logic_error &e) {
885 Ptr<Vector<Real>> ijv_dual = (gs.
get_1())->clone();
891 catch (
const std::logic_error &e) {
932 Ptr<Vector<Real>> jv2 = jv.
clone();
947 Ptr<Vector<Real>> ajv1 = (ajvs.
get_1())->clone();
950 Ptr<Vector<Real>> ajv2 = (ajvs.
get_2())->clone();
968 Ptr<Vector<Real>> C11 = (ahwvs.
get_1())->clone();
969 Ptr<Vector<Real>> C21 = (ahwvs.
get_1())->clone();
975 Ptr<Vector<Real>> C12 = (ahwvs.
get_2())->clone();
976 Ptr<Vector<Real>> C22 = (ahwvs.
get_2())->clone();
988 const bool printToStream =
true,
989 std::ostream & outStream = std::cout) {
991 Real tol = ROL_EPSILON<Real>();
992 Ptr<Vector<Real>> r = c.
clone();
993 Ptr<Vector<Real>> s = u.
clone();
996 Ptr<Vector<Real>> cs = c.
clone();
1000 Real rnorm = r->norm();
1001 Real cnorm = cs->norm();
1002 if ( printToStream ) {
1003 std::stringstream hist;
1004 hist << std::scientific << std::setprecision(8);
1005 hist <<
"\nTest SimOpt solve at feasible (u,z):\n";
1006 hist <<
" Solver Residual = " << rnorm <<
"\n";
1007 hist <<
" ||c(u,z)|| = " << cnorm <<
"\n";
1008 outStream << hist.str();
1030 const bool printToStream =
true,
1031 std::ostream & outStream = std::cout) {
1057 const bool printToStream =
true,
1058 std::ostream & outStream = std::cout) {
1059 Real tol = ROL_EPSILON<Real>();
1060 Ptr<Vector<Real>> Jv = dualw.
clone();
1064 Real wJv = w.
apply(*Jv);
1065 Ptr<Vector<Real>> Jw = dualv.
clone();
1069 Real vJw = v.
apply(*Jw);
1070 Real diff = std::abs(wJv-vJw);
1071 if ( printToStream ) {
1072 std::stringstream hist;
1073 hist << std::scientific << std::setprecision(8);
1074 hist <<
"\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1076 hist <<
" |<w,Jv>| = " << std::abs(wJv) <<
"\n";
1077 hist <<
" Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) <<
"\n";
1078 outStream << hist.str();
1100 const bool printToStream =
true,
1101 std::ostream & outStream = std::cout) {
1126 const bool printToStream =
true,
1127 std::ostream & outStream = std::cout) {
1128 Real tol = ROL_EPSILON<Real>();
1129 Ptr<Vector<Real>> Jv = dualw.
clone();
1133 Real wJv = w.
apply(*Jv);
1134 Ptr<Vector<Real>> Jw = dualv.
clone();
1138 Real vJw = v.
apply(*Jw);
1139 Real diff = std::abs(wJv-vJw);
1140 if ( printToStream ) {
1141 std::stringstream hist;
1142 hist << std::scientific << std::setprecision(8);
1143 hist <<
"\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = "
1145 hist <<
" |<w,Jv>| = " << std::abs(wJv) <<
"\n";
1146 hist <<
" Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) <<
"\n";
1147 outStream << hist.str();
1156 const bool printToStream =
true,
1157 std::ostream & outStream = std::cout) {
1158 Real tol = ROL_EPSILON<Real>();
1159 Ptr<Vector<Real>> Jv = jv.
clone();
1162 Ptr<Vector<Real>> iJJv = u.
clone();
1165 Ptr<Vector<Real>> diff = v.
clone();
1167 diff->axpy(-1.0,*iJJv);
1168 Real dnorm = diff->norm();
1169 Real vnorm = v.
norm();
1170 if ( printToStream ) {
1171 std::stringstream hist;
1172 hist << std::scientific << std::setprecision(8);
1173 hist <<
"\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = "
1175 hist <<
" ||v|| = " << vnorm <<
"\n";
1176 hist <<
" Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) <<
"\n";
1177 outStream << hist.str();
1186 const bool printToStream =
true,
1187 std::ostream & outStream = std::cout) {
1188 Real tol = ROL_EPSILON<Real>();
1189 Ptr<Vector<Real>> Jv = jv.
clone();
1192 Ptr<Vector<Real>> iJJv = v.
clone();
1195 Ptr<Vector<Real>> diff = v.
clone();
1197 diff->axpy(-1.0,*iJJv);
1198 Real dnorm = diff->norm();
1199 Real vnorm = v.
norm();
1200 if ( printToStream ) {
1201 std::stringstream hist;
1202 hist << std::scientific << std::setprecision(8);
1203 hist <<
"\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = "
1205 hist <<
" ||v|| = " << vnorm <<
"\n";
1206 hist <<
" Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) <<
"\n";
1207 outStream << hist.str();
1218 const bool printToStream =
true,
1219 std::ostream & outStream = std::cout,
1221 const int order = 1) {
1222 std::vector<Real> steps(numSteps);
1223 for(
int i=0;i<numSteps;++i) {
1224 steps[i] = pow(10,-i);
1237 const std::vector<Real> &steps,
1238 const bool printToStream =
true,
1239 std::ostream & outStream = std::cout,
1240 const int order = 1) {
1242 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1243 "Error: finite difference order must be 1,2,3, or 4" );
1250 Real tol = std::sqrt(ROL_EPSILON<Real>());
1252 int numSteps = steps.size();
1254 std::vector<Real> tmp(numVals);
1255 std::vector<std::vector<Real>> jvCheck(numSteps, tmp);
1258 nullstream oldFormatState;
1259 oldFormatState.copyfmt(outStream);
1262 Ptr<Vector<Real>> c = jv.
clone();
1264 this->
value(*c, u, z, tol);
1267 Ptr<Vector<Real>> Jv = jv.
clone();
1269 Real normJv = Jv->norm();
1272 Ptr<Vector<Real>> cdif = jv.
clone();
1273 Ptr<Vector<Real>> cnew = jv.
clone();
1274 Ptr<Vector<Real>> unew = u.
clone();
1276 for (
int i=0; i<numSteps; i++) {
1278 Real eta = steps[i];
1283 cdif->scale(weights[order-1][0]);
1285 for(
int j=0; j<order; ++j) {
1287 unew->axpy(eta*shifts[order-1][j], v);
1289 if( weights[order-1][j+1] != 0 ) {
1291 this->
value(*cnew,*unew,z,tol);
1292 cdif->axpy(weights[order-1][j+1],*cnew);
1297 cdif->scale(one/eta);
1300 jvCheck[i][0] = eta;
1301 jvCheck[i][1] = normJv;
1302 jvCheck[i][2] = cdif->norm();
1303 cdif->axpy(-one, *Jv);
1304 jvCheck[i][3] = cdif->norm();
1306 if (printToStream) {
1307 std::stringstream hist;
1310 << std::setw(20) <<
"Step size"
1311 << std::setw(20) <<
"norm(Jac*vec)"
1312 << std::setw(20) <<
"norm(FD approx)"
1313 << std::setw(20) <<
"norm(abs error)"
1315 << std::setw(20) <<
"---------"
1316 << std::setw(20) <<
"-------------"
1317 << std::setw(20) <<
"---------------"
1318 << std::setw(20) <<
"---------------"
1321 hist << std::scientific << std::setprecision(11) << std::right
1322 << std::setw(20) << jvCheck[i][0]
1323 << std::setw(20) << jvCheck[i][1]
1324 << std::setw(20) << jvCheck[i][2]
1325 << std::setw(20) << jvCheck[i][3]
1327 outStream << hist.str();
1333 outStream.copyfmt(oldFormatState);
1343 const bool printToStream =
true,
1344 std::ostream & outStream = std::cout,
1346 const int order = 1) {
1347 std::vector<Real> steps(numSteps);
1348 for(
int i=0;i<numSteps;++i) {
1349 steps[i] = pow(10,-i);
1362 const std::vector<Real> &steps,
1363 const bool printToStream =
true,
1364 std::ostream & outStream = std::cout,
1365 const int order = 1) {
1367 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
1368 "Error: finite difference order must be 1,2,3, or 4" );
1375 Real tol = std::sqrt(ROL_EPSILON<Real>());
1377 int numSteps = steps.size();
1379 std::vector<Real> tmp(numVals);
1380 std::vector<std::vector<Real>> jvCheck(numSteps, tmp);
1383 nullstream oldFormatState;
1384 oldFormatState.copyfmt(outStream);
1387 Ptr<Vector<Real>> c = jv.
clone();
1389 this->
value(*c, u, z, tol);
1392 Ptr<Vector<Real>> Jv = jv.
clone();
1394 Real normJv = Jv->norm();
1397 Ptr<Vector<Real>> cdif = jv.
clone();
1398 Ptr<Vector<Real>> cnew = jv.
clone();
1399 Ptr<Vector<Real>> znew = z.
clone();
1401 for (
int i=0; i<numSteps; i++) {
1403 Real eta = steps[i];
1408 cdif->scale(weights[order-1][0]);
1410 for(
int j=0; j<order; ++j) {
1412 znew->axpy(eta*shifts[order-1][j], v);
1414 if( weights[order-1][j+1] != 0 ) {
1416 this->
value(*cnew,u,*znew,tol);
1417 cdif->axpy(weights[order-1][j+1],*cnew);
1422 cdif->scale(one/eta);
1425 jvCheck[i][0] = eta;
1426 jvCheck[i][1] = normJv;
1427 jvCheck[i][2] = cdif->norm();
1428 cdif->axpy(-one, *Jv);
1429 jvCheck[i][3] = cdif->norm();
1431 if (printToStream) {
1432 std::stringstream hist;
1435 << std::setw(20) <<
"Step size"
1436 << std::setw(20) <<
"norm(Jac*vec)"
1437 << std::setw(20) <<
"norm(FD approx)"
1438 << std::setw(20) <<
"norm(abs error)"
1440 << std::setw(20) <<
"---------"
1441 << std::setw(20) <<
"-------------"
1442 << std::setw(20) <<
"---------------"
1443 << std::setw(20) <<
"---------------"
1446 hist << std::scientific << std::setprecision(11) << std::right
1447 << std::setw(20) << jvCheck[i][0]
1448 << std::setw(20) << jvCheck[i][1]
1449 << std::setw(20) << jvCheck[i][2]
1450 << std::setw(20) << jvCheck[i][3]
1452 outStream << hist.str();
1458 outStream.copyfmt(oldFormatState);
1470 const bool printToStream =
true,
1471 std::ostream & outStream = std::cout,
1473 const int order = 1 ) {
1474 std::vector<Real> steps(numSteps);
1475 for(
int i=0;i<numSteps;++i) {
1476 steps[i] = pow(10,-i);
1488 const std::vector<Real> &steps,
1489 const bool printToStream =
true,
1490 std::ostream & outStream = std::cout,
1491 const int order = 1 ) {
1497 Real tol = std::sqrt(ROL_EPSILON<Real>());
1499 int numSteps = steps.size();
1501 std::vector<Real> tmp(numVals);
1502 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1505 Ptr<Vector<Real>> AJdif = hv.
clone();
1506 Ptr<Vector<Real>> AJp = hv.
clone();
1507 Ptr<Vector<Real>> AHpv = hv.
clone();
1508 Ptr<Vector<Real>> AJnew = hv.
clone();
1509 Ptr<Vector<Real>> unew = u.
clone();
1512 nullstream oldFormatState;
1513 oldFormatState.copyfmt(outStream);
1521 Real normAHpv = AHpv->norm();
1523 for (
int i=0; i<numSteps; i++) {
1525 Real eta = steps[i];
1531 AJdif->scale(weights[order-1][0]);
1533 for(
int j=0; j<order; ++j) {
1535 unew->axpy(eta*shifts[order-1][j],v);
1537 if( weights[order-1][j+1] != 0 ) {
1540 AJdif->axpy(weights[order-1][j+1],*AJnew);
1544 AJdif->scale(one/eta);
1547 ahpvCheck[i][0] = eta;
1548 ahpvCheck[i][1] = normAHpv;
1549 ahpvCheck[i][2] = AJdif->norm();
1550 AJdif->axpy(-one, *AHpv);
1551 ahpvCheck[i][3] = AJdif->norm();
1553 if (printToStream) {
1554 std::stringstream hist;
1557 << std::setw(20) <<
"Step size"
1558 << std::setw(20) <<
"norm(adj(H)(u,v))"
1559 << std::setw(20) <<
"norm(FD approx)"
1560 << std::setw(20) <<
"norm(abs error)"
1562 << std::setw(20) <<
"---------"
1563 << std::setw(20) <<
"-----------------"
1564 << std::setw(20) <<
"---------------"
1565 << std::setw(20) <<
"---------------"
1568 hist << std::scientific << std::setprecision(11) << std::right
1569 << std::setw(20) << ahpvCheck[i][0]
1570 << std::setw(20) << ahpvCheck[i][1]
1571 << std::setw(20) << ahpvCheck[i][2]
1572 << std::setw(20) << ahpvCheck[i][3]
1574 outStream << hist.str();
1580 outStream.copyfmt(oldFormatState);
1593 const bool printToStream =
true,
1594 std::ostream & outStream = std::cout,
1596 const int order = 1 ) {
1597 std::vector<Real> steps(numSteps);
1598 for(
int i=0;i<numSteps;++i) {
1599 steps[i] = pow(10,-i);
1614 const std::vector<Real> &steps,
1615 const bool printToStream =
true,
1616 std::ostream & outStream = std::cout,
1617 const int order = 1 ) {
1623 Real tol = std::sqrt(ROL_EPSILON<Real>());
1625 int numSteps = steps.size();
1627 std::vector<Real> tmp(numVals);
1628 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1631 Ptr<Vector<Real>> AJdif = hv.
clone();
1632 Ptr<Vector<Real>> AJp = hv.
clone();
1633 Ptr<Vector<Real>> AHpv = hv.
clone();
1634 Ptr<Vector<Real>> AJnew = hv.
clone();
1635 Ptr<Vector<Real>> znew = z.
clone();
1638 nullstream oldFormatState;
1639 oldFormatState.copyfmt(outStream);
1647 Real normAHpv = AHpv->norm();
1649 for (
int i=0; i<numSteps; i++) {
1651 Real eta = steps[i];
1657 AJdif->scale(weights[order-1][0]);
1659 for(
int j=0; j<order; ++j) {
1661 znew->axpy(eta*shifts[order-1][j],v);
1663 if( weights[order-1][j+1] != 0 ) {
1666 AJdif->axpy(weights[order-1][j+1],*AJnew);
1670 AJdif->scale(one/eta);
1673 ahpvCheck[i][0] = eta;
1674 ahpvCheck[i][1] = normAHpv;
1675 ahpvCheck[i][2] = AJdif->norm();
1676 AJdif->axpy(-one, *AHpv);
1677 ahpvCheck[i][3] = AJdif->norm();
1679 if (printToStream) {
1680 std::stringstream hist;
1683 << std::setw(20) <<
"Step size"
1684 << std::setw(20) <<
"norm(adj(H)(u,v))"
1685 << std::setw(20) <<
"norm(FD approx)"
1686 << std::setw(20) <<
"norm(abs error)"
1688 << std::setw(20) <<
"---------"
1689 << std::setw(20) <<
"-----------------"
1690 << std::setw(20) <<
"---------------"
1691 << std::setw(20) <<
"---------------"
1694 hist << std::scientific << std::setprecision(11) << std::right
1695 << std::setw(20) << ahpvCheck[i][0]
1696 << std::setw(20) << ahpvCheck[i][1]
1697 << std::setw(20) << ahpvCheck[i][2]
1698 << std::setw(20) << ahpvCheck[i][3]
1700 outStream << hist.str();
1706 outStream.copyfmt(oldFormatState);
1719 const bool printToStream =
true,
1720 std::ostream & outStream = std::cout,
1722 const int order = 1 ) {
1723 std::vector<Real> steps(numSteps);
1724 for(
int i=0;i<numSteps;++i) {
1725 steps[i] = pow(10,-i);
1738 const std::vector<Real> &steps,
1739 const bool printToStream =
true,
1740 std::ostream & outStream = std::cout,
1741 const int order = 1 ) {
1747 Real tol = std::sqrt(ROL_EPSILON<Real>());
1749 int numSteps = steps.size();
1751 std::vector<Real> tmp(numVals);
1752 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1755 Ptr<Vector<Real>> AJdif = hv.
clone();
1756 Ptr<Vector<Real>> AJp = hv.
clone();
1757 Ptr<Vector<Real>> AHpv = hv.
clone();
1758 Ptr<Vector<Real>> AJnew = hv.
clone();
1759 Ptr<Vector<Real>> unew = u.
clone();
1762 nullstream oldFormatState;
1763 oldFormatState.copyfmt(outStream);
1771 Real normAHpv = AHpv->norm();
1773 for (
int i=0; i<numSteps; i++) {
1775 Real eta = steps[i];
1781 AJdif->scale(weights[order-1][0]);
1783 for(
int j=0; j<order; ++j) {
1785 unew->axpy(eta*shifts[order-1][j],v);
1787 if( weights[order-1][j+1] != 0 ) {
1790 AJdif->axpy(weights[order-1][j+1],*AJnew);
1794 AJdif->scale(one/eta);
1797 ahpvCheck[i][0] = eta;
1798 ahpvCheck[i][1] = normAHpv;
1799 ahpvCheck[i][2] = AJdif->norm();
1800 AJdif->axpy(-one, *AHpv);
1801 ahpvCheck[i][3] = AJdif->norm();
1803 if (printToStream) {
1804 std::stringstream hist;
1807 << std::setw(20) <<
"Step size"
1808 << std::setw(20) <<
"norm(adj(H)(u,v))"
1809 << std::setw(20) <<
"norm(FD approx)"
1810 << std::setw(20) <<
"norm(abs error)"
1812 << std::setw(20) <<
"---------"
1813 << std::setw(20) <<
"-----------------"
1814 << std::setw(20) <<
"---------------"
1815 << std::setw(20) <<
"---------------"
1818 hist << std::scientific << std::setprecision(11) << std::right
1819 << std::setw(20) << ahpvCheck[i][0]
1820 << std::setw(20) << ahpvCheck[i][1]
1821 << std::setw(20) << ahpvCheck[i][2]
1822 << std::setw(20) << ahpvCheck[i][3]
1824 outStream << hist.str();
1830 outStream.copyfmt(oldFormatState);
1840 const bool printToStream =
true,
1841 std::ostream & outStream = std::cout,
1843 const int order = 1 ) {
1844 std::vector<Real> steps(numSteps);
1845 for(
int i=0;i<numSteps;++i) {
1846 steps[i] = pow(10,-i);
1858 const std::vector<Real> &steps,
1859 const bool printToStream =
true,
1860 std::ostream & outStream = std::cout,
1861 const int order = 1 ) {
1867 Real tol = std::sqrt(ROL_EPSILON<Real>());
1869 int numSteps = steps.size();
1871 std::vector<Real> tmp(numVals);
1872 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp);
1875 Ptr<Vector<Real>> AJdif = hv.
clone();
1876 Ptr<Vector<Real>> AJp = hv.
clone();
1877 Ptr<Vector<Real>> AHpv = hv.
clone();
1878 Ptr<Vector<Real>> AJnew = hv.
clone();
1879 Ptr<Vector<Real>> znew = z.
clone();
1882 nullstream oldFormatState;
1883 oldFormatState.copyfmt(outStream);
1891 Real normAHpv = AHpv->norm();
1893 for (
int i=0; i<numSteps; i++) {
1895 Real eta = steps[i];
1901 AJdif->scale(weights[order-1][0]);
1903 for(
int j=0; j<order; ++j) {
1905 znew->axpy(eta*shifts[order-1][j],v);
1907 if( weights[order-1][j+1] != 0 ) {
1910 AJdif->axpy(weights[order-1][j+1],*AJnew);
1914 AJdif->scale(one/eta);
1917 ahpvCheck[i][0] = eta;
1918 ahpvCheck[i][1] = normAHpv;
1919 ahpvCheck[i][2] = AJdif->norm();
1920 AJdif->axpy(-one, *AHpv);
1921 ahpvCheck[i][3] = AJdif->norm();
1923 if (printToStream) {
1924 std::stringstream hist;
1927 << std::setw(20) <<
"Step size"
1928 << std::setw(20) <<
"norm(adj(H)(u,v))"
1929 << std::setw(20) <<
"norm(FD approx)"
1930 << std::setw(20) <<
"norm(abs error)"
1932 << std::setw(20) <<
"---------"
1933 << std::setw(20) <<
"-----------------"
1934 << std::setw(20) <<
"---------------"
1935 << std::setw(20) <<
"---------------"
1938 hist << std::scientific << std::setprecision(11) << std::right
1939 << std::setw(20) << ahpvCheck[i][0]
1940 << std::setw(20) << ahpvCheck[i][1]
1941 << std::setw(20) << ahpvCheck[i][2]
1942 << std::setw(20) << ahpvCheck[i][3]
1944 outStream << hist.str();
1950 outStream.copyfmt(oldFormatState);
Contains definitions of custom data types in ROL.
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
Defines the constraint operator interface for simulation-based optimization.
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 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 update_2(const Vector< Real > &z, UpdateType type, int iter=-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 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.
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)
Ptr< Vector< Real > > unew_
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 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,...
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 .
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 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...
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 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,...
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,...
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 ...
virtual void value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol)=0
Evaluate the constraint operator at .
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 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...
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)
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)
, , , ,
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,...
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update constraint function.
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...
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 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 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.
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
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...
const int DEFAULT_solverType_
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)
virtual void setSolveParameters(ParameterList &parlist)
Set solve parameters.
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)
const Real DEFAULT_factor_
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)
, , , ,
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 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 update_1(const Vector< Real > &u, UpdateType type, int iter=-1)
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)
const bool DEFAULT_print_
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 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 ,...
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)
virtual void update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-1)
Ptr< Vector< Real > > jv_
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)
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
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 void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
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...
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)
virtual void solve_update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-1)
Update SimOpt constraint during solve (disconnected from optimization updates).
virtual Real checkSolve(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &c, const bool printToStream=true, std::ostream &outStream=std::cout)
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
Defines the general constraint operator interface.
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 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
Defines the linear algebra or vector space interface for simulation-based optimization.
ROL::Ptr< const Vector< Real > > get_2() const
ROL::Ptr< const Vector< Real > > get_1() const
void set_1(const Vector< Real > &vec)
void set_2(const Vector< Real > &vec)
const Vector< Real > & dual(void) const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
Defines the linear algebra or vector space interface.
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
virtual void scale(const Real alpha)=0
Compute where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual void plus(const Vector &x)=0
Compute , where .
virtual void zero()
Set to zero vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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 axpy(const Real alpha, const Vector &x)
Compute where .
const double weights[4][5]
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.