2013-07-20 67 views
2

我遇到了一個簡單的代碼下面的問題。我正在嘗試與GFortran一起使用OpenMP。 x的代碼結果應該與AND沒有!$OMP語句相同,因爲並行代碼和串行代碼應輸出相同的結果。OpenMP工作精度

program test 
implicit none 
!INCLUDE 'omp_lib.h' 
integer i,j 
Real(8) :: x,t1,t2 

x=0.0d0 
!$OMP PARALLEL DO PRIVATE(i,j) shared(X) 
Do i=1,3 
    Write(*,*) I 
    !pause 
    Do j=1,10000000 
    !$OMP ATOMIC 
    X=X+2.d0*Cos(i*j*1.0d0) 
    end do 
end do 
!$OMP END PARALLEL Do 

write(*,*) x  
end program test 

但奇怪的是我得到以下結果x

並行:-3.17822355415XXXXX

編號:-3.1782235541569084

其中XXXXX是一些隨機數字。每次運行串行代碼時,我都會得到相同的結果(-3.1782235541569084)。我該如何解決它?這是由於某些OpenMP工作精度選項造成的問題嗎?

回答

4

浮點算術不嚴格聯想。在f-p算術中,既不是a+(b+c)==(a+b)+c也不是a*(b*c)==(a*b)*c總是如此,因爲它們都是實數算術。這是衆所周知的,並且廣泛解釋了對SO和其他網絡上其他有聲地方的其他問題的回答。這裏我不會詳細闡述這一點。

由於您已經編寫了程序,因此計算最終值X的操作順序是非確定性的,也就是說它可能(也可能確實)因執行而異。 atomic指令一次只允許一個線程更新X,但它不會對到達指令的線程施加任何排序約束。

考慮到程序中計算的性質,我相信你在串行執行和並行執行之間所看到的差異可能完全被這種非確定性所解釋。

在你考慮'解決'這個問題之前,你應該首先確定它是一個問題。是什麼讓你認爲這個串行代碼的答案是一個真正的答案?如果你要向後(依次)運行循環,並得到不同的答案(很可能),那麼你正在尋找哪個答案?在許多科學計算中,這可能是OpenMP的核心領域,可用的數據和所使用的數字方法不支持超出少數有效數字的程序結果準確性的斷言。

如果您仍然認爲這是一個需要解決的問題,最簡單的方法是簡單地取出OpenMP指令。

+0

我意識到你在說什麼。我在我的問題中展示的代碼僅僅是我在真實代碼中面臨的一個示例,這些代碼執行了更多計算。我正在編寫一個代碼來爲我計算流體動力學(CFD)開發的新方法學進行基準測試。所以這就是爲什麼我關心這個問題,因爲我會盡可能地使用最大精度數字來比較和基準方法。 – Eleteroboltz

1

要增加高性能標記所說的,另一個差異來源是編譯器可能發出了x87 FPU指令來完成數學運算。 x87使用80位內部精度,優化後的串行代碼在實際將最終值寫入X的內存位置之前僅使用寄存器算術。在並行的情況下,由於X是一個共享變量,因此在每次迭代時都會更新內存位置。這意味着80位x87 FPU寄存器被刷新到64位存儲單元,然後回讀,因此在每次迭代時都會丟失一些精度位,這會增加觀察到的差異。

如果現代64位CPU與發出SIMD指令的編譯器一起使用,則不會出現這種效果,例如, SSE2 +或AVX。這些只能用於64位內部精度,然後僅使用寄存器尋址不會比在每次迭代中刷新和重新加載內存值時的精度更高。在這種情況下,差異來自高性能標記解釋的非關聯性。

這些影響是相當多的預期和通常考慮。他們被充分研究和理解,如果你的CFD算法在平行運行時發生故障,那麼算法在數值上是非常不穩定的,即使在串行情況下,我也絕不會相信它給出的結果。

順便說一句,一個更好的方式來實現你的循環是使用還原:

!$OMP PARALLEL DO PRIVATE(j) REDUCTION(+:X) 
Do i=1,3 
    Write(*,*) I 
    !pause 
    Do j=1,10000000 
     X=X+2.d0*Cos(i*j*1.0d0) 
    end do 
end do 

這將允許編譯器爲每個線程註冊優化的代碼,然後精度損失不僅會在線程總和它們的局部部分值以便獲得最終值X時在最後出現。

+0

謝謝,Hristo Iliev。我將用'reduction'子句實現算法。我的算法沒有崩潰,而且非常穩定。我注意到的唯一問題是我在問題中指出的問題。我認爲這是由於高性能標記所說的F-P算術非相關性。我試圖用'Real(16)'而不是'Real(8)',並且我得到的結果對於串行和並行實現是相同的。 – Eleteroboltz

+0

其實,我不能使用'Reduction',因爲在我的真實代碼中變量是一個非常大的數組。我在「使用OpenMP編寫Fortran 95中的並行編程」一書中讀到,「變量''x'受'Reduction'子句影響,必須具有標量性質和內部類型。儘管'x'可能是標量項一個非 延遲形狀或假定大小的數組如果是這種情況,如果'x'是一個大數組,則由於變量的私有副本的創建和初始化而導致的計算開銷可能非常大。 – Eleteroboltz

+0

我會測試你引用的語句的有效性。*使用OpenMP *的Fortran 95中的並行編程,它根據編譯器(版本),平臺,實現細節以及各種各樣的東西來判斷「真相」。 –