44#define OPTIMIZATION_PROBLEM_REFACTOR
46#include "Teuchos_GlobalMPISession.hpp"
50#include "ROL_NonlinearProgram.hpp"
55#include "HS_Problem_041.hpp"
73 ROL::Ptr<const std::vector<Real> > xp =
76 outStream <<
"Standard Vector" << std::endl;
77 for(
size_t i=0; i<xp->size(); ++i ) {
78 outStream << (*xp)[i] << std::endl;
81 catch(
const std::bad_cast& e ) {
82 outStream <<
"Partitioned Vector" << std::endl;
87 const PV &xpv =
dynamic_cast<const PV&
>(x);
89 for(
size_type i=0; i<xpv.numVectors(); ++i ) {
90 outStream <<
"--------------------" << std::endl;
93 outStream <<
"--------------------" << std::endl;
103 typedef std::vector<Real> vector;
107 typedef typename PV::size_type
size_type;
111 using ROL::dynamicPtrCast;
118 const PV &sol_pv =
dynamic_cast<const PV&
>(sol);
119 const vector &x = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(OPT))->getVector());
120 const vector &l = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(EQUAL))->getVector());
121 const vector &zl = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(LOWER))->getVector());
122 const vector &zu = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(UPPER))->getVector());
124 PV &c_pv =
dynamic_cast<PV&
>(c);
125 vector &cx = *(ROL::dynamicPtrCast<SV>(c_pv.get(OPT))->getVector());
126 vector &cl = *(ROL::dynamicPtrCast<SV>(c_pv.get(EQUAL))->getVector());
127 vector &czl = *(ROL::dynamicPtrCast<SV>(c_pv.get(LOWER))->getVector());
128 vector &czu = *(ROL::dynamicPtrCast<SV>(c_pv.get(UPPER))->getVector());
130 cx[0] = -x[1]*x[2] + l[0] - zl[0] + zu[0];
131 cx[1] = -x[0]*x[1] + 2*l[0] - zl[1] + zu[1];
132 cx[2] = -x[0]*x[1] + 2*l[0] - zl[2] + zu[2];
133 cx[3] = - l[0] - zl[3] + zu[3];
135 cl[0] = x[0] + 2*x[1] + 2*x[2] - x[3];
137 czl[0] = x[0]*zl[0] - mu;
138 czl[1] = x[1]*zl[1] - mu;
139 czl[2] = x[2]*zl[2] - mu;
140 czl[3] = x[3]*zl[3] - mu;
142 czu[0] = (1.0-x[0])*zu[0] - mu;
143 czu[1] = (1.0-x[1])*zu[1] - mu;
144 czu[2] = (1.0-x[2])*zl[2] - mu;
145 czu[3] = (2.0-x[3])*zl[3] - mu;
154 typedef std::vector<Real> vector;
158 typedef typename PV::size_type
size_type;
162 using ROL::dynamicPtrCast;
169 const PV &sol_pv =
dynamic_cast<const PV&
>(sol);
170 const vector &x = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(OPT))->getVector());
172 const vector &zl = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(LOWER))->getVector());
173 const vector &zu = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(UPPER))->getVector());
175 const PV &v_pv =
dynamic_cast<const PV&
>(v);
176 const vector &vx = *(ROL::dynamicPtrCast<const SV>(v_pv.get(OPT))->getVector());
177 const vector &vl = *(ROL::dynamicPtrCast<const SV>(v_pv.get(EQUAL))->getVector());
178 const vector &vzl = *(ROL::dynamicPtrCast<const SV>(v_pv.get(LOWER))->getVector());
179 const vector &vzu = *(ROL::dynamicPtrCast<const SV>(v_pv.get(UPPER))->getVector());
181 PV &jv_pv =
dynamic_cast<PV&
>(jv);
182 vector &jvx = *(ROL::dynamicPtrCast<SV>(jv_pv.get(OPT))->getVector());
183 vector &jvl = *(ROL::dynamicPtrCast<SV>(jv_pv.get(EQUAL))->getVector());
184 vector &jvzl = *(ROL::dynamicPtrCast<SV>(jv_pv.get(LOWER))->getVector());
185 vector &jvzu = *(ROL::dynamicPtrCast<SV>(jv_pv.get(UPPER))->getVector());
187 jvx[0] = -x[1]*vx[2] - x[2]*vx[1] + vl[0] - vzl[0] + vzu[0];
188 jvx[1] = -x[0]*vx[2] - x[2]*vx[0] + 2*vl[0] - vzl[1] + vzu[1];
189 jvx[2] = -x[0]*vx[1] - x[1]*vx[0] + 2*vl[0] - vzl[2] + vzu[2];
190 jvx[3] = - vl[0] - vzl[3] + vzu[3];
192 jvl[0] = vx[0] + 2*vx[1] + 2*vx[2] - vx[3];
194 jvzl[0] = zl[0]*vx[0] + vzl[0]*x[0];
195 jvzl[1] = zl[1]*vx[1] + vzl[1]*x[1];
196 jvzl[2] = zl[2]*vx[2] + vzl[2]*x[2];
197 jvzl[3] = zl[3]*vx[3] + vzl[3]*x[3];
199 jvzu[0] = -zu[0]*vx[0] + vzu[0]*(1.0-x[0]);
200 jvzu[1] = -zu[1]*vx[1] + vzu[1]*(1.0-x[1]);
201 jvzu[2] = -zu[2]*vx[2] + vzu[2]*(1.0-x[2]);
202 jvzu[3] = -zu[3]*vx[3] + vzu[3]*(2.0-x[3]);
209int main(
int argc,
char *argv[]) {
213 typedef ROL::ParameterList PL;
221 typedef ROL::NonlinearProgram<RealT> NLP;
228 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
230 int iprint = argc - 1;
231 ROL::Ptr<std::ostream> outStream;
234 outStream = ROL::makePtrFromRef(std::cout);
236 outStream = ROL::makePtrFromRef(bhs);
244 RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
248 PL &iplist = parlist.sublist(
"Step").sublist(
"Primal Dual Interior Point");
249 PL &lblist = iplist.sublist(
"Barrier Objective");
251 iplist.set(
"Symmetrize Primal Dual System",
false);
253 lblist.set(
"Use Linear Damping",
true);
254 lblist.set(
"Linear Damping Coefficient",1.e-4);
255 lblist.set(
"Initial Barrier Parameter", mu);
261 ROL::Ptr<NLP> nlp = ROL::makePtr<HS::Problem_041<RealT>>();
262 ROL::Ptr<OPT> opt = nlp->getOptimizationProblem();
264 ROL::Ptr<V> x = opt->getSolutionVector();
265 ROL::Ptr<V> l = opt->getMultiplierVector();
266 ROL::Ptr<V> zl = x->clone();
267 ROL::Ptr<V> zu = x->clone();
269 ROL::Ptr<V> scratch = x->clone();
271 ROL::Ptr<PV> x_pv = ROL::dynamicPtrCast<PV>(x);
273 ROL::Ptr<V> sol = CreatePartitionedVector(x,l,zl,zu);
275 ROL::Ptr<V> c = sol->clone();
276 ROL::Ptr<V> v = sol->clone();
277 ROL::Ptr<V> jv = sol->clone();
279 ROL::Ptr<V> c_exact = c->clone();
280 ROL::Ptr<V> jv_exact = jv->clone();
286 ROL::Ptr<OBJ> obj = opt->getObjective();
287 ROL::Ptr<CON> con = opt->getConstraint();
288 ROL::Ptr<BND> bnd = opt->getBoundConstraint();
290 PENALTY penalty(obj,bnd,parlist);
292 ROL::Ptr<const V> maskL = penalty.getLowerMask();
293 ROL::Ptr<const V> maskU = penalty.getUpperMask();
298 ROL::Ptr<CON> res = ROL::makePtr<RESIDUAL>(obj,con,bnd,*sol,maskL,maskU,scratch,mu,
false);
301 *outStream <<
"\n[x|lambda|zl|zu]" << std::endl;
304 res->value(*c,*sol,tol);
305 value(*c_exact,*sol,mu);
307 *outStream <<
"\n[cx|clambda|czl|czu]" << std::endl;
311 c->axpy(-1.0,*c_exact);
313 RealT cerror = c->norm();
319 *outStream <<
"\n\n||c-c_exact|| = " << cerror << std::endl;
321 res->applyJacobian(*jv,*v,*sol,tol);
324 *outStream <<
"\n[vx|vlambda|vzl|vzu]" << std::endl;
328 *outStream <<
"\n[jvx|jvlambda|jvzl|jvzu]" << std::endl;
332 jv->axpy(-1.0,*jv_exact);
334 RealT jverror = jv->norm();
336 if( jverror > tol ) {
340 *outStream <<
"\n\n||jv-jv_exact|| = " << jverror << std::endl;
342 *outStream <<
"Residual Jacobian Finite Difference Check" << std::endl;
343 res->checkApplyJacobian(*sol,*v,*v);
348 catch (std::logic_error& err) {
349 *outStream << err.what() << std::endl;
354 std::cout <<
"End Result: TEST FAILED" << std::endl;
356 std::cout <<
"End Result: TEST PASSED" << std::endl;
PartitionedVector< Real > PV
typename PV< Real >::size_type size_type
Provides the interface to apply upper and lower bound constraints.
Defines the general constraint operator interface.
Provides the interface to evaluate the Interior Pointy log barrier penalty function with upper and lo...
Provides the interface to evaluate objective functions.
Defines the linear algebra of vector space on a generic partitioned vector.
Symmetrized form of the KKT operator for the Type-EB problem with equality and bound multipliers.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
void RandomizeVector(Vector< Real > &x, const Real &lower=0.0, const Real &upper=1.0)
Fill a ROL::Vector with uniformly-distributed random numbers in the interval [lower,...
int main(int argc, char *argv[])
void applyJacobian(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &sol)
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &sol, const Real &mu)
void printVector(const ROL::Vector< Real > &x, std::ostream &outStream)