Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral Class Reference

#include <NumericalIndefiniteIntegral.h>

Collaboration diagram for Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral:

Detailed Description

The simple class for numerical integrations.

It allows to evaluate following indefinite integrals:

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_i , x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{x_i} 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 , x_{i+1}, \dots , x_n \right) = \int\limits_{x_i}^{a} 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 , x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{x_i} 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 , x_{i+1}, \dots , x_n \right) = \int\limits_{x_i}^{+\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 76 of file NumericalIndefiniteIntegral.h.


Public Types
typedef std::vector< double >	Points
	typedef for vector of singular points
Public Member Functions
	FUNCTION_OBJECT_DEF (NumericalIndefiniteIntegral)
	From CLHEP/GenericFunctions.
	NumericalIndefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const GaudiMath::Integration::Limit limit=GaudiMath::Integration::VariableHighLimit, const GaudiMath::Integration::Type type=GaudiMath::Integration::Adaptive, const GaudiMath::Integration::KronrodRule rule=GaudiMath::Integration::Default, const double epsabs=1e-10, const double epsrel=1.e-7, const size_t size=1000)
	Standard constructor The function, created with this constructor evaluates following indefinite integral:.
	NumericalIndefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const Points &points, const GaudiMath::Integration::Limit limit=GaudiMath::Integration::VariableHighLimit, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
	standard constructor
	NumericalIndefiniteIntegral (const AbsFunction &function, const size_t index, const GaudiMath::Integration::Limit limit=GaudiMath::Integration::VariableHighLimit, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
	standard constructor The function, created with this constructor evaluates following indefinite integral:
	NumericalIndefiniteIntegral (const NumericalIndefiniteIntegral &)
	copy constructor
virtual	~NumericalIndefiniteIntegral ()
	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
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::Limit	limit () const
	integration limit
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
	NumericalIndefiniteIntegral ()
NumericalIndefiniteIntegral &	operator= (const NumericalIndefiniteIntegral &)
Private Attributes
const AbsFunction *	m_function
size_t	m_DIM
size_t	m_index
double	m_a
GaudiMath::Integration::Limit	m_limit
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
Classes
struct	_Function
struct	_Workspace

Member Typedef Documentation

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

typedef for vector of singular points

Definition at line 80 of file NumericalIndefiniteIntegral.h.

Constructor & Destructor Documentation

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

Standard constructor The function, created with this constructor evaluates following indefinite integral:.

Standard constructor.

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

for value of limit = VariableHighLimit and the integral

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

for value of limit = VariableLowLimit

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 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	the base function
	index	the variable index
	a	integration limit
	limit	flag to distinguisch low variable limit from high variable limit
	type	the integration type (adaptive, non-adaptive or adaptive for singular functions
	key	Gauss-Kronrad integration rule
	epsabs	absolute precision for integration
	epsrel	relative precision for integration
	lim	bisection limit
	function	the base function
	index	the variable index
	a	integration limit
	limit	flag to distinguisch low variable limit from high variable limit
	type	the integration type (adaptive, non-adaptive or adaptive for singular functions
	key	Gauss-Kronrad integration rule
	epsabs	absolute precision for integration
	epsrel	relative precision for integration
	lim	bisection limit

Definition at line 71 of file NumericalIndefiniteIntegral.cpp.

00080       : AbsFunction () 
00081       , m_function  ( function.clone()                    )
00082       , m_DIM       ( function.dimensionality()           )
00083       , m_index     ( index                               )
00084       , m_a         ( a                                   ) 
00085       , m_limit     ( limit                               )
00086       , m_type      ( type                                ) 
00087       , m_category  ( GaudiMath::Integration::Finite      ) 
00088       , m_rule      ( rule                                ) 
00089       //
00090       , m_points    (                                     ) 
00091       , m_pdata     ( 0                                   )
00092       //
00093       , m_epsabs    ( epsabs                              ) 
00094       , m_epsrel    ( epsrel                              ) 
00095       //
00096       , m_result    ( GSL_NEGINF                          ) 
00097       , m_error     ( GSL_POSINF                          ) 
00098       //
00099       , m_size      ( size                                ) 
00100       , m_ws        ( 0                                   ) 
00101       , m_argument  ( function.dimensionality()           ) 
00102     {
00103       if ( GaudiMath::Integration::Fixed == m_rule ) 
00104         { m_rule = GaudiMath::Integration::Default ; }
00105       if ( m_index >= m_DIM  ) 
00106         { Exception("::constructor: invalid variable index") ; }
00107     };

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::NumericalIndefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const double	a,
		const Points &	points,
		const GaudiMath::Integration::Limit	limit = `GaudiMath::Integration::VariableHighLimit`,
		const double	epsabs = `1e-9`,
		const double	epsrel = `1.e-6`,
		const size_t	size = `1000`
	)

standard constructor

The function, created with this constructor evaluates following indefinite integral:

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

for value of limit = VariableHighLimit and the integral

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

for value of limit = VariableLowLimit

The integrand is assumed to have a known discontinuities

The GSL routine gsl_integration_qagp is used for integration

Parameters:

	function	the base function
	index	the variable index
	a	integration limit
	limit	flag to distinguisch low variable limit from high variable limit
	points	list of known function singularities
	epsabs	absolute precision for integration
	epsrel	relative precision for integration
	function	the base function
	index	the variable index
	a	integration limit
	limit	flag to distinguisch low variable limit from high variable limit
	points	list of known function singularities
	epsabs	absolute precision for integration
	epsrel	relative precision for integration

Definition at line 121 of file NumericalIndefiniteIntegral.cpp.

00129       : AbsFunction () 
00130       , m_function  ( function.clone()          )
00131       , m_DIM       ( function.dimensionality() ) 
00132       , m_index     ( index                     ) 
00133       , m_a         ( a                         )
00134       , m_limit     ( limit                     )
00135       , m_type      ( GaudiMath::Integration::    Other )
00136       , m_category  ( GaudiMath::Integration:: Singular )
00137       , m_rule      ( GaudiMath::Integration::    Fixed )
00138       , m_points    ( points  ) 
00139       , m_pdata     ( 0       ) 
00140       , m_epsabs    ( epsabs  )
00141       , m_epsrel    ( epsrel  ) 
00142       //
00143       , m_result    ( GSL_NEGINF                          ) 
00144       , m_error     ( GSL_POSINF                          ) 
00145       //
00146       , m_size      ( size                                ) 
00147       , m_ws        ( 0                                   ) 
00148       , m_argument  ( function.dimensionality()           ) 
00149     {
00150       if ( m_index >= m_DIM ) 
00151         { Exception("::constructor: invalid variable index") ; } 
00152       m_pdata = new double[ 2 + m_points.size() ] ;
00153       m_points.push_back( a ) ;
00154       std::sort( m_points.begin() , m_points.end() ) ;
00155       m_points.erase ( std::unique( m_points.begin () , 
00156                                     m_points.end   () ) , m_points.end() );
00157     };

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::NumericalIndefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const GaudiMath::Integration::Limit	limit = `GaudiMath::Integration::VariableHighLimit`,
		const double	epsabs = `1e-9`,
		const double	epsrel = `1.e-6`,
		const size_t	size = `1000`
	)

standard constructor The function, created with this constructor evaluates following indefinite integral:

standard constructor the integral limt is assumed to be infinity

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

for value of limit = VariableHighLimit and the integral

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

for value of limit = VariableLowLimit

The GSL routines gsl_integration_qagil and gsl_integration_qagiu are used for adapive integration

Parameters:

	function	the base function
	index	the variable index
	limit	flag to distinguisch low variable limit from high variable limit
	singularities	list of known function singularities
	function	the base function
	index	the variable index
	limit	flag to distinguisch low variable limit from high variable limit

Definition at line 169 of file NumericalIndefiniteIntegral.cpp.

00175       : AbsFunction () 
00176       , m_function  ( function.clone() ) 
00177       , m_DIM       ( function.dimensionality() )
00178       , m_index     ( index            )
00179       , m_a         ( GSL_NEGINF       ) // should not be used! 
00180       , m_limit     ( limit            )
00181       , m_type      ( GaudiMath::Integration::    Other ) 
00182       , m_category  ( GaudiMath::Integration:: Infinite )
00183       , m_rule      ( GaudiMath::Integration::    Fixed )
00184       , m_points    (            ) 
00185       , m_pdata     ( 0          )
00186       , m_epsabs    ( epsabs     ) 
00187       , m_epsrel    ( epsrel     )
00188       , m_result    ( GSL_NEGINF ) 
00189       , m_error     ( GSL_POSINF ) 
00190       , m_size      ( size       ) 
00191       , m_ws        ( 0          ) 
00192       , m_argument  ( function.dimensionality()           ) 
00193     {
00194       if ( m_index >= m_DIM ) 
00195         { Exception("::constructor: invalid variable index") ; }
00196     };

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

copy constructor

Definition at line 204 of file NumericalIndefiniteIntegral.cpp.

00205       : AbsFunction () 
00206       , m_function  ( right.m_function->clone() ) 
00207       , m_DIM       ( right.m_DIM      ) 
00208       , m_index     ( right.m_index    ) 
00209       , m_a         ( right.m_a        ) 
00210       , m_limit     ( right.m_limit    ) 
00211       , m_type      ( right.m_type     ) 
00212       , m_category  ( right.m_category )
00213       , m_rule      ( right.m_rule     )
00214       , m_points    ( right.m_points   ) 
00215       , m_pdata     ( 0 )           // attention 
00216       , m_epsabs    ( right.m_epsabs   ) 
00217       , m_epsrel    ( right.m_epsrel   ) 
00218       , m_result    ( GSL_NEGINF       ) 
00219       , m_error     ( GSL_POSINF       )
00220       , m_size      ( right.m_size     ) 
00221       , m_ws        ( 0                )
00222       , m_argument  ( right.m_argument )
00223     {
00224       m_pdata = new double[ 2 + m_points.size() ] ; // attention! 
00225     };

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

destructor

Definition at line 232 of file NumericalIndefiniteIntegral.cpp.

00233     {
00234       if( 0 != m_ws ) 
00235         { 
00236           gsl_integration_workspace_free ( m_ws->ws ) ;
00237           delete m_ws ;
00238           m_ws = 0 ;
00239         }
00240       if ( 0 != m_pdata    ) { delete m_pdata    ; m_pdata    = 0 ; }
00241       if ( 0 != m_function ) { delete m_function ; m_function = 0 ; }
00242     };

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

Member Function Documentation

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::FUNCTION_OBJECT_DEF ( NumericalIndefiniteIntegral )

From CLHEP/GenericFunctions.

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

dimensionality of the problem

Definition at line 278 of file NumericalIndefiniteIntegral.h.

00278 { return m_DIM ; }

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

Function value.

evaluate the function

Definition at line 262 of file NumericalIndefiniteIntegral.cpp.

00263     {
00264       // reset the result and the error  
00265       m_result = GSL_NEGINF ;
00266       m_error  = GSL_POSINF ;
00267       
00268       // check the argument 
00269       if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
00270       
00271       m_argument[0] = argument ;
00272       return (*this) ( m_argument );
00273     };

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

Function value.

evaluate the function

Definition at line 302 of file NumericalIndefiniteIntegral.cpp.

00303     {
00304       // reset the result and the error  
00305       m_result = GSL_NEGINF ;
00306       m_error  = GSL_POSINF ;
00307       
00308       // check the argument 
00309       if( argument.dimension() != m_DIM ) 
00310         { Exception ( "operator(): invalid argument size " ) ; };
00311       
00312       // copy the argument 
00313       {for( size_t i  = 0 ; i < m_DIM ; ++i ){ m_argument[i] = argument[i];}}
00314       
00315       // create the helper object 
00316       GSL_Helper helper( *m_function , m_argument , m_index );
00317       
00318       // use GSL to evaluate the numerical derivative 
00319       gsl_function F ;
00320       F.function = &GSL_Adaptor ;
00321       F.params   = &helper                 ;
00322       _Function F1    ;
00323       F1.fn      = &F ;
00324       
00325       if        (  GaudiMath::Integration::Infinite         == category () )
00326         { return   QAGI ( &F1 ) ; }                                // RETURN
00327       else if   (  GaudiMath::Integration::Singular         == category () )
00328         { return   QAGP ( &F1 ) ; }                                // RETURN
00329       else if   (  GaudiMath::Integration::Finite           == category () )
00330         if      (  GaudiMath::Integration::NonAdaptive      == type     () ) 
00331           { return QNG  ( &F1 ) ; }                                // RETURN
00332         else if (  GaudiMath::Integration::Adaptive         == type     () ) 
00333           { return QAG  ( &F1 ) ; }                                // RETURN
00334         else if (  GaudiMath::Integration::AdaptiveSingular == type     () ) 
00335           { return QAGS ( &F1 ) ; }                                // RETURN
00336         else 
00337           { Exception ( "::operator(): invalid type "  ); }
00338       else 
00339         { Exception ( "::operator(): invalid category "  ); }
00340       
00341       return 0 ;
00342     };

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

Does this function have an analytic derivative?

Definition at line 286 of file NumericalIndefiniteIntegral.h.

00286 { return true ;}

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

Derivatives.

Definition at line 280 of file NumericalIndefiniteIntegral.cpp.

00281     {
00282       if      ( idx >= m_DIM   )
00283         { Exception ( "::partial(i): invalid variable index " ) ; };
00284       if      ( idx != m_index )
00285         {
00286           const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
00287           return Genfun::FunctionNoop( &aux ) ;
00288         }
00289       else if ( GaudiMath::Integration::VariableLowLimit == limit () ) 
00290         { 
00291           const AbsFunction& aux = -1 * function() ;
00292           return Genfun::FunctionNoop( &aux ) ;
00293         }
00294       const AbsFunction& aux = function() ;
00295       return Genfun::FunctionNoop( &aux ) ;
00296     };

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

accessor to the function itself

Definition at line 294 of file NumericalIndefiniteIntegral.h.

00294 { return *m_function ; }

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

integration limit

Definition at line 296 of file NumericalIndefiniteIntegral.h.

00296 { return         m_a ; }

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

known singularities

Definition at line 298 of file NumericalIndefiniteIntegral.h.

00298 { return    m_points ; }

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

absolute precision

Definition at line 300 of file NumericalIndefiniteIntegral.h.

00300 { return    m_epsabs ; }

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

relatiove precision

Definition at line 302 of file NumericalIndefiniteIntegral.h.

00302 { return    m_epsrel ; }

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

previous result

Definition at line 305 of file NumericalIndefiniteIntegral.h.

00305 { return    m_result ; }

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

evaluate of previous error

Definition at line 307 of file NumericalIndefiniteIntegral.h.

00307 { return     m_error ; }

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

Definition at line 310 of file NumericalIndefiniteIntegral.h.

00310 { return      m_size ; }

GaudiMath::Integration::Limit Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::limit ( ) const [inline]

integration limit

Definition at line 314 of file NumericalIndefiniteIntegral.h.

00314 { return     m_limit ; }

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

integration type

Definition at line 317 of file NumericalIndefiniteIntegral.h.

00317 { return      m_type ; }

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

integration category

Definition at line 320 of file NumericalIndefiniteIntegral.h.

00320 { return  m_category ; }

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

integration rule

Definition at line 323 of file NumericalIndefiniteIntegral.h.

00323 { return      m_rule ; }

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

Definition at line 364 of file NumericalIndefiniteIntegral.cpp.

00365     {
00366       // check the argument 
00367       if( 0 == F ) { Exception("::QAGI: invalid function"); }
00368       
00369       const double x = m_argument[m_index] ;
00370       
00371       // allocate workspace 
00372       if( 0 == ws() ) { allocate() ; }
00373       
00374       int ierror = 0 ;
00375       switch ( limit() ) 
00376         {
00377         case GaudiMath::Integration::VariableLowLimit  : 
00378           ierror = gsl_integration_qagiu ( F->fn     , x        ,  
00379                                            m_epsabs  , m_epsrel , 
00380                                            size ()   , ws()->ws , 
00381                                            &m_result , &m_error ) ; break ;
00382         case GaudiMath::Integration::VariableHighLimit :
00383           ierror = gsl_integration_qagil ( F->fn     , x        ,
00384                                            m_epsabs  , m_epsrel , 
00385                                            size ()   , ws()->ws , 
00386                                            &m_result , &m_error ) ; break ;
00387         default :
00388           Exception ( "::QAGI: invalid mode" ) ;
00389         };
00390       
00391       if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI" , 
00392                                 __FILE__ , __LINE__ , ierror ) ;}
00393       
00394       return m_result ;
00395     };

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

Definition at line 401 of file NumericalIndefiniteIntegral.cpp.

00402     {
00403       if( 0 == F ) { Exception("QAGP::invalid function") ; }
00404       
00405       const double x = m_argument[m_index] ;
00406       if ( m_a == x  ) 
00407         { 
00408           m_result = 0    ;
00409           m_error  = 0    ;   // EXACT !
00410           return m_result ; 
00411         }
00412       
00413       // no known singular points ?
00414       if( points().empty() ) { return QAGS( F ) ; }
00415       
00416       // integration limits 
00417       const double a = std::min ( m_a , x ) ;
00418       const double b = std::max ( m_a , x ) ;
00419       
00420       // "active" singular points
00421       Points::const_iterator lower = 
00422         std::lower_bound ( points().begin() , points().end() , a ) ;
00423       Points::const_iterator upper = 
00424         std::upper_bound ( points().begin() , points().end() , b ) ;
00425       
00426       Points pnts ( upper - lower ) ;
00427       std::copy( lower , upper , pnts.begin() );
00428       if ( *lower != a       ) { pnts.insert( pnts.begin () , a ) ; }
00429       if ( *upper != b       ) { pnts.insert( pnts.end   () , b ) ; }
00430       std::copy( pnts.begin() , pnts.end() , m_pdata );
00431       const size_t npts = pnts.size() ;
00432       
00433       // use GSL 
00434       int ierror = 
00435         gsl_integration_qagp ( F->fn                , 
00436                                m_pdata   , npts     , 
00437                                m_epsabs  , m_epsrel ,
00438                                size ()   , ws()->ws , 
00439                                &m_result , &m_error ) ;
00440       
00441       if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI " , 
00442                                 __FILE__ , __LINE__ , ierror ) ; }
00443       
00444       // sign
00445       if      ( GaudiMath::Integration::VariableHighLimit == limit() 
00446                 &&  x < m_a  ) { m_result *= -1 ; }
00447       else if ( GaudiMath::Integration::VariableLowLimit  == limit() 
00448                 &&  x > m_a  ) { m_result *= -1 ; }
00449       
00450       return m_result ;
00451     };

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

Definition at line 457 of file NumericalIndefiniteIntegral.cpp.

00458     {
00459       if( 0 == F ) { Exception("QNG::invalid function") ; }
00460       
00461       const double x = m_argument[m_index] ;
00462       if ( m_a == x  ) 
00463         { 
00464           m_result = 0    ;
00465           m_error  = 0    ;   // EXACT !
00466           return m_result ; 
00467         }
00468       
00469       // integration limits 
00470       const double a = std::min ( m_a , x ) ;
00471       const double b = std::max ( m_a , x ) ;
00472       
00473       size_t neval = 0 ;
00474       int ierror = 
00475         gsl_integration_qng ( F->fn                 ,  
00476                               a         ,         b , 
00477                               m_epsabs  ,  m_epsrel , 
00478                               &m_result , &m_error  , &neval  ) ;
00479       
00480       if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " , 
00481                                 __FILE__ , __LINE__ , ierror ) ; }
00482       
00483       // sign
00484       if      ( GaudiMath::Integration::VariableHighLimit == limit() 
00485                 &&  x < m_a  ) { m_result *= -1 ; }
00486       else if ( GaudiMath::Integration::VariableLowLimit  == limit() 
00487                 &&  x > m_a  ) { m_result *= -1 ; }
00488 
00489       return m_result ;
00490     };

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

Definition at line 495 of file NumericalIndefiniteIntegral.cpp.

00496     {
00497       if( 0 == F ) { Exception("QAG::invalid function") ; }
00498       
00499       const double x = m_argument[m_index] ;
00500       if ( m_a == x  ) 
00501         { 
00502           m_result = 0    ;
00503           m_error  = 0    ;   // EXACT !
00504           return m_result ; 
00505         }
00506       
00507       // allocate workspace 
00508       if( 0 == ws () ) { allocate () ; }
00509       
00510       // integration limits 
00511       const double a = std::min ( m_a , x ) ;
00512       const double b = std::max ( m_a , x ) ;
00513       
00514       int ierror = 
00515         gsl_integration_qag ( F->fn                    , 
00516                               a         ,            b , 
00517                               m_epsabs  ,     m_epsrel , 
00518                               size   () , (int) rule() , ws ()->ws , 
00519                               &m_result ,     &m_error             );
00520       
00521       if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAG " , 
00522                                 __FILE__ , __LINE__ , ierror ) ; }
00523       
00524       // sign
00525       if      ( GaudiMath::Integration::VariableHighLimit == limit() 
00526                 &&  x < m_a  ) { m_result *= -1 ; }
00527       else if ( GaudiMath::Integration::VariableLowLimit  == limit() 
00528                 &&  x > m_a  ) { m_result *= -1 ; }
00529       
00530       return m_result ;
00531     };

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

Definition at line 537 of file NumericalIndefiniteIntegral.cpp.

00538     {
00539       if( 0 == F ) { Exception("QAG::invalid function") ; }
00540       
00541       const double x = m_argument[m_index] ;
00542       if ( m_a == x  ) 
00543         { 
00544           m_result = 0    ;
00545           m_error  = 0    ;   // EXACT !
00546           return m_result ; 
00547         }
00548       
00549       // allocate workspace 
00550       if( 0 == ws () ) { allocate () ; }
00551       
00552       // integration limits 
00553       const double a = std::min ( m_a , x ) ;
00554       const double b = std::max ( m_a , x ) ;
00555       
00556       int ierror = 
00557         gsl_integration_qags ( F->fn                , 
00558                                a         , b        , 
00559                                m_epsabs  , m_epsrel , 
00560                                size   () , ws()->ws , 
00561                                &m_result , &m_error );
00562       
00563       if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGS " , 
00564                                 __FILE__ , __LINE__ , ierror ) ; }
00565       
00566       // sign
00567       if      ( GaudiMath::Integration::VariableHighLimit == limit() 
00568                 &&  x < m_a  ) { m_result *= -1 ; }
00569       else if ( GaudiMath::Integration::VariableLowLimit  == limit() 
00570                 &&  x > m_a  ) { m_result *= -1 ; }
00571       
00572       return m_result ;
00573     };

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

allocate the integration workspace

Definition at line 349 of file NumericalIndefiniteIntegral.cpp.

00350     {
00351       if ( 0 != m_ws ) { return m_ws; }
00352       gsl_integration_workspace* aux = 
00353         gsl_integration_workspace_alloc( size () );
00354       if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; };
00355       m_ws = new _Workspace() ;
00356       m_ws->ws = aux ;
00357       return m_ws ;
00358     };

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

Definition at line 341 of file NumericalIndefiniteIntegral.h.

00342       { return m_ws ; };

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

Definition at line 249 of file NumericalIndefiniteIntegral.cpp.

00251     {
00252       throw GaudiException( "NumericalIndefiniteIntegral::" + message , 
00253                             "*GaudiMath*" , sc ) ;
00254       return sc ;
00255     };

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

Member Data Documentation

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

Definition at line 359 of file NumericalIndefiniteIntegral.h.

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

Definition at line 360 of file NumericalIndefiniteIntegral.h.

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

Definition at line 361 of file NumericalIndefiniteIntegral.h.

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

Definition at line 363 of file NumericalIndefiniteIntegral.h.

GaudiMath::Integration::Limit Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_limit [private]

Definition at line 365 of file NumericalIndefiniteIntegral.h.

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

Definition at line 366 of file NumericalIndefiniteIntegral.h.

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

Definition at line 367 of file NumericalIndefiniteIntegral.h.

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

Definition at line 368 of file NumericalIndefiniteIntegral.h.

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

Definition at line 370 of file NumericalIndefiniteIntegral.h.

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

Definition at line 371 of file NumericalIndefiniteIntegral.h.

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

Definition at line 373 of file NumericalIndefiniteIntegral.h.

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

Definition at line 374 of file NumericalIndefiniteIntegral.h.

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

Definition at line 376 of file NumericalIndefiniteIntegral.h.

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

Definition at line 377 of file NumericalIndefiniteIntegral.h.

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

Definition at line 379 of file NumericalIndefiniteIntegral.h.

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

Definition at line 380 of file NumericalIndefiniteIntegral.h.

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

Definition at line 382 of file NumericalIndefiniteIntegral.h.

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

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

Gaudi Framework, version v21r4

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral 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