可分配数组的f2py错误

Question

我有一个Fortran子程序,我想在Python中使用.

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,然后相应地分配用于阵列中的空间bins和hist.

虽然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不能容忍在运行时为数组分配空间(不知道为什么,因为这似乎是一种自然的科学需求).

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

Answer 1

您可以使用一些有效的选项/工作轮.

1.最简单的可能是安排python调用该函数n_bins,然后使用其中的结果调用(稍微修改过的版本)函数,并使用create_hist正确大小的不可分配输出数组.

即在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

与Fortran接口create_hist现在定义为

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

这只适用于您可以n_bins便宜地从外面打电话的情况create_hist.我知道有两种模拟可分配数组的方法,适用于不适用的情况(即计算数组大小的代码很昂贵而且不容易分开).

2.第一种是使用模块级可分配数组(在此处的文档中描述).这本质上是一个可分配的全局变量 - 你调用你的函数,它将数据保存到全局变量中,然后从Python访问它.缺点是它不是线程安全的(即如果你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.我非常喜欢的另一种方法是仅在函数内分配数组(即不将它们作为参数传入/传出).但是,你传递的是一个Python回调函数,它在最后被调用并保存数组.这是线程安全的(我相信).

然后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(这是为了构建一个名为的Python模块fortran_module- 根据需要更改名称),然后修改它的相关行以阐明回调函数的维度:

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,因为没有多线程陷阱要避免(但实际上这是一个你可能永远看不到的角落情况).选项2更容易实现.

Answer 2

据我所知， 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 位。