ROL
step/interiorpoint/test_03.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#define OPTIMIZATION_PROBLEM_REFACTOR
45
46#include "Teuchos_GlobalMPISession.hpp"
47
48#include "ROL_RandomVector.hpp"
49#include "ROL_StdVector.hpp"
50#include "ROL_NonlinearProgram.hpp"
55#include "ROL_KrylovFactory.hpp"
56
57#include "HS_ProblemFactory.hpp"
58
59#include <iomanip>
60
71template<class Real>
72void printVector( const ROL::Vector<Real> &x, std::ostream &outStream ) {
73
74 try {
75 ROL::Ptr<const std::vector<Real> > xp =
76 dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
77
78 outStream << "Standard Vector" << std::endl;
79 for( size_t i=0; i<xp->size(); ++i ) {
80 outStream << (*xp)[i] << std::endl;
81 }
82 }
83 catch( const std::bad_cast& e ) {
84 outStream << "Partitioned Vector" << std::endl;
85
87 typedef typename PV::size_type size_type;
88
89 const PV &xpv = dynamic_cast<const PV&>(x);
90
91 for( size_type i=0; i<xpv.numVectors(); ++i ) {
92 outStream << "--------------------" << std::endl;
93 printVector( *(xpv.get(i)), outStream );
94 }
95 outStream << "--------------------" << std::endl;
96 }
97}
98
99template<class Real>
100void printMatrix( const std::vector<ROL::Ptr<ROL::Vector<Real> > > &A,
101 const std::vector<ROL::Ptr<ROL::Vector<Real> > > &I,
102 std::ostream &outStream ) {
103 typedef typename std::vector<Real>::size_type uint;
104 uint dim = A.size();
105
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]));
109 }
110 outStream << std::endl;
111 }
112}
113
114
115template<class Real>
117public:
118 void apply( ROL::Vector<Real> &Hv, const ROL::Vector<Real> &v, Real &tol ) const {
119 Hv.set(v);
120 }
121}; // IdentityOperator
122
123
124typedef double RealT;
125
126int main(int argc, char *argv[]) {
127
128// typedef std::vector<RealT> vector;
129
130 typedef ROL::ParameterList PL;
131
132 typedef ROL::Vector<RealT> V;
134 typedef ROL::Objective<RealT> OBJ;
135 typedef ROL::Constraint<RealT> CON;
136 typedef ROL::BoundConstraint<RealT> BND;
138 typedef ROL::NonlinearProgram<RealT> NLP;
139 typedef ROL::LinearOperator<RealT> LOP;
141 typedef ROL::Krylov<RealT> KRYLOV;
142
143
144 typedef ROL::InteriorPointPenalty<RealT> PENALTY;
146
147
148
149 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
150
151 int iprint = argc - 1;
152 ROL::Ptr<std::ostream> outStream;
153 ROL::nullstream bhs;
154 if( iprint > 0 )
155 outStream = ROL::makePtrFromRef(std::cout);
156 else
157 outStream = ROL::makePtrFromRef(bhs);
158
159 int errorFlag = 0;
160
161 try {
162
163 RealT mu = 0.1;
164
165 RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
166
167 PL parlist;
168
169 PL &iplist = parlist.sublist("Step").sublist("Primal Dual Interior Point");
170 PL &lblist = iplist.sublist("Barrier Objective");
171
172 lblist.set("Use Linear Damping", true); // Not used in this test
173 lblist.set("Linear Damping Coefficient",1.e-4); // Not used in this test
174 lblist.set("Initial Barrier Parameter", mu);
175
176 PL &krylist = parlist.sublist("General").sublist("Krylov");
177
178 krylist.set("Absolute Tolerance", 1.e-6);
179 krylist.set("Relative Tolerance", 1.e-6);
180 krylist.set("Iteration Limit", 50);
181
182 // Create a Conjugate Gradients solver
183 krylist.set("Type","Conjugate Gradients");
184 ROL::Ptr<KRYLOV> cg = ROL::KrylovFactory<RealT>(parlist);
185 HS::ProblemFactory<RealT> problemFactory;
186
187 // Choose an example problem with inequality constraints and
188 // a mixture of finite and infinite bounds
189 ROL::Ptr<NLP> nlp = problemFactory.getProblem(16);
190 ROL::Ptr<OPT> opt = nlp->getOptimizationProblem();
191
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();
196
197 ROL::Ptr<V> scratch = x->clone();
198
199 ROL::Ptr<PV> x_pv = ROL::dynamicPtrCast<PV>(x);
200 // New slack variable initialization does not guarantee strict feasibility.
201 // This ensures that the slack variables are the same as the previous
202 // implementation.
203 (*ROL::dynamicPtrCast<ROL::StdVector<RealT> >(x_pv->get(1))->getVector())[0] = 1.0;
204
205 ROL::Ptr<V> sol = CreatePartitionedVector(x,l,zl,zu);
206
207 std::vector< ROL::Ptr<V> > I;
208 std::vector< ROL::Ptr<V> > J;
209
210 for( int k=0; k<sol->dimension(); ++k ) {
211 I.push_back(sol->basis(k));
212 J.push_back(sol->clone());
213 }
214
215 ROL::Ptr<V> u = sol->clone();
216 ROL::Ptr<V> v = sol->clone();
217
218 ROL::Ptr<V> rhs = sol->clone();
219 ROL::Ptr<V> symrhs = sol->clone();
220
221 ROL::Ptr<V> gmres_sol = sol->clone(); gmres_sol->set(*sol);
222 ROL::Ptr<V> cg_sol = sol->clone(); cg_sol->set(*sol);
223
225
226 RandomizeVector(*u,-1.0,1.0);
227 RandomizeVector(*v,-1.0,1.0);
228
229 ROL::Ptr<OBJ> obj = opt->getObjective();
230 ROL::Ptr<CON> con = opt->getConstraint();
231 ROL::Ptr<BND> bnd = opt->getBoundConstraint();
232
233 PENALTY penalty(obj,bnd,parlist);
234
235 ROL::Ptr<const V> maskL = penalty.getLowerMask();
236 ROL::Ptr<const V> maskU = penalty.getUpperMask();
237
238 zl->set(*maskL);
239 zu->set(*maskU);
240
241 /********************************************************************************/
242 /* Nonsymmetric representation test */
243 /********************************************************************************/
244
245 int gmres_iter = 0;
246 int gmres_flag = 0;
247
248 // Form the residual's Jacobian operator
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 );
251
252 // Evaluate the right-hand-side
253 res->value(*rhs,*sol,tol);
254
255 // Create a GMRES solver
256 krylist.set("Type","GMRES");
257 ROL::Ptr<KRYLOV> gmres = ROL::KrylovFactory<RealT>(parlist);
258
259 for( int k=0; k<sol->dimension(); ++k ) {
260 res->applyJacobian(*(J[k]),*(I[k]),*sol,tol);
261 }
262
263 *outStream << "Nonsymmetric Jacobian" << std::endl;
264 printMatrix(J,I,*outStream);
265
266 // Solve the system
267 gmres->run( *gmres_sol, *lop, *rhs, identity, gmres_iter, gmres_flag );
268
269 errorFlag += gmres_flag;
270
271 *outStream << "GMRES terminated after " << gmres_iter << " iterations "
272 << "with exit flag " << gmres_flag << std::endl;
273
274
275 /********************************************************************************/
276 /* Symmetric representation test */
277 /********************************************************************************/
278
279 int cg_iter = 0;
280 int cg_flag = 0;
281
282 ROL::Ptr<V> jv = v->clone();
283 ROL::Ptr<V> ju = u->clone();
284
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);
289
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;
294
295 cg->run( *cg_sol, *symlop, *symrhs, identity, cg_iter, cg_flag );
296
297 *outStream << "CG terminated after " << cg_iter << " iterations "
298 << "with exit flag " << cg_flag << std::endl;
299
300 *outStream << "Check that GMRES solution also is a solution to the symmetrized system"
301 << std::endl;
302
303 symres->applyJacobian(*ju,*gmres_sol,*sol,tol);
304 ju->axpy(-1.0,*symrhs);
305 RealT mismatch = ju->norm();
306 if(mismatch>tol) {
307 errorFlag++;
308 }
309 *outStream << "||J_sym*sol_nonsym-rhs_sym|| = " << mismatch << std::endl;
310
311 }
312 catch (std::logic_error& err) {
313 *outStream << err.what() << std::endl;
314 errorFlag = -1000;
315 }
316
317 if (errorFlag != 0)
318 std::cout << "End Result: TEST FAILED" << std::endl;
319 else
320 std::cout << "End Result: TEST PASSED" << std::endl;
321
322 return 0;
323}
324
Vector< Real > V
PartitionedVector< Real > PV
typename PV< Real >::size_type size_type
void apply(ROL::Vector< Real > &Hv, const ROL::Vector< Real > &v, Real &tol) const
Apply linear operator.
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 definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:58
A simple wrapper which allows application of constraint Jacobians through the LinearOperator interfac...
Provides the interface to apply a linear operator.
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.
Definition: ROL_Vector.hpp:84
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
int main(int argc, char *argv[])
void printVector(const ROL::Vector< Real > &x, std::ostream &outStream)
void printMatrix(const std::vector< ROL::Ptr< ROL::Vector< Real > > > &A, const std::vector< ROL::Ptr< ROL::Vector< Real > > > &I, std::ostream &outStream)
constexpr auto dim