python – F2Py con matrices de formas asignables y asumidas

Pregunta:

Me gustaría usar f2py con el Fortran moderno. En particular, estoy tratando de que funcione el siguiente ejemplo básico. Este es el ejemplo útil más pequeño que pude generar.

! alloc_test.f90
subroutine f(x, z)
  implicit none

! Argument Declarations !
  real*8, intent(in) ::  x(:)
  real*8, intent(out) :: z(:)

! Variable Declarations !
  real*8, allocatable :: y(:)
  integer :: n

! Variable Initializations !
  n = size(x)
  allocate(y(n))

! Statements !
  y(:) = 1.0
  z = x + y

  deallocate(y)
  return
end subroutine f

Tenga en cuenta que n se infiere de la forma del parámetro de entrada x . Tenga en cuenta que y se asigna y se desasigna dentro del cuerpo de la subrutina.

Cuando compilo esto con f2py

f2py -c alloc_test.f90 -m alloc

Y luego ejecuta en Python

from alloc import f
from numpy import ones
x = ones(5)
print f(x)

Obtuve el siguiente error

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (-1,)

Así que voy y creo y edito el archivo pyf manualmente

f2py -h alloc_test.pyf -m alloc alloc_test.f90

Original

python module alloc ! in 
    interface  ! in :alloc
        subroutine f(x,z) ! in :alloc:alloc_test.f90
            real*8 dimension(:),intent(in) :: x
            real*8 dimension(:),intent(out) :: z
        end subroutine f
    end interface 
end python module alloc

Modificado

python module alloc ! in 
    interface  ! in :alloc
        subroutine f(x,z,n) ! in :alloc:alloc_test.f90
            integer, intent(in) :: n
            real*8 dimension(n),intent(in) :: x
            real*8 dimension(n),intent(out) :: z
        end subroutine f
    end interface 
end python module alloc

Ahora se ejecuta, pero los valores de la salida z son siempre 0 . Algunas impresiones de depuración revelan que n tiene el valor 0 dentro de la subrutina f . Supongo que me falta algo de magia de encabezado f2py para manejar esta situación correctamente.

De manera más general, ¿cuál es la mejor manera de vincular la subrutina anterior a Python? Preferiría encarecidamente no tener que modificar la subrutina en sí.

Respuesta:

No estoy muy familiarizado con las funciones internas de f2py, pero estoy muy familiarizado con la envoltura de Fortran. F2py simplemente automatiza algunas o todas las cosas a continuación.

Primero debe exportar a C usando el módulo iso_c_binding, como se describe, por ejemplo, aquí:

http://fortran90.org/src/best-practices.html#interfacing-with-c

Descargo de responsabilidad: soy el autor principal de las páginas de fortran90.org. Esta es la única forma independiente de la plataforma y el compilador de llamar a Fortran desde C. Esto es F2003, por lo que en estos días no hay razón para usar otra forma.
Solo puede exportar / llamar matrices con la longitud completa especificada (forma explícita), es decir:
```
 integer(c_int), intent(in) :: N real(c_double), intent(out) :: mesh(N)
```
pero no asumir forma:
```
 real(c_double), intent(out) :: mesh(:)
```
Eso se debe a que el lenguaje C no admite tales matrices en sí. Se habla de incluir dicho soporte en F2008 o posterior (no estoy seguro), y la forma en que funcionaría es a través de algunas estructuras de datos C de soporte, ya que necesita llevar información de forma sobre la matriz.

En Fortran, debe usar principalmente la forma de asumir, solo en casos especiales debe usar la forma explícita, como se describe aquí:

http://fortran90.org/src/best-practices.html#arrays

Eso significa que necesita escribir un contenedor simple alrededor de su subrutina de forma asumida, que envolverá las cosas en matrices de formas explícitas, según mi primer enlace anterior.
Una vez que tenga una firma C, simplemente llámela desde Python de la forma que desee, yo uso Cython, pero puede usar ctype o C / API a mano.
La instrucción deallocate(y) no es necesaria, Fortran desasigna automáticamente.

http://fortran90.org/src/best-practices.html#allocatable-arrays
real*8 no debe usarse, sino real(dp) :

http://fortran90.org/src/best-practices.html#floating-point-numbers
La declaración y(:) = 1.0 está asignando 1.0 con precisión simple, por lo que el resto de dígitos serán aleatorios. Este es un error común:

http://fortran90.org/src/gotchas.html#floating-point-numbers

Necesita usar y(:) = 1.0_dp .
En lugar de escribir y(:) = 1.0_dp , puede escribir y = 1 , eso es todo. Puede asignar un número entero a un número de punto flotante sin perder precisión, y no necesita poner el redundante (:) allí. Mucho más simple.
En vez de
```
 y = 1 z = x + y
```
Solo usa
```
 z = x + 1
```
y no se moleste en absoluto con la matriz y .
No necesita la instrucción "return" al final de la subrutina.
Finalmente, probablemente debería usar módulos, y simplemente poner implicit none en el nivel del módulo y no necesita repetirlo en cada subrutina.

De lo contrario, me parece bien. Aquí está el código siguiendo las sugerencias 1-10 anteriores:
```
 module test use iso_c_binding, only: c_double, c_int implicit none integer, parameter :: dp=kind(0.d0) contains subroutine f(x, z) real(dp), intent(in) :: x(:) real(dp), intent(out) :: z(:) z = x + 1 end subroutine subroutine c_f(n, x, z) bind(c) integer(c_int), intent(in) :: n real(c_double), intent(in) :: x(n) real(c_double), intent(out) :: z(n) call f(x, z) end subroutine end module
```
Muestra la subrutina simplificada así como un contenedor C.

En cuanto a f2py, probablemente intente escribir este contenedor por usted y falle. Tampoco estoy seguro de si está utilizando el módulo iso_c_binding . Entonces, por todas estas razones, prefiero envolver las cosas a mano. Entonces está exactamente claro lo que está sucediendo.

Pregunta:

Respuesta:

Leave a Comment Cancel Reply