gsl_quadrature_adaptive_gauss_kronrod.h

#ifndef GSL_QUADRATURE_ADAPTIVE_GAUSS_KRONROD_H
#define GSL_QUADRATURE_ADAPTIVE_GAUSS_KRONROD_H
 
#include <gsl/gsl_integration.h>
#include <omp.h>
 
#include <cassert>
 
class GslQuadratureAdaptiveGaussKronrod {
 public:
  using fT = double (*)(double x, void* params);
 
  GslQuadratureAdaptiveGaussKronrod(double epsabs = 0.00, double epsrel = 1e-10)
      : epsabs(epsabs), epsrel(epsrel) {
    thread_num = (size_t)omp_get_max_threads();
    workspace = new gsl_integration_workspace*[thread_num];
    for (size_t i = 0; i < thread_num; ++i) {
      workspace[i] = gsl_integration_workspace_alloc(1000);
    }
  }
 
  ~GslQuadratureAdaptiveGaussKronrod() {
    if (workspace) {
      for (size_t i = 0; i < thread_num; ++i) {
        gsl_integration_workspace_free(workspace[i]);
      }
      delete[] workspace;
    }
  }
 
  inline bool integrate(fT f, void* params, double lowerLimit,
                        double upperLimit, double* result,
                        double* abserr) const {
    gsl_function F;
    F.function = f;
    F.params = params;
    int tid = omp_get_thread_num();
    int status =
        gsl_integration_qag(&F, lowerLimit, upperLimit, epsabs, epsrel, 1000,
                            GSL_INTEG_GAUSS31, workspace[tid], result, abserr);
    return status == GSL_SUCCESS;
  }
 
 private:
  gsl_integration_workspace** workspace;
  size_t thread_num;
  double epsabs;
  double epsrel;
};
 
#endif  // GSL_QUADRATURE_ADAPTIVE_GAUSS_KRONROD_H