Efficient double integration in C with GSL
<p>Consider the double integral</p>

<pre><code>I = int int [(a^k)*b] da db
</code></pre>

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

<p>My code is as follows</p>

<pre><code>#include &lt;stdlib.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;math.h&gt;
#include &lt;gsl/gsl_integration.h&gt;

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

    return pow(a,k)*b;
}
</code></pre>

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

<pre><code>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 = &amp;innerIntegrand;
    F.params = params;

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

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

    // Free the integration workspace
    gsl_integration_workspace_free(giw);

    // Return result
    return result
}
</code></pre>

<p>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</p>

<pre><code>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 = &amp;outerIntegrand;
    F.params = params;

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

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

    // Free the integration workspace
    gsl_integration_workspace_free(giw);


    // Free memory
    free(params);

    // Return result
    return result
}
</code></pre>

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

<p>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 <code>params</code> argument although it then starts to get quite messy.</p>