F2Py с размещаемыми и предполагаемыми массивами форм

18

Я хотел бы использовать f2pyс современным Fortran. В частности, я пытаюсь заставить работать следующий базовый пример. Это самый маленький полезный пример, который я мог привести.

! 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

Обратите внимание, что nвыводится из формы входного параметра x. Обратите внимание, что yэто распределяется и освобождается в теле подпрограммы.

Когда я собираю это с f2py

f2py -c alloc_test.f90 -m alloc

А затем запустить в Python

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

Я получаю следующую ошибку

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

Поэтому я иду и создаю и редактирую pyfфайл вручную

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

оригинал

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

модифицированный

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

Теперь он работает, но значения вывода zвсегда 0. Некоторая отладочная печать показывает, что nимеет значение 0в подпрограмме f. Я предполагаю, что мне не хватает некоторой f2pyмагии заголовка, чтобы правильно управлять этой ситуацией.

В более общем смысле, как лучше связать вышеуказанную подпрограмму с Python? Я бы предпочел не изменять саму подпрограмму.

MRocklin
источник
Мэтт, вы знакомы с руководством по передовым практикам Ондрея Certik, в частности, раздел Интерфейс с Python ? Мы обсуждали аналогичную проблему с интерфейсом PyClaw и пока не решили ее :)
Aron Ahmadia

Ответы:

23

Я не очень знаком с внутренностями f2py, но я очень хорошо знаком с упаковкой Фортрана. F2py просто автоматизирует некоторые или все вещи ниже.

  1. Сначала вам нужно экспортировать в C, используя модуль iso_c_binding, как описано, например, здесь:

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

    Отказ от ответственности: я являюсь основным автором страниц fortran90.org. Это единственный независимый от платформы и компилятора способ вызова Fortran из C. Это F2003, поэтому в наши дни нет никаких причин использовать какой-либо другой способ.

  2. Вы можете экспортировать / вызывать только массивы с полной длиной (явная форма), то есть:

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

    но не принимать форму:

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

    Это потому, что язык C не поддерживает такие массивы сам по себе. Говорят о том, чтобы включить такую ​​поддержку в F2008 или более позднюю версию (я не уверен), и способ, которым она будет работать, заключается в поддержке структур данных C, так как вам нужно нести информацию о форме массива.

    В Фортране вы должны в основном использовать предположительную форму, только в особых случаях вы должны использовать явную форму, как описано здесь:

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

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

  3. Если у вас есть подпись C, просто вызывайте ее из Python любым удобным для вас способом, я использую Cython, но вы можете использовать ctype или C / API вручную.

  4. deallocate(y)Заявление не требуется, Fortran автоматически освобождает.

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

  5. real*8не следует использовать, а real(dp):

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

  6. Оператор y(:) = 1.0присваивает 1.0 с одинарной точностью, поэтому остальные цифры будут случайными! Это распространенная ошибка:

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

    Вам нужно использовать y(:) = 1.0_dp.

  7. Вместо того, чтобы писать y(:) = 1.0_dp, вы можете просто написать y = 1, вот и все. Вы можете назначить целое число на число с плавающей запятой без потери точности, и вам не нужно помещать (:)туда избыточное число . Гораздо проще

  8. Вместо того

    y = 1
    z = x + y

    просто используйте

    z = x + 1

    и не беспокойтесь с yмассивом вообще.

  9. Вам не нужно выражение "return" в конце подпрограммы.

  10. Наконец, вы, вероятно, должны использовать модули и просто поставить implicit noneна уровень модулей, и вам не нужно повторять это в каждой подпрограмме.

    В противном случае это выглядит хорошо для меня. Вот код, следующий за предложениями 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

    Он показывает упрощенную подпрограмму, а также оболочку на языке C.

    Что касается f2py, он, вероятно, пытается написать эту обертку для вас и терпит неудачу. Я также не уверен, использует ли он iso_c_bindingмодуль. Поэтому по всем этим причинам я предпочитаю оборачивать вещи вручную. Тогда точно понятно, что происходит.

Ондржей Чертик
источник
Насколько я знаю, f2py не использует привязки ISO C (его основной целью является код Fortran 77 и код Fortran 90).
Арон Ахмадиа
Я знал, что был немного глуп, yно я хотел, чтобы что-то было выделено (мой реальный код имеет нетривиальное распределение). Я не знал о многих других моментах, хотя. Похоже, мне стоит заглянуть в какое-нибудь руководство по лучшим практикам Fortran90 .... Спасибо за подробный ответ!
MRocklin
Обратите внимание, что используя современные компиляторы Fortran, вы оборачиваете F77 точно таким же образом - записывая простую оболочку iso_c_binding и вызывая из нее устаревшую подпрограмму f77.
Ондржей Чертик
6

Все, что вам нужно сделать, это следующее:

!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

Хотя размер массива x и z теперь передается в качестве явного аргумента, f2py делает аргумент n необязательным. Ниже приведена строка документации функции в том виде, как она выглядит в 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)

Импортировать и вызывать его из python:

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

дает следующий вывод:

[ 2.  2.  2.  2.  2.]
Prometheous
источник
Есть ли способ использовать какое-то нетривиальное выражение в качестве размера? Например, я передаю nи хочу получить массив размера 2 ** n. До сих пор я должен также передать 2 ** n в качестве отдельного аргумента.
Alleo