f2py，将Python函数传递给Fortran时出现问题

Question

f2py，将Python函数传递给Fortran时出现问题

3

I have a this simple Fortran code (stack.f90):

      subroutine fortran_sum(f,xs,nf,nxs)
      integer nf,nxs
      double precision xs,result
      dimension xs(nxs),result(nf)
      external f 
      result = 0.0
      do I = 1,nxs
         result = result + f(xs(I))
         print *,xs(I),f(xs(I))
      enddo
      return
      end

我正在使用以下方式进行编译：

f2py -c --compiler=mingw32 -m stack2 stack2.f90

然后使用这个 Python 脚本（stack.py）进行测试：

import numpy as np
from numpy import cos, sin , exp
from stack import fortran_sum
def func(x):
    return x**2

if __name__ == '__main__':
    xs = np.linspace(0.,10.,10)
    ans =  fortran_sum(func,xs,1)
    print 'Fortran:',ans
    print 'Python:',func(xs).sum()

当我使用 "python stack.py" 运行时，会得到以下结果：

   0.0000000000000000        0.00000000
   1.1111111111111112              Infinity
   2.2222222222222223              Infinity
   3.3333333333333335        9.19089638E-26
   4.4444444444444446              Infinity
   5.5555555555555554        0.00000000
   6.6666666666666670        9.19089638E-26
   7.7777777777777786        1.60398502E+09
   8.8888888888888893              Infinity
   10.000000000000000        0.00000000
Fortran: None
Python: 351.851851852

我的问题是：

为什么函数没有被正确评估？
如何将result返回给Python？
在Fortran中是否可能一次性评估数组xs？

谢谢！

编辑：通过@SethMMorton的极好提示，我最终得到了以下结果：

      subroutine fortran_sum(f,xs,nxs,result)
      implicit none
      integer :: I
      integer, intent(in) :: nxs
      double precision, intent(in) :: xs(nxs)
      double precision, intent(out) :: result
      double precision :: f
      external :: f 
      ! "FIX" will be added here
      result = 0.0
      do I = 1,nxs
        result = result + f(xs(I))
        print *,xs(I),f(xs(I))
      enddo
      return
      end

使用修改后的命令来运行 stack.py： ans = fortran_sum(func,xs)，结果如下：

   0.0000000000000000        0.0000000000000000
   1.1111111111111112        3.8883934247189009E+060
   2.2222222222222223        3.8883934247189009E+060
   3.3333333333333335        9.1908962428537221E-026
   4.4444444444444446        3.8883934247189009E+060
   5.5555555555555554        5.1935286092977251E-060
   6.6666666666666670        9.1908962428537221E-026
   7.7777777777777786        1603984978.1728516
   8.8888888888888893        3.8883934247189009E+060
   10.000000000000000        0.0000000000000000
Fortran: 1.55535736989e+61
Python: 351.851851852

这是错误的。如果我添加一个中间变量 x=x(I) 并使用这个变量 f(x) 调用函数，就不会出现这种奇怪的行为。有趣的是，如果我调用一次 f(x)，所需的调用 f(x(I)) 也可以正常工作。应用此“修复”后：

double precision :: x, f_dummy
x = xs(I)
f_dummy = f(x)

编译和运行后，可以得到正确的结果：

   0.0000000000000000        0.0000000000000000
   1.1111111111111112        1.2345679012345681
   2.2222222222222223        4.9382716049382722
   3.3333333333333335        11.111111111111112
   4.4444444444444446        19.753086419753089
   5.5555555555555554        30.864197530864196
   6.6666666666666670        44.444444444444450
   7.7777777777777786        60.493827160493836
   8.8888888888888893        79.012345679012356
   10.000000000000000        100.00000000000000
Fortran: 351.851851852
Python: 351.851851852

希望有人能够解释一下为什么？

- Saullo G. P. Castro

如果您在Python函数中添加了一个打印语句，它是否显示您所期望的内容？我唯一可能建议的是使用调试标志进行编译，以查看是否发生了某些奇怪的事情，这是Fortran会通过正确的标志捕获到的。 - SethMMorton

谢谢你的建议。使用 f2py 进行编译没有显示任何警告... 我现在会接受这个“FIX”，这可能是由于 f2py 导致的，太过神秘了... - Saullo G. P. Castro

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- SethMMorton · Accepted Answer

为什么函数的计算结果不正确？

你直接将回调函数 f 的结果发送给了 print 方法。由于Fortran不知道 f 将返回的数据类型（因为它没有声明，也不是一个Fortran函数），print 无法正确地解释数据以正确打印。如果你将结果赋值给一个变量，然后打印该变量，你应该会看到预期的结果：

double precision x
....
x = f(xs(1))
print *, x

如何将 result 返回给 Python？

您需要告诉 f2py 变量的意图。这实际上可以使用标准的 Fortran 90 语法来完成，我强烈建议这样做，因为这可以使您的代码更清晰（我知道您正在使用 Fortran 90，因为您能够打印整个数组而不需要循环遍历它）。

subroutine stack2(f,xs,result,nf,nxs)
    implicit none
    integer,          intent(in)  :: nf, nxs
    double precision, intent(in)  :: xs(nxs)
    double precision, intent(out) :: result(nf)
    external f
    ...
end subroutine stack2

现在，哪些变量进入例程，哪些变量出去已经很清楚了。请注意，为了传递出去，result必须成为子程序的一个参数。 f2py现在将理解应该从stack2函数返回result。

是否可能一次性在Fortran中评估数组xs？

我不确定您确切的意思，但我假设您想知道是否可以一次对整个数组进行操作，而不是在do循环中。

可能可以。

由于xs是一个数组，您应该只需将整个数组传递给f，它应该对整个数组进行操作。这可能无法正常工作，因为Fortran要求函数是ELEMENTAL才能这样做，而这可能不在f2py的范围内。如果f2py没有使函数成为ELEMENTAL，则必须使用循环来完成此操作。

您可能需要处理的一个问题

行result = result + f(xs(I))将相同的值分配给result的每个元素。如果这不是您想要的，那么就需要解决这个问题。

基于代码的概述

下面是使用这些新建议的版本的示例代码。

subroutine stack2(f,xs,result,nf,nxs)
    implicit none
    integer,          intent(in)  :: nf, nxs
    double precision, intent(in)  :: xs(nxs)
    double precision, intent(out) :: result(nf)
    double precision :: x
    integer          :: I
    external f
    result = 0.0d0 ! Make this a double constant since result is double
    do I = 1,nxs
        x = f(xs(I))
        result = result + x
        print *, x
    end do
    print *, xs
    return ! Not needed in Fortran 90
end subroutine stack2

请注意，f2py会正确计算输入数组的长度，因此当您调用它时，只需要定义result的长度即可：

import numpy as np
from stack2 import stack2
def func(x):
    return x**2

if __name__ == '__main__':
    xs = np.linspace(0.,10.,10)
    ans =  stack2(func,xs,5) # This was changed
    print 'ans:',ans