ROL
zakharov/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_Algorithm.hpp"
51 #include "ROL_LineSearchStep.hpp"
52 #include "ROL_RandomVector.hpp"
53 #include "ROL_StatusTest.hpp"
54 #include "ROL_StdVector.hpp"
55 #include "ROL_Zakharov.hpp"
56 #include "ROL_ParameterList.hpp"
58 #include "ROL_Stream.hpp"
59 #include "Teuchos_GlobalMPISession.hpp"
60 
61 #include <iostream>
62 
63 typedef double RealT;
64 
65 int main(int argc, char *argv[]) {
66 
67  using namespace Teuchos;
68 
69  typedef std::vector<RealT> vector;
70  typedef ROL::Vector<RealT> V; // Abstract vector
71  typedef ROL::StdVector<RealT> SV; // Concrete vector containing std::vector data
72 
73  GlobalMPISession mpiSession(&argc, &argv);
74 
75  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
76  auto outStream = ROL::makeStreamPtr( std::cout, argc > 1 );
77 
78  int errorFlag = 0;
79 
80  // *** Example body.
81 
82  try {
83 
84  int dim = 10; // Set problem dimension.
85 
86  std::string paramfile = "parameters.xml";
87  auto parlist = ROL::getParametersFromXmlFile( paramfile );
88 
89  // Define algorithm.
90  ROL::Algorithm<RealT> algo("Trust-Region",*parlist);
91 
92  // Iteration vector.
93  ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(dim, 0.0);
94 
95  // Vector of natural numbers.
96  ROL::Ptr<vector> k_ptr = ROL::makePtr<vector>(dim, 0.0);
97 
98  // For gradient and Hessian checks.
99  ROL::Ptr<vector> xtest_ptr = ROL::makePtr<vector>(dim, 0.0);
100  ROL::Ptr<vector> d_ptr = ROL::makePtr<vector>(dim, 0.0);
101  ROL::Ptr<vector> v_ptr = ROL::makePtr<vector>(dim, 0.0);
102  ROL::Ptr<vector> hv_ptr = ROL::makePtr<vector>(dim, 0.0);
103  ROL::Ptr<vector> ihhv_ptr = ROL::makePtr<vector>(dim, 0.0);
104 
105 
106  RealT left = -1e0, right = 1e0;
107  for (int i=0; i<dim; i++) {
108  (*x_ptr)[i] = 2;
109  (*k_ptr)[i] = i+1.0;
110  }
111 
112  ROL::Ptr<V> k = ROL::makePtr<SV>(k_ptr);
113  SV x(x_ptr);
114 
115  // Check gradient and Hessian.
116  SV xtest(xtest_ptr);
117  SV d(d_ptr);
118  SV v(v_ptr);
119  SV hv(hv_ptr);
120  SV ihhv(ihhv_ptr);
121 
122  ROL::RandomizeVector( xtest, left, right );
123  ROL::RandomizeVector( d, left, right );
124  ROL::RandomizeVector( v, left, right );
125 
127 
128  obj.checkGradient(xtest, d, true, *outStream); *outStream << "\n";
129  obj.checkHessVec(xtest, v, true, *outStream); *outStream << "\n";
130  obj.checkHessSym(xtest, d, v, true, *outStream); *outStream << "\n";
131 
132  // Check inverse Hessian.
133  RealT tol=0;
134  obj.hessVec(hv,v,xtest,tol);
135  obj.invHessVec(ihhv,hv,xtest,tol);
136  ihhv.axpy(-1,v);
137  *outStream << "Checking inverse Hessian" << std::endl;
138  *outStream << "||H^{-1}Hv-v|| = " << ihhv.norm() << std::endl;
139 
140 
141  // Run algorithm.
142  algo.run(x, obj, true, *outStream);
143 
144  // Get True Solution
145  ROL::Ptr<vector> xtrue_ptr = ROL::makePtr<vector>(dim, 0.0);
146  SV xtrue(xtrue_ptr);
147 
148 
149  // Compute Error
150  x.axpy(-1.0, xtrue);
151  RealT abserr = x.norm();
152  *outStream << std::scientific << "\n Absolute Error: " << abserr << std::endl;
153  if ( abserr > sqrt(ROL::ROL_EPSILON<RealT>()) ) {
154  errorFlag += 1;
155  }
156  }
157  catch (std::logic_error err) {
158  *outStream << err.what() << "\n";
159  errorFlag = -1000;
160  }; // end try
161 
162  if (errorFlag != 0)
163  std::cout << "End Result: TEST FAILED\n";
164  else
165  std::cout << "End Result: TEST PASSED\n";
166 
167  return 0;
168 
169 }
170 
void invHessVec(Vector< Real > &ihv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply inverse Hessian approximation to vector.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual std::vector< std::string > run(Vector< Real > &x, Objective< Real > &obj, bool print=false, std::ostream &outStream=std::cout, bool printVectors=false, std::ostream &vectorStream=std::cout)
Run algorithm on unconstrained problems (Type-U). This is the primary Type-U interface.
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].
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Ptr< ostream > makeStreamPtr(ostream &os, bool noSuppressOutput=true)
Definition: ROL_Stream.hpp:75
virtual std::vector< std::vector< Real > > checkGradient(const Vector< Real > &x, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference gradient check.
Vector< Real > V
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Provides an interface to run optimization algorithms.
Contains definitions for the Zakharov function as evaluated using only the ROL::Vector interface...
int main(int argc, char *argv[])
virtual std::vector< std::vector< Real > > checkHessVec(const Vector< Real > &x, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference Hessian-applied-to-vector check.
virtual std::vector< Real > checkHessSym(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &w, const bool printToStream=true, std::ostream &outStream=std::cout)
Hessian symmetry check.