Genfun::GaudiMathImplementation::NumericalDefiniteIntegral Class Reference

#include <NumericalDefiniteIntegral.h>

Collaboration diagram for Genfun::GaudiMathImplementation::NumericalDefiniteIntegral:

Detailed Description

This class allows the numerical evaluation of the following functions:.

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

Author:: Vanya BELYAEV Ivan.Belyaev@itep.ru

Date:: 2003-08-31

Definition at line 72 of file NumericalDefiniteIntegral.h.


Public Types
typedef std::vector< double >	Points
	typedef for vector of singular points
Public Member Functions
	FUNCTION_OBJECT_DEF (NumericalDefiniteIntegral)
	From CLHEP/GenericFunctions.
	NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const double b, const GaudiMath::Integration::Type type=GaudiMath::Integration::Adaptive, const GaudiMath::Integration::KronrodRule rule=GaudiMath::Integration::Default, const double epsabs=1.e-10, const double epsrel=1.e-7, const size_t size=1000)
	Standard constructor The function created with this constructor compute the following integral:.
	NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const double b, const Points &points, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
	Standard constructor The function created with this constructor compute the following integral:.
	NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const GaudiMath::Integration::Inf b=GaudiMath::Integration::Infinity, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
	Standard constructor The function created with this constructor compute the following integral:.
	NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const GaudiMath::Integration::Inf a, const double b, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
	Standard constructor The function created with this constructor compute the following integral:.
	NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const float epsabs=1e-9, const float epsrel=1.e-6, const size_t size=1000)
	Standard constructor The function created with this constructor compute the following integral:.
	NumericalDefiniteIntegral (const NumericalDefiniteIntegral &)
	copy constructor
virtual	~NumericalDefiniteIntegral ()
	destructor
virtual unsigned int	dimensionality () const
	dimensionality of the problem
virtual double	operator() (double argument) const
	Function value.
virtual double	operator() (const Argument &argument) const
	Function value.
virtual bool	hasAnalyticDerivative () const
	Does this function have an analytic derivative?
virtual Genfun::Derivative	partial (unsigned int index) const
	Derivatives.
const AbsFunction &	function () const
	accessor to the function itself
double	a () const
	integration limit
double	b () const
const Points &	points () const
	known singularities
double	epsabs () const
	absolute precision
double	epsrel () const
	relatiove precision
double	result () const
	previous result
double	error () const
	evaluate of previous error
size_t	size () const
GaudiMath::Integration::Type	type () const
	integration type
GaudiMath::Integration::Category	category () const
	integration category
GaudiMath::Integration::KronrodRule	rule () const
	integration rule
Protected Member Functions
double	QAGI (_Function *fun) const
double	QAGP (_Function *fun) const
double	QNG (_Function *fun) const
double	QAG (_Function *fun) const
double	QAGS (_Function *fun) const
_Workspace *	allocate () const
	allocate the integration workspace
_Workspace *	ws () const
StatusCode	Exception (const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const
Private Member Functions
	NumericalDefiniteIntegral ()
NumericalDefiniteIntegral &	operator= (const NumericalDefiniteIntegral &)
Private Attributes
const AbsFunction *	m_function
size_t	m_DIM
size_t	m_index
double	m_a
double	m_b
bool	m_ia
bool	m_ib
GaudiMath::Integration::Type	m_type
GaudiMath::Integration::Category	m_category
GaudiMath::Integration::KronrodRule	m_rule
Points	m_points
double *	m_pdata
double	m_epsabs
double	m_epsrel
double	m_result
double	m_error
size_t	m_size
_Workspace *	m_ws
Argument	m_argument
Argument	m_argF
Classes
struct	_Function
struct	_Workspace

Member Typedef Documentation

typedef std::vector<double> Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Points

typedef for vector of singular points

Definition at line 76 of file NumericalDefiniteIntegral.h.

Constructor & Destructor Documentation

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const double	a,
		const double	b,
		const GaudiMath::Integration::Type	type = `GaudiMath::Integration::Adaptive`,
		const GaudiMath::Integration::KronrodRule	rule = `GaudiMath::Integration::Default`,
		const double	epsabs = `1.e-10`,
		const double	epsrel = `1.e-7`,
		const size_t	size = `1000`
	)

Standard constructor The function created with this constructor compute the following integral:.

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

If function contains singularities, the type = Type::AdaptiveSingular need to be used

For faster integration of smooth function non-adaptive integration can be used: type = Type::NonAdaptive need to be used

For adaptive integration type = Type>>Adaptive one can specify the order of Gauss-Kronrad integration rule rule = KronrodRule::Gauss15 The higher-order rules give better accuracy for smooth functions, while lower-order rules save the time when the function contains local difficulties, such as discontinuites.

The GSL routine gsl_integration_qng is used for type = Type:NonAdaptive :
The GSL routine gsl_integration_qag is used for type = Type:Adaptive :
The GSL routine gsl_integration_qags is used for type = Type:AdaptiveSingular :

Parameters:

	function	funtion to be integrated
	index	variable index
	a	lower intgeration limit
	b	high integration limit
	type	integration type
	rule	integration rule (for adaptive integration)
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

If function contains singularities, the type = Type::AdaptiveSingular need to be used

For faster integration of smooth function non-adaptive integration can be used: type = Type::NonAdaptive need to be used

For adaptive integration type = Type>>Adaptive one can specify the order of Gauss-Kronrad integration rule rule = KronrodRule::Gauss15 The higher-order rules give better accuracy for smooth functions, while lower-order rules save the time when the function contains local difficulties, such as discontinuites.

The GSL routine gsl_integration_qng is used for type = Type:NonAdaptive :
The GSL routine gsl_integration_qag is used for type = Type:Adaptive :
The GSL routine gsl_integration_qags is used for type = Type:AdaptiveSingular :

Parameters:

	function	funtion to be integrated
	index	variable index
	a	lower intgeration limit
	b	high integration limit
	type	integration type
	rule	integration rule (for adaptive integration)
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration

Definition at line 99 of file NumericalDefiniteIntegral.cpp.

00108       : AbsFunction () 
00109       , m_function  ( function.clone          ()    ) 
00110       , m_DIM       ( 0                             )
00111       , m_index     ( index                         )
00112       , m_a         ( a      ) 
00113       , m_b         ( b      )
00114       , m_ia        ( false  ) 
00115       , m_ib        ( false  )
00116       , m_type      ( type   ) 
00117       , m_category  ( GaudiMath::Integration::Finite )
00118       , m_rule      ( rule       ) 
00119       , m_points    (            ) 
00120       , m_pdata     ( 0          ) 
00121       , m_epsabs    ( epsabs     )
00122       , m_epsrel    ( epsrel     )
00123       , m_result    ( GSL_NEGINF ) 
00124       , m_error     ( GSL_POSINF ) 
00125       , m_size      ( size       ) 
00126       , m_ws        ()
00127       , m_argument  ()
00128       , m_argF      () 
00129     {
00130       if ( GaudiMath::Integration::Fixed == m_rule ) 
00131         { m_rule = GaudiMath::Integration::Default ; }
00132       if ( function.dimensionality() < 2          ) 
00133         { Exception("::constructor: invalid dimensionality ") ; }        
00134       if ( m_index >= function.dimensionality()   ) 
00135         { Exception("::constructor: invalid variable index") ; }
00136 
00137       m_DIM = function.dimensionality() - 1 ;
00138       m_argument = Argument( m_DIM      ) ;
00139       m_argF     = Argument( m_DIM + 1  ) ;
00140 
00141     };

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const double	a,
		const double	b,
		const Points &	points,
		const double	epsabs = `1e-9`,
		const double	epsrel = `1.e-6`,
		const size_t	size = `1000`
	)

Standard constructor The function created with this constructor compute the following integral:.

standard constructor

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

where function have known singulatities

The GSL routine gsl_integration_qagp is used for integration

Parameters:

	function	funtion to be integrated
	index	variable index
	a	lower intgeration limit
	b	high integration limit
	pointns	vector of know singular points
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration
	function	the base function
	index	the variable index
	a	integration limit
	b	integration limit
	points	list of known function singularities
	epsabs	absolute precision for integration
	epsrel	relative precision for integration

Definition at line 154 of file NumericalDefiniteIntegral.cpp.

00162       : AbsFunction () 
00163       , m_function  ( function.clone()             )
00164       , m_DIM       ( 0                            ) 
00165       , m_index     ( index                        ) 
00166       , m_a         ( a                            )
00167       , m_b         ( b                            )
00168       , m_ia        ( false  ) 
00169       , m_ib        ( false  )
00170       , m_type      ( GaudiMath::Integration::    Other )
00171       , m_category  ( GaudiMath::Integration:: Singular )
00172       , m_rule      ( GaudiMath::Integration::    Fixed )
00173       , m_points    ( points  ) 
00174       , m_pdata     ( 0       ) 
00175       , m_epsabs    ( epsabs  )
00176       , m_epsrel    ( epsrel  ) 
00177       //
00178       , m_result    ( GSL_NEGINF                          ) 
00179       , m_error     ( GSL_POSINF                          ) 
00180       //
00181       , m_size      ( size                                ) 
00182       , m_ws        ( 0                                   ) 
00183       , m_argument  () 
00184       , m_argF      () 
00185     {
00186       if ( function.dimensionality() < 2          ) 
00187         { Exception("::constructor: invalid dimensionality ") ; }        
00188       if ( m_index >= function.dimensionality()   ) 
00189         { Exception("::constructor: invalid variable index") ; }
00190 
00191       m_DIM = function.dimensionality() - 1 ;
00192       m_argument = Argument( m_DIM      ) ;
00193       m_argF     = Argument( m_DIM + 1  ) ;
00194 
00195       const double l1 = std::min ( a , b ) ;
00196       const double l2 = std::max ( a , b ) ;
00197       m_points.push_back ( l1 ) ;
00198       m_points.push_back ( l2 ) ;
00199       std::sort ( m_points.begin() , m_points.end() ) ;
00200       m_points.erase( std::unique( m_points.begin () , 
00201                                    m_points.end   () ) , 
00202                       m_points.end() );
00203       
00204       Points::iterator lower = 
00205         std::lower_bound ( m_points.begin () , m_points.end () , l1 ) ;
00206       m_points.erase     ( m_points.begin () , lower                ) ;
00207       Points::iterator upper = 
00208         std::upper_bound ( m_points.begin () , m_points.end () , l2 ) ;
00209       m_points.erase     ( upper             , m_points.end ()      ) ;
00210       
00211       m_pdata = new double[ m_points.size() ] ;
00212       std::copy( m_points.begin() , m_points.end() , m_pdata );
00213     };

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const double	a,
		const GaudiMath::Integration::Inf	b = `GaudiMath::Integration::Infinity`,
		const double	epsabs = `1e-9`,
		const double	epsrel = `1.e-6`,
		const size_t	size = `1000`
	)

Standard constructor The function created with this constructor compute the following integral:.

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

The GSL routine gsl_integration_qagiu is used for integration

Parameters:

	function	funtion to be integrated
	index	variable index
	a	lower intgeration limit
	b	indication for infinity
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

The GSL routine gsl_integration_qagiu is used for integration

Parameters:

	function	funtion to be integrated
	index	variable index
	a	lower intgeration limit
	b	indication for infinity
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration

Definition at line 241 of file NumericalDefiniteIntegral.cpp.

00248       : AbsFunction() 
00249       , m_function  ( function.clone()             )
00250       , m_DIM       ( 0                            ) 
00251       , m_index     ( index                        ) 
00252       , m_a         ( a                            )
00253       , m_b         ( GSL_POSINF                   )
00254       , m_ia        ( false  ) 
00255       , m_ib        ( true   )
00256       , m_type      ( GaudiMath::Integration::    Other )
00257       , m_category  ( GaudiMath::Integration:: Infinite )
00258       , m_rule      ( GaudiMath::Integration::    Fixed )
00259       , m_points    (         ) 
00260       , m_pdata     ( 0       ) 
00261       , m_epsabs    ( epsabs  )
00262       , m_epsrel    ( epsrel  ) 
00263       //
00264       , m_result    ( GSL_NEGINF                          ) 
00265       , m_error     ( GSL_POSINF                          ) 
00266       //
00267       , m_size      ( size                                ) 
00268       , m_ws        ( 0                                   ) 
00269       , m_argument  () 
00270       , m_argF      () 
00271     {
00272       if ( function.dimensionality() < 2          ) 
00273         { Exception("::constructor: invalid dimensionality ") ; }        
00274       if ( m_index >= function.dimensionality()   ) 
00275         { Exception("::constructor: invalid variable index") ; }
00276 
00277       m_DIM = function.dimensionality() - 1 ;
00278       m_argument = Argument( m_DIM      ) ;
00279       m_argF     = Argument( m_DIM + 1  ) ;
00280 
00281     };

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const GaudiMath::Integration::Inf	a,
		const double	b,
		const double	epsabs = `1e-9`,
		const double	epsrel = `1.e-6`,
		const size_t	size = `1000`
	)

Standard constructor The function created with this constructor compute the following integral:.

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

The GSL routine gsl_integration_qagiu is used for integration

Parameters:

	function	funtion to be integrated
	index	variable index
	a	lower intgeration limit
	b	indication for infinity
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

The GSL routine gsl_integration_qagiu is used for integration

Parameters:

	function	funtion to be integrated
	index	variable index
	a	lower intgeration limit
	b	indication for infinity
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration

Definition at line 310 of file NumericalDefiniteIntegral.cpp.

00317       : AbsFunction() 
00318       , m_function  ( function.clone()             )
00319       , m_DIM       ( 0                            ) 
00320       , m_index     ( index                        ) 
00321       , m_a         ( GSL_NEGINF                   )
00322       , m_b         ( b                            )
00323       , m_ia        ( true   ) 
00324       , m_ib        ( false  )
00325       , m_type      ( GaudiMath::Integration::    Other )
00326       , m_category  ( GaudiMath::Integration:: Infinite )
00327       , m_rule      ( GaudiMath::Integration::    Fixed )
00328       , m_points    (         ) 
00329       , m_pdata     ( 0       ) 
00330       , m_epsabs    ( epsabs  )
00331       , m_epsrel    ( epsrel  ) 
00332       //
00333       , m_result    ( GSL_NEGINF                          ) 
00334       , m_error     ( GSL_POSINF                          ) 
00335       //
00336       , m_size      ( size                                ) 
00337       , m_ws        ( 0                                   ) 
00338       , m_argument  () 
00339       , m_argF      () 
00340     {
00341       if ( function.dimensionality() < 2          ) 
00342         { Exception("::constructor: invalid dimensionality ") ; }        
00343       if ( m_index >= function.dimensionality()   ) 
00344         { Exception("::constructor: invalid variable index") ; }
00345 
00346       m_DIM = function.dimensionality() - 1 ;
00347       m_argument = Argument( m_DIM      ) ;
00348       m_argF     = Argument( m_DIM + 1  ) ;
00349     };

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const float	epsabs = `1e-9`,
		const float	epsrel = `1.e-6`,
		const size_t	size = `1000`
	)

Standard constructor The function created with this constructor compute the following integral:.

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

The GSL routine gsl_integration_qagi is used for integration

Parameters:

	function	funtion to be integrated
	index	variable index
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i$

The GSL routine gsl_integration_qagi is used for integration

Parameters:

	function	funtion to be integrated
	index	variable index
	epsabs	required absolute precision
	epsrel	required relative precision
	size	maximal number of bisections for adaptive integration

Definition at line 375 of file NumericalDefiniteIntegral.cpp.

00380       : AbsFunction() 
00381       , m_function  ( function.clone()             )
00382       , m_DIM       ( 0      ) 
00383       , m_index     ( index                        ) 
00384       , m_a         ( GSL_NEGINF                   )
00385       , m_b         ( GSL_POSINF                   )
00386       , m_ia        ( true   ) 
00387       , m_ib        ( true   )
00388       , m_type      ( GaudiMath::Integration::    Other )
00389       , m_category  ( GaudiMath::Integration:: Infinite )
00390       , m_rule      ( GaudiMath::Integration::    Fixed )
00391       , m_points    (         ) 
00392       , m_pdata     ( 0       ) 
00393       , m_epsabs    ( epsabs  )
00394       , m_epsrel    ( epsrel  ) 
00395       //
00396       , m_result    ( GSL_NEGINF                          ) 
00397       , m_error     ( GSL_POSINF                          ) 
00398       //
00399       , m_size      ( size                                ) 
00400       , m_ws        ( 0                                   ) 
00401       , m_argument  () 
00402       , m_argF      () 
00403     {
00404       if ( function.dimensionality() < 2          ) 
00405         { Exception("::constructor: invalid dimensionality ") ; }        
00406       if ( m_index >= function.dimensionality()   ) 
00407         { Exception("::constructor: invalid variable index") ; }
00408       
00409       m_DIM = function.dimensionality() - 1 ;
00410       m_argument = Argument( m_DIM      ) ;
00411       m_argF     = Argument( m_DIM + 1  ) ;
00412       
00413     };

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral ( const NumericalDefiniteIntegral & right )

copy constructor

Definition at line 417 of file NumericalDefiniteIntegral.cpp.

00418       : AbsFunction () 
00419       , m_function  ( right.m_function->clone() ) 
00420       , m_DIM       ( right.m_DIM      )
00421       , m_index     ( right.m_index    )
00422       , m_a         ( right.m_a        ) 
00423       , m_b         ( right.m_b        )
00424       , m_ia        ( right.m_ia       ) 
00425       , m_ib        ( right.m_ib       )
00426       , m_type      ( right.m_type     ) 
00427       , m_category  ( right.m_category )
00428       , m_rule      ( right.m_rule     ) 
00429       , m_points    ( right.m_points   ) 
00430       , m_pdata     ( 0                ) 
00431       , m_epsabs    ( right.m_epsabs   )
00432       , m_epsrel    ( right.m_epsrel   )
00433       , m_result    ( GSL_NEGINF       ) 
00434       , m_error     ( GSL_POSINF       ) 
00435       , m_size      ( right.m_size     ) 
00436       , m_ws        ( 0                )
00437       , m_argument  ( right.m_argument )
00438       , m_argF      ( right.m_argF     )
00439     {
00440       if( 0 != right.m_pdata ) 
00441         {
00442           m_pdata = new double[m_points.size()] ;
00443           std::copy( m_points.begin() , m_points.end() , m_pdata );
00444         } 
00445     };

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::~NumericalDefiniteIntegral ( ) [virtual]

destructor

Definition at line 447 of file NumericalDefiniteIntegral.cpp.

00448     {
00449       if( 0 != m_function ) { delete m_function ; m_function = 0 ; }
00450       if( 0 != m_pdata    ) { delete m_pdata    ; m_pdata    = 0 ; }      
00451     };

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral ( ) [private]

Member Function Documentation

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::FUNCTION_OBJECT_DEF ( NumericalDefiniteIntegral )

From CLHEP/GenericFunctions.

virtual unsigned int Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::dimensionality ( ) const [inline, virtual]

dimensionality of the problem

Definition at line 304 of file NumericalDefiniteIntegral.h.

00304 { return m_DIM ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::operator() ( double argument ) const [virtual]

Function value.

evaluate the function

Definition at line 470 of file NumericalDefiniteIntegral.cpp.

00471     {
00472       // reset the result and the error  
00473       m_result = GSL_NEGINF ;
00474       m_error  = GSL_POSINF ;
00475       
00476       // check the argument 
00477       if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
00478       
00479       m_argument[0] = argument ;
00480       return (*this) ( m_argument );
00481     };

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::operator() ( const Argument & argument ) const [virtual]

Function value.

evaluate the function

Definition at line 503 of file NumericalDefiniteIntegral.cpp.

00504     {
00505       // reset the result and the error  
00506       m_result = GSL_NEGINF ;
00507       m_error  = GSL_POSINF ;
00508       
00509       // check the argument 
00510       if( argument.dimension() != m_DIM ) 
00511         { Exception ( "operator(): invalid argument size " ) ; };
00512       
00513       // copy the argument 
00514       
00515       for( size_t i  = 0 ; i < m_DIM ; ++i )
00516         { 
00517           m_argument [i] = argument [i] ;
00518           const size_t j =  i < m_index ? i : i + 1 ;
00519           m_argF     [j] = argument [i] ;
00520         };
00521       
00522       // create the helper object 
00523       GSL_Helper helper( *m_function , m_argF , m_index );
00524       
00525       // use GSL to evaluate the numerical derivative 
00526       gsl_function F ;
00527       F.function = &GSL_Adaptor ;
00528       F.params   = &helper                 ;
00529       _Function F1    ;
00530       F1.fn      = &F ;
00531       
00532       if        (  GaudiMath::Integration::Infinite         == category () )
00533         { return   QAGI ( &F1 ) ; }                                // RETURN
00534       
00535       if ( m_a == m_b ) 
00536         { 
00537           m_result = 0    ; m_error  = 0    ;                      // EXACT 
00538           return m_result ;                                       // RETURN 
00539         }
00540       
00541       if        (  GaudiMath::Integration::Singular         == category () )
00542          { return   QAGP ( &F1 ) ; }                                // RETURN
00543       else if   (  GaudiMath::Integration::Finite           == category () )
00544         if      (  GaudiMath::Integration::NonAdaptive      == type     () ) 
00545           { return QNG  ( &F1 ) ; }                                // RETURN
00546         else if (  GaudiMath::Integration::Adaptive         == type     () ) 
00547           { return QAG  ( &F1 ) ; }                                // RETURN
00548         else if (  GaudiMath::Integration::AdaptiveSingular == type     () ) 
00549           { return QAGS ( &F1 ) ; }                                // RETURN
00550         else 
00551           { Exception ( "::operator(): invalid type "  ); }
00552       else 
00553         { Exception ( "::operator(): invalid category "  ); }
00554       
00555       return 0 ;
00556     };

virtual bool Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::hasAnalyticDerivative ( ) const [inline, virtual]

Does this function have an analytic derivative?

Definition at line 312 of file NumericalDefiniteIntegral.h.

00312 { return true ;}

Genfun::Derivative Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::partial ( unsigned int index ) const [virtual]

Derivatives.

Definition at line 488 of file NumericalDefiniteIntegral.cpp.

00489     {
00490       if      ( idx >= m_DIM   )
00491         { Exception ( "::partial(i): invalid variable index " ) ; };
00492       //
00493       const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
00494       return Genfun::FunctionNoop( &aux ) ;
00495     };

const AbsFunction& Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::function ( ) const [inline]

accessor to the function itself

Definition at line 320 of file NumericalDefiniteIntegral.h.

00320 { return *m_function ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::a ( ) const [inline]

integration limit

Definition at line 322 of file NumericalDefiniteIntegral.h.

00322 { return         m_a ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::b ( ) const [inline]

Definition at line 323 of file NumericalDefiniteIntegral.h.

00323 { return         m_b ; }

const Points& Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::points ( ) const [inline]

known singularities

Definition at line 325 of file NumericalDefiniteIntegral.h.

00325 { return    m_points ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::epsabs ( ) const [inline]

absolute precision

Definition at line 327 of file NumericalDefiniteIntegral.h.

00327 { return    m_epsabs ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::epsrel ( ) const [inline]

relatiove precision

Definition at line 329 of file NumericalDefiniteIntegral.h.

00329 { return    m_epsrel ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::result ( ) const [inline]

previous result

Definition at line 332 of file NumericalDefiniteIntegral.h.

00332 { return    m_result ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::error ( ) const [inline]

evaluate of previous error

Definition at line 334 of file NumericalDefiniteIntegral.h.

00334 { return     m_error ; }

size_t Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::size ( void ) const [inline]

Definition at line 337 of file NumericalDefiniteIntegral.h.

00337 { return      m_size ; }

GaudiMath::Integration::Type Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::type ( ) const [inline]

integration type

Definition at line 341 of file NumericalDefiniteIntegral.h.

00341 { return      m_type ; }

GaudiMath::Integration::Category Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::category ( ) const [inline]

integration category

Definition at line 344 of file NumericalDefiniteIntegral.h.

00344 { return  m_category ; }

GaudiMath::Integration::KronrodRule Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::rule ( ) const [inline]

integration rule

Definition at line 347 of file NumericalDefiniteIntegral.h.

00347 { return      m_rule ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QAGI ( _Function * fun ) const [protected]

Definition at line 578 of file NumericalDefiniteIntegral.cpp.

00579     {
00580       // check the argument 
00581       if ( 0 == F    ) { Exception("::QAGI: invalid function"); }
00582       
00583       // allocate workspace 
00584       if ( 0 == ws() ) { allocate() ; }
00585       
00586       int ierror = 0 ;
00587       
00588       if ( m_ia && m_ib ) 
00589         {
00590           ierror = gsl_integration_qagi  ( F->fn                , 
00591                                            m_epsabs  , m_epsrel , 
00592                                            size ()   , ws()->ws , 
00593                                            &m_result , &m_error ) ;
00594         }
00595       else if ( m_ia )
00596         {
00597           ierror = gsl_integration_qagil ( F->fn     , m_b      , 
00598                                            m_epsabs  , m_epsrel , 
00599                                            size ()   , ws()->ws , 
00600                                            &m_result , &m_error ) ;
00601         }
00602       else if ( m_ib )
00603         {
00604           ierror = gsl_integration_qagiu ( F->fn     , m_a      , 
00605                                            m_epsabs  , m_epsrel , 
00606                                            size ()   , ws()->ws , 
00607                                            &m_result , &m_error ) ;
00608         }
00609       else 
00610         { Exception ( "::QAGI: invalid mode" ) ; };
00611       
00612       if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" , 
00613                                 __FILE__ , __LINE__ , ierror ) ;}
00614       
00615       return m_result ;
00616     };

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QAGP ( _Function * fun ) const [protected]

Definition at line 622 of file NumericalDefiniteIntegral.cpp.

00623     {
00624       if( 0 == F ) { Exception("QAGP::invalid function") ; }
00625       
00626       // no known singular points ?
00627       if( points().empty() || 0 == m_pdata ) { return QAGS( F ) ; }
00628       
00629       const size_t npts = points().size();
00630       
00631       // use GSL 
00632       int ierror =
00633         gsl_integration_qagp ( F->fn                , 
00634                                m_pdata   , npts     , 
00635                                m_epsabs  , m_epsrel ,
00636                                size ()   , ws()->ws , 
00637                                &m_result , &m_error ) ;
00638       
00639       if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " , 
00640                                 __FILE__ , __LINE__ , ierror ) ; }
00641       
00642       // the sign
00643       if ( m_a > m_b ) { m_result *= -1 ; }
00644       
00645       return m_result ;
00646     };

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QNG ( _Function * fun ) const [protected]

Definition at line 652 of file NumericalDefiniteIntegral.cpp.

00653     {
00654       if( 0 == F ) { Exception("QNG::invalid function") ; }
00655       
00656       // integration limits 
00657       const double low  = std::min ( m_a , m_b ) ;
00658       const double high = std::max ( m_a , m_b ) ;
00659       
00660       size_t neval = 0 ;
00661       int ierror = 
00662         gsl_integration_qng ( F->fn                 ,  
00663                               low       ,      high , 
00664                               m_epsabs  ,  m_epsrel , 
00665                               &m_result , &m_error  , &neval  ) ;
00666       
00667       if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " , 
00668                                 __FILE__ , __LINE__ , ierror ) ; }
00669       
00670       // sign
00671       if ( m_a > m_b ) { m_result *= -1 ; }
00672       
00673       return m_result ;
00674     };

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QAG ( _Function * fun ) const [protected]

Definition at line 681 of file NumericalDefiniteIntegral.cpp.

00682     {
00683       if( 0 == F ) { Exception("QAG::invalid function") ; }
00684       
00685       // allocate workspace 
00686       if( 0 == ws () ) { allocate () ; }
00687       
00688       // integration limits 
00689       const double low  = std::min ( m_a , m_b ) ;
00690       const double high = std::max ( m_a , m_b ) ;
00691       
00692       int ierror = 
00693         gsl_integration_qag ( F->fn                    , 
00694                               low       ,         high , 
00695                               m_epsabs  ,     m_epsrel , 
00696                               size   () , (int) rule() , ws ()->ws , 
00697                               &m_result ,     &m_error             );
00698       
00699       if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " , 
00700                                 __FILE__ , __LINE__ , ierror ) ; }
00701       
00702       // sign
00703       if ( m_a > m_b ) { m_result *= -1 ; }
00704       
00705       return m_result ;
00706     };

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QAGS ( _Function * fun ) const [protected]

Definition at line 712 of file NumericalDefiniteIntegral.cpp.

00713     {
00714       if( 0 == F ) { Exception("QAG::invalid function") ; }
00715       
00716       if ( m_a == m_b ) 
00717         { 
00718           m_result = 0    ;
00719           m_error  = 0    ;   // EXACT !
00720           return m_result ; 
00721         }
00722       
00723       // allocate workspace 
00724       if( 0 == ws () ) { allocate () ; }
00725       
00726       // integration limits 
00727       const double low  = std::min ( m_a , m_b ) ;
00728       const double high = std::max ( m_a , m_b ) ;
00729       
00730       int ierror = 
00731         gsl_integration_qags ( F->fn                , 
00732                                low       , high     , 
00733                                m_epsabs  , m_epsrel , 
00734                                size   () , ws()->ws , 
00735                                &m_result , &m_error );
00736       
00737       if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " , 
00738                                 __FILE__ , __LINE__ , ierror ) ; }
00739       
00740       // sign
00741       if ( m_a > m_b ) { m_result *= -1 ; }
00742       
00743       return m_result ;
00744     };

NumericalDefiniteIntegral::_Workspace * Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::allocate ( ) const [protected]

allocate the integration workspace

Definition at line 563 of file NumericalDefiniteIntegral.cpp.

00564     {
00565       if ( 0 != m_ws ) { return m_ws; }
00566       gsl_integration_workspace* aux = 
00567         gsl_integration_workspace_alloc( size () );
00568       if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; };
00569       m_ws = new _Workspace() ;
00570       m_ws->ws = aux ;
00571       return m_ws ;
00572     };

_Workspace* Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::ws ( ) const [inline, protected]

Definition at line 365 of file NumericalDefiniteIntegral.h.

00366       { return m_ws ; };

StatusCode Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Exception	(	const std::string &	message,
		const StatusCode &	sc = `StatusCode::FAILURE`
	)			const `[protected]`

Definition at line 457 of file NumericalDefiniteIntegral.cpp.

00459     {
00460       throw GaudiException( "NumericalDefiniteIntegral::" + message , 
00461                             "*GaudiMath*" , sc ) ;
00462       return sc ;
00463     };

NumericalDefiniteIntegral& Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::operator= ( const NumericalDefiniteIntegral & ) [private]

Member Data Documentation

const AbsFunction* Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_function [private]

Definition at line 382 of file NumericalDefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_DIM [private]

Definition at line 383 of file NumericalDefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_index [private]

Definition at line 384 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_a [private]

Definition at line 386 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_b [private]

Definition at line 387 of file NumericalDefiniteIntegral.h.

bool Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_ia [private]

Definition at line 388 of file NumericalDefiniteIntegral.h.

bool Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_ib [private]

Definition at line 389 of file NumericalDefiniteIntegral.h.

GaudiMath::Integration::Type Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_type [private]

Definition at line 391 of file NumericalDefiniteIntegral.h.

GaudiMath::Integration::Category Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_category [private]

Definition at line 392 of file NumericalDefiniteIntegral.h.

GaudiMath::Integration::KronrodRule Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_rule [private]

Definition at line 393 of file NumericalDefiniteIntegral.h.

Points Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_points [private]

Definition at line 395 of file NumericalDefiniteIntegral.h.

double* Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_pdata [private]

Definition at line 396 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_epsabs [private]

Definition at line 398 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_epsrel [private]

Definition at line 399 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_result [mutable, private]

Definition at line 401 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_error [mutable, private]

Definition at line 402 of file NumericalDefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_size [private]

Definition at line 404 of file NumericalDefiniteIntegral.h.

_Workspace* Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_ws [mutable, private]

Definition at line 405 of file NumericalDefiniteIntegral.h.

Argument Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_argument [mutable, private]

Definition at line 407 of file NumericalDefiniteIntegral.h.

Argument Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_argF [mutable, private]

Definition at line 408 of file NumericalDefiniteIntegral.h.

The documentation for this class was generated from the following files:

/afs/cern.ch/sw/Gaudi/releases/GAUDI/GAUDI_v21r4/GaudiGSL/GaudiMath/NumericalDefiniteIntegral.h
/afs/cern.ch/sw/Gaudi/releases/GAUDI/GAUDI_v21r4/GaudiGSL/src/Lib/NumericalDefiniteIntegral.cpp

Gaudi Framework, version v21r4

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral Class Reference

Detailed Description

Public Types

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Classes

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation