I'm trying to numerically calculate some tricky integrals with GSL, and I'd like the implementation to separate out the integrands from the GSL implementation to make my code more readable (and debuggable). These are the functions that actually evaluates my models, what's in them isn't important:
double Integrand1(double *k, size_t dim, void *params) {
double result = k[0]*k[1]*k[2];
return result;
}
double Integrand2(double *k, size_t dim, void *params) {
double result = k[0]-k[1]-k[2];
return result;
}
And this is where I'm running into trouble, trying to abstract the calls to GSL. Note the type constraints for the gsl_monte_function G
:
void EvalIntegral(double &func, int_info ¶ms,double *xl,double *xu,size_t dim,
double &result,double &error){
const gsl_rng_type *T;
gsl_rng *r;
gsl_monte_function G;
G.f = &func; //This must be double (*)(double *,size_t,void *)
G.dim = dim; //This must be size_t
G.params = ¶ms; //This is a struct holding some variables
size_t calls = 5e5;
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc(T);
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(dim);
gsl_monte_vegas_integrate(&G,xl,xu,dim,1e5,r,s, &result,&error);
gsl_monte_vegas_free (s);
gsl_rng_free (r);
}
I set up the above function to reduce some clutter in another function:
void Evaluate(double x, double y){
int_info params;
params.x = x;
params.y = y;
//Lower limits to the variables
double xl[3] = {0.,-100.,0.};
//Upper limits to the variables
double xu[3] = {100.,100.,1.};
double result1, result2, error1, error2;
EvalIntegral(&Integrand1,¶ms,xl,xu,3,&result1,&error1);
EvalIntegral(&Integrand2,¶ms,xl,xu,3,&result2,&error2);
}
When I compile, though, it seems to be refusing my passing the integrand functions:
error: cannot convert ‘double*’ to ‘double (*)(double*, size_t, void*)’ {aka ‘double (*)(double*, long unsigned int, void*)’} in assignment
G.f = &func;
^~~~
I'm a little confused, all of this works when I copy the contents of EvalIntegral
into Evaluate
twice and replace func
with Integrand1
and Integrand2
. How can I convince the compiler that func
is a function with the right prototype?
The compiler gives you the answer: instead of double&
you need a function pointer type:
void EvalIntegral(double (*func)(double*, std::size_t, void*), ...)