From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me X-Spam-Level: X-Spam-Status: No, score=-1.9 required=5.0 tests=BAYES_00 autolearn=unavailable autolearn_force=no version=3.4.4 X-Received: by 10.99.115.66 with SMTP id d2mr11151964pgn.144.1490750960660; Tue, 28 Mar 2017 18:29:20 -0700 (PDT) X-Received: by 10.157.11.150 with SMTP id 22mr1256673oth.18.1490750960611; Tue, 28 Mar 2017 18:29:20 -0700 (PDT) Path: eternal-september.org!reader01.eternal-september.org!reader02.eternal-september.org!news.eternal-september.org!news.eternal-september.org!feeder.eternal-september.org!news.glorb.com!y18no1327343itc.0!news-out.google.com!m191ni13359itc.0!nntp.google.com!y18no1327340itc.0!postnews.google.com!glegroupsg2000goo.googlegroups.com!not-for-mail Newsgroups: comp.lang.ada Date: Tue, 28 Mar 2017 18:29:20 -0700 (PDT) Complaints-To: groups-abuse@google.com Injection-Info: glegroupsg2000goo.googlegroups.com; posting-host=118.209.150.170; posting-account=l8GBMwoAAADCbqdOJSbg4dBRqkD14dJd NNTP-Posting-Host: 118.209.150.170 User-Agent: G2/1.0 MIME-Version: 1.0 Message-ID: <2b035819-e55e-4805-9ff1-a2cec09f13e0@googlegroups.com> Subject: Example Ada calling Gnu Scientific Library (GSL) From: Leo Brewin Injection-Date: Wed, 29 Mar 2017 01:29:20 +0000 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: quoted-printable Xref: news.eternal-september.org comp.lang.ada:46497 Date: 2017-03-28T18:29:20-07:00 List-Id: Hi Folks, I've been following the discussion about the thin bindings to the GSL libra= ry (and thanks to all for the useful information). I decided to follow the suggestion by Per Sandberg (Mar 19) to use gcc to c= reate the header files. After a small amount of fiddling I got things worki= ng. 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 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 th= e 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 :=3D 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 :=3D 0.0; -- lower limit b :=3D 1.0; -- upper limit epsabs :=3D 1.0e-25; -- target actual error epsrel :=3D 1.0e-13; -- target relative error, 1.0e-14 gives a GSL e= rror limit :=3D 1000; -- maximum number of sub-intervals key :=3D 5; -- which integration method to use, range from = 1 to 6 integrand.c_function :=3D fun.f'access; -- the function to integrate integrand.params :=3D alpha'address; -- parameters for the functio= n workspace :=3D gsl_gsl_integration_h.gsl_integration_workspace_alloc (li= mit); return_code :=3D gsl_gsl_integration_h.gsl_integration_qag (integrand'access,a,b,epsabs,epsrel,limit,key,workspace= ,result'access,abserr'access); put ("Intgeral (qag) =3D "); put (Interfaces.C.double'Image(result)); new_line; return_code :=3D gsl_gsl_integration_h.gsl_integration_qags (integrand'access,a,b,epsabs,epsrel,limit,workspace,res= ult'access,abserr'access); put ("Intgeral (qags) =3D "); put (Interfaces.C.double'Image(result)); new_line; size :=3D workspace.size; level :=3D workspace.level; put ("Workspace size =3D "); put (sys_utypes_usize_t_h.size_t'Image(size)); new_line; put ("Workspace level =3D "); 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 Int= erfaces.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 (Inte= rfaces.C.double); use GSL_Maths; function f (x : Interfaces.C.double; params : System.Address) return Int= erfaces.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 h= as 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) =3D -1.22741127776022E+00 Intgeral (qags) =3D -1.22741127776022E+00 Workspace size =3D 12 Workspace level =3D 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 crashe= s with a PROGRAM_ERROR. I tried to catch this with an exception block (arou= nd 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 docume= ntation states that even when the code fails, it returns the best answer. S= o I'd rather have an answer even if it didn't meet my too stringent error r= equest. One further note. On experiments with the GSL Bessel library I found that t= he GSL code uses the case sensitive nature of C to provide different versio= ns of the Bessel functions. This is a problem for Ada which is cases insens= tive. The solution is to add some unique descriptor to selected function na= mes in the Ada headers to avoid the ambiguities. Cheers, Leo