Calling GSL routine CQUAD from Fortran

Question

I'm trying to call the GSL rountine CQUAD from Fortran. My idea was to write a .c subroutine that calls the gsl rountine and depends on a function and bounds. Two problems: I have only very little idea about c and fortrans iso_c_binding. My attempt is as follows:

A simple calling program (similar to M.S.B's post in No output from a Fortran program using the Gnu Scientific Library via a c wrapper ):

program test

   use iso_c_binding

   interface
      function  my_cquad (f,a,b)  bind(c)
         import
         real (kind=c_double) :: my_cquad

         interface
            function f(x) bind(c)
               import
               real(kind=c_double) :: f,x 
            end function
         end interface

         real (kind=c_double) :: a,b
      end function my_cquad
   end interface

   real (kind=c_double) :: y,a,b

   a=0. ; b=1.
   y=my_cquad(g,a,b)
   print*,y

   stop

contains

   function g(x) bind(C)
      real(kind=c_double) :: g,x
      g=sin(x)/x
      return
   end function g

end program test

The .c subroutine (basically taken from the example given by the author of CQUAD in https://scicomp.stackexchange.com/questions/4730/numerical-integration-handling-nans-c-fortran):

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

double my_cquad ( double my_f() , double a , double b  )
{
    gsl_function f;
    gsl_integration_cquad_workspace *ws = NULL;
    double res, abserr;
    size_t neval;

    /* Prepare the function. */
    f.function = my_f;
    f.params = NULL;

    /* Initialize the workspace. */
    if ( ( ws = gsl_integration_cquad_workspace_alloc( 200 ) ) == NULL ) {
        printf( "call to gsl_integration_cquad_workspace_alloc failed.\n" );
        abort();
        }

    /* Call the integrator. */
    if ( gsl_integration_cquad( &f, a , b , 1.0e-10 , 1.0e-10 , ws , &res , &abserr , &neval ) != 0 ) {
        printf( "call to gsl_integration_cquad failed.\n" );
        abort();
        }

    /* Free the workspace. */
    gsl_integration_cquad_workspace_free( ws );

    /* Bye. */
    return res;
}

The .c subroutine alone seems to work fine. This can be tested with:

double g (double x)
{
return sin(x)/x;
}

int main () {
double y;
y=my_cquad(g,0.,1.);
printf("y:    %2.18f\n", y);

return 0;
}

But together with the .f90 calling program, at the moment it compiles but at runtime I get a segmentation fault that I don't quite get.

Additionally, it would of course be good to have some kind of wrapper that creates a c-type function depending on a fortran type function. I'm thinking about something like:

function f_to_c(f) bind(c)
   real(kind=c_double) :: f_to_c
   real(kind=8) :: f
   f_to_c=real(f,kind=c_double)
end function

But this desn't cover the dummy variables.

Thanks in advance and very sorry for the amount of code.

Answer 1

Beware, according to the Fortran standard internal functions shall not have the bind(C) attribute. I moved the function to a module.

The a and b must be passed by value to my_cquad and x to the integrated function:

module functions_to_integrate

   use iso_c_binding

contains
   function g(x) bind(C)
      real(kind=c_double) :: g
      real(kind=c_double), value :: x
      g = sin(x)/x
   end function g
end module

program test

   use iso_c_binding

   use functions_to_integrate

   interface
      function  my_cquad (f,a,b)  bind(c)
         import
         real (kind=c_double) :: my_cquad

         interface
            function f(x) bind(c)
               import
               real(kind=c_double) :: f
               real(kind=c_double), value :: x 
            end function
         end interface

         real (kind=c_double), value :: a,b
      end function my_cquad
   end interface

   real (kind=c_double) :: y,a,b

   a = 0 ; b = 1
   y = my_cquad(g,a,b)

   print *,y

end program test

test:

> gfortran my_cquad.c test_cquad.f90 -lgsl -lopenblas
> ./a.out 
  0.94608307036718275  

This is an example of a wrapper to a normal Fortran function (please do not use kind=8 for many reasons explained in another questions):

module functions_to_integrate
   use iso_fortran_env
   use iso_c_binding


   integer, parameter :: wp = real64
contains

    pure function g(x)
      real(kind=wp) :: g
      real(kind=wp), intent(in) :: x
      g = sin(x)/x
    end function g

    function g_to_c(x) bind(C)
      real(kind=c_double) :: g_to_c
      real(kind=c_double), value :: x
      g_to_c = real(g(x),kind=c_double)
    end function

end module

program test

    use iso_c_binding

    use functions_to_integrate

    interface
      function  my_cquad (f,a,b)  bind(c)
          import
          real (kind=c_double) :: my_cquad

          interface
            function f(x) bind(c)
                import
                real(kind=c_double) :: f
                real(kind=c_double), value :: x 
            end function
          end interface

          real (kind=c_double), value :: a,b
      end function my_cquad
    end interface

    real (kind=c_double) :: y,a,b

    a = 0 ; b = 1
    y = my_cquad(g_to_c,a,b)

    print *,y

end program test

P.S. I also deleted your stop and return statement before the end. Somehow it is always driving me mad, but that may be just my OCD. I was too used to see it in old programs coming from ancient times.

P.P.S: You may wish to see the FGSL interface package linked by Vincius Miranda http://www.lrz.de/services/software/mathematik/gsl/fortran/ . I knew about that one, but I tried mainly to point out the errors so that you can make similar interfaces yourself, where no ready-made package is available.