Here's a minimum working example.
I used gfortran and wrote the compile commands directly into the setup file.
gfunc.f90
module gfunc_module
implicit none
contains
subroutine 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 do
end subroutine
end module
pygfunc.f90
module gfunc1_interface
use iso_c_binding, only: c_double, c_int
use gfunc_module, only: gfunc
implicit none
contains
subroutine 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 subroutine
end 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, empty
from numpy cimport ndarray as ar
cdef 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 setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
# This line only needed if building with NumPy in Cython file.
from numpy import get_include
from os import system
# compile the fortran modules without linking
fortran_mod_comp = 'gfortran gfunc.f90 -c -o gfunc.o -O3 -fPIC'
print fortran_mod_comp
system(fortran_mod_comp)
shared_obj_comp = 'gfortran pygfunc.f90 -c -o pygfunc.o -O3 -fPIC'
print shared_obj_comp
system(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 correctness
from pygfunc import f
print f(1., a=-1., b=1., n=4)
import numpy as np
a = np.linspace(-1, 1, 4)**2
A, 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.
