Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
GSLIntegrator.cxx
Go to the documentation of this file.
1// @(#)root/mathmore:$Id$
2// Authors: L. Moneta, A. Zsenei 08/2005
3
4 /**********************************************************************
5 * *
6 * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7 * *
8 * This library is free software; you can redistribute it and/or *
9 * modify it under the terms of the GNU General Public License *
10 * as published by the Free Software Foundation; either version 2 *
11 * of the License, or (at your option) any later version. *
12 * *
13 * This library is distributed in the hope that it will be useful, *
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16 * General Public License for more details. *
17 * *
18 * You should have received a copy of the GNU General Public License *
19 * along with this library (see file COPYING); if not, write *
20 * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21 * 330, Boston, MA 02111-1307 USA, or contact the author. *
22 * *
23 **********************************************************************/
24
25// Implementation file for class GSLIntegrator
26//
27// Created by: Lorenzo Moneta at Thu Nov 11 14:22:32 2004
28//
29// Last update: Thu Nov 11 14:22:32 2004
30//
31
32#include "Math/GSLIntegrator.h"
33
34#include "gsl/gsl_integration.h"
35
36#include "Math/IFunction.h"
38#include "GSLFunctionWrapper.h"
39
40// for toupper
41#include <algorithm>
42#include <functional>
43#include <ctype.h> // need to use c version of tolower defined here
44
45
46#include <iostream>
47
48
49
50namespace ROOT {
51namespace Math {
52
53
54
55
56GSLIntegrator::GSLIntegrator(const Integration::Type type , const Integration::GKRule rule, double absTol, double relTol, size_t size) :
57 fType(type),
58 fRule(rule),
59 fAbsTol(absTol),
60 fRelTol(relTol),
61 fSize(size),
62 fMaxIntervals(size),
63 fResult(0),fError(0),fStatus(-1),fNEval(-1),
64 fFunction(0),
65 fWorkspace(0)
66{
67 // constructor for all types of integrations
68 // allocate workspace (only if not adaptive algorithm)
71
72
73}
74
75
76
77GSLIntegrator::GSLIntegrator(double absTol, double relTol, size_t size) :
78 fType(Integration::kADAPTIVESINGULAR),
79 fRule(Integration::kGAUSS31),
80 fAbsTol(absTol),
81 fRelTol(relTol),
82 fSize(size),
83 fMaxIntervals(size),
84 fResult(0),fError(0),fStatus(-1),fNEval(-1),
85 fFunction(0),
86 fWorkspace(0)
87{
88 // constructor with default type (ADaptiveSingular) , rule is not needed
90
91}
92
93
94
95GSLIntegrator::GSLIntegrator(const Integration::Type type , double absTol, double relTol, size_t size) :
96 fType(type),
97 fRule(Integration::kGAUSS31),
98 fAbsTol(absTol),
99 fRelTol(relTol),
100 fSize(size),
101 fMaxIntervals(size),
102 fResult(0),fError(0),fStatus(-1),fNEval(-1),
103 fFunction(0),
104 fWorkspace(0)
105{
106
107 // constructor with default rule (gauss31) passing the type
108 // allocate workspace (only if not adaptive algorithm)
111
112}
113
114 GSLIntegrator::GSLIntegrator(const char * type , int rule, double absTol, double relTol, size_t size) :
115 fRule(Integration::kGAUSS31),
116 fAbsTol(absTol),
117 fRelTol(relTol),
118 fSize(size),
119 fMaxIntervals(size),
120 fResult(0),fError(0),fStatus(-1),fNEval(-1),
121 fFunction(0),
122 fWorkspace(0)
123{
124 //std::cout << type << std::endl;
125
127 if (type != 0) { // use this dafault
128 std::string typeName(type);
129 std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );
130 if (typeName == "NONADAPTIVE")
132 else if (typeName == "ADAPTIVE")
134 else {
135 if (typeName != "ADAPTIVESINGULAR")
136 MATH_WARN_MSG("GSLIntegrator","Use default type: AdaptiveSingular");
137 }
138 }
139
140
141 // constructor with default rule (gauss31) passing the type
142 // allocate workspace (only if not adaptive algorithm)
145
147
148}
149
150
152{
153 // delete workspace and function
154 if (fFunction) delete fFunction;
155 if (fWorkspace) delete fWorkspace;
156}
157
160{
161 // dummy copy ctr
162}
163
165{
166 // dummy operator=
167 if (this == &rhs) return *this; // time saving self-test
168
169 return *this;
170}
171
172
173
174
176 // fill GSLFunctionWrapper with the pointer to the function
177 if (fFunction ==0) fFunction = new GSLFunctionWrapper();
179 fFunction->SetParams ( p );
180}
181
183 // set function (make a copy of it)
184 if (fFunction ==0) fFunction = new GSLFunctionWrapper();
186}
187
188// evaluation methods
189
190double GSLIntegrator::Integral(double a, double b) {
191 // defined integral evaluation
192 // need here look at all types of algorithms
193 // find more elegant solution ? Use template OK, but need to chose algorithm statically , t.b.i.
194
195 if (!CheckFunction()) return 0;
196
198 size_t neval = 0; // need to export this ?
200 fNEval = neval;
201 }
202 else if (fType == Integration::kADAPTIVE) {
204 const int npts[6] = {15,21,31,41,51,61};
205 assert(fRule>=1 && fRule <=6);
206 fNEval = (fWorkspace->GetWS()->size)*npts[fRule-1]; // get size of workspace (number of iterations)
207 }
209
210 // singular integration - look if we know about singular points
211
212
214 fNEval = (fWorkspace->GetWS()->size) * 21; //since 21 point rule is used in qags
215 }
216 else {
217
218 fResult = 0;
219 fError = 0;
220 fStatus = -1;
221 std::cerr << "GSLIntegrator - Error: Unknown integration type" << std::endl;
222 throw std::exception(); //"Unknown integration type");
223 }
224
225 return fResult;
226
227}
228
229//=============================
230double GSLIntegrator::IntegralCauchy(double a, double b, double c) {
231 //eval integral with Cauchy principal value defined at the value c
232 if (!CheckFunction()) return 0;
233
234 fStatus = gsl_integration_qawc( fFunction->GetFunc(), a, b , c, fAbsTol, fRelTol, fMaxIntervals, fWorkspace->GetWS(), &fResult, &fError);
235 fNEval = (fWorkspace->GetWS()->size) * 15; // 15 point rule is used ?
236
237 return fResult;
238
239}
240
241double GSLIntegrator::IntegralCauchy(const IGenFunction & f, double a, double b, double c) {
242 //eval integral with Cauchy principal value defined at the value c
243
244 if (!CheckFunction()) return 0;
245 SetFunction(f);
246 return IntegralCauchy(a, b, c);
247
248}
249
250//==============================
251
252double GSLIntegrator::Integral( const std::vector<double> & pts) {
253 // integral eval with singularities
254
255 if (!CheckFunction()) return 0;
256
257 if (fType == Integration::kADAPTIVESINGULAR && pts.size() >= 2 ) {
258 // remove constness ( should be const in GSL ? )
259 double * p = const_cast<double *>(&pts.front() );
260 fStatus = gsl_integration_qagp( fFunction->GetFunc(), p, pts.size() , fAbsTol, fRelTol, fMaxIntervals, fWorkspace->GetWS(), &fResult, &fError);
261 fNEval = (fWorkspace->GetWS()->size) * 15; // 15 point rule is used ?
262 }
263 else {
264 fResult = 0;
265 fError = 0;
266 fStatus = -1;
267 std::cerr << "GSLIntegrator - Error: Unknown integration type or not enough singular points defined" << std::endl;
268 return 0;
269 }
270 return fResult;
271}
272
273
275 // Eval for indefined integrals: use QAGI method
276 // if method was chosen NO_ADAPTIVE WS does not exist create it
277
278 if (!CheckFunction()) return 0;
279
281
283 fNEval = (fWorkspace->GetWS()->size) * 15; // 15 point rule is used ?
284
285 return fResult;
286}
287
288
289
290double GSLIntegrator::IntegralUp( double a ) {
291 // Integral between [a, + inf]
292 // if method was chosen NO_ADAPTIVE WS does not exist create it
293
294 if (!CheckFunction()) return 0;
295
297
299 fNEval = (fWorkspace->GetWS()->size) * 21; // 21 point rule is used ?
300
301 return fResult;
302}
303
304
305
307 // Integral between [-inf, + b]
308 // if method was chosen NO_ADAPTIVE WS does not exist create it
309
310 if (!CheckFunction()) return 0;
311
313
315 fNEval = (fWorkspace->GetWS()->size) * 21; // 21 point rule is used ?
316
317 return fResult;
318}
319
320
321
322
323double GSLIntegrator::Integral(const IGenFunction & f, double a, double b) {
324 // use generic function interface
325 SetFunction(f);
326 return Integral(a,b);
327}
328
330 // use generic function interface
331 SetFunction(f);
332 return Integral();
333}
334
335double GSLIntegrator::IntegralUp(const IGenFunction & f, double a) {
336 // use generic function interface
337 SetFunction(f);
338 return IntegralUp(a);
339}
340
342 // use generic function interface
343 SetFunction(f);
344 return IntegralLow(b);
345}
346
347double GSLIntegrator::Integral(const IGenFunction & f, const std::vector<double> & pts) {
348 // use generic function interface
349 SetFunction(f);
350 return Integral(pts);
351}
352
353
354
355
356double GSLIntegrator::Integral( GSLFuncPointer f , void * p, double a, double b) {
357 // use c free function pointer
358 SetFunction( f, p);
359 return Integral(a,b);
360}
361
363 // use c free function pointer
364 SetFunction( f, p);
365 return Integral();
366}
367
368double GSLIntegrator::IntegralUp( GSLFuncPointer f, void * p, double a ) {
369 // use c free function pointer
370 SetFunction( f, p);
371 return IntegralUp(a);
372}
373
374double GSLIntegrator::IntegralLow( GSLFuncPointer f, void * p, double b ) {
375 // use c free function pointer
376 SetFunction( f, p);
377 return IntegralLow(b);
378}
379
380double GSLIntegrator::Integral( GSLFuncPointer f, void * p, const std::vector<double> & pts ) {
381 // use c free function pointer
382 SetFunction( f, p);
383 return Integral(pts);
384}
385
386
387
388double GSLIntegrator::Result() const { return fResult; }
389
390double GSLIntegrator::Error() const { return fError; }
391
392int GSLIntegrator::Status() const { return fStatus; }
393
394
395// get and setter methods
396
397// double GSLIntegrator::getAbsTolerance() const { return fAbsTol; }
398
399void GSLIntegrator::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }
400
401// double GSLIntegrator::getRelTolerance() const { return fRelTol; }
402
403void GSLIntegrator::SetRelTolerance(double relTol){ this->fRelTol = relTol; }
404
405
407
409 // check if a function has been previously set.
410 if (fFunction && fFunction->IsValid()) return true;
411 fStatus = -1; fResult = 0; fError = 0;
412 std::cerr << "GSLIntegrator - Error : Function has not been specified " << std::endl;
413 return false;
414}
415
417{
418 // set integration options
419 fType = opt.IntegratorType();
424 MATH_WARN_MSG("GSLIntegrator::SetOptions","Invalid rule options - use default ADAPTIVESINGULAR");
426 }
429 fSize = opt.WKSize();
432 int npts = opt.NPoints();
433 if ( npts >= Integration::kGAUSS15 && npts <= Integration::kGAUSS61)
434 fRule = (Integration::GKRule) npts;
435 else {
436 MATH_WARN_MSG("GSLIntegrator::SetOptions","Invalid rule options - use default GAUSS31");
438 }
439 }
440}
441
446 opt.SetWKSize(fSize);
448
450 opt.SetNPoints(fRule);
452 opt.SetNPoints( Integration::kGAUSS31 ); // fixed rule for adaptive singular
453 else
454 opt.SetNPoints( 0 ); // not available for the rest
455
456 return opt;
457}
458
459const char * GSLIntegrator::GetTypeName() const {
460 if (fType == IntegrationOneDim::kADAPTIVE) return "Adaptive";
461 if (fType == IntegrationOneDim::kADAPTIVESINGULAR) return "AdaptiveSingular";
462 if (fType == IntegrationOneDim::kNONADAPTIVE) return "NonAdaptive";
463 return "Undefined";
464}
465
466
467} // namespace Math
468} // namespace ROOT
size_t fSize
#define MATH_WARN_MSG(loc, str)
Definition Error.h:80
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
int gsl_integration_qagil(gsl_function *f, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
int gsl_integration_qag(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, int key, gsl_integration_workspace *workspace, double *result, double *abserr)
int gsl_integration_qagi(gsl_function *f, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
int gsl_integration_qagiu(gsl_function *f, double a, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
int gsl_integration_qags(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
int gsl_integration_qng(const gsl_function *f, double a, double b, double epsabs, double epsrel, double *result, double *abserr, size_t *neval)
int type
Definition TGX11.cxx:121
void SetAbsTolerance(double tol)
non-static methods for setting options
double RelTolerance() const
absolute tolerance
void SetRelTolerance(double tol)
set the relative tolerance
unsigned int WKSize() const
size of the workspace
double AbsTolerance() const
non-static methods for retrivieng options
void SetWKSize(unsigned int size)
set workspace size
Wrapper class to the gsl_function C structure.
void SetFunction(const FuncType &f)
fill the GSL C struct from a generic C++ callable object implementing operator()
void SetFuncPointer(GSLFuncPointer f)
set in the GSL C struct the pointer to the function evaluation
void SetParams(void *p)
set in the GSL C struct the extra-object pointer
bool IsValid()
check if function is valid (has been set)
Class for performing numerical integration of a function in one dimension.
virtual ROOT::Math::IntegratorOneDimOptions Options() const
get the option used for the integration
void SetAbsTolerance(double absTolerance)
set the desired absolute Error
virtual void SetOptions(const ROOT::Math::IntegratorOneDimOptions &opt)
set the options
Integration::GKRule fRule
void SetIntegrationRule(Integration::GKRule)
set the integration rule (Gauss-Kronrod rule).
GSLIntegrator & operator=(const GSLIntegrator &)
double IntegralCauchy(double a, double b, double c)
evaluate the Cauchy principal value of the integral of a previously defined function f over the defin...
double IntegralUp(const IGenFunction &f, double a)
evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
double IntegralLow(const IGenFunction &f, double b)
evaluate the Integral of a function f over the over the semi-infinite interval (-inf,...
GSLIntegrationWorkspace * fWorkspace
GSLFunctionWrapper * fFunction
int Status() const
return the Error Status of the last Integral calculation
double Error() const
return the estimate of the absolute Error of the last Integral calculation
void SetRelTolerance(double relTolerance)
set the desired relative Error
GSLIntegrator(double absTol=1.E-9, double relTol=1E-6, size_t size=1000)
Default constructor of GSL Integrator for Adaptive Singular integration.
const char * GetTypeName() const
return the name
Integration::Type fType
double Integral()
evaluate the Integral over the infinite interval (-inf,+inf) using the function previously set with G...
double Result() const
return the Result of the last Integral calculation
void SetFunction(const IGenFunction &f)
method to set the a generic integration function
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition IFunction.h:135
Numerical one dimensional integration options.
IntegrationOneDim::Type IntegratorType() const
type of the integrator (return the enumeration type)
void SetIntegrator(const char *name)
set 1D integrator name
void SetNPoints(unsigned int n)
Set number of points for active integration rule.
unsigned int NPoints() const
Number of points used by current integration rule.
Interface (abstract) class for 1D numerical integration It must be implemented by the concrate Integr...
GKRule
enumeration specifying the Gauss-KronRod integration rule for ADAPTIVE integration type
Type
enumeration specifying the integration types.
Namespace for new Math classes and functions.
double(* GSLFuncPointer)(double, void *)
Function pointer corresponding to gsl_function signature.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...