3

我有一个简单的 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))也有效。应用此“修复”后：

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

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

4

1 回答 1

6

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

您正在将回调函数的结果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
于 2013-07-03T22:47:03.847 回答