# Calling Fortran from Python, Julia and Matlab

I continue this series with illustrating some more ways of calling the Mandelbrot function implemented in Fortran from other languages. What unifies the approaches presented here is that they all use a foreign function interface (FFI), to call a function in a shared library.

Again, the code I use here can be accessed with

git clone https://maurow@bitbucket.org/maurow/mauro_learning_fortran.git cd mauro_learning_fortran git checkout v0.5 cd mandel

## Shared library

First I need to make a shared library out of the `mandel_wrap.f90`

gfortran -c -fPIC mandel.f90 -o mandel.o gfortran -shared -fPIC mandel.o mandel_wrap.f90 -o libmandel.so

I also made a header file `mandel_wrap.h`

which contains the line

void c_calc_num_iter(int nre, double re[nre], int nim, double im[nim], int itermax, double escape, int out[nre][nim]);

This shared library can be used from C using the existing
`run_mandel_from_c.c`

by compiling it like so:

gfortran -L. -lmandel run_mandel_from_c.c -o cout_so

or from Fortran (without the wrapper)

gfortran -shared -fPIC mandel.f90 -o libmandel_fortran.so gfortran -L. -lmandel_fortran run_mandel_from_fortran.f90 -o fout_so

## Python-cffi

This is Python's take on FFI (note that I used version 0.7.2).

`cffi_test.py`

:

from cffi import FFI import numpy as np itermax = 20 escape = 2. nn = 5000 ni = 1001 ffi = FFI() # this defines the function signature: ffi.cdef("void c_calc_num_iter(int nre, double *re, int nim, double *im, int itermax, double escape, int out[%d][%d]);" % (nn, ni)) lib = ffi.dlopen('./libmandel.so') #lib = ffi.dlopen('./libmandel_openMP.so') # using openMP does not really work with cffi 0.7.2 ... re = ffi.new("double[]", np.linspace(-2.1,1,nn).tolist()) im = ffi.new("double[]", np.linspace(-1.3,1.3,ni).tolist()) res = ffi.new("int[%d][%d]" % (nn, ni)) lib.c_calc_num_iter(nn, re, ni, im, itermax, escape, res)

The `ffi.cdef("void ...`

defines the function signature, like in a C
header file, `ffi.dlopen`

loads the shared library, `ffi.new`

creates
a new C-variable and `lib.c_calc_num_iter`

then calls the function.

This is quite a bit shorter than the Cython version and requires no
build step. However, one would still want to write a Python wrapper
function and thus the total amount of code would be similar. In
particular the `res`

variable is not a Numpy array.

The advantage of this approach is that one does not need to learn a glue language (Cython) but can do all in Python + C. But Cython's extra power will probably make it better suited to more complicated wrapping jobs, in particular when using Numpy (e.g.).

Note that the openMP versions of the shared library did not run faster in python-cffi but in all other bindings it did…

## Julia

Julia is a new cool programming language (seriously, check it out).
It has a built-in FFI with the `ccall`

function. I updated the
wrapper a bit by dropping the `transpose`

, because Julia (like Fortran
and Matlab) uses colum-major array storage.

`mandel_wrap_colmajor.f90`

:

module mandel_wrap ! To wrap calc_num_iter for use in Julia with column-major arrays. ! Still using iso_c_binding. use iso_c_binding, only: c_double, c_int use mandel, only: calc_num_iter implicit none contains ! need to make a subroutine as only scalars can be returned subroutine f_calc_num_iter(nre, re, nim, im, itermax, escape, out) bind(c) real(c_double), intent(in):: re(nre), im(nim) real(c_double), intent(in), value:: escape integer(c_int), intent(in), value:: itermax, nre, nim integer(c_int), intent(out):: out(nre, nim) out = calc_num_iter(re, im, itermax, escape) end subroutine f_calc_num_iter end module mandel_wrap

Compile to a shared library and call from Julia:

`julia_test.jl`

nn = 10000 ni = 3000 re = linspace(-2.1,1,nn) im = linspace(-1.3,1.3,ni) itermax = 50 escape = 2.0 result = zeros(Cint, nn, ni) ccall((:f_calc_num_iter, "libmandel_colmajor"), Void, (Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cdouble, Ptr{Cint}), nn, re, ni, im, itermax, escape, result)

(The general syntax of Julia is very similar to Matlab's.) All the
FFI magic is in the `ccall`

function with arguments: ```
(:f_calc_num_iter,
"libmandel_colmajor")
```

gives the function name and name of the shared
library; `Void`

is the return type of the function; ```
(Cint,
Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cdouble, Ptr{Cint})
```

is the
function signature; and `nn, re, ni, im, itermax, escape, result`

are
the arguments to the function call.

Note that, unlike with python-cffi, `result`

is usable as an ordinary
Julia-array. However, a wrapper function would still be good.

## Matlab

Matlab also has a FFI usable through the functions `loadlibrary`

and
`calllib`

. It needs a header file to figure out the calling signature
(or a complicated Matlab file):

`mandel_wrap_colmajor.h`

:

void f_calc_num_iter(int nre, double re[], int nim, double im[], int itermax, double escape, int *out);

Then the shared library can be called like so (again using the library for the colum-major wrapper):

`matlab_test.m`

:

libname = 'libmandel_colmajor'; [notfound,warnings]=loadlibrary(libname, 'mandel_wrap_colmajor4matlab.h') nn = 1000; ni = 1100; re = linspace(-2.1,1,nn); im = linspace(-1.3,1.3,ni); itermax = 50; escape = 2.0; result = zeros(nn, ni)-5; [re, im, result] = calllib(libname, 'f_calc_num_iter', nn, re, ni, im, itermax, escape, result);

So, `loadlibrary(libname, 'mandel_wrap_colmajor4matlab.h')`

loads the
library using both the library and above header file. `calllib`

then
calls the function `f_calc_num_iter`

with ```
nn, re, ni, im, itermax,
escape, result
```

as arguments. Matlab does strictly pass-by-value thus
`calllib`

does some magic and returns all the arguments which are
pointer types in case they get modified and are needed as output.

## Asides

### The function signature

The type needed for the last argument in the function signature varies:

C and python-cffi are happy with:

void c_calc_num_iter(int nre, double re[nre], int nim, double im[nim], int itermax, double escape, int out[nre][nim]);

(python-cffi actually needs the actual numbers for `nre`

and `nim`

)

whereas Cython, Julia and Matlab want:

void f_calc_num_iter(int nre, double re[], int nim, double im[], int itermax, double escape, int *out);

I.e. note the difference in the last argument, the 2D array. I think this is due to the slight oddness of C 2D arrays (see).

`iso_c_binding`

It is actually possible to do the wrapper without using the
`iso_c_binding`

module.

`mandel_wrap_f90.f90`

:

module mandel_wrap_f90 ! To wrap calc_num_iter into a subroutine for use in Julia without ! the use of iso_c_binding use mandel, only: calc_num_iter, dp implicit none contains ! need to make a subroutine as only scalars can be returned subroutine f90_calc_num_iter(nre, re, nim, im, itermax, escape, out) bind(c) ! the bind(c) avoids name mangling in shared lib real(dp), intent(in):: re(nre), im(nim) real(dp), intent(in):: escape integer, intent(in):: itermax, nre, nim integer, intent(out):: out(nre, nim) out = calc_num_iter(re, im, itermax, escape) end subroutine f90_calc_num_iter end module mandel_wrap_f90

Note the `bind(c)`

avoids the name-mangling applied by Fortran compilers.

After compilation to a shared library, this can be called in, e.g., Julia with:

ccall((:f90_calc_num_iter, "libmandel_f90"), Void, (Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}), &nn, re, &ni, im, &itermax, &escape, result)

Note that the signature is now all pointers as Fortran does call by reference.

I think the advantage of using the `iso_c_binding`

is that the types
of C are more clearly defined than Fortran's and cannot change with
compile options. But I'm not sure.