Calling Fortran from Python

[2013-04-19 Fri] (edited [2013-09-12 Thu])

The way I use here to call Fortran from Python is to make the Fortran code callable from C and then bind that with Cython (version 0.19.1). This is what is recommended in "Fortran best practices". There are other ways, like F2Py and fwrap, but (I think) the former struggles with modern Fortran and the latter's development has stalled. Maybe another avenue to pursue could be CFFI. In fact, all this linking Python to compiled code seems to be quite in a state of flux at the moment and this could be outdated soon (or is already!).

I build on the code of Calling Fortran from C. The idea is to call the subroutine c_calc_num_iter in mandel\_wrap.f90 from Cython. Once that is done, then a call from Python is trivial. 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.2.1

and is in the folder mandel/.


The cython wrapper program wraps the same mandel\_wrap.f90 as I used in Calling Fortran from C.


import numpy as np

# this defines the external function's interface.  Why the out needs
# to be defined as 'int *out' and not 'int **out' I do not know...
cdef extern:
    void c_calc_num_iter(int nre, double *re, int nim, double *im,
                        int itermax, double escape, int *out)

def calc_num_iter(double[::1] re not None,
                   double[::1] im not None,
                     int itermax=20, double escape=2.):
    cdef int nre, nim
    nre = len(re)
    nim = len(im)
    # initialize the output array with cython so Python manages the
    # memory:
    out = np.empty((nre, nim), dtype=np.int32)
    cdef int [:,::1] res = out

    c_calc_num_iter(nre, &re[0], nim, &im[0],
                     itermax, escape, &res[0,0])
    return out


  • the cdef extern defines the function signature of the "c" function. In this case the signature is provided by using Fortran's iso_c_binding module in mandel_wrap.f90. I don't know why the last argument, a 2D array in Fortran, is just a simple pointer (int *out) as opposed to a pointer to a pointer (int **out).
  • This implementation uses memoryviews which makes the code nicer than the older numpy buffer style.
  • out is a Python variable, so Python manages the memory, which is good.
  • cdef int [:,::1] res makes a C-style array, which is good as that is whats returned from mandel_wrap.f90.

Compiling the Cython extension module

Compilation was the trickiest part! First make the the *.o files of the Fortran code:

gfortran -fPIC -c mandel.f90 mandel_wrap.f90

The -fPIC is needed.

Then, for cython compilation, it is usually best to make a setup.py file:

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy
npy_include_dir = numpy.get_include()

ext_modules = [Extension("mandel", ["pymandel.pyx"],
                         include_dirs = [npy_include_dir],
                         extra_objects=["mandel.o", "mandel_wrap.o"])]

setup(name = 'Mandelbrot fractals',
      cmdclass = {'build_ext': build_ext},
      ext_modules = ext_modules)

The none-standard parts is to include the extra *.o files. I think depending on the Fortran code there might be some libraries missing which can be included with libraries["somelib"]=. To figure out which other libraries are needed, run the last gcc command which running python setup.py build_ext --inplace spits out but replace gcc with (1) gcc -v and (2) gfortran -v and see what the difference are.

Running in Python

Now Python can be used to plot the Mandelbrot set:


import numpy as np
import pylab as plt
import mandel as md
nn = 1000
re = np.linspace(-2.1,1,nn)
im = np.linspace(-1.3,1.3,nn+1)

out = md.calc_num_iter(re, im) # calling Fortran

plt.imshow(out.T, extent=(re[0],re[-1],im[0],im[-1]))


Summary of steps

  1. write Fortran module
  2. write wrapper in Fortran to make function available in C
  3. make a Cython wrapper, wrapping between the Fortran/C function and Python
  4. call the thing from Python

Seems quite convoluted…