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_Algorithm.hpp"
52 #include "ROL_TrustRegionStep.hpp"
53 #include "ROL_StatusTest.hpp"
54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
56 #include <iostream>
57 
58 typedef double RealT;
59 
60 
61 /*** Declare two vector spaces. ***/
62 
63 // Forward declarations:
64 
65 template <class Real, class Element=Real>
66 class OptStdVector; // Optimization space.
67 
68 template <class Real, class Element=Real>
69 class OptDualStdVector; // Dual optimization space.
70 
71 
72 // Vector space definitions:
73 
74 // Optimization space.
75 template <class Real, class Element>
76 class OptStdVector : public ROL::Vector<Real> {
77 
78  typedef std::vector<Element> vector;
80 
81  typedef typename vector::size_type uint;
82 
83 private:
84 Teuchos::RCP<std::vector<Element> > std_vec_;
85 mutable Teuchos::RCP<OptDualStdVector<Real> > dual_vec_;
86 
87 public:
88 
89 OptStdVector(const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
90 
91 void plus( const ROL::Vector<Real> &x ) {
92  using Teuchos::RCP; using Teuchos::dyn_cast;
93  RCP<const vector> xvalptr = dyn_cast<const OptStdVector>(x).getVector();
94  uint dimension = std_vec_->size();
95  for (uint i=0; i<dimension; i++) {
96  (*std_vec_)[i] += (*xvalptr)[i];
97  }
98 }
99 
100 void scale( const Real alpha ) {
101  uint dimension = std_vec_->size();
102  for (uint i=0; i<dimension; i++) {
103  (*std_vec_)[i] *= alpha;
104  }
105 }
106 
107 Real dot( const ROL::Vector<Real> &x ) const {
108  Real val = 0;
109  using Teuchos::RCP; using Teuchos::dyn_cast;
110  RCP<const vector> xvalptr = dyn_cast<const OptStdVector>(x).getVector();
111  uint dimension = std_vec_->size();
112  for (uint i=0; i<dimension; i++) {
113  val += (*std_vec_)[i]*(*xvalptr)[i];
114  }
115  return val;
116 }
117 
118 Real norm() const {
119  Real val = 0;
120  val = std::sqrt( dot(*this) );
121  return val;
122 }
123 
124 Teuchos::RCP<ROL::Vector<Real> > clone() const {
125  return Teuchos::rcp( new OptStdVector( Teuchos::rcp( new std::vector<Element>(std_vec_->size()) ) ) );
126 }
127 
128 Teuchos::RCP<const std::vector<Element> > getVector() const {
129  return std_vec_;
130 }
131 
132 Teuchos::RCP<std::vector<Element> > getVector() {
133  return std_vec_;
134 }
135 
136 Teuchos::RCP<ROL::Vector<Real> > basis( const int i ) const {
137  using Teuchos::RCP; using Teuchos::rcp;
138  RCP<vector> e_rcp = rcp( new vector(std_vec_->size(),0.0) );
139  (*e_rcp)[i] = 1.0;
140 
141  RCP<V> e = rcp( new OptStdVector( e_rcp ) );
142 
143  return e;
144 }
145 
146 int dimension() const {return static_cast<int>(std_vec_->size());}
147 
148 const ROL::Vector<Real> & dual() const {
149  dual_vec_ = Teuchos::rcp( new OptDualStdVector<Real>( Teuchos::rcp( new std::vector<Element>(*std_vec_) ) ) );
150  return *dual_vec_;
151 }
152 
153 }; // class OptStdVector
154 
155 
156 // Dual optimization space.
157 template <class Real, class Element>
158 class OptDualStdVector : public ROL::Vector<Real> {
159 
160  typedef std::vector<Element> vector;
162 
163  typedef typename vector::size_type uint;
164 
165 private:
166 Teuchos::RCP<std::vector<Element> > std_vec_;
167 mutable Teuchos::RCP<OptStdVector<Real> > dual_vec_;
168 
169 public:
170 
171 OptDualStdVector(const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
172 
173 void plus( const ROL::Vector<Real> &x ) {
174  using Teuchos::RCP; using Teuchos::dyn_cast;
175  RCP<const vector> xvalptr = dyn_cast<const OptDualStdVector>(x).getVector();
176 
177  uint dimension = std_vec_->size();
178  for (uint i=0; i<dimension; i++) {
179  (*std_vec_)[i] += (*xvalptr)[i];
180  }
181 }
182 
183 void scale( const Real alpha ) {
184  uint dimension = std_vec_->size();
185  for (uint i=0; i<dimension; i++) {
186  (*std_vec_)[i] *= alpha;
187  }
188 }
189 
190 Real dot( const ROL::Vector<Real> &x ) const {
191  Real val = 0;
192  using Teuchos::RCP; using Teuchos::dyn_cast;
193  RCP<const vector> xvalptr = dyn_cast<const OptDualStdVector>(x).getVector();
194  uint dimension = std_vec_->size();
195  for (uint i=0; i<dimension; i++) {
196  val += (*std_vec_)[i]*(*xvalptr)[i];
197  }
198  return val;
199 }
200 
201 Real norm() const {
202  Real val = 0;
203  val = std::sqrt( dot(*this) );
204  return val;
205 }
206 
207 Teuchos::RCP<ROL::Vector<Real> > clone() const {
208  return Teuchos::rcp( new OptDualStdVector( Teuchos::rcp( new std::vector<Element>(std_vec_->size()) ) ) );
209 }
210 
211 Teuchos::RCP<const std::vector<Element> > getVector() const {
212  return std_vec_;
213 }
214 
215 Teuchos::RCP<std::vector<Element> > getVector() {
216  return std_vec_;
217 }
218 
219 Teuchos::RCP<ROL::Vector<Real> > basis( const int i ) const {
220  using Teuchos::RCP; using Teuchos::rcp;
221  RCP<vector> e_rcp = rcp( new vector(std_vec_->size(),0.0) );
222  (*e_rcp)[i] = 1.0;
223 
224  RCP<V> e = rcp( new OptDualStdVector( e_rcp ) );
225  return e;
226 }
227 
228 int dimension() const {return std_vec_->size();}
229 
230 const ROL::Vector<Real> & dual() const {
231  dual_vec_ = Teuchos::rcp( new OptStdVector<Real>( Teuchos::rcp( new std::vector<Element>(*std_vec_) ) ) );
232  return *dual_vec_;
233 }
234 
235 }; // class OptDualStdVector
236 
237 
238 /*** End of declaration of two vector spaces. ***/
239 
240 
241 
242 
243 
244 
245 int main(int argc, char *argv[]) {
246 
247  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
248 
249  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
250  int iprint = argc - 1;
251  Teuchos::RCP<std::ostream> outStream;
252  Teuchos::oblackholestream bhs; // outputs nothing
253  if (iprint > 0)
254  outStream = Teuchos::rcp(&std::cout, false);
255  else
256  outStream = Teuchos::rcp(&bhs, false);
257 
258  int errorFlag = 0;
259 
260  // *** Example body.
261 
262  try {
263 
265  int dim = 100; // Set problem dimension. Must be even.
266 
267  // Define algorithm.
268  Teuchos::ParameterList parlist;
269  std::string stepname = "Trust Region";
270  parlist.sublist("Step").sublist(stepname).set("Subproblem Solver", "Truncated CG");
271  parlist.sublist("General").sublist("Krylov").set("Iteration Limit",10);
272  parlist.sublist("General").sublist("Krylov").set("Relative Tolerance",1e-2);
273  parlist.sublist("General").sublist("Krylov").set("Absolute Tolerance",1e-4);
274  parlist.sublist("General").sublist("Secant").set("Use as Hessian",true);
275  parlist.sublist("Status Test").set("Gradient Tolerance",1.e-12);
276  parlist.sublist("Status Test").set("Step Tolerance",1.e-14);
277  parlist.sublist("Status Test").set("Iteration Limit",100);
278  ROL::Algorithm<RealT> algo(stepname,parlist);
279 
280  // Iteration Vector
281  Teuchos::RCP<std::vector<RealT> > x_rcp = Teuchos::rcp( new std::vector<RealT> (dim, 0.0) );
282  Teuchos::RCP<std::vector<RealT> > g_rcp = Teuchos::rcp( new std::vector<RealT> (dim, 0.0) );
283  // Set Initial Guess
284  for (int i=0; i<dim/2; i++) {
285  (*x_rcp)[2*i] = -1.2;
286  (*x_rcp)[2*i+1] = 1.0;
287  (*g_rcp)[2*i] = 0;
288  (*g_rcp)[2*i+1] = 0;
289  }
290 
291  OptStdVector<RealT> x(x_rcp); // Iteration Vector
292  OptDualStdVector<RealT> g(g_rcp); // zeroed gradient vector in dual space
293 
294  Teuchos::RCP<std::vector<RealT> > aa_rcp = Teuchos::rcp( new std::vector<RealT> (1, 1.0) );
295  OptDualStdVector<RealT> av(aa_rcp);
296  Teuchos::RCP<std::vector<RealT> > bb_rcp = Teuchos::rcp( new std::vector<RealT> (1, 2.0) );
297  OptDualStdVector<RealT> bv(bb_rcp);
298  Teuchos::RCP<std::vector<RealT> > cc_rcp = Teuchos::rcp( new std::vector<RealT> (1, 3.0) );
299  OptDualStdVector<RealT> cv(cc_rcp);
300  std::vector<RealT> std_vec_err = av.checkVector(bv,cv,true,*outStream);
301 
302  // Run Algorithm
303  std::vector<std::string> output = algo.run(x,g, obj, true, *outStream);
304 
305  // Get True Solution
306  Teuchos::RCP<std::vector<RealT> > xtrue_rcp = Teuchos::rcp( new std::vector<RealT> (dim, 1.0) );
307  OptStdVector<RealT> xtrue(xtrue_rcp);
308 
309  // Compute Errors
310  x.axpy(-1.0, xtrue);
311  RealT abserr = x.norm();
312  RealT relerr = abserr/xtrue.norm();
313  *outStream << std::scientific << "\n Absolute solution error: " << abserr;
314  *outStream << std::scientific << "\n Relative solution error: " << relerr;
315  if ( relerr > sqrt(ROL::ROL_EPSILON<RealT>()) ) {
316  errorFlag += 1;
317  }
318  Teuchos::RCP<std::vector<RealT> > vec_err_rcp = Teuchos::rcp( new std::vector<RealT> (std_vec_err) );
319  ROL::StdVector<RealT> vec_err(vec_err_rcp);
320  *outStream << std::scientific << "\n Linear algebra error: " << vec_err.norm() << std::endl;
321  if ( vec_err.norm() > 1e2*ROL::ROL_EPSILON<RealT>() ) {
322  errorFlag += 1;
323  }
324  }
325  catch (std::logic_error err) {
326  *outStream << err.what() << "\n";
327  errorFlag = -1000;
328  }; // end try
329 
330  if (errorFlag != 0)
331  std::cout << "End Result: TEST FAILED\n";
332  else
333  std::cout << "End Result: TEST PASSED\n";
334 
335  return 0;
336 
337 }
338 
Teuchos::RCP< std::vector< Element > > std_vec_
Teuchos::RCP< OptDualStdVector< Real > > dual_vec_
Rosenbrock&#39;s function.
void scale(const Real alpha)
Compute where .
OptStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec)
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Real norm() const
Returns where .
Teuchos::RCP< std::vector< Element > > getVector()
Teuchos::RCP< std::vector< Element > > getVector()
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
int main(int argc, char *argv[])
Teuchos::RCP< OptStdVector< Real > > dual_vec_
Contains definitions for Rosenbrock&#39;s function.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
int dimension() const
Return dimension of the vector space.
OptDualStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec)
Provides the std::vector implementation of the ROL::Vector interface.
void scale(const Real alpha)
Compute where .
Teuchos::RCP< const std::vector< Element > > getVector() const
Provides an interface to run optimization algorithms.
Real norm() const
Returns where .
int dimension() const
Return dimension of the vector space.
virtual std::vector< std::string > run(Vector< Real > &x, Objective< Real > &obj, bool print=false, std::ostream &outStream=std::cout)
Run algorithm on unconstrained problems (Type-U). This is the primary Type-U interface.
double RealT
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
virtual std::vector< Real > checkVector(const Vector< Real > &x, const Vector< Real > &y, const bool printToStream=true, std::ostream &outStream=std::cout) const
Verify vector-space methods.
Definition: ROL_Vector.hpp:259
Real dot(const ROL::Vector< Real > &x) const
Compute where .
Teuchos::RCP< std::vector< Element > > std_vec_
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
Teuchos::RCP< const std::vector< Element > > getVector() const
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Real norm() const
Returns where .