126int 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;
205 ROL::Ptr<V> sol = CreatePartitionedVector(x,l,zl,zu);
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);
226 RandomizeVector(*u,-1.0,1.0);
227 RandomizeVector(*v,-1.0,1.0);
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;