Fortran - Cython Workflow

Here's a minimum working example.I used gfortran and wrote the compile commands directly into the setup file.

gfunc.f90

module gfunc_moduleimplicit nonecontainssubroutine gfunc(x, n, m, a, b, c)    double precision, intent(in) :: x    integer, intent(in) :: n, m    double precision, dimension(n), intent(in) :: a    double precision, dimension(m), intent(in) :: b    double precision, dimension(n, m), intent(out) :: c    integer :: i, j    do j=1,m        do i=1,n             c(i,j) = exp(-x * (a(i)**2 + b(j)**2))        end do    end doend subroutineend module

pygfunc.f90

module gfunc1_interfaceuse iso_c_binding, only: c_double, c_intuse gfunc_module, only: gfuncimplicit nonecontainssubroutine c_gfunc(x, n, m, a, b, c) bind(c)    real(c_double), intent(in) :: x    integer(c_int), intent(in) ::  n, m    real(c_double), dimension(n), intent(in) :: a    real(c_double), dimension(m), intent(in) :: b    real(c_double), dimension(n, m), intent(out) :: c    call gfunc(x, n, m, a, b, c)end subroutineend module

pygfunc.h

extern void c_gfunc(double* x, int* n, int* m, double* a, double* b, double* c);

pygfunc.pyx

from numpy import linspace, emptyfrom numpy cimport ndarray as arcdef extern from "pygfunc.h":    void c_gfunc(double* a, int* n, int* m, double* a, double* b, double* c)def f(double x, double a=-10.0, double b=10.0, int n=100):    cdef:        ar[double] ax = linspace(a, b, n)        ar[double,ndim=2] c = empty((n, n), order='F')    c_gfunc(&x, &n, &n, <double*> ax.data, <double*> ax.data, <double*> c.data)    return c

setup.py

from distutils.core import setupfrom distutils.extension import Extensionfrom Cython.Distutils import build_ext# This line only needed if building with NumPy in Cython file.from numpy import get_includefrom os import system# compile the fortran modules without linkingfortran_mod_comp = 'gfortran gfunc.f90 -c -o gfunc.o -O3 -fPIC'print fortran_mod_compsystem(fortran_mod_comp)shared_obj_comp = 'gfortran pygfunc.f90 -c -o pygfunc.o -O3 -fPIC'print shared_obj_compsystem(shared_obj_comp)ext_modules = [Extension(# module name:                         'pygfunc',                         # source file:                         ['pygfunc.pyx'],                         # other compile args for gcc                         extra_compile_args=['-fPIC', '-O3'],                         # other files to link to                         extra_link_args=['gfunc.o', 'pygfunc.o'])]setup(name = 'pygfunc',      cmdclass = {'build_ext': build_ext},      # Needed if building with NumPy.      # This includes the NumPy headers when compiling.      include_dirs = [get_include()],      ext_modules = ext_modules)

test.py

# A script to verify correctnessfrom pygfunc import fprint f(1., a=-1., b=1., n=4)import numpy as npa = np.linspace(-1, 1, 4)**2A, B = np.meshgrid(a, a, copy=False)print np.exp(-(A + B))

Most of the changes I made aren't terribly fundamental. Here are the important ones.

• You were mixing double precision and single precision floating point numbers. Don't do that. Use real (Fortran), float (Cython), and float32 (NumPy) together and use double precision (Fortran), double (Cyton), and float64 (NumPy) together. Try not to mix them unintentionally. I assumed you wanted doubles in my example.

• You should pass all variables to Fortran as pointers. It does not match the C calling convention in that regard. The iso_c_binding module in Fortran only matches the C naming convention. Pass arrays as pointers with their size as a separate value. There may be other ways of doing this, but I don't know any.

I also added some stuff in the setup file to show where you can add some of the more useful extra arguments when building.

To compile, run python setup.py build_ext --inplace. To verify that it works, run the test script.

Here is the example shown on fortran90.org: mesh_exp

Here are two more that I put together some time ago: ftridiag, fssor I'm certainly not an expert at this, but these examples may be a good place to start.