# What Ever Happened to Fortran?

## 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`

.