55 #include "Teuchos_GlobalMPISession.hpp" 57 template<
typename Real>
60 void value(std::vector<Real> &c,
const std::vector<Real> &x, Real &tol) {
61 c[0] = x[0]+x[1]-
static_cast<Real
>(1);
63 void applyJacobian(std::vector<Real> &jv,
const std::vector<Real> &v,
const std::vector<Real> &x, Real &tol) {
66 void applyAdjointJacobian(std::vector<Real> &ajv,
const std::vector<Real> &v,
const std::vector<Real> &x, Real &tol) {
74 int main(
int argc,
char *argv[]) {
76 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
80 int iprint = argc - 1;
81 ROL::Ptr<std::ostream> outStream;
84 outStream = ROL::makePtrFromRef(std::cout);
86 outStream = ROL::makePtrFromRef(bhs);
92 RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
94 ROL::Ptr<con2d<RealT>> con = ROL::makePtr<con2d<RealT>>();
96 ROL::ParameterList list;
97 list.sublist(
"General").set(
"Output Level",2);
98 list.sublist(
"General").sublist(
"Polyhedral Projection").set(
"Type",
"Dai-Fletcher");
103 ROL::Ptr<std::vector<RealT>> yptr = ROL::makePtr<std::vector<RealT>>(2);
104 (*yptr)[0] =
static_cast<RealT>(10)*(static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX)-half);
105 (*yptr)[1] =
static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
109 ROL::Ptr<std::vector<RealT>> xptr = ROL::makePtr<std::vector<RealT>>(2);
110 (*xptr)[0] = (*yptr)[0];
111 (*xptr)[1] = (*yptr)[1];
114 ROL::Ptr<std::vector<RealT>> Pxptr = ROL::makePtr<std::vector<RealT>>(2,0.0);
115 (*Pxptr)[0] = (*yptr)[0];
116 (*Pxptr)[1] = (*yptr)[1];
119 ROL::Ptr<ROL::Vector<RealT>> l0 = x.
clone(); l0->setScalar(static_cast<RealT>(0));
120 ROL::Ptr<ROL::Vector<RealT>> u0 = x.
clone(); u0->setScalar(static_cast<RealT>(1));
121 ROL::Ptr<ROL::Bounds<RealT>> bnd0 = ROL::makePtr<ROL::Bounds<RealT>>(l0,u0);
123 ROL::Ptr<ROL::PolyhedralProjection<RealT>> pp0 = ROL::PolyhedralProjectionFactory<RealT>(x,x.
dual(),bnd0,con,r,r.dual(),list);
124 pp0->project(Px,*outStream);
126 ROL::Ptr<std::vector<RealT>> x0ptr = ROL::makePtr<std::vector<RealT>>(2);
127 RealT k0 = std::max(
zero,std::min(one,half*(one+(*yptr)[0]-(*yptr)[1])));
129 (*x0ptr)[1] = one-k0;
134 *outStream << std::setprecision(6) << std::scientific << std::endl;
135 *outStream <<
" x[0] = " << (*xptr)[0] <<
" x[1] = " << (*xptr)[1] << std::endl;
136 *outStream <<
" Px[0] = " << (*Pxptr)[0] <<
" Px[1] = " << (*Pxptr)[1] << std::endl;
137 *outStream <<
" x*[0] = " << (*x0ptr)[0] <<
" x*[1] = " << (*x0ptr)[1] << std::endl;
139 e0.
set(x0); e0.
axpy(static_cast<RealT>(-1),Px);
141 *outStream <<
" Error in Euclidean Projection: " << err << std::endl;
143 e0.
set(x); e0.
axpy(static_cast<RealT>(-1),x0);
144 *outStream <<
" ||x*-x||^2 = " << e0.
norm() << std::endl;
146 e0.
set(x); e0.
axpy(static_cast<RealT>(-1),Px);
147 *outStream <<
" ||Px-x||^2 = " << e0.
norm() << std::endl << std::endl;
149 errorFlag += (err > tol);
151 ROL::Ptr<std::vector<RealT>> dptr = ROL::makePtr<std::vector<RealT>>(2);
152 (*dptr)[0] =
static_cast<RealT>(1)+static_cast<RealT>(2)*
static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
153 (*dptr)[1] =
static_cast<RealT>(1)+static_cast<RealT>(5)*
static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
155 ROL::Ptr<std::vector<RealT>> x1ptr = ROL::makePtr<std::vector<RealT>>(2);
156 RealT k1 = std::max(
zero,std::min(one,((*dptr)[1]*(one-(*yptr)[1])+(*dptr)[0]*(*yptr)[0])/((*dptr)[0]+(*dptr)[1])));
158 (*x1ptr)[1] = one-k1;
161 ROL::Ptr<std::vector<RealT>> zptr = ROL::makePtr<std::vector<RealT>>(2);
162 (*zptr)[0] = (*yptr)[0];
163 (*zptr)[1] = (*yptr)[1];
166 ROL::Ptr<std::vector<RealT>> Pzptr = ROL::makePtr<std::vector<RealT>>(2,0.0);
167 (*Pzptr)[0] = (*yptr)[0];
168 (*Pzptr)[1] = (*yptr)[1];
171 ROL::Ptr<ROL::Vector<RealT>> l1 = z.
clone(); l1->setScalar(static_cast<RealT>(0));
172 ROL::Ptr<ROL::Vector<RealT>> u1 = z.
clone(); u1->setScalar(static_cast<RealT>(1));
173 ROL::Ptr<ROL::Bounds<RealT>> bnd1 = ROL::makePtr<ROL::Bounds<RealT>>(l1,u1);
175 ROL::Ptr<ROL::PolyhedralProjection<RealT>> pp1 = ROL::PolyhedralProjectionFactory<RealT>(z,z.
dual(),bnd1,con,r,r.dual(),list);
176 pp1->project(Pz,*outStream);
178 ROL::Ptr<std::vector<RealT>> e1ptr = ROL::makePtr<std::vector<RealT>>(2);
181 *outStream << std::endl;
182 *outStream <<
" x[0] = " << (*zptr)[0] <<
" x[1] = " << (*zptr)[1] << std::endl;
183 *outStream <<
" Px[0] = " << (*Pzptr)[0] <<
" Px[1] = " << (*Pzptr)[1] << std::endl;
184 *outStream <<
" x*[0] = " << (*x1ptr)[0] <<
" x*[1] = " << (*x1ptr)[1] << std::endl;
186 e1.set(x1); e1.axpy(static_cast<RealT>(-1),Pz);
188 *outStream <<
" Error in Scaled Projection: " << err << std::endl;
190 e1.set(z); e1.axpy(static_cast<RealT>(-1),x1);
191 *outStream <<
" ||x*-x||^2 = " << e1.norm() << std::endl;
193 e1.set(z); e1.axpy(static_cast<RealT>(-1),Pz);
194 *outStream <<
" ||Px-x||^2 = " << e1.norm() << std::endl << std::endl;
196 errorFlag += (err > tol);
199 catch (std::logic_error& err) {
200 *outStream << err.what() <<
"\n";
205 std::cout <<
"End Result: TEST FAILED\n";
207 std::cout <<
"End Result: TEST PASSED\n";
void axpy(const Real alpha, const Vector< Real > &x)
Compute where .
basic_nullstream< char, char_traits< char > > nullstream
int main(int argc, char *argv[])
Defines the equality constraint operator interface for StdVectors.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
void applyJacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
Provides the std::vector implementation of the ROL::Vector interface that handles scalings in the inn...
void set(const Vector< Real > &x)
Set where .
virtual Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
const Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void applyAdjointJacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
void value(std::vector< Real > &c, const std::vector< Real > &x, Real &tol)
Real norm() const
Returns where .