F2PY - Calling Fortran routines from Python

Last edited: February 11th, 2017


This tutorial gives a quick introduction to the F2PY package. Check out F2PY's website for a more complete reference and installation procedures. Note that F2PY is often included when installing NumPy and SciPy. From the documentation:

The purpose of the F2PY --Fortran to Python interface generator-- project is to provide a connection between Python and Fortran languages. F2PY is a Python package (with a command line tool f2py and a module f2py2e) that facilitates creating/building Python C/API extension modules that make it possible

How does F2PY work?

F2PY works by creating an extension module that can be imported in Python using the import keyword. The module contains automatically generated wrapper functions that can be called from Python, acting as an interface between Python and the compiled Fortran routines. First, F2PY reads the Fortran source file and creates a so-called signature file that contains all the necessary information about the Fortran routines needed to make the wrapper functions. The signature file is then read and the source code of the extension module is generated in C, using the Python C API. In the last step, F2PY compiles all the source code and builds the extension module containing the wrappers and the compiled Fortran routines.


Why should you use F2PY?

The choice of programming language can be challenging at times, especially when it comes to finding a balance between computational efficiency and implementation time and effort. While scripting languages like MATLAB and Python may provide intuitive code which is fast to implement, compiled languages like C/C++ and Fortran yield superior computational speed. By wrapping a compiled code for Python, we can get the best of both worlds! Our notebook Calling Fortran(95) routines from a Python Script shows an example of the usage and the gain in computational time.

When should you use F2PY?

This is perhaps the ultimate question, and unfortunately, there is no definite answer. A good rule of thumb however, is to use F2PY, or compiled languages in general, when performing multiple operations/- computations within (nested) loops. Possibly, the most typical example would be operations on elements in multidimensional matrices. That is, linear algebra in general. Other good examples could be programs calculating integrals or conducting Monte Carlo simulations. At this point, you might wonder if anyone has already made F2PY-modules fitting your particular problem. The answer is most likely yes! Most of the functions and routines found in NumPy and SciPy are actually compiled Fortran (or C/C++) routines which provide highly efficient and fast solvers for multiple problems. We thus advice you to always check if one of these two packages/libraries already provide a routine in which may be suitable for your problem. If not, you should first implement your solver in a pure Python script to investigate whether or not computational efficiency really is an issue. If it is, then F2PY may possibly provide the best solution strategy for your problem.

Getting started

We will be considering a simple example in which the sieve of Eratosthenes algorithm is used to compute prime numbers. The Fortran code is saved in primes.f95. Note that the code only includes subroutines and that the variables are defined with a new keyword intent. The latter is explained in the next section.

subroutine sieve(is_prime, n_max)
! =====================================================
! Uses the sieve of Eratosthenes to compute a logical
! array of size n_max, where .true. in element i
! indicates that i is a prime.
! =====================================================
    integer, intent(in)   :: n_max
    logical, intent(out)  :: is_prime(n_max)
    integer :: i
    is_prime = .true.
    is_prime(1) = .false.
    do i = 2, int(sqrt(real(n_max)))
        if (is_prime (i)) is_prime (i * i : n_max : i) = .false.
    end do
    return
end subroutine

subroutine logical_to_integer(prime_numbers, is_prime, num_primes, n)
! =====================================================
! Translates the logical array from sieve to an array
! of size num_primes of prime numbers.
! =====================================================
    integer                 :: i, j=0
    integer, intent(in)     :: n
    logical, intent(in)     :: is_prime(n)
    integer, intent(in)     :: num_primes
    integer, intent(out)    :: prime_numbers(num_primes)
    do i = 1, size(is_prime)
        if (is_prime(i)) then
            j = j + 1
            prime_numbers(j) = i
        end if
    end do
end subroutine

The simplest way to wrap this subroutine to python is to run

f2py -c primes.f95 -m primes

Note that you might need to run f2py3 to use Python 3! This command builds (-c flag) an extension module primes.so to the current directory. If the -m flag is excluded, the extension module will be named untitled.so.

We can now access these subroutines from Python:

>>> import primes
>>> print(primes.__doc__)
This module 'primes' is auto-generated with f2py (version:2).
Functions:
  is_prime = sieve(n_max)
  prime_numbers = logical_to_integer(is_prime,num_primes,n=len(is_prime))
.
>>> print(primes.logical_to_integer.__doc__)
prime_numbers = logical_to_integer(is_prime,num_primes,[n])

Wrapper for ``logical_to_integer``.

Parameters
----------
is_prime : input rank-1 array('i') with bounds (n)
num_primes : input int

Other Parameters
----------------
n : input int, optional
    Default: len(is_prime)

Returns
-------
prime_numbers : rank-1 array('i') with bounds (num_primes)

>>> sieve_array = primes.sieve(100)
>>> prime_numbers = primes.logical_to_integer(sieve_array, sum(sieve_array))
>>> print(prime_numbers)
[ 2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97]

Note that F2PY automatically found that the last argument (n) of the logical_to_integer subroutine was the dimension of the input array is_prime. F2PY concluded that n can be optional, with the default value len(is_prime)! One can use different values for the optional argument n. However, an exception is raised when it is incompatible with is_prime.

Specifying input and output arguments

In the example above, the different arguments of the subroutine were defined as input or output using the intent() attribute. The three most useful are:

If intent is excluded, the arguments become input-only arguments (same as using intent(inout)) by default. It is considered good practice to specify all arguments using the intent attribute. It is also preferred to use intent(out) (and not intent(inout)) to have a returned value.

The intent and optionality of the arguments can also be edited manually in the signature file primes.pyf generated by running

f2py primes.f95 -m primes -h primes.pyf

The final module is built from the signature file by running

f2py -c primes.pyf primes.f95

The attributes can also be specified as comments, which is done in our Calling Fortran(95) routines from a Python Script notebook.

Check out the documentation for the signature file for more options.

Pitfalls

integer, allocatable, intent(out) :: array1(:)  ! Not valid
integer, intent(out)              :: array2(:)  ! Not valid
integer, intent(out)              :: array3(10) ! Valid
integer, intent(in)               :: array4(:)  ! Valid

Custom docstrings

As we have seen, F2PY creates a default documentation for the module and functions which can be reached using e.g. help() or .__doc__. As far as we know, there are no options in F2PY in which we can modify this documentation. However, it can be changed upon import (<module>.__doc__=<string>) or one create a python function with its own (custom) docstring which calls the module. Many of SciPy's modules are built using F2PY, and their docstring are created using the latter method.