Herwig/Utilities/GSLBisection.tcc

0001 // -*- C++ -*-
0002 //
0003 // GSLBisection.tcc is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 //
0010 // This is the implementation of the non-inlined templated member
0011 // functions of the GSLBisection class.
0012 //
0013 using namespace Herwig;
0014 using namespace ThePEG;
0015
0016 namespace {
0017
0018   template<class T> double func(double x , void * p) {
0019     //Units of the argument and return type
0020     const typename T::ValType ValUnit = ((struct GSLparam<T> *)p)->function.vUnit();
0021     const typename T::ArgType ArgUnit = ((struct GSLparam<T> *)p)->function.aUnit();
0022
0023     const T & f = ((struct GSLparam<T> *)p)->function;
0024     return f(x * ArgUnit ) / ValUnit;
0025   }
0026
0027 }
0028
0029 namespace Herwig {
0030 using namespace ThePEG;
0031
0032 template <class T> inline typename T::ArgType
0033 GSLBisection::value(const T & fn,
0034                     const typename T::ArgType lower,
0035                     const typename T::ArgType upper) const {
0036
0037   typedef typename T::ArgType ArgType;
0038   const ArgType ArgUnit = fn.aUnit();
0039
0040   //use own error handler
0041   gsl_error_handler_t *old_handler =
0042     gsl_set_error_handler(& GSLsubstHandler);
0043
0044
0045   int status(0), iter(0);
0046   const gsl_root_fsolver_type *solverType;
0047   gsl_root_fsolver *solver;
0048   double result(0);
0049   double x_lo(lower/ArgUnit), x_hi(upper/ArgUnit);
0050
0051   GSLparam<T> parameters = { fn };
0052   gsl_function F;
0053   F.function = & func<T>;
0054   F.params = &parameters;
0055
0056   solverType = gsl_root_fsolver_brent;
0057   solver = gsl_root_fsolver_alloc (solverType);
0058
0059   try{
0060     gsl_root_fsolver_set (solver, &F, x_lo, x_hi);
0061   }catch(GSLerror){
0062     //cerr << "GSLBisection: initial interval does not contain zero\n";
0063     throw IntervalError();
0064   }
0065
0066 /*
0067   printf ("Root finding is using %s method\n",
0068           gsl_root_fsolver_name (solver));
0069   printf ("%5s [%9s, %9s] %9s %10s\n",
0070           "iter", "lower", "upper", "root",
0071           "err");
0072 */
0073   do{
0074     iter++;
0075     status = gsl_root_fsolver_iterate (solver);
0076     result = gsl_root_fsolver_root (solver);
0077     x_lo = gsl_root_fsolver_x_lower (solver);
0078     x_hi = gsl_root_fsolver_x_upper (solver);
0079     status = gsl_root_test_interval (x_lo, x_hi, abserr_, relerr_);
0080
0081 /*
0082     if (status == GSL_SUCCESS)
0083       printf ("Converged:\n");
0084
0085     printf ("%5d [%.7f, %.7f] %.7f %.7f\n",
0086             iter, x_lo, x_hi,
0087             result, x_hi - x_lo);
0088 */
0089   }
0090   while (status == GSL_CONTINUE && iter < maxPoints_);
0091
0092   gsl_root_fsolver_free (solver);
0093   //use default GSL error handler again
0094   gsl_set_error_handler(old_handler);
0095
0096   //fix units and return
0097   return result * ArgUnit;
0098 }
0099
0100 }