How to Call Fortran from Python

There are basically two basic ways to call Fortran from Python:

• Use an automatic wrapper such as F2py or fmodpy. This is the simplest way, as the required wrappers are produced automatically.
• Combine iso_c_bindings on the Fortran side, with ctypes or Cython on the Python side. This is the most general way of using Fortran in Python.

For a thorough explanation of what is going on under the hood in the different alternatives, a good reference is the Numpy user manual section on ‘using Python as glue’, which you can find here.

F2py is the most established wrapper-generation library, as it is actually a part of Numpy, but regrettably supports a limited subset of Fortran — it is great for older Fortran code, though. You might also want to try out the newer fmodpy, that supports modern Fortran.

In this tutorial, we will go through some quick examples. In particular, we will test the different methods of passing a complex Numpy array to Fortran.

F2PY Example

This method is absolutely straightforward. Let’s say we have the following Fortran function,

!
!   my_fortran_lib.f90:
!
function foo(a) result(b)
    implicit none
    
    real(kind=8), intent(in)    :: a(:,:)
    complex(kind=8)             :: b(size(a,1),size(a,2))
    
    b = exp((0,1)*a)
    
end function foo

To compile it and create a Python extension module, use the following command: (pay attention to the “3” in F2py3 in case you intend to use Python 3).

f2py3 -c -m myflib my_fortran_lib.f90

That’s it, now the function is available to import it as a Python module. Simply call it by

import myflib
import numpy as np

a = np.array([[1,2,3,4], [5,6,7,8]], order='F')
myflib.foo(a)

Also, it is important to remember that C and Fortran employ different array orderings. Python normally employs the c-style row-major-order, while in Fortran the column-major-order is employed. By setting the Numpy order property to ‘F’ we are complying with the Fortran ordering.

Fmodpy example

An alternative method to automatically produce wrappers is to rely on fmodpy.

Let’s test the above example with fmodpy.

First, we have to install fmodpy, which is a Python package.

python3 -m pip install fmodpy

And then we can use it from within Python directly:

import fmodpy
import numpy as np

# This will compile and import the Fortran code.
# (will also recompile if Fortran source changed since last access!)
myflib = fmodpy.fimport("my_fortran_lib.f90")

a = np.array([[1,2,3,4], [5,6,7,8]], order='F')
myflib.foo(a)

That should give the same result as above.

A few advatages of fmodpy over f2py, besides this nice live recompilation feature, is that as it is built using iso_c_bindings on the fortran side and ctypes on the Python side, which are more modern features to what F2PY , which since the beginning has empasized compatibility with old F77 codebases, is using. (to the best of my knowledge, as there are no mentions to iso_c_bindings in the f2py source code).

For example, f2py is known to have issues under windows (see official docs on F2PY and Windows). I’m not a windows user, so I can’t really tell you.

How iso_c_bindings and ctypes work

Besides using automated libraries like f2py or fmodpy, sometimes wrappers need to be written manually.

In this case, the best way to do this (for enhanced reliability and compatibility), is by relying on iso_c_bindings on the Fortran side, and ctpyes on the Python side.

iso_c_bindings is a feature introduced in the Fortran 2003 standard to improve interoperability with C.

On the Python side, ctypes is a part of the stdlib which implements interfaces to foreign functions.

These two tools can be combined to create an interface that enables Fortran and Python compatibility.

Let’s see an example of this in action.

2.1 Getting started

As this is significantly more complicated than F2py we will start with a simple (yet meaningful) example, to get familiar with the basics. We will build a function that takes a double-precision float and adds the number 2.

As a first step, we need to declare the C types of each argument of our Fortran functions or subroutines. A good reference on how to translate each variable type can be found in the GCC docs here

In a file named my_fortran_lib.f90 we write:

function sum2(a) result(b) bind(c, name='sum2')
    use iso_c_binding
    implicit none
    
    real(c_double), intent(in)  :: a
    real(c_double)              :: b
    
    b = a + 2.d0
    
end function sum2

We now compile this code as a shared library, with the command

gfortran -shared my_fortran_lib.f90 -o myflib.so

Now we have a shared library and we need to import it into Python. We use ctypes for this matter. It is also important to remember that, by default, Fortran passes values by reference.

import ctypes as ct

# import the shared library
fortlib = ct.CDLL('myflib.so') 

# Specify arguments and result types
fortlib.sum2.argtypes = [ct.POINTER(ct.c_double)]
fortlib.sum2.restype = ct.c_double

# Create a double and pass it to Fotran (by reference)
a = ct.c_double(5)
b = fortlib.sum2(ct.byref(a))
print(b)

2.2 Dealing with complex numbers

Things start getting trickier here. Even though Fortran interoperability with C supports a ‘c_double_complex’ type which is part of the C99 standard, this is not supported by Python’s ctypes, at least at the time of writing this article. Follow this issue for updates on this matter. The workaround, for the moment, is to rely on wrapper functions, which can be easily coded in Fortran.

2.3 Dealing with multidimensional arrays

It is usually the case that we want to pass a large Numpy array to Fortran, as in our F2PY example, which also returns an array.

An important thing to notice is that the return type of a C-interoperable Fortran function can’t be an array. That’s because Fortran arrays are quite different from C arrays. The workaround is to rely on subroutines instead.

Our Fortran subroutine should look something like the following:

subroutine double_array(x,N) bind(C, name="double_array")
    use iso_c_binding
    implicit none
    
    integer(c_int), intent(in), value   :: N
    real(c_double), intent(inout)       :: x(N,N)
    
    x  = exp(x)
    
end subroutine double_array

After compiling it as a shared library in the same way as before, we must build a Python interface. For this matter, we can resort to Numpy.ctypeslib to help us in building the required interface.

import ctypes as ct
import numpy as np

# import the shared library
fortlib = ct.CDLL('myflib.so') 
f = fortlib.double_array

# Specify arguments type as a pointer to double and an integer
f.argtypes=[ct.POINTER(ct.c_double), ct.c_int]

# Create a double array, pass it to Fotran as a pointer
x = np.ones((3,3), order="F")
x_ptr = x.ctypes.data_as(ct.POINTER(ct.c_double))

# Call function
rint = f(x_ptr, ct.c_int(3))
print(x)

Do you need help with Fortran?

You might be also interested in my article on How to Work with Legacy Fortran Code: A Short Guideline.

If you are looking for general advice, or have a small bug, etc, a great resource to ask for help is the Fortran Forum.

You can also hire a freelancer! Make sure to contact them and ask them a question on how they can help you.