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

#include <NumericalDefiniteIntegral.h>

Collaboration diagram for Genfun::GaudiMathImplementation::NumericalDefiniteIntegral:

Classes
struct	_Function
struct	_Workspace
Public Types
typedef std::vector< double >	Points
	typedef for vector of singular points
Public Member Functions
	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)
	From CLHEP/GenericFunctions.
	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

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.

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`
	)

From CLHEP/GenericFunctions.

from CLHEP/GenericFunctions

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

Definition at line 110 of file NumericalDefiniteIntegral.cpp.

      : AbsFunction ()
      , m_function  ( function.clone          ()    )
      , m_DIM       ( 0                             )
      , m_index     ( index                         )
      , m_a         ( a      )
      , m_b         ( b      )
      , m_ia        ( false  )
      , m_ib        ( false  )
      , m_type      ( type   )
      , m_category  ( GaudiMath::Integration::Finite )
      , m_rule      ( rule       )
      , m_points    (            )
      , m_pdata     ( 0          )
      , m_epsabs    ( epsabs     )
      , m_epsrel    ( epsrel     )
      , m_result    ( GSL_NEGINF )
      , m_error     ( GSL_POSINF )
      , m_size      ( size       )
      , m_ws        ()
      , m_argument  ()
      , m_argF      ()
    {
      if ( GaudiMath::Integration::Fixed == m_rule )
        { m_rule = GaudiMath::Integration::Default ; }
      if ( function.dimensionality() < 2          )
        { Exception("::constructor: invalid dimensionality ") ; }
      if ( m_index >= function.dimensionality()   )
        { Exception("::constructor: invalid variable index") ; }

      m_DIM = function.dimensionality() - 1 ;
      m_argument = Argument( m_DIM      ) ;
      m_argF     = Argument( m_DIM + 1  ) ;

    }

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 165 of file NumericalDefiniteIntegral.cpp.

      : AbsFunction ()
      , m_function  ( function.clone()             )
      , m_DIM       ( 0                            )
      , m_index     ( index                        )
      , m_a         ( a                            )
      , m_b         ( b                            )
      , m_ia        ( false  )
      , m_ib        ( false  )
      , m_type      ( GaudiMath::Integration::    Other )
      , m_category  ( GaudiMath::Integration:: Singular )
      , m_rule      ( GaudiMath::Integration::    Fixed )
      , m_points    ( points  )
      , m_pdata     ( 0       )
      , m_epsabs    ( epsabs  )
      , m_epsrel    ( epsrel  )
      //
      , m_result    ( GSL_NEGINF                          )
      , m_error     ( GSL_POSINF                          )
      //
      , m_size      ( size                                )
      , m_ws        ( 0                                   )
      , m_argument  ()
      , m_argF      ()
    {
      if ( function.dimensionality() < 2          )
        { Exception("::constructor: invalid dimensionality ") ; }
      if ( m_index >= function.dimensionality()   )
        { Exception("::constructor: invalid variable index") ; }

      m_DIM = function.dimensionality() - 1 ;
      m_argument = Argument( m_DIM      ) ;
      m_argF     = Argument( m_DIM + 1  ) ;

      const double l1 = std::min ( a , b ) ;
      const double l2 = std::max ( a , b ) ;
      m_points.push_back ( l1 ) ;
      m_points.push_back ( l2 ) ;
      std::sort ( m_points.begin() , m_points.end() ) ;
      m_points.erase( std::unique( m_points.begin () ,
                                   m_points.end   () ) ,
                      m_points.end() );

      Points::iterator lower =
        std::lower_bound ( m_points.begin () , m_points.end () , l1 ) ;
      m_points.erase     ( m_points.begin () , lower                ) ;
      Points::iterator upper =
        std::upper_bound ( m_points.begin () , m_points.end () , l2 ) ;
      m_points.erase     ( upper             , m_points.end ()      ) ;

      m_pdata = new double[ m_points.size() ] ;
      std::copy( m_points.begin() , m_points.end() , m_pdata );
    }

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

Definition at line 252 of file NumericalDefiniteIntegral.cpp.

      : AbsFunction()
      , m_function  ( function.clone()             )
      , m_DIM       ( 0                            )
      , m_index     ( index                        )
      , m_a         ( a                            )
      , m_b         ( GSL_POSINF                   )
      , m_ia        ( false  )
      , m_ib        ( true   )
      , m_type      ( GaudiMath::Integration::    Other )
      , m_category  ( GaudiMath::Integration:: Infinite )
      , m_rule      ( GaudiMath::Integration::    Fixed )
      , m_points    (         )
      , m_pdata     ( 0       )
      , m_epsabs    ( epsabs  )
      , m_epsrel    ( epsrel  )
      //
      , m_result    ( GSL_NEGINF                          )
      , m_error     ( GSL_POSINF                          )
      //
      , m_size      ( size                                )
      , m_ws        ( 0                                   )
      , m_argument  ()
      , m_argF      ()
    {
      if ( function.dimensionality() < 2          )
        { Exception("::constructor: invalid dimensionality ") ; }
      if ( m_index >= function.dimensionality()   )
        { Exception("::constructor: invalid variable index") ; }

      m_DIM = function.dimensionality() - 1 ;
      m_argument = Argument( m_DIM      ) ;
      m_argF     = Argument( m_DIM + 1  ) ;

    }

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

Definition at line 321 of file NumericalDefiniteIntegral.cpp.

      : AbsFunction()
      , m_function  ( function.clone()             )
      , m_DIM       ( 0                            )
      , m_index     ( index                        )
      , m_a         ( GSL_NEGINF                   )
      , m_b         ( b                            )
      , m_ia        ( true   )
      , m_ib        ( false  )
      , m_type      ( GaudiMath::Integration::    Other )
      , m_category  ( GaudiMath::Integration:: Infinite )
      , m_rule      ( GaudiMath::Integration::    Fixed )
      , m_points    (         )
      , m_pdata     ( 0       )
      , m_epsabs    ( epsabs  )
      , m_epsrel    ( epsrel  )
      //
      , m_result    ( GSL_NEGINF                          )
      , m_error     ( GSL_POSINF                          )
      //
      , m_size      ( size                                )
      , m_ws        ( 0                                   )
      , m_argument  ()
      , m_argF      ()
    {
      if ( function.dimensionality() < 2          )
        { Exception("::constructor: invalid dimensionality ") ; }
      if ( m_index >= function.dimensionality()   )
        { Exception("::constructor: invalid variable index") ; }

      m_DIM = function.dimensionality() - 1 ;
      m_argument = Argument( m_DIM      ) ;
      m_argF     = Argument( m_DIM + 1  ) ;
    }

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

Definition at line 386 of file NumericalDefiniteIntegral.cpp.

      : AbsFunction()
      , m_function  ( function.clone()             )
      , m_DIM       ( 0      )
      , m_index     ( index                        )
      , m_a         ( GSL_NEGINF                   )
      , m_b         ( GSL_POSINF                   )
      , m_ia        ( true   )
      , m_ib        ( true   )
      , m_type      ( GaudiMath::Integration::    Other )
      , m_category  ( GaudiMath::Integration:: Infinite )
      , m_rule      ( GaudiMath::Integration::    Fixed )
      , m_points    (         )
      , m_pdata     ( 0       )
      , m_epsabs    ( epsabs  )
      , m_epsrel    ( epsrel  )
      //
      , m_result    ( GSL_NEGINF                          )
      , m_error     ( GSL_POSINF                          )
      //
      , m_size      ( size                                )
      , m_ws        ( 0                                   )
      , m_argument  ()
      , m_argF      ()
    {
      if ( function.dimensionality() < 2          )
        { Exception("::constructor: invalid dimensionality ") ; }
      if ( m_index >= function.dimensionality()   )
        { Exception("::constructor: invalid variable index") ; }

      m_DIM = function.dimensionality() - 1 ;
      m_argument = Argument( m_DIM      ) ;
      m_argF     = Argument( m_DIM + 1  ) ;

    }

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

copy constructor

Definition at line 428 of file NumericalDefiniteIntegral.cpp.

      : AbsFunction ()
      , m_function  ( right.m_function->clone() )
      , m_DIM       ( right.m_DIM      )
      , m_index     ( right.m_index    )
      , m_a         ( right.m_a        )
      , m_b         ( right.m_b        )
      , m_ia        ( right.m_ia       )
      , m_ib        ( right.m_ib       )
      , m_type      ( right.m_type     )
      , m_category  ( right.m_category )
      , m_rule      ( right.m_rule     )
      , m_points    ( right.m_points   )
      , m_pdata     ( 0                )
      , m_epsabs    ( right.m_epsabs   )
      , m_epsrel    ( right.m_epsrel   )
      , m_result    ( GSL_NEGINF       )
      , m_error     ( GSL_POSINF       )
      , m_size      ( right.m_size     )
      , m_ws        ( 0                )
      , m_argument  ( right.m_argument )
      , m_argF      ( right.m_argF     )
    {
      if( 0 != right.m_pdata )
        {
          m_pdata = new double[m_points.size()] ;
          std::copy( m_points.begin() , m_points.end() , m_pdata );
        }
    }

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

destructor

Definition at line 458 of file NumericalDefiniteIntegral.cpp.

    {
      if( 0 != m_function ) { delete   m_function ; m_function = 0 ; }
      if( 0 != m_pdata    ) { delete[] m_pdata    ; m_pdata    = 0 ; }
    }

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

Member Function Documentation

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

integration limit

Definition at line 322 of file NumericalDefiniteIntegral.h.

{ return         m_a ; }

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

allocate the integration workspace

Definition at line 574 of file NumericalDefiniteIntegral.cpp.

    {
      if ( 0 != m_ws ) { return m_ws; }
      gsl_integration_workspace* aux =
        gsl_integration_workspace_alloc( size () );
      if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; };
      m_ws = new _Workspace() ;
      m_ws->ws = aux ;
      return m_ws ;
    }

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

Definition at line 323 of file NumericalDefiniteIntegral.h.

{ return         m_b ; }

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

integration category

Definition at line 344 of file NumericalDefiniteIntegral.h.

{ return  m_category ; }

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

dimensionality of the problem

Definition at line 304 of file NumericalDefiniteIntegral.h.

{ return m_DIM ; }

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

absolute precision

Definition at line 327 of file NumericalDefiniteIntegral.h.

{ return    m_epsabs ; }

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

relatiove precision

Definition at line 329 of file NumericalDefiniteIntegral.h.

{ return    m_epsrel ; }

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

evaluate of previous error

Definition at line 334 of file NumericalDefiniteIntegral.h.

{ return     m_error ; }

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

Definition at line 468 of file NumericalDefiniteIntegral.cpp.

    {
      throw GaudiException( "NumericalDefiniteIntegral::" + message ,
                            "*GaudiMath*" , sc ) ;
      return sc ;
    }

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

accessor to the function itself

Definition at line 320 of file NumericalDefiniteIntegral.h.

{ return *m_function ; }

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

Does this function have an analytic derivative?

Definition at line 312 of file NumericalDefiniteIntegral.h.

{ return true ;}

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

Function value.

evaluate the function

Definition at line 481 of file NumericalDefiniteIntegral.cpp.

    {
      // reset the result and the error
      m_result = GSL_NEGINF ;
      m_error  = GSL_POSINF ;

      // check the argument
      if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };

      m_argument[0] = argument ;
      return (*this) ( m_argument );
    }

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

Function value.

evaluate the function

Definition at line 514 of file NumericalDefiniteIntegral.cpp.

    {
      // reset the result and the error
      m_result = GSL_NEGINF ;
      m_error  = GSL_POSINF ;

      // check the argument
      if( argument.dimension() != m_DIM )
        { Exception ( "operator(): invalid argument size " ) ; };

      // copy the argument

      for( size_t i  = 0 ; i < m_DIM ; ++i )
        {
          m_argument [i] = argument [i] ;
          const size_t j =  i < m_index ? i : i + 1 ;
          m_argF     [j] = argument [i] ;
        };

      // create the helper object
      GSL_Helper helper( *m_function , m_argF , m_index );

      // use GSL to evaluate the numerical derivative
      gsl_function F ;
      F.function = &GSL_Adaptor ;
      F.params   = &helper                 ;
      _Function F1    ;
      F1.fn      = &F ;

      if        (  GaudiMath::Integration::Infinite         == category () )
        { return   QAGI ( &F1 ) ; }                                // RETURN

      if ( m_a == m_b )
        {
          m_result = 0    ; m_error  = 0    ;                      // EXACT
          return m_result ;                                       // RETURN
        }

      if        (  GaudiMath::Integration::Singular         == category () )
         { return   QAGP ( &F1 ) ; }                                // RETURN
      else if   (  GaudiMath::Integration::Finite           == category () )
        if      (  GaudiMath::Integration::NonAdaptive      == type     () )
          { return QNG  ( &F1 ) ; }                                // RETURN
        else if (  GaudiMath::Integration::Adaptive         == type     () )
          { return QAG  ( &F1 ) ; }                                // RETURN
        else if (  GaudiMath::Integration::AdaptiveSingular == type     () )
          { return QAGS ( &F1 ) ; }                                // RETURN
        else
          { Exception ( "::operator(): invalid type "  ); }
      else
        { Exception ( "::operator(): invalid category "  ); }

      return 0 ;
    }

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

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

Derivatives.

Definition at line 499 of file NumericalDefiniteIntegral.cpp.

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

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

known singularities

Definition at line 325 of file NumericalDefiniteIntegral.h.

{ return    m_points ; }

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

Definition at line 692 of file NumericalDefiniteIntegral.cpp.

    {
      if( 0 == F ) { Exception("QAG::invalid function") ; }

      // allocate workspace
      if( 0 == ws () ) { allocate () ; }

      // integration limits
      const double low  = std::min ( m_a , m_b ) ;
      const double high = std::max ( m_a , m_b ) ;

      int ierror =
        gsl_integration_qag ( F->fn                    ,
                              low       ,         high ,
                              m_epsabs  ,     m_epsrel ,
                              size   () , (int) rule() , ws ()->ws ,
                              &m_result ,     &m_error             );

      if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " ,
                                __FILE__ , __LINE__ , ierror ) ; }

      // sign
      if ( m_a > m_b ) { m_result *= -1 ; }

      return m_result ;
    }

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

Definition at line 589 of file NumericalDefiniteIntegral.cpp.

    {
      // check the argument
      if ( 0 == F    ) { Exception("::QAGI: invalid function"); }

      // allocate workspace
      if ( 0 == ws() ) { allocate() ; }

      int ierror = 0 ;

      if ( m_ia && m_ib )
        {
          ierror = gsl_integration_qagi  ( F->fn                ,
                                           m_epsabs  , m_epsrel ,
                                           size ()   , ws()->ws ,
                                           &m_result , &m_error ) ;
        }
      else if ( m_ia )
        {
          ierror = gsl_integration_qagil ( F->fn     , m_b      ,
                                           m_epsabs  , m_epsrel ,
                                           size ()   , ws()->ws ,
                                           &m_result , &m_error ) ;
        }
      else if ( m_ib )
        {
          ierror = gsl_integration_qagiu ( F->fn     , m_a      ,
                                           m_epsabs  , m_epsrel ,
                                           size ()   , ws()->ws ,
                                           &m_result , &m_error ) ;
        }
      else
        { Exception ( "::QAGI: invalid mode" ) ; };

      if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" ,
                                __FILE__ , __LINE__ , ierror ) ;}

      return m_result ;
    }

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

Definition at line 633 of file NumericalDefiniteIntegral.cpp.

    {
      if( 0 == F ) { Exception("QAGP::invalid function") ; }

      // no known singular points ?
      if( points().empty() || 0 == m_pdata ) { return QAGS( F ) ; }

      const size_t npts = points().size();

      // use GSL
      int ierror =
        gsl_integration_qagp ( F->fn                ,
                               m_pdata   , npts     ,
                               m_epsabs  , m_epsrel ,
                               size ()   , ws()->ws ,
                               &m_result , &m_error ) ;

      if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " ,
                                __FILE__ , __LINE__ , ierror ) ; }

      // the sign
      if ( m_a > m_b ) { m_result *= -1 ; }

      return m_result ;
    }

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

Definition at line 723 of file NumericalDefiniteIntegral.cpp.

    {
      if( 0 == F ) { Exception("QAG::invalid function") ; }

      if ( m_a == m_b )
        {
          m_result = 0    ;
          m_error  = 0    ;   // EXACT !
          return m_result ;
        }

      // allocate workspace
      if( 0 == ws () ) { allocate () ; }

      // integration limits
      const double low  = std::min ( m_a , m_b ) ;
      const double high = std::max ( m_a , m_b ) ;

      int ierror =
        gsl_integration_qags ( F->fn                ,
                               low       , high     ,
                               m_epsabs  , m_epsrel ,
                               size   () , ws()->ws ,
                               &m_result , &m_error );

      if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " ,
                                __FILE__ , __LINE__ , ierror ) ; }

      // sign
      if ( m_a > m_b ) { m_result *= -1 ; }

      return m_result ;
    }

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

Definition at line 663 of file NumericalDefiniteIntegral.cpp.

    {
      if( 0 == F ) { Exception("QNG::invalid function") ; }

      // integration limits
      const double low  = std::min ( m_a , m_b ) ;
      const double high = std::max ( m_a , m_b ) ;

      size_t neval = 0 ;
      int ierror =
        gsl_integration_qng ( F->fn                 ,
                              low       ,      high ,
                              m_epsabs  ,  m_epsrel ,
                              &m_result , &m_error  , &neval  ) ;

      if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
                                __FILE__ , __LINE__ , ierror ) ; }

      // sign
      if ( m_a > m_b ) { m_result *= -1 ; }

      return m_result ;
    }

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

previous result

Definition at line 332 of file NumericalDefiniteIntegral.h.

{ return    m_result ; }

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

integration rule

Definition at line 347 of file NumericalDefiniteIntegral.h.

{ return      m_rule ; }

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

Definition at line 337 of file NumericalDefiniteIntegral.h.

{ return      m_size ; }

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

integration type

Definition at line 341 of file NumericalDefiniteIntegral.h.

{ return      m_type ; }

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

Definition at line 365 of file NumericalDefiniteIntegral.h.

      { return m_ws ; };

Member Data Documentation

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

Definition at line 386 of file NumericalDefiniteIntegral.h.

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

Definition at line 408 of file NumericalDefiniteIntegral.h.

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

Definition at line 407 of file NumericalDefiniteIntegral.h.

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

Definition at line 387 of file NumericalDefiniteIntegral.h.

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

Definition at line 392 of file NumericalDefiniteIntegral.h.

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

Definition at line 383 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_error [mutable, private]

Definition at line 402 of file NumericalDefiniteIntegral.h.

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

Definition at line 382 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.

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

Definition at line 384 of file NumericalDefiniteIntegral.h.

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

Definition at line 396 of file NumericalDefiniteIntegral.h.

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

Definition at line 395 of file NumericalDefiniteIntegral.h.

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

Definition at line 401 of file NumericalDefiniteIntegral.h.

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

Definition at line 393 of file NumericalDefiniteIntegral.h.

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

Definition at line 404 of file NumericalDefiniteIntegral.h.

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

Definition at line 391 of file NumericalDefiniteIntegral.h.

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

Definition at line 405 of file NumericalDefiniteIntegral.h.

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

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

Gaudi Framework, version v23r4

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral Class Reference

Classes

Public Types

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation