Spicing up your Python with a little bit of Fortran - part 1
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
- choosing a compiler used by
f2pyand configuring its flags.
The above should be enough to get you started with some simple routines.
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:
- 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:
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
compute_colors: basically run
colorfor 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
Step 2: wrapping Fortran routines with f2py
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:
-margument controls how the created Python module will be named. It is the name that you will import to use your Fortran functions.
f2pyto 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:
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
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
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
--f90flagscommand line arguments
- by specifying
So for instance, assuming that you are using
gfortran as your Fortran compiler, you might specify
-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:
--fcompilercommand line argument
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.
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_colorsto 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.