Gaudi Framework, version v22r0
Home   Generated: 9 Feb 2011

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral Class Reference

The simple class for numerical integrations. More...

#include <NumericalIndefiniteIntegral.h>

Collaboration diagram for Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral:
Collaboration graph
[legend]

List of all members.

Classes

struct  _Function
struct  _Workspace

Public Types

typedef std::vector< double > Points
 typedef for vector of singular points

Public Member Functions

 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)
 From CLHEP/GenericFunctions.
 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

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.

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 
)

From CLHEP/GenericFunctions.

from CLHEP/GenericFunctions

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

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

Standard constructor

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

Definition at line 77 of file NumericalIndefiniteIntegral.cpp.

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

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 127 of file NumericalIndefiniteIntegral.cpp.

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

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 175 of file NumericalIndefiniteIntegral.cpp.

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

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

copy constructor

Definition at line 210 of file NumericalIndefiniteIntegral.cpp.

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

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

destructor

Definition at line 238 of file NumericalIndefiniteIntegral.cpp.

00239     {
00240       if( 0 != m_ws ) 
00241         { 
00242           gsl_integration_workspace_free ( m_ws->ws ) ;
00243           delete m_ws ;
00244           m_ws = 0 ;
00245         }
00246       if ( 0 != m_pdata    ) { delete m_pdata    ; m_pdata    = 0 ; }
00247       if ( 0 != m_function ) { delete m_function ; m_function = 0 ; }
00248     }

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

Member Function Documentation

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

integration limit

Definition at line 296 of file NumericalIndefiniteIntegral.h.

00296 { return         m_a ; }

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

allocate the integration workspace

Definition at line 355 of file NumericalIndefiniteIntegral.cpp.

00356     {
00357       if ( 0 != m_ws ) { return m_ws; }
00358       gsl_integration_workspace* aux = 
00359         gsl_integration_workspace_alloc( size () );
00360       if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; };
00361       m_ws = new _Workspace() ;
00362       m_ws->ws = aux ;
00363       return m_ws ;
00364     }

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

integration category

Definition at line 320 of file NumericalIndefiniteIntegral.h.

00320 { return  m_category ; }

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::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::error (  )  const [inline]

evaluate of previous error

Definition at line 307 of file NumericalIndefiniteIntegral.h.

00307 { return     m_error ; }

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

Definition at line 255 of file NumericalIndefiniteIntegral.cpp.

00257     {
00258       throw GaudiException( "NumericalIndefiniteIntegral::" + message , 
00259                             "*GaudiMath*" , sc ) ;
00260       return sc ;
00261     }

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 ; }

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 ;}

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

integration limit

Definition at line 314 of file NumericalIndefiniteIntegral.h.

00314 { return     m_limit ; }

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

Function value.

evaluate the function

Definition at line 308 of file NumericalIndefiniteIntegral.cpp.

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

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

Function value.

evaluate the function

Definition at line 268 of file NumericalIndefiniteIntegral.cpp.

00269     {
00270       // reset the result and the error  
00271       m_result = GSL_NEGINF ;
00272       m_error  = GSL_POSINF ;
00273       
00274       // check the argument 
00275       if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
00276       
00277       m_argument[0] = argument ;
00278       return (*this) ( m_argument );
00279     }

NumericalIndefiniteIntegral& Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::operator= ( const NumericalIndefiniteIntegral &   )  [private]
Genfun::Derivative Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::partial ( unsigned int  index  )  const [virtual]

Derivatives.

Definition at line 286 of file NumericalIndefiniteIntegral.cpp.

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

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::QAG ( _Function *  fun  )  const [protected]

Definition at line 501 of file NumericalIndefiniteIntegral.cpp.

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

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

Definition at line 370 of file NumericalIndefiniteIntegral.cpp.

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

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

Definition at line 407 of file NumericalIndefiniteIntegral.cpp.

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

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

Definition at line 543 of file NumericalIndefiniteIntegral.cpp.

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

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

Definition at line 463 of file NumericalIndefiniteIntegral.cpp.

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

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

previous result

Definition at line 305 of file NumericalIndefiniteIntegral.h.

00305 { return    m_result ; }

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

integration rule

Definition at line 323 of file NumericalIndefiniteIntegral.h.

00323 { return      m_rule ; }

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

Definition at line 310 of file NumericalIndefiniteIntegral.h.

00310 { return      m_size ; }

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

integration type

Definition at line 317 of file NumericalIndefiniteIntegral.h.

00317 { return      m_type ; }

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

Definition at line 341 of file NumericalIndefiniteIntegral.h.

00342       { return m_ws ; };

Member Data Documentation

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

Definition at line 363 of file NumericalIndefiniteIntegral.h.

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

Definition at line 382 of file NumericalIndefiniteIntegral.h.

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

Definition at line 367 of file NumericalIndefiniteIntegral.h.

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

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

Definition at line 377 of file NumericalIndefiniteIntegral.h.

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

Definition at line 359 of file NumericalIndefiniteIntegral.h.

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

Definition at line 361 of file NumericalIndefiniteIntegral.h.

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

Definition at line 365 of file NumericalIndefiniteIntegral.h.

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

Definition at line 371 of file NumericalIndefiniteIntegral.h.

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

Definition at line 370 of file NumericalIndefiniteIntegral.h.

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

Definition at line 376 of file NumericalIndefiniteIntegral.h.

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

Definition at line 368 of file NumericalIndefiniteIntegral.h.

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

Definition at line 379 of file NumericalIndefiniteIntegral.h.

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

Definition at line 366 of file NumericalIndefiniteIntegral.h.

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

Definition at line 380 of file NumericalIndefiniteIntegral.h.

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines
Generated at Wed Feb 9 16:33:43 2011 for Gaudi Framework, version v22r0 by Doxygen version 1.6.2 written by Dimitri van Heesch, © 1997-2004