Example Ada calling Gnu Scientific Library (GSL)

[Leo Brewin <leo.brewin@internode.on.net>]

Hi Folks,

I've been following the discussion about the thin bindings to the GSL library (and thanks to all for the useful information).

I decided to follow the suggestion by Per Sandberg (Mar 19) to use gcc to create the header files. After a small amount of fiddling I got things working. And I'm very impressed with just how painless the whole process was. So here is a short summary of what I did.

1. Created a simple C-file, "integral.c"

#include <gsl/gsl_integration.h>
int main (void) {}

2. Run "gcc -fdump-ada-spec integral.c"
3. Replace "limited with" by "with" on line 6 of gsl_gsl_integration_h.ads
4. Used the generated headers to glean the required Ada syntax to create the
   thin binding needed for a simple GSL example. Here is the main program

with Ada.Text_IO;                 use Ada.Text_IO;
with Interfaces.C;                use Interfaces.C;

with gsl_gsl_math_h;
with gsl_gsl_integration_h;
with sys_utypes_usize_t_h;

with fun; -- function to integerate, must be in a separate pacakge

procedure GSL_Test is

   a,b         : Interfaces.C.double;
   epsabs      : Interfaces.C.double;
   epsrel      : Interfaces.C.double;
   abserr      : aliased Interfaces.C.double;
   result      : aliased Interfaces.C.double;
   limit       : sys_utypes_usize_t_h.size_t;
   key         : Interfaces.C.int;
   return_code : Interfaces.C.int;

   integrand   : aliased gsl_gsl_math_h.gsl_function;
   alpha       : Interfaces.C.double := 4.0;

   workspace   : access  gsl_gsl_integration_h.gsl_integration_workspace;
   size        : aliased sys_utypes_usize_t_h.size_t;
   level       : access  sys_utypes_usize_t_h.size_t;

begin

   a := 0.0;              -- lower limit
   b := 1.0;              -- upper limit

   epsabs := 1.0e-25;     -- target actual error
   epsrel := 1.0e-13;     -- target relative error, 1.0e-14 gives a GSL error
   limit  := 1000;        -- maximum number of sub-intervals
   key    := 5;           -- which integration method to use, range from 1 to 6

   integrand.c_function := fun.f'access;    -- the function to integrate
   integrand.params     := alpha'address;   -- parameters for the function

   workspace := gsl_gsl_integration_h.gsl_integration_workspace_alloc (limit);

   return_code := gsl_gsl_integration_h.gsl_integration_qag
                    (integrand'access,a,b,epsabs,epsrel,limit,key,workspace,result'access,abserr'access);

   put ("Intgeral (qag)  = ");
   put (Interfaces.C.double'Image(result));
   new_line;

   return_code := gsl_gsl_integration_h.gsl_integration_qags
                    (integrand'access,a,b,epsabs,epsrel,limit,workspace,result'access,abserr'access);

   put ("Intgeral (qags) = ");
   put (Interfaces.C.double'Image(result));
   new_line;

   size  := workspace.size;
   level := workspace.level;

   put ("Workspace size  = ");
   put (sys_utypes_usize_t_h.size_t'Image(size));
   new_line;

   put ("Workspace level = ");
   put (sys_utypes_usize_t_h.size_t'Image(level.all));
   new_line;

   gsl_gsl_integration_h.gsl_integration_workspace_free (workspace);

end GSL_Test;

5. The integrand is in a separate package, here is the spec file (fun.ads)

with Interfaces.C; use Interfaces.C;
with System;

package fun is

   function f (x : Interfaces.C.double; params : System.Address) return Interfaces.C.double;
      pragma Convention (C, f);

end fun;

6. And here is the body, fun.adb

with Ada.Numerics.Generic_Elementary_Functions;

package body fun is

   package GSL_Maths is new Ada.Numerics.Generic_Elementary_Functions (Interfaces.C.double); use GSL_Maths;

   function f (x : Interfaces.C.double; params : System.Address) return Interfaces.C.double is
      alpha : Interfaces.C.double;
      for alpha'address use params;
   begin
      return log(alpha*x) / sqrt(x);
   end f;

end fun;

7. The code can be compiled using the "-lgsl" linker option (assuming GSL has been installed
   in the standard locations. I used "homebrew" on my macOS, "brew install gsl").

8. The output from the above code was

Intgeral (qag)  = -1.22741127776022E+00
Intgeral (qags) = -1.22741127776022E+00
Workspace size  =  12
Workspace level =  11

I hope this helps (sorry for the long code listing).

But I do have a question. If I chose the target relative error (epsrel) to be 1.0e-14 the GSL code fails to meet that accuracy and the Ada code crashes with a PROGRAM_ERROR. I tried to catch this with an exception block (around the call to the gsl code) but that didn't work. The exception is raised within the GSL library -- is there any way I can catch this? The GSL documentation states that even when the code fails, it returns the best answer. So I'd rather have an answer even if it didn't meet my too stringent error request.

One further note. On experiments with the GSL Bessel library I found that the GSL code uses the case sensitive nature of C to provide different versions of the Bessel functions. This is a problem for Ada which is cases insenstive. The solution is to add some unique descriptor to selected function names in the Ada headers to avoid the ambiguities.

Cheers,
Leo