Fortran - Cython 工作流程

Question

我想设置一个工作流程，使用Cython在Windows机器上从Python调用Fortran例程。

经过一番搜索，我找到了以下资料：http://www.fortran90.org/src/best-practices.html#interfacing-with-c 和 https://stackoverflow.com/tags/fortran-iso-c-binding/info

还有一些代码片段：

Fortran端：

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端：

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

以及安装文件：

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 )

所有文件都在一个目录中。

Fortran文件编译（使用NAG Fortran Builder），pygfunc编译。

但将它们链接起来会出现以下错误：

error LNK2019：未解析的外部符号_c_gfunc，该符号在函数___pyx_pf_7pygfunc_f中被引用

当然还有：

fatal error LNK1120：1个无法解析的外部

我漏掉了什么？或者说这种Python和Fortran之间的工作流程设置从一开始就注定要失败吗？

谢谢 Martin

- martburg

---

这很奇怪。你难道不是在Fortran中没有明确接口的情况下，从某个地方调用c_gfunc，或者以Fortran过程而不是C过程的方式调用它吗？ - Vladimir F Героям слава

抱歉，我不明白这个问题。 - martburg

不是解决你的问题，但如果你不知道，还有一个备选方案可以考虑：f2py。 - steabert

没关系，这没有任何进展。显然Cython根本没有使用你的Fortran对象文件。我不太了解Cython，无法帮助你。我能够使用“ctypes”运行你的代码。 - Vladimir F Героям слава

@steabert，我认为ctypes更接近于Cython，虽然这并不是解决问题的方法。Fortran代码可以保持原样。 - Vladimir F Героям слава

Answer 1

这是一个最小工作示例。我使用了gfortran，并将编译命令直接写入设置文件中。

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

我所做的更改并不是非常基础，以下是一些重要的更改。

• 你混合了双精度和单精度浮点数。不要这样做。在代码中同时使用 real (Fortran)、float (Cython) 和 float32 (NumPy)，并同时使用 double precision (Fortran)、double (Cython) 和 float64 (NumPy)。尽量不要意外混用它们。我在我的示例中假设您想使用双精度。

• 应该将所有变量作为指针传递给 Fortran。在这方面，它与 C 调用约定不匹配。Fortran 中的 iso_c_binding 模块只与 C 命名约定匹配。将数组作为指针传递，并将它们的大小作为单独的值传递。可能有其他方法来解决这个问题，但我不知道。

我还在设置文件中添加了一些内容，以显示在构建时可以添加一些更有用的额外参数。

要进行编译，请运行python setup.py build_ext --inplace。要验证它是否有效，请运行测试脚本。

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

以下是我一段时间前编写的两个示例：ftridiag、fssor。虽然我肯定不是这方面的专家，但这些示例可能是一个好的起点。

- IanH

---

谢谢，明确的setup.py是一个启示。 - martburg

我有一篇Cython帖子，你可能能够提供见解。 - ballade4op52