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 

我的問題是：

• 爲什麼函數沒有被正確評估？

• 如何將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.2345679
    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.
    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 

如何返回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函數返回。

是否有可能在陣列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