ROOT  6.06/09
Reference Guide
testGSLRootFinder.cxx
Go to the documentation of this file.
1 #include "Math/Polynomial.h"
2 #include "Math/Functor.h"
3 #include "Math/RootFinder.h"
5 
6 #ifdef HAVE_ROOTLIBS
7 #include "TStopwatch.h"
8 #else
9 struct TStopwatch {
10  void Start(){}
11  void Stop(){}
12  void Reset(){}
13  double RealTime() { return 0; }
14  double CpuTime() { return 0; }
15 };
16 #endif
17 
18 #include <iostream>
19 
20 
21 typedef double ( * FP ) ( double, void * );
22 
23 const int iterTest = 10000;
24 int myfuncCalls = 0;
25 
26 const double absTol = 1E-3;
27 //const double relTol = 1E-6;
28 
29 double myfunc ( double x) {
30  myfuncCalls += 1;
31  return x*x - 5;
32 }
33 // required for GSL type func
34 double myfunc_gsl(double x, void *) {
35  return myfunc(x);
36 }
37 
38 double myfunc_deriv ( double x ) {
39  return 2.0*x;
40 }
41 double myfunc_deriv_gsl ( double x, void * /*params*/) {
42  return myfunc_deriv(x);
43 }
44 
45 void myfunc_fdf( double x, void * /*params*/, double *y, double *dy) {
46  *y = x*x - 5;
47  *dy = 2.0*x;
48 }
49 
50 template<class RF>
51 int findRoot( RF * r ) {
52  //int returnCode = r->Solve( 100, absTol, relTol);
53  bool ret = r->Solve();
54  int returnCode = !ret; //r->Status();
55 
56  return returnCode;
57 }
58 
59 template<class RF>
60 int printStats( RF * r, int returnCode, TStopwatch& timer ) {
61  std::cout << "\nTest " << r->Name() << " algorithm " << std::endl;
62 
63  double root = r->Root();
64 
65  std::cout << "Return code: " << returnCode << std::endl;
66  std::cout << "Result: " << root << " n iters = " << r->Iterations() << std::endl;
67  std::cout << "Exact result: " << sqrt(5.0) << " difference: " << root - sqrt(5.0) << std::endl;
68  std::cout << "Time: " << timer.RealTime()/(double) iterTest << std::endl;
69  std::cout << "Number of calls to function: " << myfuncCalls/iterTest << std::endl;
70 
71  if ( fabs(root - sqrt(5.0)) > absTol ) {
72  std::cerr << "Test Root finder with " << r->Name() << " failed " << std::endl;
73  return 1;
74  }
75  return 0;
76 }
77 
78 
80 
81  int returnCode = 0;
82  int status = 0;
84 
85  ROOT::Math::Polynomial polyf(2);
86  std::vector<double> p(3);
87  p[0] = -5;
88  p[1] = 0;
89  p[2] = 1;
90 
91  polyf.SetParameters(&p[0]);
92 
93  //ROOT::Math::IGenFunction *func = &polyf;
95 
97  timer.Reset(); timer.Start(); myfuncCalls = 0;
98  for (int i = 0; i < iterTest; ++i)
99  {
100  rf1->SetFunction( func, 0, 5);
101  returnCode = findRoot(rf1);
102  }
103  timer.Stop();
104  status += printStats(rf1, returnCode, timer);
105 
106 
108  timer.Reset(); timer.Start(); myfuncCalls = 0;
109  for (int i = 0; i < iterTest; ++i)
110  {
111  rf2->SetFunction( func, 0, 5);
112  returnCode = findRoot(rf2);
113  }
114  timer.Stop();
115  status += printStats(rf2, returnCode, timer);
116 
117  // methods using derivatives
118 
119 
122  timer.Reset(); timer.Start(); myfuncCalls = 0;
123  for (int i = 0; i < iterTest; ++i)
124  {
125  rf3->SetFunction( gfunc, 1);
126  returnCode = findRoot(rf3);
127  }
128  timer.Stop();
129  status += printStats(rf3, returnCode, timer);
130 
131 
133  timer.Reset(); timer.Start(); myfuncCalls = 0;
134  for (int i = 0; i < iterTest; ++i)
135  {
136  rf4->SetFunction( gfunc, 1);
137  returnCode = findRoot(rf4);
138  }
139  timer.Stop();
140  status += printStats(rf4, returnCode, timer);
141 
142 
144  void * ptr2 = 0;
145  timer.Reset(); timer.Start(); myfuncCalls = 0;
146  for (int i = 0; i < iterTest; ++i)
147  {
149  returnCode = findRoot(rf5);
150  }
151  timer.Stop();
152  status += printStats(rf5, returnCode, timer);
153 
154 
155 
156  // the following two examples won't work when interpreted CINT
157  //const FP funcPtr = &myfunc;
159  void * ptr1 = 0;
161  //ROOT::Math::RootFinder<ROOT::Math::Roots::Brent> *rf6 = new ROOT::Math::RootFinder<ROOT::Math::Roots::Brent>;
162  timer.Reset(); timer.Start(); myfuncCalls = 0;
163  for (int i = 0; i < iterTest; ++i)
164  {
165  rf6->SetFunction( funcPtr, ptr1, 0.0, 5.0);
166  returnCode = findRoot(rf6);
167  }
168  timer.Stop();
169  status += printStats(rf6, returnCode, timer);
170 
171  return status;
172 
173 }
174 
175 int main() {
176 
177  int status = 0;
178 
179  status += testGSLRootFinder();
180 
181  return status;
182 }
User Class to find the Root of one dimensional functions.
Definition: RootFinder.h:88
double myfunc(double x)
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
const double absTol
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
GradFunctor1D class for one-dimensional gradient functions.
Definition: Functor.h:705
static byte * ptr1
Definition: gifdecode.c:16
double(* FP)(double, void *)
int myfuncCalls
double myfunc_deriv_gsl(double x, void *)
TAlienJobStatus * status
Definition: TAlienJob.cxx:51
double myfunc_gsl(double x, void *)
double sqrt(double)
TStopwatch timer
Definition: pirndm.C:37
Double_t x[n]
Definition: legend1.C:17
double(* GSLFuncPointer)(double, void *)
a Newton algorithm, which computes the derivative at each iteration See the GSL manual for more infor...
bool SetFunction(const IGradFunction &f, double xstart)
Sets the function for algorithms using derivatives.
void myfunc_fdf(double x, void *, double *y, double *dy)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
virtual void SetParameters(const double *p)
Set the parameter values.
ROOT::R::TRInterface & r
Definition: Object.C:4
static byte * ptr2
Definition: gifdecode.c:16
Double_t E()
Definition: TMath.h:54
int findRoot(RF *r)
double myfunc_deriv(double x)
Parametric Function class describing polynomials of order n.
Definition: Polynomial.h:63
bool SetFunction(const IGenFunction &f, double xlow, double xup)
Provide to the solver the function and the initial search interval [xlow, xup] for algorithms not usi...
Definition: RootFinder.h:124
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
int testGSLRootFinder()
Brent-Dekker algorithm which combines an interpolation strategy with the bisection algorithm See the ...
int printStats(RF *r, int returnCode, TStopwatch &timer)
const int iterTest
Functor1D class for one-dimensional functions.
Definition: Functor.h:492
int main()
bool SetFunction(const IGenFunction &f, double xlow, double xup)
Sets the function for the rest of the algorithms.
Stopwatch class.
Definition: TStopwatch.h:30