Fortran, GSL and MPI

As a minor side-issue in my MPI course, this page contains a few notes on calling a GSL routine from Fortran.

In C/C++, function prototypes etc exist in header files. These are compiler-independent, and rarely read by Humans. Fortran modules are a similar concept, but not quite the same. A Fortran module file is a compiler-dependent binary file, perhaps even compiler version dependent, and as such it is unlikely to be provided with a third-party library. It can be generated from a compiler-independent Human-readable source file containing an interface block.

In simple cases, these are simple, but the example routine here is not entirely simple, requiring a structure to be defined, and a function to be passed as an argument. The Fortran is, in my opinion, rather verbose, but one must not overlook the fact that the corresponding C is a similar number of lines. The relevant lines from the C include files are:

    GSL_INTEG_GAUSS15 = 1,      /* 15 point Gauss-Kronrod rule */

struct gsl_function_struct 
  double (* function) (double x, void * params);
  void * params;

typedef struct gsl_function_struct gsl_function ;

typedef struct
    size_t limit;
    size_t size;
    size_t nrmax;
    size_t i;
    size_t maximum_level;
    double *alist;
    double *blist;
    double *rlist;
    double *elist;
    size_t *order;
    size_t *level;

gsl_integration_workspace *
  gsl_integration_workspace_alloc (const size_t n);

  gsl_integration_workspace_free (gsl_integration_workspace * w);

int gsl_integration_qag (const gsl_function * f,
                         double a, double b,
                         double epsabs, double epsrel, size_t limit,
                         int key,
                         gsl_integration_workspace * workspace,
                         double *result, double *abserr);

Our Fortran is slightly lazier, and uses a pointer to an unspecified type in place of one to the defined gsl_integration_workspace structure. So the Fortran translation of the above became

module gsl_int
  use, intrinsic :: iso_c_binding
! Value taken from /usr/include/gsl/gsl_integration.h
  integer, parameter :: GSL_INTEG_GAUSS15 = 1

  type, bind(c) :: gsl_function
     type (c_funptr) :: func
     type (c_ptr) :: params
  end type gsl_function

     function gsl_integration_workspace_alloc(wlen) bind(c)
       use, intrinsic :: iso_c_binding
       integer(c_size_t),value :: wlen
       type(c_ptr) :: gsl_integration_workspace_alloc
     end function gsl_integration_workspace_alloc
  end interface
     subroutine gsl_integration_workspace_free(work) bind(c)
       use, intrinsic :: iso_c_binding
       type(c_ptr),value :: work
     end subroutine gsl_integration_workspace_free
  end interface

     function gsl_integration_qag(F,xmin,xmax,epsabs,epsrel,limit, &
       key,work,out,abserr) bind(c)
       use, intrinsic :: iso_c_binding
       import :: gsl_function
       real(c_double), value :: xmin,xmax,epsabs,epsrel
       real(c_double) :: out,abserr
       integer(c_size_t), value :: limit
       integer(c_int),value :: key
       integer(c_int) :: gsl_integration_qag
       type(c_ptr),value :: work
       type(gsl_function) :: F
     end function gsl_integration_qag
  end interface
end module gsl_int

And finally the serial example, gsl.f90, and the example with MPI, gsl_mpi.f90. The serial version should compile as follows:

  $ gfortran gsl.f90 -lgsl
  $ ./a.out
   Integral is   0.88622692545275761