Wrapping Fortran with F2py

Background information

  • Python can really only call C – it can’t call Fortran directly. So f2py exists to automatically build C bindings so that you can call an automatically generated C file which will in turn call Fortran

  • If you have a pure Fortran file called foo.f, we typically create a new file called p_foo.f that (also written in Fortran) will contain f2py pragmas’ so that f2py knows how to build the C bindings

Guide

Let’s suppose you’ve got a subroutine like this:

SUBROUTINE add(array1,array2,sums)
  include "b02.h"
  REAL *8 array1(NJM), array2(NJM), sums(NJM)
  INTEGER i
  do i=1,NJT
     sums(i) = array1(i) + array2(i)
  enddo
END
RETURN

You’d “wrap” it by creating p_Add that just calls Add, but sets up some appropriate Cf2py hints. Note that in Fortran, any line that starts with a C is a comment; so the Cf2py directives are ignored by the Fortran compiler, and used only by f2py while processing the file.

SUBROUTINE p_add(array1,array2,sums, array1x, array2x, sumsx)
  IMPLICIT NONE
  INTEGER array1x, array2x, sumsx
  REAL *8 array1(array1x), array2(array2x), sums(sumsx)
Cf2py intent(in,out) sums
  INCLUDE '../include/b03.h'
  CALL add(array1, array2, sums)
RETURN
END

If you want to create a Cf2py intent(in) for array1 and array2 you can, but apparently those are optional. array1x, array2x, sumsx are automatically passed by f2py’s generated C code based on what the input data contains. So if array1 were of length 1000, then array1x would automatically be equal to 1000, etc.

After you compile this, you could call it from python like this:

import myfortranmodule
import numpy

array1 = numpy.ones(1000, numpy.float64)
array2 = numpy.ones(1000, numpy.float64)
sums = numpy.ones(1000, numpy.float64)

sums = myfortranmodule.p_Add(array1, array2, sums)

This was also assuming that NJM and NJT exist in common block b02 and were already setup to be the correct values (1000 in this case). You can set common blocks by doing something like this (from python):

myfortranmodule.b02.NJM = 1000
myfortranmodule.b02.NJT = 1000

To compile this, you probably just want to add it to an existing module. Let’s say you wanted to add it to fem_meshmodule.so. You’d edit the Makefile.in in <continuity>/src/mesh. You’d add p_add.f to MESHMODSRC, and Add.f to MESHSRC. You’d probably put p_add.f in <continuity>/src/mesh and Add.f in <continuity>/src/fsrc. If they were electrophysiology files, then you’d put them in ep and electrophys respectively. I think fsrc was originally intended to contain files used by multiple modules (mesh, nodes, bm, etc.). And it does. But it also has files that are exclusive to mesh, which it probably shouldn’t.

Writing Callbacks

It is possible to pass a pointer to a python Function to Fortran so that you can “call-back” to python. This is quite tricky and we have generally moved away from this strategy over the last few years. Instead you are probably better off refactoring your Fortran code such that a callback is no longer necessary

Row vs Column ordering

Row-major order is used in C and python; column-major order is used in Fortran. This means that if you take a numpy array and send it directly to Fortran it will need to convert it from row-major to column major. If you do this in the middle of a critical loop you will see a major performance penalty. Instead what you can do is use numpy.asfortranarray to convert your data once. For example, this code:

array1 = numpy.ones(1000, numpy.float64)
array2 = numpy.ones(1000, numpy.float64)
sums = numpy.ones(1000, numpy.float64)

for i in range(10000000000):
    sums = myfortranmodule.p_add(array1, array2, sums)

Would run much more quickly like this:

array1 = numpy.ones(1000, numpy.float64)
array2 = numpy.ones(1000, numpy.float64)
sums = numpy.ones(1000, numpy.float64)

array1 = numpy.asfortranarray(array1)
array2 = numpy.asfortranarray(array2)
sums = numpy.asfortranarray(sums)

for i in range(10000000000):
    sums = myfortranmodule.p_add(array1, array2, sums)