Fortran-Cython工作流程

Fortran-Cython工作流程

这是一个最小的工作示例。我使用gfortran并将编译命令直接写入安装文件。

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))

我所做的大多数更改都不是根本性的。这是重要的。

• 您正在混合双精度和单精度浮点数。 不要那样做 一起使用real（Fortran），float（Cython）和float32（NumPy），并一起使用double precision（Fortran），double（Cyton）和float64（NumPy）。尽量不要无意间将它们混合在一起。我以为您想在我的例子中加倍。

• 您应该将所有变量作为指针传递给Fortran。在这方面，它与C调用约定不匹配。Fortran中的iso_c_binding模块仅匹配C命名约定。将数组作为指针传递，其大小作为单独的值。可能还有其他方法可以执行此 *** 作，但我不知道。

我还在安装文件中添加了一些东西，以显示在构建时可以在其中添加一些更有用的额外参数的地方。

要编译，请运行

python setup.py build_ext --inplace

。要验证它是否有效，请运行测试脚本。

这是在fortran90.org上显示的示例：mesh_exp

这里还有两个，我放在一起前段时间：ftridiag，fssor
我当然不会在这方面的专家，但这些例子可能是一个良好的开端。