53 #include "Teuchos_oblackholestream.hpp" 54 #include "Teuchos_GlobalMPISession.hpp" 64 template <
class Real,
class Element=Real>
67 template <
class Real,
class Element=Real>
70 template <
class Real,
class Element=Real>
73 template <
class Real,
class Element=Real>
79 template <
class Real,
class Element>
84 typedef typename vector::size_type
uint;
87 Teuchos::RCP<std::vector<Element> >
std_vec_;
88 mutable Teuchos::RCP<OptDualStdVector<Real> >
dual_vec_;
92 OptStdVector(
const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
95 using Teuchos::RCP;
using Teuchos::dyn_cast;
101 (*std_vec_)[i] += (*xp)[i];
108 (*std_vec_)[i] *= alpha;
114 using Teuchos::RCP;
using Teuchos::dyn_cast;
121 val += (*std_vec_)[i]*(*xp)[i];
128 val = std::sqrt(
dot(*
this) );
132 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
133 return Teuchos::rcp(
new OptStdVector( Teuchos::rcp(
new std::vector<Element>(std_vec_->size()) ) ) );
136 Teuchos::RCP<const std::vector<Element> >
getVector()
const {
144 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
145 using Teuchos::RCP;
using Teuchos::rcp;
146 RCP<vector> e_rcp = rcp(
new vector(std_vec_->size(),0.0) );
152 int dimension()
const {
return static_cast<int>(std_vec_->size());}
155 dual_vec_ = Teuchos::rcp(
new OptDualStdVector<Real>( Teuchos::rcp(
new std::vector<Element>(*std_vec_) ) ) );
163 template <
class Real,
class Element>
168 typedef typename vector::size_type
uint;
171 Teuchos::RCP<std::vector<Element> >
std_vec_;
172 mutable Teuchos::RCP<OptStdVector<Real> >
dual_vec_;
176 OptDualStdVector(
const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
179 using Teuchos::RCP;
using Teuchos::dyn_cast;
183 (*std_vec_)[i] += (*xp)[i];
190 (*std_vec_)[i] *= alpha;
195 using Teuchos::RCP;
using Teuchos::dyn_cast;
200 val += (*std_vec_)[i]*(*xp)[i];
207 val = std::sqrt(
dot(*
this) );
211 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
212 return Teuchos::rcp(
new OptDualStdVector( Teuchos::rcp(
new std::vector<Element>(std_vec_->size()) ) ) );
215 Teuchos::RCP<const std::vector<Element> >
getVector()
const {
223 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
224 using Teuchos::RCP;
using Teuchos::rcp;
225 RCP<vector> e_rcp = rcp(
new vector( std_vec_->size(), 0.0 ) );
231 int dimension()
const {
return static_cast<int>(std_vec_->size());}
234 dual_vec_ = Teuchos::rcp(
new OptStdVector<Real>( Teuchos::rcp(
new std::vector<Element>(*std_vec_) ) ) );
242 template <
class Real,
class Element>
247 typedef typename vector::size_type
uint;
251 mutable Teuchos::RCP<ConDualStdVector<Real> >
dual_vec_;
255 ConStdVector(
const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
258 using Teuchos::RCP;
using Teuchos::dyn_cast;
262 (*std_vec_)[i] += (*xp)[i];
269 (*std_vec_)[i] *= alpha;
274 using Teuchos::RCP;
using Teuchos::dyn_cast;
279 val += (*std_vec_)[i]*(*xp)[i];
286 val = std::sqrt(
dot(*
this) );
290 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
291 return Teuchos::rcp(
new ConStdVector( Teuchos::rcp(
new std::vector<Element>(std_vec_->size())) ) );
294 Teuchos::RCP<const std::vector<Element> >
getVector()
const {
302 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
303 using Teuchos::RCP;
using Teuchos::rcp;
304 RCP<vector> e_rcp = rcp(
new vector(std_vec_->size(), 0.0) );
310 int dimension()
const {
return static_cast<int>(std_vec_->size());}
313 dual_vec_ = Teuchos::rcp(
new ConDualStdVector<Real>( Teuchos::rcp(
new std::vector<Element>(*std_vec_) ) ) );
321 template <
class Real,
class Element>
326 typedef typename vector::size_type
uint;
335 ConDualStdVector(
const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
338 using Teuchos::RCP;
using Teuchos::dyn_cast;
342 (*std_vec_)[i] += (*xp)[i];
349 (*std_vec_)[i] *= alpha;
354 using Teuchos::RCP;
using Teuchos::dyn_cast;
359 val += (*std_vec_)[i]*(*xp)[i];
366 val = std::sqrt(
dot(*
this) );
370 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
371 return Teuchos::rcp(
new ConDualStdVector( Teuchos::rcp(
new std::vector<Element>(std_vec_->size())) ) );
374 Teuchos::RCP<const std::vector<Element> >
getVector()
const {
382 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
383 using Teuchos::RCP;
using Teuchos::rcp;
384 RCP<vector> e_rcp = rcp(
new vector(std_vec_->size(),0.0) );
390 int dimension()
const {
return static_cast<int>(std_vec_->size());}
393 dual_vec_ = Teuchos::rcp(
new ConStdVector<Real>( Teuchos::rcp(
new std::vector<Element>(*std_vec_) ) ) );
402 int main(
int argc,
char *argv[]) {
404 typedef std::vector<RealT>
vector;
405 typedef vector::size_type
uint;
407 using Teuchos::RCP;
using Teuchos::rcp;
409 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
412 int iprint = argc - 1;
413 Teuchos::RCP<std::ostream> outStream;
414 Teuchos::oblackholestream bhs;
416 outStream = Teuchos::rcp(&std::cout,
false);
418 outStream = Teuchos::rcp(&bhs,
false);
426 RCP<ROL::Objective<RealT> > obj;
427 RCP<ROL::EqualityConstraint<RealT> > constr;
428 RCP<vector> x_rcp = rcp(
new vector(0, 0.0) );
429 RCP<vector> sol_rcp = rcp(
new vector(0, 0.0) );
439 RealT left = -1e0, right = 1e0;
440 RCP<vector> xtest_rcp = rcp(
new vector(dim, 0.0) );
441 RCP<vector> g_rcp = rcp(
new vector(dim, 0.0) );
442 RCP<vector> d_rcp = rcp(
new vector(dim, 0.0) );
443 RCP<vector> gd_rcp = rcp(
new vector(dim, 0.0) );
444 RCP<vector> v_rcp = rcp(
new vector(dim, 0.0) );
445 RCP<vector> vc_rcp = rcp(
new vector(nc, 0.0) );
446 RCP<vector> vl_rcp = rcp(
new vector(nc, 0.0) );
448 OptDualStdVector<RealT> g(g_rcp);
450 OptDualStdVector<RealT> gd(gd_rcp);
452 ConStdVector<RealT> vc(vc_rcp);
455 for (uint i=0; i<dim; i++) {
456 (*xtest_rcp)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
457 (*d_rcp)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
458 (*gd_rcp)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
459 (*v_rcp)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
462 for (uint i=0; i<nc; i++) {
463 (*vc_rcp)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
464 (*vl_rcp)[i] = ( (
RealT)rand() / (
RealT)RAND_MAX ) * (right - left) + left;
466 obj->checkGradient(xtest, g, d,
true, *outStream); *outStream <<
"\n";
467 obj->checkHessVec(xtest, g, v,
true, *outStream); *outStream <<
"\n";
468 obj->checkHessSym(xtest, g, d, v,
true, *outStream); *outStream <<
"\n";
469 constr->checkApplyJacobian(xtest, v, vc,
true, *outStream); *outStream <<
"\n";
470 constr->checkApplyAdjointJacobian(xtest, vl, vc, g,
true, *outStream); *outStream <<
"\n";
471 constr->checkApplyAdjointHessian(xtest, vl, d, g,
true, *outStream); *outStream <<
"\n";
473 RCP<vector> v1_rcp = rcp(
new vector(dim, 0.0) );
474 RCP<vector> v2_rcp = rcp(
new vector(nc, 0.0) );
478 constr->solveAugmentedSystem(v1, v2, gd, vc, xtest, augtol);
482 Teuchos::ParameterList parlist;
483 std::string stepname =
"Composite Step";
484 parlist.sublist(
"Step").sublist(stepname).sublist(
"Optimality System Solver").set(
"Nominal Relative Tolerance",1e-4);
485 parlist.sublist(
"Step").sublist(stepname).sublist(
"Optimality System Solver").set(
"Fix Tolerance",
true);
486 parlist.sublist(
"Step").sublist(stepname).sublist(
"Tangential Subproblem Solver").set(
"Iteration Limit",20);
487 parlist.sublist(
"Step").sublist(stepname).sublist(
"Tangential Subproblem Solver").set(
"Relative Tolerance",1e-2);
488 parlist.sublist(
"Step").sublist(stepname).set(
"Output Level",0);
489 parlist.sublist(
"Status Test").set(
"Gradient Tolerance",1.e-12);
490 parlist.sublist(
"Status Test").set(
"Constraint Tolerance",1.e-12);
491 parlist.sublist(
"Status Test").set(
"Step Tolerance",1.e-18);
492 parlist.sublist(
"Status Test").set(
"Iteration Limit",100);
499 algo.
run(x, g, vl, vc, *obj, *constr,
true, *outStream);
502 *outStream <<
"\nReference solution x_r =\n";
503 *outStream << std::scientific <<
" " << (*sol_rcp)[0] <<
"\n";
504 *outStream << std::scientific <<
" " << (*sol_rcp)[1] <<
"\n";
505 *outStream << std::scientific <<
" " << (*sol_rcp)[2] <<
"\n";
506 *outStream << std::scientific <<
" " << (*sol_rcp)[3] <<
"\n";
507 *outStream << std::scientific <<
" " << (*sol_rcp)[4] <<
"\n";
508 *outStream <<
"\nOptimal solution x =\n";
509 *outStream << std::scientific <<
" " << (*x_rcp)[0] <<
"\n";
510 *outStream << std::scientific <<
" " << (*x_rcp)[1] <<
"\n";
511 *outStream << std::scientific <<
" " << (*x_rcp)[2] <<
"\n";
512 *outStream << std::scientific <<
" " << (*x_rcp)[3] <<
"\n";
513 *outStream << std::scientific <<
" " << (*x_rcp)[4] <<
"\n";
515 RealT abserr = x.norm();
517 *outStream << std::scientific <<
"\n Absolute Error: " << abserr;
518 *outStream << std::scientific <<
"\n Relative Error: " << relerr <<
"\n";
523 catch (std::logic_error err) {
524 *outStream << err.what() <<
"\n";
529 std::cout <<
"End Result: TEST FAILED\n";
531 std::cout <<
"End Result: TEST PASSED\n";
std::vector< Element > vector
Teuchos::RCP< std::vector< Element > > std_vec_
Teuchos::RCP< OptDualStdVector< Real > > dual_vec_
void scale(const Real alpha)
Compute where .
OptStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec)
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
std::vector< Element > vector
ConStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec)
Real norm() const
Returns where .
Teuchos::RCP< std::vector< Element > > getVector()
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Teuchos::RCP< std::vector< Element > > getVector()
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Teuchos::RCP< const std::vector< Element > > getVector() const
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
int dimension() const
Return dimension of the vector space.
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Teuchos::RCP< std::vector< Element > > std_vec_
virtual void zero()
Set to zero vector.
Teuchos::RCP< const std::vector< Element > > getVector() const
std::vector< Element > vector
Defines the linear algebra or vector space interface.
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void plus(const ROL::Vector< Real > &x)
Compute , where .
Teuchos::RCP< ConDualStdVector< Real > > dual_vec_
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...
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
OptDualStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec)
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void scale(const Real alpha)
Compute where .
void scale(const Real alpha)
Compute where .
Teuchos::RCP< const std::vector< Element > > getVector() const
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.
Provides an interface to run optimization algorithms.
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Real norm() const
Returns where .
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
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.
Teuchos::RCP< std::vector< Element > > std_vec_
Real norm() const
Returns where .
ConDualStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec)
virtual std::vector< std::string > run(Vector< Real > &x, Objective< Real > &obj, bool print=false, std::ostream &outStream=std::cout)
Run algorithm on unconstrained problems (Type-U). This is the primary Type-U interface.
Teuchos::RCP< ConStdVector< Real > > dual_vec_
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Teuchos::RCP< std::vector< Element > > getVector()
Teuchos::RCP< const std::vector< Element > > getVector() const
void plus(const ROL::Vector< Real > &x)
Compute , where .
int main(int argc, char *argv[])
void plus(const ROL::Vector< Real > &x)
Compute , where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
Teuchos::RCP< std::vector< Element > > getVector()
Real norm() const
Returns where .
void scale(const Real alpha)
Compute where .
static const double ROL_EPSILON
Platform-dependent machine epsilon.