ROL
dual-spaces/rosenbrock-1/example_01.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
48#define USE_HESSVEC 1
49
50#include "ROL_Rosenbrock.hpp"
51#include "ROL_Solver.hpp"
52#include "ROL_Stream.hpp"
53#include "Teuchos_GlobalMPISession.hpp"
54#include <iostream>
55
56typedef double RealT;
57
58
59/*** Declare two vector spaces. ***/
60
61// Forward declarations:
62
63template <class Real, class Element=Real>
64class OptStdVector; // Optimization space.
65
66template <class Real, class Element=Real>
67class OptDualStdVector; // Dual optimization space.
68
69
70// Vector space definitions:
71
72// Optimization space.
73template <class Real, class Element>
74class OptStdVector : public ROL::Vector<Real> {
75
76 typedef std::vector<Element> vector;
78
79 typedef typename vector::size_type uint;
80
81private:
82ROL::Ptr<std::vector<Element> > std_vec_;
83mutable ROL::Ptr<OptDualStdVector<Real> > dual_vec_;
84
85public:
86
87OptStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}
88
89void plus( const ROL::Vector<Real> &x ) {
90
91 ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector&>(x).getVector();
92 uint dimension = std_vec_->size();
93 for (uint i=0; i<dimension; i++) {
94 (*std_vec_)[i] += (*xvalptr)[i];
95 }
96}
97
98void scale( const Real alpha ) {
99 uint dimension = std_vec_->size();
100 for (uint i=0; i<dimension; i++) {
101 (*std_vec_)[i] *= alpha;
102 }
103}
104
105Real dot( const ROL::Vector<Real> &x ) const {
106 Real val = 0;
107
108 ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector&>(x).getVector();
109 uint dimension = std_vec_->size();
110 for (uint i=0; i<dimension; i++) {
111 val += (*std_vec_)[i]*(*xvalptr)[i];
112 }
113 return val;
114}
115
116Real norm() const {
117 Real val = 0;
118 val = std::sqrt( dot(*this) );
119 return val;
120}
121
122ROL::Ptr<ROL::Vector<Real> > clone() const {
123 return ROL::makePtr<OptStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()) );
124}
125
126ROL::Ptr<const std::vector<Element> > getVector() const {
127 return std_vec_;
128}
129
130ROL::Ptr<std::vector<Element> > getVector() {
131 return std_vec_;
132}
133
134ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
135
136 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
137 (*e_ptr)[i] = 1.0;
138
139 ROL::Ptr<V> e = ROL::makePtr<OptStdVector>( e_ptr );
140
141 return e;
142}
143
144int dimension() const {return static_cast<int>(std_vec_->size());}
145
146const ROL::Vector<Real> & dual() const {
147 dual_vec_ = ROL::makePtr<OptDualStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
148 return *dual_vec_;
149}
150
151Real apply( const ROL::Vector<Real> &x ) const {
152 Real val = 0;
153
154 ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector<Real,Element>&>(x).getVector();
155 uint dimension = std_vec_->size();
156 for (uint i=0; i<dimension; i++) {
157 val += (*std_vec_)[i]*(*xvalptr)[i];
158 }
159 return val;
160}
161
162}; // class OptStdVector
163
164
165// Dual optimization space.
166template <class Real, class Element>
167class OptDualStdVector : public ROL::Vector<Real> {
168
169 typedef std::vector<Element> vector;
171
172 typedef typename vector::size_type uint;
173
174private:
175ROL::Ptr<std::vector<Element> > std_vec_;
176mutable ROL::Ptr<OptStdVector<Real> > dual_vec_;
177
178public:
179
180OptDualStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}
181
182void plus( const ROL::Vector<Real> &x ) {
183
184 ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector&>(x).getVector();
185
186 uint dimension = std_vec_->size();
187 for (uint i=0; i<dimension; i++) {
188 (*std_vec_)[i] += (*xvalptr)[i];
189 }
190}
191
192void scale( const Real alpha ) {
193 uint dimension = std_vec_->size();
194 for (uint i=0; i<dimension; i++) {
195 (*std_vec_)[i] *= alpha;
196 }
197}
198
199Real dot( const ROL::Vector<Real> &x ) const {
200 Real val = 0;
201
202 ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector&>(x).getVector();
203 uint dimension = std_vec_->size();
204 for (uint i=0; i<dimension; i++) {
205 val += (*std_vec_)[i]*(*xvalptr)[i];
206 }
207 return val;
208}
209
210Real norm() const {
211 Real val = 0;
212 val = std::sqrt( dot(*this) );
213 return val;
214}
215
216ROL::Ptr<ROL::Vector<Real> > clone() const {
217 return ROL::makePtr<OptDualStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()) );
218}
219
220ROL::Ptr<const std::vector<Element> > getVector() const {
221 return std_vec_;
222}
223
224ROL::Ptr<std::vector<Element> > getVector() {
225 return std_vec_;
226}
227
228ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
229
230 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
231 (*e_ptr)[i] = 1.0;
232
233 ROL::Ptr<V> e = ROL::makePtr<OptDualStdVector>( e_ptr );
234 return e;
235}
236
237int dimension() const {return std_vec_->size();}
238
239const ROL::Vector<Real> & dual() const {
240 dual_vec_ = ROL::makePtr<OptStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
241 return *dual_vec_;
242}
243
244Real apply( const ROL::Vector<Real> &x ) const {
245 Real val = 0;
246
247 ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector<Real,Element>&>(x).getVector();
248 uint dimension = std_vec_->size();
249 for (uint i=0; i<dimension; i++) {
250 val += (*std_vec_)[i]*(*xvalptr)[i];
251 }
252 return val;
253}
254
255}; // class OptDualStdVector
256
257
258/*** End of declaration of two vector spaces. ***/
259
260
261
262
263
264
265int main(int argc, char *argv[]) {
266
267 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
268
269 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
270 int iprint = argc - 1;
271 ROL::Ptr<std::ostream> outStream;
272 ROL::nullstream bhs; // outputs nothing
273 if (iprint > 0)
274 outStream = ROL::makePtrFromRef(std::cout);
275 else
276 outStream = ROL::makePtrFromRef(bhs);
277
278 int errorFlag = 0;
279
280 // *** Example body.
281
282 try {
283
284 // Define objective function.
285 ROL::Ptr<ROL::ZOO::Objective_Rosenbrock<RealT, OptStdVector<RealT>, OptDualStdVector<RealT>>> obj
286 = ROL::makePtr<ROL::ZOO::Objective_Rosenbrock<RealT, OptStdVector<RealT>, OptDualStdVector<RealT>>>();
287 int dim = 100; // Set problem dimension. Must be even.
288
289 // Iteration Vector
290 ROL::Ptr<std::vector<RealT>> x_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
291 ROL::Ptr<std::vector<RealT>> g_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
292 // Set Initial Guess
293 for (int i=0; i<dim/2; i++) {
294 (*x_ptr)[2*i] = -1.2;
295 (*x_ptr)[2*i+1] = 1.0;
296 (*g_ptr)[2*i] = 0;
297 (*g_ptr)[2*i+1] = 0;
298 }
299 ROL::Ptr<OptStdVector<RealT>> x = ROL::makePtr<OptStdVector<RealT>>(x_ptr); // Iteration Vector
300 ROL::Ptr<OptDualStdVector<RealT>> g = ROL::makePtr<OptDualStdVector<RealT>>(g_ptr); // zeroed gradient vector in dual space
301
302 // Check vector.
303 ROL::Ptr<std::vector<RealT> > aa_ptr = ROL::makePtr<std::vector<RealT>>(1, 1.0);
304 OptDualStdVector<RealT> av(aa_ptr);
305 ROL::Ptr<std::vector<RealT> > bb_ptr = ROL::makePtr<std::vector<RealT>>(1, 2.0);
306 OptDualStdVector<RealT> bv(bb_ptr);
307 ROL::Ptr<std::vector<RealT> > cc_ptr = ROL::makePtr<std::vector<RealT>>(1, 3.0);
308 OptDualStdVector<RealT> cv(cc_ptr);
309 std::vector<RealT> std_vec_err = av.checkVector(bv,cv,true,*outStream);
310
311 // Build optimization problem.
312 ROL::Ptr<ROL::Problem<RealT>> problem
313 = ROL::makePtr<ROL::Problem<RealT>>(obj,x,g);
314 problem->finalize(false,true,*outStream);
315
316 // Define algorithm.
317 ROL::ParameterList parlist;
318 std::string stepname = "Trust Region";
319 parlist.sublist("Step").set("Type",stepname);
320 //parlist.sublist("Step").sublist(stepname).sublist("Descent Method").set("Type","Quasi-Newton Method");
321 parlist.sublist("Step").sublist(stepname).set("Subproblem Solver", "Truncated CG");
322 parlist.sublist("General").set("Output Level",1);
323 parlist.sublist("General").sublist("Krylov").set("Relative Tolerance",1e-2);
324 parlist.sublist("General").sublist("Krylov").set("Iteration Limit",10);
325 parlist.sublist("General").sublist("Krylov").set("Absolute Tolerance",1e-4);
326 parlist.sublist("General").sublist("Secant").set("Type","Limited-Memory BFGS");
327 parlist.sublist("General").sublist("Secant").set("Use as Hessian",true);
328 parlist.sublist("Status Test").set("Gradient Tolerance",1.e-12);
329 parlist.sublist("Status Test").set("Step Tolerance",1.e-14);
330 parlist.sublist("Status Test").set("Iteration Limit",100);
331 ROL::Ptr<ROL::Solver<RealT>> solver
332 = ROL::makePtr<ROL::Solver<RealT>>(problem,parlist);
333
334 // Run Algorithm
335 solver->solve(*outStream);
336
337 // Get True Solution
338 ROL::Ptr<std::vector<RealT> > xtrue_ptr = ROL::makePtr<std::vector<RealT>>(dim, 1.0);
339 OptStdVector<RealT> xtrue(xtrue_ptr);
340
341 // Compute Errors
342 x->axpy(-1.0, xtrue);
343 RealT abserr = x->norm();
344 RealT relerr = abserr/xtrue.norm();
345 *outStream << std::scientific << "\n Absolute solution error: " << abserr;
346 *outStream << std::scientific << "\n Relative solution error: " << relerr;
347 if ( relerr > sqrt(ROL::ROL_EPSILON<RealT>()) ) {
348 errorFlag += 1;
349 }
350 ROL::Ptr<std::vector<RealT> > vec_err_ptr = ROL::makePtr<std::vector<RealT>>(std_vec_err);
351 ROL::StdVector<RealT> vec_err(vec_err_ptr);
352 *outStream << std::scientific << "\n Linear algebra error: " << vec_err.norm() << std::endl;
353 if ( vec_err.norm() > 1e2*ROL::ROL_EPSILON<RealT>() ) {
354 errorFlag += 1;
355 }
356 }
357 catch (std::logic_error& err) {
358 *outStream << err.what() << "\n";
359 errorFlag = -1000;
360 }; // end try
361
362 if (errorFlag != 0)
363 std::cout << "End Result: TEST FAILED\n";
364 else
365 std::cout << "End Result: TEST PASSED\n";
366
367 return 0;
368
369}
370
Contains definitions for Rosenbrock's function.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Real norm() const
Returns where .
int dimension() const
Return dimension of the vector space.
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void scale(const Real alpha)
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
ROL::Ptr< std::vector< Element > > std_vec_
ROL::Ptr< const std::vector< Element > > getVector() const
ROL::Ptr< std::vector< Element > > getVector()
ROL::Ptr< OptStdVector< Real > > dual_vec_
OptDualStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
ROL::Ptr< const std::vector< Element > > getVector() const
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void plus(const ROL::Vector< Real > &x)
Compute , where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
ROL::Ptr< std::vector< Element > > getVector()
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
int dimension() const
Return dimension of the vector space.
OptStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Real norm() const
Returns where .
ROL::Ptr< OptDualStdVector< Real > > dual_vec_
ROL::Ptr< std::vector< Element > > std_vec_
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void scale(const Real alpha)
Compute where .
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
int main(int argc, char *argv[])
constexpr auto dim