title: Efficient double integration in C with GSL
body: <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>