F2Py z alokowanymi i zakładanymi tablicami kształtów

Chciałbym używać f2pyz nowoczesnym Fortranem. W szczególności próbuję uzyskać następujący podstawowy przykład do działania. To jest najmniejszy użyteczny przykład, jaki mogłem wygenerować.

! 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

Zauważ, że nwynika to z kształtu parametru wejściowego x. Należy pamiętać, że yjest on przydzielany i zwalniany w treści podprogramu.

Kiedy to skompiluję f2py

f2py -c alloc_test.f90 -m alloc

A następnie uruchom w Pythonie

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

Pojawia się następujący błąd

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

Więc idę pyfręcznie tworzyć i edytować plik

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

Oryginalny

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

Zmodyfikowano

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

Teraz działa, ale wartości wyjściowe zsą zawsze 0. Niektóre drukowanie debugowania ujawnia, że nma wartość 0w podprogramie f. Zakładam, że brakuje mi f2pymagii nagłówka, aby właściwie zarządzać tą sytuacją.

Mówiąc bardziej ogólnie, jaki jest najlepszy sposób połączenia powyższego podprogramu z Pythonem? Zdecydowanie wolałbym nie modyfikować samego podprogramu.

Nie jestem zbyt obeznany z wewnętrznymi elementami f2py, ale doskonale znam się na pakowaniu Fortrana. F2py po prostu automatyzuje niektóre lub wszystkie rzeczy poniżej.

1. Najpierw musisz wyeksportować do C za pomocą modułu iso_c_binding, jak opisano na przykład tutaj:

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

   Oświadczenie: Jestem głównym autorem stron fortran90.org. Jest to jedyny niezależny od platformy i kompilatora sposób wywoływania Fortran z C. To jest F2003, więc w dzisiejszych czasach nie ma powodu, aby używać innego sposobu.

2. Możesz eksportować / wywoływać tylko tablice o pełnej długości (wyraźny kształt), to znaczy:

   integer(c_int), intent(in) :: N
   real(c_double), intent(out) :: mesh(N)

   ale nie przyjmują kształtu:

   real(c_double), intent(out) :: mesh(:)

   Wynika to z faktu, że język C nie obsługuje takich tablic. Mówi się o włączeniu takiego wsparcia w F2008 lub nowszym (nie jestem pewien), a sposób, w jaki działałby, to poprzez niektóre wspierające struktury danych C, ponieważ musisz przenosić informacje o kształcie o tablicy.

   W Fortranie powinieneś używać głównie zakładania kształtu, tylko w szczególnych przypadkach powinieneś używać wyraźnego kształtu, jak opisano tutaj:

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

   Oznacza to, że musisz napisać prostą owijkę wokół podprogramu „Przyjmij kształt”, który zawinie rzeczy w jawne tablice kształtów, za pomocą mojego pierwszego linku powyżej.

3. Gdy masz już podpis C, po prostu wywołaj go z Pythona w dowolny sposób, używam Cython, ale możesz ręcznie używać ctype lub C / API.

4. deallocate(y)Instrukcja nie jest potrzebna, Fortran dealokuje automatycznie.

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

5. real*8nie należy go używać, ale real(dp):

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

6. Instrukcja y(:) = 1.0przypisuje 1.0 z pojedynczą precyzją, więc reszta cyfr będzie losowa! Jest to częsta pułapka:

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

   Musisz użyć y(:) = 1.0_dp.

7. Zamiast pisać y(:) = 1.0_dp, możesz po prostu pisać y = 1, to wszystko. Możesz przypisać liczbę całkowitą do liczby zmiennoprzecinkowej bez utraty dokładności i nie musisz umieszczać tam nadmiaru (:). O wiele prostsze.

8. Zamiast

   y = 1
   z = x + y

   po prostu użyj

   z = x + 1

   i nie zawracaj sobie głowy ytablicą.

9. Nie potrzebujesz instrukcji „return” na końcu podprogramu.

10. Wreszcie, prawdopodobnie powinieneś używać modułów i po prostu umieścić implicit none na poziomie modułu i nie musisz powtarzać go w każdym podprogramie.

    W przeciwnym razie wygląda mi to dobrze. Oto kod zgodny z powyższymi sugestiami 1-10:

    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

    Pokazuje uproszczony podprogram, a także opakowanie typu C.

    Jeśli chodzi o f2py, prawdopodobnie próbuje napisać to opakowanie dla ciebie i kończy się niepowodzeniem. Nie jestem też pewien, czy używa iso_c_bindingmodułu. Z tych wszystkich powodów wolę owijać rzeczy ręcznie. Wtedy jest dokładnie jasne, co się dzieje.

Wszystko, co musisz zrobić, to:

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

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

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

! Variable Initializations !
  allocate(y(n))

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

  deallocate(y)
  return
end subroutine f

Chociaż rozmiar tablicy xiz jest teraz przekazywany jako jawny argument, f2py sprawia, że ​​argument n jest opcjonalny. Poniżej znajduje się dokumentacja funkcji, tak jak wygląda na python:

Type:       fortran
String Form:<fortran object>
Docstring:
f - Function signature:
  z = f(x,[n])
Required arguments:
  x : input rank-1 array('d') with bounds (n)
Optional arguments:
  n := len(x) input int
Return objects:
  z : rank-1 array('d') with bounds (n)

Importowanie i wywoływanie go z Pythona:

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

daje następujący wynik:

[ 2.  2.  2.  2.  2.]