python - Fortran-Cython Workflow

fortran-iso-c-binding

Question

Me gustaría configurar un flujo de trabajo para llegar a las rutinas de Fortran desde Python usando Cython en una máquina con Windows

después de buscar, encontré: http://www.fortran90.org/src/best-practices.html#interfacing-with-c y https://stackoverflow.com/tags/fortran-iso-c-binding/info

y algunos pices de código:

Lado de 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

Lado de 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

y el archivo de configuración:

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 )

todos los archivos están en un directorio

los compiladores fortran compilan (usando NAG Fortran Builder) compilaciones pygfunc

pero unirlos arroja un:

error LNK2019: símbolo externo no resuelto _c_gfunc al que se hace referencia en la función ___pyx_pf_7pygfunc_f

y por supuesto:

error fatal LNK1120: 1 external externo sin resolver

Qué me estoy perdiendo ? ¿O es así para establecer un flujo de trabajo entre Python y Fortran maldito desde el principio?

THX Martin

Answer 1

Aquí hay un ejemplo mínimo de trabajo. Usé gfortran y escribí los comandos de compilación directamente en el archivo de configuración.

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

La mayoría de los cambios que hice no son terriblemente fundamentales. Aquí están los importantes.

• Usted estaba mezclando números de coma flotante de doble precisión y precisión simple. No hagas eso. Use real (Fortran), float (Cython) y float32 (NumPy) juntos y use doble precisión (Fortran), doble (Cyton) y float64 (NumPy) juntos. Intenta no mezclarlos involuntariamente. Supuse que querías dobles en mi ejemplo.

• Debes pasar todas las variables a Fortran como punteros. No coincide con la convención de llamadas C a ese respecto. El módulo iso_c_binding en Fortran solo coincide con la convención de nomenclatura C. Pase las matrices como punteros con su tamaño como un valor separado. Puede haber otras formas de hacerlo, pero no conozco ninguna.

También agregué algunas cosas en el archivo de configuración para mostrar dónde puede agregar algunos de los argumentos adicionales más útiles al compilar.

Para compilar, ejecute python setup.py build_ext --inplace . Para verificar que funcione, ejecute el script de prueba.

Aquí está el ejemplo que se muestra en fortran90.org: mesh_exp

Aquí hay dos más que reuní hace un tiempo: ftridiag , fssor Ciertamente no soy un experto en esto, pero estos ejemplos pueden ser un buen lugar para comenzar.