44 #ifndef ROL_CONSTRAINT_SIMOPT_H
45 #define ROL_CONSTRAINT_SIMOPT_H
107 template <
class Real>
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>>();
301 Ptr<Vector<Real>> l = c.
dual().clone();
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);
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);
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();
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();
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();
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();
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);
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
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 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 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 ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Ptr< const Vector< Real > > get_2() const
virtual int dimension() const
Return dimension of the vector space.
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
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 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 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 void axpy(const Real alpha, const Vector &x)
Compute where .
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 .
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 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 update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update constraint function.
Defines the linear algebra or vector space interface for simulation-based optimization.
virtual void update_1(const Vector< Real > &u, UpdateType type, int iter=-1)
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 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.
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 is ...
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, flag = true if optimization variable is changed, iter is the outer algorithm iterations count.
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.
void set_1(const Vector< Real > &vec)
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 void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply 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 void update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-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 . In general, this preconditioner satisfies the fo...
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
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 is ...
Ptr< ostream > makeStreamPtr(ostream &os, bool noSuppressOutput=true)
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)
Ptr< Vector< Real > > jv_
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 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 .
const int DEFAULT_solverType_
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 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.
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 value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol)=0
Evaluate the constraint operator at .
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 > > 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)
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 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, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
const Real DEFAULT_factor_
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
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...
basic_nullstream< char, char_traits< char >> nullstream
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol)
Given , solve for .
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 update_2(const Vector< Real > &z, UpdateType type, int iter=-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 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, for use with dual spaces where the user does not define the dual() operation.
virtual void setSolveParameters(ParameterList &parlist)
Set solve parameters.
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 Real norm() const =0
Returns where .
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.
Ptr< Vector< Real > > unew_
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)
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Defines the constraint operator interface for simulation-based optimization.
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 .
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)
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)
, , , ,
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 ...
void set_2(const Vector< Real > &vec)
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)
, , , ,
Defines the general constraint operator interface.
const bool DEFAULT_print_
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 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 Vector< Real > & dual(void) const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
ROL::Ptr< const Vector< Real > > get_1() const