Calling Fortran from Python

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

pymandel.pyx

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

Notes:

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

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.

Now Python can be used to plot the Mandelbrot set:

test.py

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]))
plt.show()

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…