18 #include "Teuchos_GlobalMPISession.hpp"
29 template <
class Real,
class Element=Real>
32 template <
class Real,
class Element=Real>
35 template <
class Real,
class Element=Real>
38 template <
class Real,
class Element=Real>
44 template <
class Real,
class Element>
52 ROL::Ptr<std::vector<Element> >
std_vec_;
53 mutable ROL::Ptr<OptDualStdVector<Real> >
dual_vec_;
66 (*std_vec_)[i] += (*xp)[i];
70 void scale(
const Real alpha ) {
73 (*std_vec_)[i] *= alpha;
86 val += (*std_vec_)[i]*(*xp)[i];
93 val = std::sqrt(
dot(*
this) );
97 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
98 return ROL::makePtr<OptStdVector>( ROL::makePtr<std::vector<Element>>(
std_vec_->size()) );
101 ROL::Ptr<const std::vector<Element> >
getVector()
const {
109 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
111 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(
std_vec_->size(),0.0);
112 ROL::Ptr<V> e = ROL::makePtr<OptStdVector>( e_ptr );
120 dual_vec_ = ROL::makePtr<OptDualStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
130 val += (*std_vec_)[i]*(*xvalptr)[i];
139 template <
class Real,
class Element>
147 ROL::Ptr<std::vector<Element> >
std_vec_;
148 mutable ROL::Ptr<OptStdVector<Real> >
dual_vec_;
159 (*std_vec_)[i] += (*xp)[i];
166 (*std_vec_)[i] *= alpha;
176 val += (*std_vec_)[i]*(*xp)[i];
183 val = std::sqrt(
dot(*
this) );
187 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
188 return ROL::makePtr<OptDualStdVector>( ROL::makePtr<std::vector<Element>>(
std_vec_->size()) );
191 ROL::Ptr<const std::vector<Element> >
getVector()
const {
199 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
201 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(
std_vec_->size(), 0.0 );
202 ROL::Ptr<V> e = ROL::makePtr<OptDualStdVector>( e_ptr );
210 dual_vec_ = ROL::makePtr<OptStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
220 val += (*std_vec_)[i]*(*xvalptr)[i];
229 template <
class Real,
class Element>
249 (*std_vec_)[i] += (*xp)[i];
256 (*std_vec_)[i] *= alpha;
266 val += (*std_vec_)[i]*(*xp)[i];
273 val = std::sqrt(
dot(*
this) );
277 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
278 return ROL::makePtr<ConStdVector>( ROL::makePtr<std::vector<Element>>(
std_vec_->size()));
281 ROL::Ptr<const std::vector<Element> >
getVector()
const {
289 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
291 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(
std_vec_->size(), 0.0);
292 ROL::Ptr<V> e = ROL::makePtr<ConStdVector>( e_ptr );
300 dual_vec_ = ROL::makePtr<ConDualStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
310 val += (*std_vec_)[i]*(*xvalptr)[i];
319 template <
class Real,
class Element>
340 (*std_vec_)[i] += (*xp)[i];
347 (*std_vec_)[i] *= alpha;
357 val += (*std_vec_)[i]*(*xp)[i];
364 val = std::sqrt(
dot(*
this) );
368 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
369 return ROL::makePtr<ConDualStdVector>( ROL::makePtr<std::vector<Element>>(
std_vec_->size()));
372 ROL::Ptr<const std::vector<Element> >
getVector()
const {
380 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
382 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(
std_vec_->size(),0.0);
383 ROL::Ptr<V> e = ROL::makePtr<ConDualStdVector>(e_ptr);
391 dual_vec_ = ROL::makePtr<ConStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
401 val += (*std_vec_)[i]*(*xvalptr)[i];
411 int main(
int argc,
char *argv[]) {
414 typedef std::vector<RealT> vector;
417 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
420 int iprint = argc - 1;
421 ROL::Ptr<std::ostream> outStream;
424 outStream = ROL::makePtrFromRef(std::cout);
426 outStream = ROL::makePtrFromRef(bhs);
436 ROL::Ptr<ROL::Objective<RealT> > obj;
437 ROL::Ptr<ROL::Constraint<RealT> > constr;
438 ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(
dim, 0.0);
439 ROL::Ptr<vector> sol_ptr = ROL::makePtr<vector>(
dim, 0.0);
440 ROL::Ptr<OptStdVector<RealT>> x = ROL::makePtr<OptStdVector<RealT>>(x_ptr);
441 ROL::Ptr<OptStdVector<RealT>> sol = ROL::makePtr<OptStdVector<RealT>>(sol_ptr);
445 obj = SEC.getObjective();
446 constr = SEC.getEqualityConstraint();
447 x->
set(*SEC.getInitialGuess());
448 sol->set(*SEC.getSolution());
451 RealT left = -1e0, right = 1e0;
452 ROL::Ptr<vector> xtest_ptr = ROL::makePtr<vector>(
dim, 0.0);
453 ROL::Ptr<vector> g_ptr = ROL::makePtr<vector>(
dim, 0.0);
454 ROL::Ptr<vector> d_ptr = ROL::makePtr<vector>(
dim, 0.0);
455 ROL::Ptr<vector> gd_ptr = ROL::makePtr<vector>(
dim, 0.0);
456 ROL::Ptr<vector> v_ptr = ROL::makePtr<vector>(
dim, 0.0);
457 ROL::Ptr<vector> vc_ptr = ROL::makePtr<vector>(nc, 0.0);
458 ROL::Ptr<vector> vl_ptr = ROL::makePtr<vector>(nc, 0.0);
459 ROL::Ptr<OptStdVector<RealT>> xtest = ROL::makePtr<OptStdVector<RealT>>(xtest_ptr);
460 ROL::Ptr<OptDualStdVector<RealT>> g = ROL::makePtr<OptDualStdVector<RealT>>(g_ptr);
461 ROL::Ptr<OptStdVector<RealT>> d = ROL::makePtr<OptStdVector<RealT>>(d_ptr);
462 ROL::Ptr<OptDualStdVector<RealT>> gd = ROL::makePtr<OptDualStdVector<RealT>>(gd_ptr);
463 ROL::Ptr<OptStdVector<RealT>> v = ROL::makePtr<OptStdVector<RealT>>(v_ptr);
464 ROL::Ptr<ConStdVector<RealT>> vc = ROL::makePtr<ConStdVector<RealT>>(vc_ptr);
465 ROL::Ptr<ConDualStdVector<RealT>> vl = ROL::makePtr<ConDualStdVector<RealT>>(vl_ptr);
467 for (uint i=0; i<
dim; i++) {
468 (*xtest_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
469 (*d_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
470 (*gd_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
471 (*v_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
474 for (uint i=0; i<nc; i++) {
475 (*vc_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
476 (*vl_ptr)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
478 obj->checkGradient(*xtest, *g, *d,
true, *outStream); *outStream << std::endl;
479 obj->checkHessVec(*xtest, *g, *v,
true, *outStream); *outStream << std::endl;
480 obj->checkHessSym(*xtest, *g, *d, *v,
true, *outStream); *outStream << std::endl;
481 constr->checkApplyJacobian(*xtest, *v, *vc,
true, *outStream); *outStream << std::endl;
482 constr->checkApplyAdjointJacobian(*xtest, *vl, *vc, *g,
true, *outStream); *outStream << std::endl;
483 constr->checkApplyAdjointHessian(*xtest, *vl, *d, *g,
true, *outStream); *outStream << std::endl;
485 ROL::Ptr<vector> v1_ptr = ROL::makePtr<vector>(
dim, 0.0);
486 ROL::Ptr<vector> v2_ptr = ROL::makePtr<vector>(nc, 0.0);
487 ROL::Ptr<OptStdVector<RealT>> v1 = ROL::makePtr<OptStdVector<RealT>>(v1_ptr);
488 ROL::Ptr<ConDualStdVector<RealT>> v2 = ROL::makePtr<ConDualStdVector<RealT>>(v2_ptr);
490 constr->solveAugmentedSystem(*v1, *v2, *gd, *vc, *xtest, augtol);
493 vl->zero(); vc->zero(); g->zero();
494 ROL::Ptr<ROL::Problem<RealT>>
495 problem = ROL::makePtr<ROL::Problem<RealT>>(obj,x,g);
497 problem->finalize(
false,
true,*outStream);
500 ROL::ParameterList parlist;
501 std::string stepname =
"Augmented Lagrangian";
502 parlist.sublist(
"General").set(
"Output Level",1);
503 parlist.sublist(
"Step").set(
"Type",stepname);
504 parlist.sublist(
"Step").sublist(stepname).sublist(
"Optimality System Solver").set(
"Nominal Relative Tolerance",1e-4);
505 parlist.sublist(
"Step").sublist(stepname).sublist(
"Optimality System Solver").set(
"Fix Tolerance",
true);
506 parlist.sublist(
"Step").sublist(stepname).sublist(
"Tangential Subproblem Solver").set(
"Iteration Limit",20);
507 parlist.sublist(
"Step").sublist(stepname).sublist(
"Tangential Subproblem Solver").set(
"Relative Tolerance",1e-2);
508 parlist.sublist(
"Step").sublist(stepname).set(
"Output Level",0);
509 parlist.sublist(
"Status Test").set(
"Gradient Tolerance",1.e-12);
510 parlist.sublist(
"Status Test").set(
"Constraint Tolerance",1.e-12);
511 parlist.sublist(
"Status Test").set(
"Step Tolerance",1.e-18);
512 parlist.sublist(
"Status Test").set(
"Iteration Limit",100);
513 ROL::Ptr<ROL::Solver<RealT>> solver
514 = ROL::makePtr<ROL::Solver<RealT>>(problem,parlist);
519 solver->solve(*outStream);
522 *outStream <<
"\nReference solution x_r =\n";
523 *outStream << std::scientific <<
" " << (*sol_ptr)[0] <<
"\n";
524 *outStream << std::scientific <<
" " << (*sol_ptr)[1] <<
"\n";
525 *outStream << std::scientific <<
" " << (*sol_ptr)[2] <<
"\n";
526 *outStream << std::scientific <<
" " << (*sol_ptr)[3] <<
"\n";
527 *outStream << std::scientific <<
" " << (*sol_ptr)[4] <<
"\n";
528 *outStream <<
"\nOptimal solution x =\n";
529 *outStream << std::scientific <<
" " << (*x_ptr)[0] <<
"\n";
530 *outStream << std::scientific <<
" " << (*x_ptr)[1] <<
"\n";
531 *outStream << std::scientific <<
" " << (*x_ptr)[2] <<
"\n";
532 *outStream << std::scientific <<
" " << (*x_ptr)[3] <<
"\n";
533 *outStream << std::scientific <<
" " << (*x_ptr)[4] <<
"\n";
535 RealT abserr = x->norm();
536 RealT relerr = abserr/sol->norm();
537 *outStream << std::scientific <<
"\n Absolute Error: " << abserr;
538 *outStream << std::scientific <<
"\n Relative Error: " << relerr <<
"\n";
539 if ( relerr > sqrt(ROL::ROL_EPSILON<RealT>()) ) {
543 catch (std::logic_error& err) {
544 *outStream << err.what() <<
"\n";
549 std::cout <<
"End Result: TEST FAILED\n";
551 std::cout <<
"End Result: TEST PASSED\n";
std::vector< Element > vector
ROL::Ptr< std::vector< Element > > std_vec_
typename PV< Real >::size_type size_type
void scale(const Real alpha)
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
std::vector< Element > vector
ROL::Ptr< const std::vector< Element > > getVector() const
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Real norm() const
Returns where .
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
ROL::Ptr< std::vector< Element > > getVector()
OptDualStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
ROL::Ptr< std::vector< Element > > std_vec_
ROL::Ptr< std::vector< Element > > getVector()
int dimension() const
Return dimension of the vector space.
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
std::vector< Element > vector
Defines the linear algebra or vector space interface.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void plus(const ROL::Vector< Real > &x)
Compute , where .
ROL::Ptr< OptStdVector< Real > > dual_vec_
ROL::Ptr< std::vector< Element > > std_vec_
ROL::Ptr< std::vector< Element > > getVector()
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
std::vector< Element > vector
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void scale(const Real alpha)
Compute where .
ConDualStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
void scale(const Real alpha)
Compute where .
ROL::Ptr< const std::vector< Element > > getVector() const
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
int dimension() const
Return dimension of the vector space.
int dimension() const
Return dimension of the vector space.
ROL::Ptr< ConStdVector< Real > > dual_vec_
ROL::Ptr< OptDualStdVector< Real > > dual_vec_
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
ConStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Real norm() const
Returns where .
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Contains definitions for the equality constrained NLP from Nocedal/Wright, 2nd edition, page 574, example 18.2; note the typo in reversing the initial guess and the solution.
int dimension() const
Return dimension of the vector space.
ROL::Ptr< std::vector< Element > > std_vec_
basic_nullstream< char, char_traits< char >> nullstream
int main(int argc, char *argv[])
Real norm() const
Returns where .
ROL::Ptr< const std::vector< Element > > getVector() const
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
virtual void set(const Vector &x)
Set where .
ROL::Ptr< const std::vector< Element > > getVector() const
ROL::Ptr< ConDualStdVector< Real > > dual_vec_
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void plus(const ROL::Vector< Real > &x)
Compute , where .
ROL::Ptr< std::vector< Element > > getVector()
void plus(const ROL::Vector< Real > &x)
Compute , where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
Real norm() const
Returns where .
OptStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void scale(const Real alpha)
Compute where .