可分配数组的 f2py 错误

Question

我有一个要在 Python 中使用的 Fortran 子例程。

subroutine create_hist(a, n, dr, bins, hist)
    integer, intent(in) :: n
    real(8), intent(in) :: a(n)
    real(8), intent(in) :: dr
    integer, intent(out), allocatable :: hist(:)
    real(8), intent(out), allocatable :: bins(:)

    n_b = n_bins(a, n, dr)  ! a function calculating the number of bins
    allocate(bins(n_b+1))
    allocate(hist(n_b))
    hist = 0

    !... Create the histogram hist by putting elems of a into bins
end subroutine

这是一个简单的程序，用于获取数字数组 a 并根据给定的 bin 大小 dr 创建直方图。首先，它使用函数 n_bins 获取 bin 的数量，然后相应地为数组 bins 和 hist.

分配 space

虽然 gfortran 可以很好地编译此代码，但 f2py 会引发错误：

/tmp/tmpY5badz/src.linux-x86_64-2.6/mat_ops-f2pywrappers2.f90:32.37:
      call create_hist_gnu(a, n, dr, bins, hist)
Error: Actual argument for 'bins' must be ALLOCATABLE at (1)

据我了解，f2py 不允许在运行时为数组分配 space（不知道为什么，因为这似乎是一种自然的科学需求）。

任何人都可以建议一种在运行时间分配 Fortran 数组的方法，以便 f2py 对此感到满意吗？

Answer 1

据我所知，f2py 不支持具有属性 ALLOCATABLE 的虚拟参数（有时称为函数参数）。另一种方法是使它成为一个足够大的显式形状数组。然后您将只使用数组的特定部分。

subroutine create_hist(a, n, dr, n_b_max, bins, hist)
    integer, intent(in) :: n
    real(8), intent(in) :: a(n)
    real(8), intent(in) :: dr
    integer, intent(in) :: n_b_max
    integer, intent(out) :: hist(n_b_max)
    real(8), intent(out) :: bins(n_b_max+1)

顺便说一句，使用 real(8) 不是指定 64 位精度的可移植方式，并且在某些编译器中会失败。即使是旧的 double precision 也更好，因为它总是会编译成某种东西，而且通常会编译成 64 位。

Answer 2

您可以使用一些有效的 options/work-rounds。

1. 可能最简单的方法是安排 python 调用函数 n_bins，然后使用该函数的结果调用 ( create_hist 函数的略微修改版本，具有正确大小的 non-allocatable 输出数组。

即在 Python:

n_b = n_bins(a,dr) # don't need to pass n - inferred from a
bins,hist = create_hist(a,dr,n_b) # bins and hist are auto-allocated

create_hist 的 Fortran 接口现在定义为

subroutine create_hist(a, n, dr, n_b, bins, hist)
    integer, intent(in) :: n
    real(8), intent(in) :: a(n)
    real(8), intent(in) :: dr
    integer, intent(in) :: n_b
    integer, intent(out) :: hist(n_b)
    real(8), intent(out) :: bins(n_b+1)

    ! code follows

! you also need to specify the function n_bins...

这仅适用于您可以从外部 create_hist 廉价调用 n_bins 的情况。我知道有 2 种模拟可分配数组的方法，用于不适用的情况（即计算数组大小的代码很昂贵并且不容易分离）。

2. 第一种是使用模块级可分配数组（在文档here中有描述）。这本质上是一个可分配的全局变量——你调用你的函数，它将数据保存到全局变量中，然后你从 Python 访问它。缺点是它不是 thread-safe（即，如果您同时并行调用 create_hist 是不好的）

module something
    real(8), allocatable :: bins(:)
    integer, allocatable :: hist(:)
contains
    subroutine create_hist(a,n,dr)
        integer, intent(in) :: n
        real(8), intent(in) :: a(n)
        real(8), intent(in) :: dr
        integer :: n_b

        n_b = n_bins(a,n,dr)

        allocate(bins(n_b+1))
        allocate(hist(n_b))
        ! code follows
    end subroutine
end module

Python 调用看起来像

something.create_hist(a,n,dr)
bins = something.bins # or possible bins = copy.copy(something.bins)
hist = something.hist # or possible hist = copy.copy(something.hist)

3. 我非常喜欢的另一种方法是只在函数内分配数组（即不要将它们作为参数传递 in/out）。但是，您传递的是一个 Python 回调函数，该函数在最后调用并保存数组。这是thread-safe（我相信）。

fortran 代码看起来像

subroutine create_hist(a,n,dr,callback)
    integer, intent(in) :: n
    real(8), intent(in) :: a(n)
    real(8), intent(in) :: dr
    external callable ! note: better to specify the type with an interface block (see http://www.fortran90.org/src/best-practices.html#callbacks)
    integer :: n_b
    real(8), allocatable :: bins(:)
    integer, allocatable :: hist(:)

    n_b = n_bins(a,n,dr)

    allocate(bins(n_b+1))
    allocate(hist(n_b))
    ! code follows
    call callable(bins,hist,n_b)
end subroutine

不幸的是，它涉及得更多了。您需要使用命令 f2py -m fortran_module -h fortran_module.pyf my_fortran_file.f90 创建一个签名文件（这是构建一个名为 fortran_module 的 Python 模块 - 根据需要更改名称），然后将其中的相关行修改为阐明回调函数的维度：

python module create_hist__user__routines 
    interface create_hist_user_interface 
        subroutine callable(bins,hist,n_b) ! in :f:my_fortran_file.f90:stuff:create_hist:unknown_interface
            real(kind=8), dimension(n_b+1) :: bins
            integer, dimension(n_b) :: hist
            integer :: n_b
        end subroutine callable
    end interface create_hist_user_interface
end python module create_hist__user__routines

用 f2py -c fortran_module.pyf my_fortran_file.f90

编译

然后是一个简短的 Python 包装器（为您提供一个简单的界面），看起来像

def create_hist(a,dr):
    class SaveArraysCallable(object):
        def __call__(self,bins,hist):
            self.bins = bins.copy()
            self.hist = hist.copy()

    f = SaveArrayCallable()
    fortran_module.create_hist(a,dr,f)
    return f.bins, f.hist

总结

在很多情况下，选项 1 可能是最好的。否则选项 2 或选项 3。我更喜欢选项 3，因为没有要避免的多线程陷阱（但实际上这是一个您可能永远不会看到的 corner-case）。选项 2 更容易实施。

可分配数组的 f2py 错误

f2py error with allocatable arrays

python

fortran

f2py