f2py,将Python函数传递给Fortran的问题

Question

我有一个这个简单的Fortran代码(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

我的问题是:

• 为什么函数没有被正确评估？

• 如何返回resultPython？

• 是否有可能xs在Fortran中一次评估阵列？

谢谢!

---

编辑:有了@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))也可以.应用此"FIX"后:

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

如果有人能解释为什么会好的吗？

Answer 1

为什么没有正确评估函数？

您将回调函数的结果f直接发送到print方法.因为Fortran不知道f将返回什么数据类型(因为它未声明且不是Fortran函数),print所以无法正确解释要正确打印的数据.如果将结果分配给变量,则打印变量,您应该看到预期的结果:

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

如何返回resultPython？

您需要通知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现在将了解result应该从stack2函数返回.

是否可以xs在Fortran中一次评估数组

我不确定你的意思究竟是什么,但我想你想知道你是否可以一次在整个数组上执行此操作,而不是在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