Фортран - рабочий процесс Cython

Ответ 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))

Большинство изменений, которые я сделал, не являются чрезвычайно фундаментальными. Вот важные.

• Вы смешивали числа с двойной точностью и одинарной точностью с плавающей запятой. Не делай этого. Используйте реальные (Fortran), float (Cython) и float32 (NumPy) вместе и используйте двойную точность (Fortran), double (Cyton) и float64 (NumPy) вместе. Старайтесь не смешивать их непреднамеренно. Я предположил, что вы хотели удвоить в моем примере.

• Вы должны передать все переменные Fortran в качестве указателей. В этом отношении он не соответствует конвенции C-вызова. Модуль iso_c_binding в Fortran соответствует только соглашению о присвоении имен. Передавайте массивы как указатели с их размером как отдельное значение. Могут быть другие способы сделать это, но я не знаю.

Я также добавил некоторые вещи в установочный файл, чтобы показать, где вы можете добавить некоторые из более полезных дополнительных аргументов при создании.

Чтобы скомпилировать, запустите python setup.py build_ext --inplace. Чтобы убедиться, что он работает, запустите тест script.

Вот пример, показанный на fortran90.org: mesh_exp

Вот еще два, которые я собрал некоторое время назад: ftridiag, fssor Я, конечно, не специалист в этом, но эти примеры могут быть хорошим местом для начала.