No output from a Fortran program using the Gnu Scientific Library via a c wrapper

Question

I'm trying to write a wrapper to use the gsl library with Fortran. I have managed to get a simple wrapper to work - the example from http://www.helsinki.fi/~fyl_tlpk/luento/ohj-13-GSL-e.html

Fortran code

program gsltest
    implicit none

    real(kind=selected_real_kind(12)) :: a = 0.11, res
    external :: glsgateway

    call gslgateway(a,res)
    write(*,*) 'x', a, 'atanh(x)', res

end program gsltest

c function

#include <gsl/gsl_math.h>

void gslgateway_(double *x, double *res){
   *res = gsl_atanh(*x);
}

That's all well and good. However, I'm having problems with a more complicated wrapper. I have the following code modified from an example at http://apwillis.staff.shef.ac.uk/aco/freesoftware.html

c wrapper (rng_initialise.c)

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

static gsl_rng* r;

void rng_initialise__(int* s) {
   r = gsl_rng_alloc(gsl_rng_taus); 
   gsl_rng_set(r, (unsigned long int)(*s)); 
}

Fortran main (main.f90)

PROGRAM main
    integer seed

    call system_clock(seed)
    WRITE (*,*) 'calling rng_initialise'
    call rng_initialise(seed)

END PROGRAM main

which I then compile and link by

gcc -c rng_initialise.c
g95 -c main.f90
g95 -o main main.o rng_initialise.o -L/usr/libs -lgsl

When I run this program, I get no output. However, if I comment out the lines inside rng_initialise

...
void rng_initialise__(int* s) {
   // r = gsl_rng_alloc(gsl_rng_taus);  
   // gsl_rng_set(r, (unsigned long int)(*s)); 
}

then I get output from the Fortran code (it writes 'calling_rng_initialise' to STDOUT).

So, the problem seems to be the calls to gsl_rng_alloc and gsl_rng_set. But I don't get any error messages, and I don't know why they would prevent the Fortran code from doing anything. Any ideas?

Answer 1

As already suggested, the best way to do this is to use the Fortran ISO C Binding because it will instruct Fortran to use the C calling conventions to match the C routines of the GSL library. The ISO C Binding is part of the Fortran 2003 standard and has been available in many Fortran 95 compilers for several years. As part of the standard it makes interfacing Fortran and C, in both directions, portable and much easier than the OS and compiler dependent hacks that used to be necessary. I recommend ignoring instructions that pre-date the ISO C Binding as obsolete.

Generally you won't need to write any C wrapper code to call the GSL library, only Fortran specifications statements to describe the C routine interfaces to Fortran in Fortran syntax. Here is a simple example that calls the GSL routine gsl_cdf_chisq_Q

program cum_chisq_prob

   use iso_c_binding

   interface GSL_CummulativeChiSq_Prob_Upper   

      function gsl_cdf_chisq_Q  ( x, nu )  bind ( C, name="gsl_cdf_chisq_Q" )

         import

         real (kind=c_double) :: gsl_cdf_chisq_Q

         real (kind=c_double), VALUE, intent (in) :: x
         real (kind=c_double), VALUE, intent (in) :: nu

      end function gsl_cdf_chisq_Q

   end interface GSL_CummulativeChiSq_Prob_Upper

   real (kind=c_double) :: chisq_value
   real (kind=c_double) :: DoF
   real (kind=c_double) :: Upper_Tail_Prob    

   write ( *, '( / "Calculates cumulative upper-tail probability for Chi-Square distribution." )' )
   write ( *, '( "Input Chisq Value, Degrees of Freedom: " )', advance='no' )
   read ( *, * ) chisq_value, DoF

   Upper_Tail_Prob = gsl_cdf_chisq_Q  ( chisq_value, DoF )

   write ( *, '( "Probability is:", 1PG17.10 )' )  Upper_Tail_Prob

   stop

end program cum_chisq_prob

Even easier: you can find a pre-written library to allow you to call GSL from Fortran at http://www.lrz.de/services/software/mathematik/gsl/fortran/

Answer 2

Most likely you have the linkage between the two routines wrong in some way. If the stack isn't dealt with correctly when you hop through that interface, dang near anything can happen.

I'm not noticing any code on either the Fortran side or the C side specifying the other's calling convention. I'm not an expert with Gnu Fortran, but I know most compilers will require some kind of note that they should be using another compiler's calling convention, or Bad Things may happen.

With just a little web searching, I see that the G95 Fortran manual (PDF) has a nice long section titled "Interfacing with G95 Programs", that appears to go into this in detail. Just from skimming, it looks like you should be using the BIND(C) attribute on your Fortran function declaration for that C routine.

Answer 3

The problem is with the static gsl_rng* r; defined in your C file. But I do not know/understand why the standard does not allow this. After studying the source file of the fgsl package a little bit, I found a tweak that works. The fortran file random_f.f90

module fgsl
    use, intrinsic :: iso_c_binding
    implicit none

    type, bind(C) :: fgsl_rng_type
        type(c_ptr) :: gsl_rng_type_ptr = c_null_ptr
    end type fgsl_rng_type

    type, bind(C) :: fgsl_rng
        type(c_ptr) :: gsl_rng_ptr = c_null_ptr
    end type fgsl_rng
end module fgsl

PROGRAM call_fgsl_rndm
    use, intrinsic :: iso_c_binding
    use fgsl
    implicit none

    interface
        subroutine srndm(seed, t, r) bind(C)
            import
            integer(C_INT) :: seed
            type(fgsl_rng) :: r
            type(fgsl_rng_type) :: t
        end subroutine srndm

        function rndm(r) bind(C)
            import
            real(C_DOUBLE) :: rndm
            type(fgsl_rng) :: r
        end function rndm
    end interface

    type(fgsl_rng) :: r
    type(fgsl_rng_type) :: t

    integer(C_INT) :: seed
    real(C_DOUBLE) :: xi
    seed = 1
    call srndm(seed, t, r)
    xi = rndm(r)
    print *, xi
    xi = rndm(r)
    print *, xi
    xi = rndm(r)
    print *, xi
END PROGRAM

and the C file random_c.c

#include <gsl/gsl_rng.h>

typedef struct{
    gsl_rng *gsl_rng_ptr;
} fgsl_rng;

typedef struct{
    gsl_rng_type *gsl_rng_type_ptr;
} fgsl_rng_type;

void srndm(int *seed, fgsl_rng_type* t, fgsl_rng* r) {
    t->gsl_rng_type_ptr = (gsl_rng_type *) gsl_rng_mt19937;  // cast to remove the const qualifier
    r->gsl_rng_ptr = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(r->gsl_rng_ptr, *seed);
}

double rndm(fgsl_rng* r) {
    return gsl_rng_uniform(r->gsl_rng_ptr);
}

Although only the pointers in the structures are used, the introduction of fgsl_rng and fgsl_rng_type is necessary. Otherwise, the program will fail. Unfortunately, I still have no clear idea why it has to work this way.