Whether you like it or not, for a long time Fortran has been dominating language in scientific computing community and there are still many maintained projects with large Fortran codebases. There are also new projects developed using Fortran from the ground up. Why is it so? Is it a power of habit?

While it is certainly true that some people code in Fortran just because they are used to it, one cannot deny that code generated by most Fortran compilers cannot be outperformed (in terms of speed) by most other languages. Combine it with native support for multidimensional arrays and you can see why scientist and engineers still use this ancient programming language.

Personally, I find Fortran to be a terrible language to code in. Its array indexing starts from 1 (what kind of monster came up with this?), it has some peculiar syntactic constructs (why do you need to `call`

a subroutine but can invoke a function like in any other language?) and so on and so forth. Writing any nontrivial Fortran program, especially one with some user interface and/or IO handling is undoubtedly a daunting task. It would be cool if we could leverage Fortran’s speed in some higher level language - and fortunately we can.

This post is meant as a first part of a series on combining Python with Fortran. I found that resources on the subject are somewhat scattered out there on the Internet, making it hard to find a single guide or series of guides covering everything from beginner level and gradually through more advanced topics. Hopefully the series I’m planning to write could serve as one.

In today’s part we’ll start with the following tasks:

- wrapping a Fortran subroutine using
`f2py`

, - choosing a compiler used by
`f2py`

and configuring its flags.

The above should be enough to get you started with some simple routines.

## Prerequisites

For the examples showed below I used Python 3.6 but everything should work just fine for 2.7. Besides Python you will also need to have the following installed:

- Python packages:
`numpy`

and`imageio`

- some Fortran compiler

## What we would like to achieve

For demonstration purposes we will implement a simple program for generating images of Mandelbrot set. I’ve choosen this particular problem since it is well known and easy to understand, yet requires handling of two-dimensional arrays.

The idea is to code the computational-heavy part, i.e. computing each point’s escape iteration, in Fortran and then write some script in Python to invoke our routine.

If you don’t remember a thing about Mandelbrot here is a quick reminder: consider point $p = p_x + i p_y$ in a complex plane and define a sequence:

The point $p$ lies in the Mandelbrot set if and only if $|z_n| < 2$ for every $n$. For a point $p$ not belonging to the Mandelbrot set, the first $n$ for which the above inequality is violated can be used for computing nice color of this point. This is another whole topic though, we will just use this value of $n$ as a number on a scale of grey.

If preceeding description does not make sense to you, you still can make sense of the forecoming sections - just treat a Fortran code as a black box. You will miss some fun though.

## Step 1: writing our Fortran routines

The below is the implementation of the required functions and subroutines:

- function
`color`

: given point’s coordinates and maximum iterations compute index of iteration after which it beceomes evident that it is not in the Mandelbrot set. After reaching maximum number of iterations we simply return`maxiter`

. - subroutine
`compute_colors`

: basically run`color`

for every point in the grid and store result in output array. Note that the grid is defined by its upper left and bottom right corners as well as number of pixels to compute - which is the same as shape of the output array.

Nothing special here. Now let’s move on to wrapping our routines with `f2py`

.

## Step 2: wrapping Fortran routines with f2py

The `f2py`

script (and acoompanying Python modules) are part of `numpy`

. The script basically creates Python extension that calls your Fortran routines. There are three ways of using `f2py`

, described in its documentation. Today we will only use the simple scenario, leaving the more advanced options for another time.

Basic usage of the `f2py`

goes as follows:

This should create file `mandelbrot.cpython-36m-x86_64-linux-gnu.so`

in current directory (exact name may differ depending on your OS and Python version).
Breaking down the above command:

- the
`-m`

argument controls how the created Python module will be named. It is the name that you will import to use your Fortran functions. - the
`-c`

tells`f2py`

to compile given files.

The newly created module should be importable. Exceute the following Python code to examine its contents:

The beginning of output should be similar to the following:

So our `mandelbrot`

module has a `mandelbrot`

Fortran module that contains two functions. If you didn’t get why we get a module inside module keep in mind that we could pass multiple Fortran files to `f2py`

, each of them containing separate Fortran module. So the top level one is named following `-m`

argument and its children correspond to Fortran modules you just compiled.

Now let’s look what are the signatures of functions that `f2py`

generated. Displaying their `__doc__`

attributes should give something similar to the output below.

There are no surprises for the `color`

function but look what happend with `compute_colors`

. The `f2py`

detected that we have a single output argument and wrapped our subroutine as a function that returns an array. Behind the scenes an array will be created for us at runtime and passed to Fortran routine and then returned back wrapped nicely as `numpy.ndarray`

. Creating arrays to be filled by subroutines is OK, but having them delivered to you is so much better!

## Step 3: write a Python script that uses it

We’ll go with a simple script that will just pass through command line arguments to the `compute_colors`

.

Executing the script like this should produce nice image of Mandelbrot set in shades of gray:

## Step 4 (optional): additional seasoning

So, suppose that for whatever reason you are planning to generate a lot of big ass images of Mandelbrot set. You are really concerned about performance so you are considering tweaking Fortran compiler flags and maybe throwing some parallelization into the mix. Or maybe something is not working as expected and you wish to include some debug information so you may solve the problem. To accomplish any of those tasks you need to specify compiler flags used by `f2py`

. You can do so in one of two, mostly equivalent, ways:

- by specifying
`--f90flags`

command line arguments - by specifying
`F90FLAGS`

environmental variable

So for instance, assuming that you are using `gfortran`

as your Fortran compiler, you might specify `-Ofast`

and `-finline-functions`

optimization flags in one of the following ways:

But wait, “assuming that you are using `gfortran`

”?. How do I tell `f2py`

to use specific Fortran compiler? Again, you have basically two options:

- speficy
`--fcompiler`

command line argument - specify
`F90`

environmental variable

This time however things are a bit more tricky. The `--fcompiler`

argument accepts some predefined values representing compilers that `f2py`

can work with. For instance:

- “pg”: Portland Group Fortran compiler
- “intelem”: ifort
- “gnu95”: gfortran

On the other hand, the F90 env variables has to be set to an executable used to compile the code. So for compiling our extension using PG fortran compiler we could do one of the following:

Why isn’t using the `--fcompiler`

as intuitive as using `F90`

environmental variable? Most probably everything boils down to the fact that `f2py`

can work only with predefiend compilers it knows how to handle - you may see that `f2py`

itself adds some compiler flags which differ depending on compiler you use. Still, the names they chose for the possible values are unintuitive at least.

## Conclusions

Hopefully after reading this post you can wrap you Fortran routines into nice Python modules and alter copmiler flags used by `f2py`

to suite your needs. There are still some interesting questions that we shall answer someday else, for instance:

- if we would like to include additional libraries in the linking process, how do we do that?
- what if it was not my intention to convert
`compute_colors`

to function returning an array and I’d rather pass the target array myself? - how do I package Python project with Fortran extensions?

We will look into the above topics in the next post in the series.