Calling Fortran from Python, Julia and Matlab

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

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 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 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.