Mixed C/Fortran Programming

Introduction/Background

Assume we need to write fortran and wish to call C code that has already been written and cannot be changed (e.g. a precompiled library). There are different levels of complexity possible depending on what is needed. Eg:

  1. Fortran calls a C function that does not take any arguments
  2. Fortran calls a C function with typical variable arguments
  3. Fortran calls a C function with pointer arguments
Here, from fortran, I will call a C function that takes in real-type and function pointer as arguments. The example I use here to describe the programming is to approximate an integral by Gauss-Legendre quadrature. The main program and the function to be integrated are in fortran, while the procedure/subprogram that actually does the quadrature calculation is in C and cannot be changed.

I have picked a function that gives us \(\pi\) on integration. Out of all such integrals, one with definite limits is required. Also, absence of singularities in the integrand is desirable. One such integral is: $$ \int_{0}^{1} \frac{4}{1+x^2} dx=\pi$$ The integrand is well behaved in the region of interest and will do the job for now.

There are, of course, more superior and efficient ways to calculate \(\pi\) numerically/computationally. For instance, the Ramanujan-type formulas or the "spigot" algorithms like the Bailey–Borwein–Plouffe formula or the Gauss-Legendre algorithm (different from the Gauss-Legendre quadrature), some of which have been used for getting a world record number of digits of \(\pi\).

C Code

Obviously, we need to know how the function is to be called. Thus, we require the name of the function, its return type, the list of arguments and the corresponding data types. The function prototype for the example here is given below. I have taken this C code from the wonderful www.mymathlib.com website.


/*a=lower limit
  b=upper limit
  f=pointer to function to be integrated
*/
long double
  Gauss_Legendre_Integration_10pts( long double a, long double b, long double (*f)(long double) )
{...
}
Now we know that the function is named 'Gauss_Legendre_Integration_10pts' and returns a 'long double'. And we can also gather what arguments it expects. This is all the information we need to call this function from a fortran program. Essentially we are calling a blackbox function.

Fortran Code

The whole fortran program is given below with discussion following.


program Integrate
  use, intrinsic :: iso_c_binding,cl=>c_long_double
  implicit none

  interface
    function gl10(a, b, f) &
    bind(c,name='Gauss_Legendre_Integration_10pts') !the function name case-sensitive because C is case-sensitive
      import::cl,c_funptr
      real(cl), intent(in), value :: a, b
      type(c_funptr), intent(in), value :: f
      real(cl) gl10
    endfunction

    function f1(x) bind(c)
      import::cl
      real(cl),intent(in), value ::x
      real(cl) f1
    endfunction
  endinterface

  type(c_funptr) :: f1_cptr
  real(cl)  :: a,b

  a=0; b=1

  f1_cptr = c_funloc (f1) !gets the C address of the fortran function f1
  write(*,*)'error= ',gl10(a,b,f1_cptr)-4*atan(1.0_cl)
endprogram

function f1(x) bind(c)
    use, intrinsic :: iso_c_binding,cl=>c_long_double
    implicit none
    real(cl),intent(in), value ::x
    real(cl) f1
    f1=4.0/(1.0+x**2)
  endfunction

C function interface

We have to use the right data types in fortran corresponding to those in C. These types are present in the module 'iso_c_binding'. We will be using two of these, i.e "c_long_double" and c_funptr. "c_long_double" is a kind of real, whereas 'c_funptr' is a derived datatype. I have just aliased the 'c_long_double' to 'cl' to reduce typing.

The interface block is used to make fortran understand how to call the function. The first function in the interface, 'gl10', is prototype for the C function that does the quadrature calculation. It is declared with the 'bind(C)' attribute along with the name of the function in C. By this interface, we give this C function, 'Gauss_Legendre_Integration_10pts', a different name in fortran, 'gl10'. Now, in fortran we can call this function as g10(...). The C function name that appears within 'bind' is case-sensitive.

There is some flexibility available in interfacing with the C functions. In general, in fortran, we can match the variables, pointers and function pointers in the argument list of a C function in two different ways. I do not know if one way is preferable over the other however.

The first two arguments in 'Gauss_Legendre_Integration_10pts' are scalar variables of type real(cl) and not pointers. As the convention in C is to pass scalars by value, the interface declares these arguments with the 'value' attribute. It is required because the C function expects a value, whereas a fortran-language call typically sends out the reference/address (usually). Now, if 'a' and 'b' in the C function definition were instead '*a' and '*b', the 'value' attribute should not be present in the interface. In our case we can edit the C code('a' to '*a' and 'b' to '*b' throughout), and remove the value attribute from 'a' and 'b' in the interface and verify that this approach also works.

Alternately, in the case with pointer arguments (like '*a','*b') for the C function, we can find the C address (using C_LOC(..)) of the corresponding variables in fortran and pass this address by value. The third argument, a C-style function pointer, is also declared with the 'value' attribute. This way the C function recieves a function address for the function pointer, f . For interoperability, for pointers that are dummy arguments, the 'value' attribute is usually required because void* matches TYPE(C_PTR), VALUE, while TYPE(C_PTR) by itself, i.e. without the 'value' attribute, matches void**. The same goes for function pointers. If the value attribute was not present in the interface for 'f', then '*f' in the C function would have to be replaced by '**f' throughout.

Here our choice of using the value attribute for all the dummy-arguments of the C function is dictated by our wish to leave the C code untouched

External fortran function

The function to be integrated, f1, is defined outside the scope of the main program. So, an interface has to be provided for this function as well. Since this function is being called from C code, it needs to be interoperable with C. So, this function is defined with the bind(C) attribute. The name of this function is only needed by the main fortran program. The C function access/calls this function using a pointer and so it doesnt need the name in this case.

But, in the general case this fortran function can be called from C code as well. Like before, The 'name' argument of bind can be provided if one wishes to call this function with a different name within the C code (see the end of the page).

C_FUNLOC(..)

C_FUNLOC(..) takes the name of a fortran function and returns a C function address. The type is a derived type called 'c_funptr' that is present in the iso_c_binding module. Here I have used the variable 'f1_cptr' to store the address of the fortran function 'f1'. Now that all the arguments have been defined, the C function is called in the program as 'gl10(a,b,f1_cptr)'.

To see how close the integral gets to \(\pi\) the error is calculated and printed. That wraps it up for now. The code can be downloaded from the links below.

Download Links:

  1. Main Fortran Program
  2. Gaussian Integration C function
To compile the code the following can be done for gcc/gfortran:

$>gcc -o gauss_legendre_10pts.o -c gauss_legendre_10pts.c
$>gfortran -o integrate.o -c integrate.f90
$>gfortran -o integrate.exe integrate.o gauss_legendre_10pts.o

     

Footnote

When compiling with gfortran & gcc, I am seeing unexpected behavior while calling a fortran function from within a C function. In the C function when I printf the value returned by the fortran function, it always shows as zero, irrespective of the actual value. Intel compilers work as expected though. My gfortran version:

Using built-in specs.
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=C:/TDM-GCC-64/bin/../libexec/gcc/x86_64-w64-mingw32/5.1.0/lto-wrapper.exe
Target: x86_64-w64-mingw32
Configured with: ../../../src/gcc-5.1.0/configure --build=x86_64-w64-mingw32 --enable-targets=all --enable-languages=ada,c,c++,fortran,lto,ob
jc,obj-c++ --enable-libgomp --enable-lto --enable-graphite --enable-cxx-flags=-DWINPTHREAD_STATIC --disable-build-with-cxx --disable-build-po
ststage1-with-cxx --enable-libstdcxx-debug --enable-threads=posix --enable-version-specific-runtime-libs --enable-fully-dynamic-string --enab
le-libstdcxx-threads --enable-libstdcxx-time --with-gnu-ld --disable-werror --disable-nls --disable-win32-registry --prefix=/mingw64tdm --wit
h-local-prefix=/mingw64tdm --with-pkgversion=tdm64-1 --with-bugurl=http://tdm-gcc.tdragon.net/bugs
Thread model: posix
gcc version 5.1.0 (tdm64-1)