Or A New Coat of Paint on your Old Shed

For decades, talented scientists implemented amazing numerical algorithms in Fortran, so much so that it became known as the “lingua franca” of the computer industry. As a result, a large body of valuable software is slowly fading into oblivion. How can we make this software accessible to modern developers?

Most dynamic programming languages support a Foreign function interface that allows to programs to call routines written in a different (hence foreign) language. ctypes, which is part of the python standard library, is only one of many similar solutions. Other approaches, like SWIG and pybind11 rely on additional declarations to generate python bindings. So far, these tools generate python wrappers for C and C++ functions. Fortran functions and subroutines can be called from C so any of these should be sufficient in theory. In practice, this requires a reasonably deep understanding of the python ABI to work correctly.

This is precisely the niche that F2PY fills. Quoting from its documentation, “The purpose of the F2PY -Fortran to Python interface generator- utility is to provide a connection between Python and Fortran.” Large chunks of numpy and scipy are made possible by f2py. Now we need an interesting piece of Fortran to experiment with F2PY.

I have chosen a frightfully clever implementation of the Cooley-Tukey FFT algorithm by Dr Norman Brenner, now at The George Washington University. Dr Brenner’s is a multi-radix, arbitrary dimension, real or complex fast fourier transform. It is described in detail in a technical note published in 1967 “with the support of the U.S. Air Force under Contract AF 19(628)-5167”. A number of digital versions can be found online, often embedded in larger single file Fortran projects, as was typical at the time.

According to this reference page, FORTRAN IV arrays already “may have 1 to 7 subscripts”. However, because indices are explicit in the syntax, arbitrary numbers of dimensions require additional effort. The technical note explains that, “As usual, the first subscript varies the fastest in storage order”. In other words, if A is a three-dimension array with dimensions L × M × N, element A(I, J, K) is stored at offset I - 1 + L * ((J - 1) + M * (K - 1)), see here.

Another wonderful quirk of older (i.e. unstructured) Fortran dialects is the arithmetic IF. Arithmetic IF does not expect a boolean expression but a general expression and is followed by 3 comma-separated integers. Control flow is transferred to one of these integer labels based on the sign of the expression (first label when the expression is negative, second label when it is zero, or third label when it is positive). This allows for a really compact binary search.

As a final note, I got everything working and felt pretty good about myself until I read about the status of numpy.distutils. I consider simply ignoring it. After all, this is not a project I maintain or anything. But I realize how painful it would be to get back to it in several months once my github actions start failing miserably. So I bit the bullet, crossed the pons asinorum, and learned meson because that’s what SciPy uses. To be frank, I take exception about the nasty line about make in the f2py documentation. Call me old-fashioned but I find it hard to place a high value on the investment it took me to pick up meson.