ROL
zakharov/example_01.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
14 #define USE_HESSVEC 1
15 
16 #include "ROL_Algorithm.hpp"
17 #include "ROL_TrustRegionStep.hpp"
18 #include "ROL_RandomVector.hpp"
19 #include "ROL_StatusTest.hpp"
20 #include "ROL_StdVector.hpp"
21 #include "ROL_Zakharov.hpp"
22 #include "ROL_ParameterList.hpp"
24 #include "ROL_Stream.hpp"
25 #include "Teuchos_GlobalMPISession.hpp"
26 
27 #include <iostream>
28 
29 typedef double RealT;
30 
31 int main(int argc, char *argv[]) {
32 
33  using namespace Teuchos;
34 
35  typedef std::vector<RealT> vector;
36  typedef ROL::Vector<RealT> V; // Abstract vector
37  typedef ROL::StdVector<RealT> SV; // Concrete vector containing std::vector data
38 
39  GlobalMPISession mpiSession(&argc, &argv);
40 
41  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
42  auto outStream = ROL::makeStreamPtr( std::cout, argc > 1 );
43 
44  int errorFlag = 0;
45 
46  // *** Example body.
47 
48  try {
49 
50  int dim = 10; // Set problem dimension.
51 
52  std::string paramfile = "parameters.xml";
53  auto parlist = ROL::getParametersFromXmlFile( paramfile );
54 
55  // Define algorithm.
56  ROL::Ptr<ROL::Step<RealT>>
57  step = ROL::makePtr<ROL::TrustRegionStep<RealT>>(*parlist);
58  ROL::Ptr<ROL::StatusTest<RealT>>
59  status = ROL::makePtr<ROL::StatusTest<RealT>>(*parlist);
60  ROL::Algorithm<RealT> algo(step,status,false);
61 
62  // Iteration vector.
63  ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(dim, 0.0);
64 
65  // Vector of natural numbers.
66  ROL::Ptr<vector> k_ptr = ROL::makePtr<vector>(dim, 0.0);
67 
68  // For gradient and Hessian checks.
69  ROL::Ptr<vector> xtest_ptr = ROL::makePtr<vector>(dim, 0.0);
70  ROL::Ptr<vector> d_ptr = ROL::makePtr<vector>(dim, 0.0);
71  ROL::Ptr<vector> v_ptr = ROL::makePtr<vector>(dim, 0.0);
72  ROL::Ptr<vector> hv_ptr = ROL::makePtr<vector>(dim, 0.0);
73  ROL::Ptr<vector> ihhv_ptr = ROL::makePtr<vector>(dim, 0.0);
74 
75 
76  RealT left = -1e0, right = 1e0;
77  for (int i=0; i<dim; i++) {
78  (*x_ptr)[i] = 2;
79  (*k_ptr)[i] = i+1.0;
80  }
81 
82  ROL::Ptr<V> k = ROL::makePtr<SV>(k_ptr);
83  SV x(x_ptr);
84 
85  // Check gradient and Hessian.
86  SV xtest(xtest_ptr);
87  SV d(d_ptr);
88  SV v(v_ptr);
89  SV hv(hv_ptr);
90  SV ihhv(ihhv_ptr);
91 
92  ROL::RandomizeVector( xtest, left, right );
93  ROL::RandomizeVector( d, left, right );
94  ROL::RandomizeVector( v, left, right );
95 
97 
98  obj.checkGradient(xtest, d, true, *outStream); *outStream << "\n";
99  obj.checkHessVec(xtest, v, true, *outStream); *outStream << "\n";
100  obj.checkHessSym(xtest, d, v, true, *outStream); *outStream << "\n";
101 
102  // Check inverse Hessian.
103  RealT tol=0;
104  obj.hessVec(hv,v,xtest,tol);
105  obj.invHessVec(ihhv,hv,xtest,tol);
106  ihhv.axpy(-1,v);
107  *outStream << "Checking inverse Hessian" << std::endl;
108  *outStream << "||H^{-1}Hv-v|| = " << ihhv.norm() << std::endl;
109 
110 
111  // Run algorithm.
112  algo.run(x, obj, true, *outStream);
113 
114  // Get True Solution
115  ROL::Ptr<vector> xtrue_ptr = ROL::makePtr<vector>(dim, 0.0);
116  SV xtrue(xtrue_ptr);
117 
118 
119  // Compute Error
120  x.axpy(-1.0, xtrue);
121  RealT abserr = x.norm();
122  *outStream << std::scientific << "\n Absolute Error: " << abserr << std::endl;
123  if ( abserr > sqrt(ROL::ROL_EPSILON<RealT>()) ) {
124  errorFlag += 1;
125  }
126  }
127  catch (std::logic_error& err) {
128  *outStream << err.what() << "\n";
129  errorFlag = -1000;
130  }; // end try
131 
132  if (errorFlag != 0)
133  std::cout << "End Result: TEST FAILED\n";
134  else
135  std::cout << "End Result: TEST PASSED\n";
136 
137  return 0;
138 
139 }
140 
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:46
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:41
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.
constexpr auto dim