445int main(
int argc,
char *argv[]) {
448 typedef std::vector<RealT> vector;
449 typedef vector::size_type uint;
451 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
454 int iprint = argc - 1;
455 ROL::Ptr<std::ostream> outStream;
458 outStream = ROL::makePtrFromRef(std::cout);
460 outStream = ROL::makePtrFromRef(bhs);
470 ROL::Ptr<ROL::Objective<RealT> > obj;
471 ROL::Ptr<ROL::Constraint<RealT> > constr;
472 ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(
dim, 0.0);
473 ROL::Ptr<vector> sol_ptr = ROL::makePtr<vector>(
dim, 0.0);
474 ROL::Ptr<OptStdVector<RealT>> x = ROL::makePtr<OptStdVector<RealT>>(x_ptr);
475 ROL::Ptr<OptStdVector<RealT>> sol = ROL::makePtr<OptStdVector<RealT>>(sol_ptr);
479 obj = SEC.getObjective();
480 constr = SEC.getEqualityConstraint();
481 x->
set(*SEC.getInitialGuess());
482 sol->set(*SEC.getSolution());
485 RealT left = -1e0, right = 1e0;
486 ROL::Ptr<vector> xtest_ptr = ROL::makePtr<vector>(
dim, 0.0);
487 ROL::Ptr<vector> g_ptr = ROL::makePtr<vector>(
dim, 0.0);
488 ROL::Ptr<vector> d_ptr = ROL::makePtr<vector>(
dim, 0.0);
489 ROL::Ptr<vector> gd_ptr = ROL::makePtr<vector>(
dim, 0.0);
490 ROL::Ptr<vector> v_ptr = ROL::makePtr<vector>(
dim, 0.0);
491 ROL::Ptr<vector> vc_ptr = ROL::makePtr<vector>(nc, 0.0);
492 ROL::Ptr<vector> vl_ptr = ROL::makePtr<vector>(nc, 0.0);
493 ROL::Ptr<OptStdVector<RealT>> xtest = ROL::makePtr<OptStdVector<RealT>>(xtest_ptr);
494 ROL::Ptr<OptDualStdVector<RealT>> g = ROL::makePtr<OptDualStdVector<RealT>>(g_ptr);
495 ROL::Ptr<OptStdVector<RealT>> d = ROL::makePtr<OptStdVector<RealT>>(d_ptr);
496 ROL::Ptr<OptDualStdVector<RealT>> gd = ROL::makePtr<OptDualStdVector<RealT>>(gd_ptr);
497 ROL::Ptr<OptStdVector<RealT>> v = ROL::makePtr<OptStdVector<RealT>>(v_ptr);
498 ROL::Ptr<ConStdVector<RealT>> vc = ROL::makePtr<ConStdVector<RealT>>(vc_ptr);
499 ROL::Ptr<ConDualStdVector<RealT>> vl = ROL::makePtr<ConDualStdVector<RealT>>(vl_ptr);
501 for (uint i=0; i<
dim; i++) {
502 (*xtest_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
503 (*d_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
504 (*gd_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
505 (*v_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
508 for (uint i=0; i<nc; i++) {
509 (*vc_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
510 (*vl_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
512 obj->checkGradient(*xtest, *g, *d,
true, *outStream); *outStream << std::endl;
513 obj->checkHessVec(*xtest, *g, *v,
true, *outStream); *outStream << std::endl;
514 obj->checkHessSym(*xtest, *g, *d, *v,
true, *outStream); *outStream << std::endl;
515 constr->checkApplyJacobian(*xtest, *v, *vc,
true, *outStream); *outStream << std::endl;
516 constr->checkApplyAdjointJacobian(*xtest, *vl, *vc, *g,
true, *outStream); *outStream << std::endl;
517 constr->checkApplyAdjointHessian(*xtest, *vl, *d, *g,
true, *outStream); *outStream << std::endl;
519 ROL::Ptr<vector> v1_ptr = ROL::makePtr<vector>(
dim, 0.0);
520 ROL::Ptr<vector> v2_ptr = ROL::makePtr<vector>(nc, 0.0);
521 ROL::Ptr<OptStdVector<RealT>> v1 = ROL::makePtr<OptStdVector<RealT>>(v1_ptr);
522 ROL::Ptr<ConDualStdVector<RealT>> v2 = ROL::makePtr<ConDualStdVector<RealT>>(v2_ptr);
524 constr->solveAugmentedSystem(*v1, *v2, *gd, *vc, *xtest, augtol);
527 vl->zero(); vc->zero(); g->zero();
528 ROL::Ptr<ROL::Problem<RealT>>
529 problem = ROL::makePtr<ROL::Problem<RealT>>(obj,x,g);
530 problem->addConstraint(
"Equality",constr,ROL::staticPtrCast<
ROL::Vector<RealT>>(vl),ROL::staticPtrCast<ROL::Vector<RealT>>(vc));
531 problem->finalize(
false,
true,*outStream);
534 ROL::ParameterList parlist;
535 std::string stepname =
"Augmented Lagrangian";
536 parlist.sublist(
"General").set(
"Output Level",1);
537 parlist.sublist(
"Step").set(
"Type",stepname);
538 parlist.sublist(
"Step").sublist(stepname).sublist(
"Optimality System Solver").set(
"Nominal Relative Tolerance",1e-4);
539 parlist.sublist(
"Step").sublist(stepname).sublist(
"Optimality System Solver").set(
"Fix Tolerance",
true);
540 parlist.sublist(
"Step").sublist(stepname).sublist(
"Tangential Subproblem Solver").set(
"Iteration Limit",20);
541 parlist.sublist(
"Step").sublist(stepname).sublist(
"Tangential Subproblem Solver").set(
"Relative Tolerance",1e-2);
542 parlist.sublist(
"Step").sublist(stepname).set(
"Output Level",0);
543 parlist.sublist(
"Status Test").set(
"Gradient Tolerance",1.e-12);
544 parlist.sublist(
"Status Test").set(
"Constraint Tolerance",1.e-12);
545 parlist.sublist(
"Status Test").set(
"Step Tolerance",1.e-18);
546 parlist.sublist(
"Status Test").set(
"Iteration Limit",100);
547 ROL::Ptr<ROL::Solver<RealT>> solver
548 = ROL::makePtr<ROL::Solver<RealT>>(problem,parlist);
553 solver->solve(*outStream);
556 *outStream <<
"\nReference solution x_r =\n";
557 *outStream << std::scientific <<
" " << (*sol_ptr)[0] <<
"\n";
558 *outStream << std::scientific <<
" " << (*sol_ptr)[1] <<
"\n";
559 *outStream << std::scientific <<
" " << (*sol_ptr)[2] <<
"\n";
560 *outStream << std::scientific <<
" " << (*sol_ptr)[3] <<
"\n";
561 *outStream << std::scientific <<
" " << (*sol_ptr)[4] <<
"\n";
562 *outStream <<
"\nOptimal solution x =\n";
563 *outStream << std::scientific <<
" " << (*x_ptr)[0] <<
"\n";
564 *outStream << std::scientific <<
" " << (*x_ptr)[1] <<
"\n";
565 *outStream << std::scientific <<
" " << (*x_ptr)[2] <<
"\n";
566 *outStream << std::scientific <<
" " << (*x_ptr)[3] <<
"\n";
567 *outStream << std::scientific <<
" " << (*x_ptr)[4] <<
"\n";
569 RealT abserr = x->norm();
570 RealT relerr = abserr/sol->norm();
571 *outStream << std::scientific <<
"\n Absolute Error: " << abserr;
572 *outStream << std::scientific <<
"\n Relative Error: " << relerr <<
"\n";
573 if ( relerr > sqrt(ROL::ROL_EPSILON<RealT>()) ) {
577 catch (std::logic_error& err) {
578 *outStream << err.what() <<
"\n";
583 std::cout <<
"End Result: TEST FAILED\n";
585 std::cout <<
"End Result: TEST PASSED\n";