Fortran-Cython工作流程

发布于 2021-01-29 18:17:13

我想建立一个工作流以在Windows机器上使用Cython从Python到达fortran例程

经过一番搜索后,我发现:http : //www.fortran90.org/src/best-
practices.html#interface-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编译

但链接它们会引发:

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

而且当然:

致命错误LNK1120:1个未解决的外部零件

我想念什么?还是从一开始就在Python和Fortran之间建立这种工作流的这种方式被指责?

THX马丁

关注者
0
被浏览
51
1 个回答
  • 面试哥
    面试哥 2021-01-29
    为面试而生,有面试问题,就找面试哥。

    这是一个最小的工作示例。我使用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(Cyton)和float64(NumPy)。尽量不要无意间将它们混合在一起。我以为您想在我的例子中加倍。

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

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

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

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

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



知识点
面圈网VIP题库

面圈网VIP题库全新上线,海量真题题库资源。 90大类考试,超10万份考试真题开放下载啦

去下载看看