44 #define OPTIMIZATION_PROBLEM_REFACTOR
46 #include "Teuchos_GlobalMPISession.hpp"
50 #include "ROL_NonlinearProgram.hpp"
57 #include "HS_ProblemFactory.hpp"
75 ROL::Ptr<const std::vector<Real> > xp =
78 outStream <<
"Standard Vector" << std::endl;
79 for(
size_t i=0; i<xp->size(); ++i ) {
80 outStream << (*xp)[i] << std::endl;
83 catch(
const std::bad_cast& e ) {
84 outStream <<
"Partitioned Vector" << std::endl;
89 const PV &xpv =
dynamic_cast<const PV&
>(x);
91 for( size_type i=0; i<xpv.numVectors(); ++i ) {
92 outStream <<
"--------------------" << std::endl;
95 outStream <<
"--------------------" << std::endl;
102 std::ostream &outStream ) {
106 for( uint i=0; i<dim; ++i ) {
107 for( uint j=0; j<dim; ++j ) {
108 outStream << std::setw(6) << A[j]->dot(*(I[i]));
110 outStream << std::endl;
126 int main(
int argc,
char *argv[]) {
130 typedef ROL::ParameterList PL;
138 typedef ROL::NonlinearProgram<RealT> NLP;
149 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
151 int iprint = argc - 1;
152 ROL::Ptr<std::ostream> outStream;
155 outStream = ROL::makePtrFromRef(std::cout);
157 outStream = ROL::makePtrFromRef(bhs);
165 RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
169 PL &iplist = parlist.sublist(
"Step").sublist(
"Primal Dual Interior Point");
170 PL &lblist = iplist.sublist(
"Barrier Objective");
172 lblist.set(
"Use Linear Damping",
true);
173 lblist.set(
"Linear Damping Coefficient",1.e-4);
174 lblist.set(
"Initial Barrier Parameter", mu);
176 PL &krylist = parlist.sublist(
"General").sublist(
"Krylov");
178 krylist.set(
"Absolute Tolerance", 1.e-6);
179 krylist.set(
"Relative Tolerance", 1.e-6);
180 krylist.set(
"Iteration Limit", 50);
183 krylist.set(
"Type",
"Conjugate Gradients");
184 ROL::Ptr<KRYLOV> cg = ROL::KrylovFactory<RealT>(parlist);
185 HS::ProblemFactory<RealT> problemFactory;
189 ROL::Ptr<NLP> nlp = problemFactory.getProblem(16);
190 ROL::Ptr<OPT> opt = nlp->getOptimizationProblem();
192 ROL::Ptr<V> x = opt->getSolutionVector();
193 ROL::Ptr<V> l = opt->getMultiplierVector();
194 ROL::Ptr<V> zl = x->clone(); zl->zero();
195 ROL::Ptr<V> zu = x->clone(); zu->zero();
197 ROL::Ptr<V> scratch = x->clone();
199 ROL::Ptr<PV> x_pv = ROL::dynamicPtrCast<PV>(x);
203 (*ROL::dynamicPtrCast<ROL::StdVector<RealT> >(x_pv->get(1))->getVector())[0] = 1.0;
207 std::vector< ROL::Ptr<V> > I;
208 std::vector< ROL::Ptr<V> > J;
210 for(
int k=0; k<sol->dimension(); ++k ) {
211 I.push_back(sol->basis(k));
212 J.push_back(sol->clone());
215 ROL::Ptr<V> u = sol->clone();
216 ROL::Ptr<V> v = sol->clone();
218 ROL::Ptr<V> rhs = sol->clone();
219 ROL::Ptr<V> symrhs = sol->clone();
221 ROL::Ptr<V> gmres_sol = sol->clone(); gmres_sol->set(*sol);
222 ROL::Ptr<V> cg_sol = sol->clone(); cg_sol->set(*sol);
229 ROL::Ptr<OBJ> obj = opt->getObjective();
230 ROL::Ptr<CON> con = opt->getConstraint();
231 ROL::Ptr<BND> bnd = opt->getBoundConstraint();
233 PENALTY penalty(obj,bnd,parlist);
235 ROL::Ptr<const V> maskL = penalty.getLowerMask();
236 ROL::Ptr<const V> maskU = penalty.getUpperMask();
249 ROL::Ptr<CON> res = ROL::makePtr<RESIDUAL>(obj,con,bnd,*sol,maskL,maskU,scratch,mu,
false);
250 ROL::Ptr<LOP> lop = ROL::makePtr<LOPEC>( sol, res );
253 res->value(*rhs,*sol,tol);
256 krylist.set(
"Type",
"GMRES");
257 ROL::Ptr<KRYLOV> gmres = ROL::KrylovFactory<RealT>(parlist);
259 for(
int k=0; k<sol->dimension(); ++k ) {
260 res->applyJacobian(*(J[k]),*(I[k]),*sol,tol);
263 *outStream <<
"Nonsymmetric Jacobian" << std::endl;
267 gmres->run( *gmres_sol, *lop, *rhs, identity, gmres_iter, gmres_flag );
269 errorFlag += gmres_flag;
271 *outStream <<
"GMRES terminated after " << gmres_iter <<
" iterations "
272 <<
"with exit flag " << gmres_flag << std::endl;
282 ROL::Ptr<V> jv = v->clone();
283 ROL::Ptr<V> ju = u->clone();
285 iplist.set(
"Symmetrize Primal Dual System",
true);
286 ROL::Ptr<CON> symres = ROL::makePtr<RESIDUAL>(obj,con,bnd,*sol,maskL,maskU,scratch,mu,
true);
287 ROL::Ptr<LOP> symlop = ROL::makePtr<LOPEC>( sol, res );
288 symres->value(*symrhs,*sol,tol);
290 symres->applyJacobian(*jv,*v,*sol,tol);
291 symres->applyJacobian(*ju,*u,*sol,tol);
292 *outStream <<
"Symmetry check |u.dot(jv)-v.dot(ju)| = "
293 << std::abs(u->dot(*jv)-v->dot(*ju)) << std::endl;
295 cg->run( *cg_sol, *symlop, *symrhs, identity, cg_iter, cg_flag );
297 *outStream <<
"CG terminated after " << cg_iter <<
" iterations "
298 <<
"with exit flag " << cg_flag << std::endl;
300 *outStream <<
"Check that GMRES solution also is a solution to the symmetrized system"
303 symres->applyJacobian(*ju,*gmres_sol,*sol,tol);
304 ju->axpy(-1.0,*symrhs);
305 RealT mismatch = ju->norm();
309 *outStream <<
"||J_sym*sol_nonsym-rhs_sym|| = " << mismatch << std::endl;
312 catch (std::logic_error err) {
313 *outStream << err.what() << std::endl;
318 std::cout <<
"End Result: TEST FAILED" << std::endl;
320 std::cout <<
"End Result: TEST PASSED" << std::endl;
Provides the interface to evaluate objective functions.
PartitionedVector< Real > PV
typename PV< Real >::size_type size_type
Defines the linear algebra of vector space on a generic partitioned vector.
ROL::Ptr< Vector< Real > > CreatePartitionedVector(const ROL::Ptr< Vector< Real >> &a)
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,upper].
Symmetrized form of the KKT operator for the Type-EB problem with equality and bound multipliers...
Defines the linear algebra or vector space interface.
void printVector(const ROL::Vector< Real > &x, std::ostream &outStream)
A simple wrapper which allows application of constraint Jacobians through the LinearOperator interfac...
Provides the interface to evaluate the Interior Pointy log barrier penalty function with upper and lo...
void printMatrix(const std::vector< ROL::Ptr< ROL::Vector< Real > > > &A, const std::vector< ROL::Ptr< ROL::Vector< Real > > > &I, std::ostream &outStream)
void apply(ROL::Vector< Real > &Hv, const ROL::Vector< Real > &v, Real &tol) const
Apply linear operator.
Provides definitions for Krylov solvers.
Provides the interface to apply a linear operator.
Provides the interface to apply upper and lower bound constraints.
basic_nullstream< char, char_traits< char >> nullstream
int main(int argc, char *argv[])
virtual void set(const Vector &x)
Set where .
Defines the general constraint operator interface.