Search code examples
c++gsl

Passing function by reference problem with typing


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 &params,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 = &params; //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,&params,xl,xu,3,&result1,&error1);
    EvalIntegral(&Integrand2,&params,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?


Solution

  • The compiler gives you the answer: instead of double& you need a function pointer type:

    void EvalIntegral(double (*func)(double*, std::size_t, void*), ...)