Fortran - Cython 工作流程

问题描述

我想设置一个工作流,以便在 Windows 机器上使用 Cython 从 Python 访问 fortran 例程

I would like to set up a workflow to reach fortran routines from Python using Cython on a Windows Machine

经过一番搜索,我发现:http://www.fortran90.org/src/best-practices.html#interfacing-with-c 和 https://stackoverflow.com/tags/fortran-iso-c-binding/信息

after some searching I found : http://www.fortran90.org/src/best-practices.html#interfacing-with-c and 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

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

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 )

所有文件都在一个目录中

all the files ar in one directory

fortran 文件编译(使用 NAG Fortran Builder)pygfunc 编译

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

但链接它们会引发:

错误 LNK2019:引用了未解析的外部符号 _c_gfunc在函数___pyx_pf_7pygfunc_f

error LNK2019: unresolved external symbol _c_gfunc referenced in function ___pyx_pf_7pygfunc_f

当然还有:

致命错误 LNK1120:1 个未解决的外部问题

fatal error LNK1120: 1 unresolved externals

我错过了什么?还是这种在 Python 和 Fortran 之间建立工作流的方式从一开始就该死?

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

THX马丁


解决方案

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

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.

  • 您混合了双精度和单精度浮点数.不要那样做. 同时使用 real (Fortran)、float (Cython) 和 float32 (NumPy),同时使用双精度 (Fortran)、double (Cyton) 和 float64 (NumPy).尽量不要无意中混合它们.我假设你在我的例子中想要双打.

  • 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.

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

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.

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

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

这是 fortran90.org 上显示的示例:mesh_exp

Here is the example shown on fortran90.org: mesh_exp

这是我前段时间整理的另外两个:ftridiag,fssor我当然不是这方面的专家,但这些示例可能是一个很好的起点.

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.

相关文章