Fortran - Cython Workflow
Question: I would like to set up a workflow to reach fortran routines from Python using
Cython on a Windows Machine

after some searching I found : <http://www.fortran90.org/src/best-
practices.html#interfacing-with-c> and <http://stackoverflow.com/tags/fortran-
iso-c-binding/info>

and some code pices:

Fortran side:

pygfunc.h:

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

pygfunc.f90

    
    
    module gfunc1_interface
        use iso_c_binding
        use gfunc_module
    
        implicit none
    
    contains
        subroutine c_gfunc(x, n, m, a, b, c) bind(c)
            real(C_FLOAT), intent(in), value :: x
            integer(C_INT), intent(in), value ::  n, m
            type(C_PTR),    intent(in), value :: a, b
            type(C_PTR),                value :: c
    
            real(C_FLOAT), dimension(:), pointer :: fa, fb
            real(C_FLOAT), dimension(:,:), pointer :: fc
    
            call c_f_pointer(a, fa, (/ n /))
            call c_f_pointer(b, fb, (/ m /))
            call c_f_pointer(c, fc, (/ n, m /))
            call gfunc(x, fa, fb, fc)
         end subroutine
    
    end module
    

gfunc.f90

    
    
    module gfunc_module
    
    use iso_c_binding
    
        implicit none
        contains
            subroutine gfunc(x, a, b, c)
                real,                 intent(in) :: x
                real, dimension(:),   intent(in) :: a, b
                real, dimension(:,:), intent(out) :: c
    
                integer :: i, j, n, m
                n = size(a)
                m = size(b)
                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
    

Cython side:

pygfunc.pyx

    
    
    cimport numpy as cnp
    import numpy as np
    
    cdef extern from "./pygfunc.h":
        void c_gfunc(double, int, int, double *, double *, double *)
    
    cdef extern from "./pygfunc.h":
        pass
    
    def f(float x, a=-10.0, b=10.0, n=100):
        cdef cnp.ndarray ax, c
        ax = np.arange(a, b, (b-a)/float(n))
        n = ax.shape[0]
        c = np.ndarray((n,n), dtype=np.float64, order='F')
        c_gfunc(x, n, n, <double *> ax.data, <double *> ax.data, <double *> c.data)
        return c
    

and the setup file:

    
    
    from distutils.core import setup
    from distutils.extension import Extension
    from Cython.Distutils import build_ext
    import numpy as np
    
    ext_modules = [Extension('pygfunc', ['pygfunc.pyx'])]
    
    setup(
        name = 'pygfunc',
        include_dirs = [np.get_include()],
        cmdclass = {'build_ext': build_ext},
        ext_modules = ext_modules )
    

all the files ar in one directory

the fortran files compile ( using NAG Fortran Builder ) pygfunc compiles

but linking them throws a:

error LNK2019: unresolved external symbol _c_gfunc referenced in function
___pyx_pf_7pygfunc_f

and of course:

fatal error LNK1120: 1 unresolved externals

What am I missing ? or is this way to set up a workflow between Python and
Fortran damned from the beginning ?

THX Martin


Answer: 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](https://github.com/arrghb2012/fortran90/tree/master/fcython_mesh)

Here are two more that I put together some time ago:
[ftridiag](https://github.com/byuimpact/numerical_computing/tree/develop/Python/cython_wrapping/ftridiag),
[fssor](https://github.com/byuimpact/numerical_computing/tree/develop/Python/cython_wrapping/fssor)
I'm certainly not an expert at this, but these examples may be a good place to
start.

