gpt4 book ai didi

python - f2py,将 Python 函数传递给 Fortran 时出现问题

转载 作者:太空宇宙 更新时间:2023-11-03 14:27:41 24 4
gpt4 key购买 nike

我有一个简单的 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

我的问题是:

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

  • 如何将结果返回给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

如果有人能解释一下为什么会很好?

最佳答案

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

您正在将回调函数 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 现在将理解 result 应该从 stack2 函数返回。

是否可以在 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

关于python - f2py,将 Python 函数传递给 Fortran 时出现问题,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/17458111/

24 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com