Wrapping Fortran with F2py
- 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
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.
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)