Efficient double integration in C with GSL

Question

Consider the double integral

I = int int [(a^k)*b] da db

where we want to integrate for a between [0,1] and b between [0,1] and k is some constant. I am using the GSL numerical integration library but have a memory allocation issue.

My code is as follows

#include <stdlib.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double innerIntegrand(double a, void *params) {
    double *cast_params = (double *) params;
    double b = params[0];
    double k = params[1];

    return pow(a,k)*b;
}

I can then evaluate the inner integral for a given b (to get an outer integrand) as follows

double outerIntegrand(double b, void *params) {
    // params = {holder for double b, k}
    double *cast_params = (double *) params;
    cast_params[0] = b;

    // Allocate integration workspace
    gsl_integration_workspace *giw = gsl_integration_workspace_alloc(100);

    // Create GSL function
    gsl_function F;
    F.function = &innerIntegrand;
    F.params = params;

    // Initialise values to put the result in
    double result;
    double abserror;

    // Perform integration
    gsl_integration_qag(&F, 0, 1, 0.001, 0.001, 100, 1, giw, &result, &abserror);

    // Free the integration workspace
    gsl_integration_workspace_free(giw);

    // Return result
    return result
}

Note however I have to allocate and free the integration workspace within the function. This means it is done many times when evaluating the final integration function

double Integral(double k) {
    // Create params
    double *params = malloc(2*sizeof(double));
    params[1] = k;

    // Allocate integration workspace
    gsl_integration_workspace *giw = gsl_integration_workspace_alloc(100);

    // Create GSL function
    gsl_function F;
    F.function = &outerIntegrand;
    F.params = params;

    // Initialise values to put the result in
    double result;
    double abserror;

    // Perform integration
    gsl_integration_qag(&F, 0, 1, 0.001, 0.001, 100, 1, giw, &result, &abserror);

    // Free the integration workspace
    gsl_integration_workspace_free(giw);


    // Free memory
    free(params);

    // Return result
    return result
}

Ideally what I want is two global gsl_integration_workspace variables, one for the integral in outerIntegrand and another for the integral in Integral. However when I try to declare them as global values I receive a initializer element is not constant error.

Can anyone see a way to do this double integral without the repeated memory allocation and freeing? I was thinking we could also pass the workspace in through the params argument although it then starts to get quite messy.

Answer 1

I managed to build a decently looking program in C++ for double integration based on GSL, avoiding repeated allocations in a clean way. I used this well known function to play:

f(x,y)=exp(-x*x-y*y)

integrating it over all the plane (the result, pi, can easily be obtained by switching to polar coordinates). It is trivial to modify it and add parameters by lambda capture.

#include <iostream>
#include <gsl/gsl_integration.h>

// Simple RAII wrapper 
class IntegrationWorkspace {
  gsl_integration_workspace * wsp;

  public:
  IntegrationWorkspace(const size_t n=1000):
    wsp(gsl_integration_workspace_alloc(n)) {}
  ~IntegrationWorkspace() { gsl_integration_workspace_free(wsp); }

  operator gsl_integration_workspace*() { return wsp; }
};

// Build gsl_function from lambda
template <typename F>
class gsl_function_pp: public gsl_function {
  const F func;
  static double invoke(double x, void *params) {
    return static_cast<gsl_function_pp*>(params)->func(x);
  }
  public:
  gsl_function_pp(const F& f) : func(f) {
    function = &gsl_function_pp::invoke; //inherited from gsl_function
    params   = this;                     //inherited from gsl_function
  }
  operator gsl_function*(){return this;}
};

// Helper function for template construction
template <typename F>
gsl_function_pp<F> make_gsl_function(const F& func) {
  return gsl_function_pp<F>(func);
}

int main() {
  double epsabs = 1e-8;
  double epsrel = 1e-8;
  size_t limit = 100;
  double result, abserr, inner_result, inner_abserr;

  IntegrationWorkspace wsp1(limit);
  IntegrationWorkspace wsp2(limit);

  auto outer = make_gsl_function( [&](double x) {
    auto inner = make_gsl_function( [&](double y) {return exp(-x*x-y*y);} );
    gsl_integration_qagi(inner, epsabs, epsrel, limit, wsp1,
                         &inner_result, &inner_abserr);
    return inner_result;
  } );
  gsl_integration_qagi(outer, epsabs, epsrel, limit, wsp2, &result, &abserr);

  std::cout << result << std::endl;
}