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
希望有人能够解释一下为什么?
f2py
进行编译没有显示任何警告... 我现在会接受这个“FIX”,这可能是由于f2py
导致的,太过神秘了... - Saullo G. P. Castro