Creating source for a basic extension module

Probably the easiest way to introduce f2py is to offer a simple example. Here is one of the subroutines contained in a file named add.f:

CSUBROUTINE ZADD(A,B,C,N)
CDOUBLE COMPLEX A(*)DOUBLE COMPLEX B(*)DOUBLE COMPLEX C(*)INTEGER NDO 20 J = 1, NC(J) = A(J)+B(J)20   CONTINUEEND

This routine simply adds the elements in two contiguous arrays and places the result in a third. The memory for all three arrays must be provided by the calling routine. A very basic interface to this routine can be automatically generated by f2py:

f2py -m add add.f

You should be able to run this command assuming your search-path is set-up properly. This command will produce an extension module named add module.c in the current directory. This extension module can now be compiled and used from Python just like any other extension module.

Creating a compiled extension module

You can also get f2py to compile add.f and also compile its produced extension module leaving only a shared-library extension file that canbe imported from Python:

f2py -c -m add add.f

This command leaves a file named add.{ext} in the current directory(where {ext} is the appropriate extension for a python extension module on your platform — so, pyd,etc. ). This module may then be imported from Python. It will contain a method for each subroutine inadd (zadd, cadd, dadd, sadd). The docstring of each method contains nformation about how the module method may be called:

>>> import add
>>> print add.zadd.__doc__
zadd - Function signature:zadd(a,b,c,n)
Required arguments:a : input rank-1 array('D') with bounds (*)b : input rank-1 array('D') with bounds (*)c : input rank-1 array('D') with bounds (*)n : input int

Improving the basic interface

The default interface is a very literal translation of the fortran code into Python. The Fortran array arguments must now be NumPy arrays and the integer argument should be an integer. The interface will attempt to convert all arguments to their required types (and shapes)and issue an error if unsuccessful. However, because it knows nothing about the semantics of the arguments (such that C is an output and n should really match the array sizes), it is possible to abuse this function in ways that can cause Python to crash. For example:

>>> add.zadd([1,2,3],[1,2],[3,4],1000)

will cause a program crash on most systems. Under the covers, thelists are being converted to proper arrays but then the underlying add loop is told to cycle way beyond the borders of the allocated memory.

In order to improve the interface, directives should be provided. This is accomplished by constructing an interface definition file. It is usually best to start from the interface file that f2py can produce(where it gets its default be havior from). To get f2py to generate the interface file use the -h option:

f2py -h add.pyf -m add add.f

This command leaves the file add.pyf in the current directory. The section of this file corresponding to zadd is:

subroutine zadd(a,b,c,n) ! in :add:add.fdouble complex dimension(*) :: adouble complex dimension(*) :: bdouble complex dimension(*) :: cinteger :: n
end subroutine zadd

By placing intent directives and checking code, the interface can becleaned up quite a bit until the Python module method is both easierto use and more robust.

subroutine zadd(a,b,c,n) ! in :add:add.fdouble complex dimension(n) :: adouble complex dimension(n) :: bdouble complex intent(out),dimension(n) :: cinteger intent(hide),depend(a) :: n=len(a)
end subroutine zadd

The intent directive, intent(out) is used to tell f2py that c is an output variable and should be created by the interface before being passed to the underlying code. The intent(hide) directive tells f2pyto not allow the user to specify the variable,n, but instead toget it from the size ofa. The depend(a ) directive is necessary to tell f2py that the value of n depends on the input a (so that it won’t try to create the variable n until the variable a is created).

The new interface has docstring:

>>> print add.zadd.__doc__
zadd - Function signature:c = zadd(a,b)
Required arguments:a : input rank-1 array('D') with bounds (n)b : input rank-1 array('D') with bounds (n)
Return objects:c : rank-1 array('D') with bounds (n)

Now, the function can be called in a much more robust way:

>>> add.zadd([1,2,3],[4,5,6])
array([ 5.+0.j,  7.+0.j,  9.+0.j])

Notice the automatic conversion to the correct format that occurred.

Inserting directives in Fortran source

The nice interface can also be generated automatically by placing the variable directives as special comments in the original fortran code.Thus, if I modify the source code to contain:

CSUBROUTINE ZADD(A,B,C,N)
C
CF2PY INTENT(OUT) :: C
CF2PY INTENT(HIDE) :: N
CF2PY DOUBLE COMPLEX :: A(N)
CF2PY DOUBLE COMPLEX :: B(N)
CF2PY DOUBLE COMPLEX :: C(N)DOUBLE COMPLEX A(*)DOUBLE COMPLEX B(*)DOUBLE COMPLEX C(*)INTEGER NDO 20 J = 1, NC(J) = A(J) + B(J)20   CONTINUEEND

Then, I can compile the extension module using:

f2py -c -m add add.f

The resulting signature for the function add.zadd is exactly the same one that was created previously. If the original source code had contained A(N) instead of A(*) and so forth with B and C, then Icould obtain (nearly) the same interface simply by placing the INTENT(OUT) :: C comment line in the source code. The only differenceis that N would be an optional input that would default to the length of A.

A filtering example

For comparison with the other methods to be discussed. Here is another example of a function that filters a two-dimensional array of double precision floating-point numbers using a fixed averaging filter. The advantage of using Fortran to index into multi-dimensional arrays should be clear from this example.

      SUBROUTINE DFILTER2D(A,B,M,N)
CDOUBLE PRECISION A(M,N)DOUBLE PRECISION B(M,N)INTEGER N, M
CF2PY INTENT(OUT) :: B
CF2PY INTENT(HIDE) :: N
CF2PY INTENT(HIDE) :: MDO 20 I = 2,M-1DO 40 J=2,N-1B(I,J) = A(I,J) +$           (A(I-1,J)+A(I+1,J) +$            A(I,J-1)+A(I,J+1) )*0.5D0 +$           (A(I-1,J-1) + A(I-1,J+1) +$            A(I+1,J-1) + A(I+1,J+1))*0.25D040      CONTINUE20   CONTINUEEND

This code can be compiled and linked into an extension module named filter using:

f2py -c -m filter filter.f

This will produce an extension module named filter.so in the current directory with a method named dfilter2d that returns a filteredversion of the input.

Calling f2py from Python¶

The f2py program is written in Python and can be run from inside your module. This provides a facility that is somewhat similar to the use of weave.ext_tools described below. An example of the final interface executed using Python code is:

import numpy.f2py as f2py
fid = open('add.f')
source = fid.read()
fid.close()
f2py.compile(source, modulename='add')
import add

The source string can be any valid Fortran code. If you want to save the extension-module source code then a suitable file-name can be provided by the source_fn keyword to the compile function.

Automatic extension module generation

If you want to distribute your f2py extension module, then you only need to include the .pyf file and the Fortran code. The distutils extensions in NumPy allow you to define an extension module entirely in terms of this interface file. A valid setup.py file allowing distribution of the add.f module (as part of the package f2py_examples so that it would be loaded as f2py_examples.add) is:

def configuration(parent_package='', top_path=None):from numpy.distutils.misc_util import Configurationconfig = Configuration('f2py_examples',parent_package, top_path)config.add_extension('add', sources=['add.pyf','add.f'])return configif __name__ == '__main__':from numpy.distutils.core import setupsetup(**configuration(top_path='').todict())

Installation of the new package is easy using:

python setup.py install

assuming you have the proper permissions to write to the main site-packages directory for the version of Python you are using. For theresulting package to work, you need to create a file named __init__.py(in the same directory as add.pyf). Notice the extension module is defined entirely in terms of the “add.pyf” and “add.f” files. The conversion of the .pyf file to a .c file is handled by numpy.disutils.

Conclusion

The interface definition file (.pyf) is how you can fine-tune theinterface between Python and Fortran. There is decent documentationfor f2py found in the numpy/f2py/docs directory where-ever NumPy isinstalled on your system (usually under site-packages). There is alsomore information on using f2py (including how to use it to wrap Ccodes) athttp://www.scipy.org/Cookbook under the “Using NumPy withOther Languages” heading.

The f2py method of linking compiled code is currently the mostsophisticated and integrated approach. It allows clean separation ofPython with compiled code while still allowing for separatedistribution of the extension module. The only draw-back is that itrequires the existence of a Fortran compiler in order for a user toinstall the code. However, with the existence of the free-compilersg77, gfortran, and g95, as well as high-quality commerical compilers,this restriction is not particularly onerous. In my opinion, Fortranis still the easiest way to write fast and clear code for scientificcomputing. It handles complex numbers, and multi-dimensional indexingin the most straightforward way. Be aware, however, that some Fortrancompilers will not be able to optimize code as well as good hand-written C-code.