有限差分方法、偽譜法的CUDA C實現以及計算效率比較

4 人贊了文章

在FWI、RTM、以及LSRTM等等所謂的高端地震數據成像演算法的計算核心其實是三維波動方程的數值計算問題，這部分計算是這些演算法的性能熱點，它們通常會採用CUDA C來實現，利用GPU的超高浮點性能來提高演算法的計算效率。

在2008年到2012年，高性能計算領域對三維波動方程FDTD的GPU實現進行了比較多的研究，這種PDE的數值離散求解在高性能計算領域又被稱為stencil computation。有一個實現方式可以在NVIDIA提供的CUDA samples里找到，也就是：

1. 採用二維的線程網格，在x-o-y平面進行區域劃分，每個線程block負責一個二維區域的計算，每個線程的計算沿著z方向一層一層的進行。
2. 每個線程block利用共享內存存放計算x-方向和y-方向二階偏導數所需要的網格點。
3. 計算z方向的偏導數所需要的網格點會通過一個寄存器隊列進行預取。

這也就是所謂的2.5D blocking的實現方式。最近有一些比較新的文章會採用3.5D blocking的方式來實現三維的FDTD，思路就是2.5D blocking + 時間軸的blocking，這樣會減少CUDA kernel函數的調用開銷，同時會提高block內的指令級別的並行度(ILP)。

提高ILP（instruction level parallelism）是使得CUDA kernel的性能接近峰值的非常重要的技巧。這個在Vasili Volkov的報告上http://www.nvidia.com/content/GTC-2010/pdfs/2238_GTC2010.pdf有提到。當然這個是非常基本的優化技術啊，高手們就不要吐槽了。

所以我想在CUDA samples提供的例子的基礎上，可以通過提高ILP的技巧進一步的提高kernel函數的性能。也就是說呢，在原來的kernel函數裡面，每個線程block計算的區域大小是和線程block的大小完全一樣的。但是，可以讓一個線程block計算更大的區域，這樣的話，在同一個線程block內，會有更多的獨立的訪存指令和浮點運算指令，指令的延遲(latency)能夠被隱藏，kernel的性能也就會提高。

現在我把CUDA samples里的實現方式寫成了C++的模板函數，每個線程block計算區域的大小BX和BY作為模板函數的參數，同理呢，線程塊的大小TX和TY也作為模板函數的參數，然後利用grid search搜索(BX, BY, TX, TY)的最佳配置，就能使kernel函數的性能優化了。經過一系列的實驗，優化後的kernel確實比原先的kernel快了不少。512×512×512的網格，計算100步，時間大概是4秒鐘，在某塊陳舊的Fermi架構的卡上達到了接近4000Mcells/s的性能，原先的實現的話大概在1500Mcell/s左右，所以還是有一定效果的。

接下來說一下偽譜法，偽譜法就是計算波動方程的拉普拉斯運算元採用FFT來實現，而不是做差分，所以實現和優化起來會比差分要簡單，因為我們的FFT是通常會用庫來做。除開FFT，偽譜法就很簡單了，只需要實現一個point wise相乘的CUDA kernel，這樣的kernel要寫得比較好還是挺簡單的。演算法完成之後，進行了測試，和FDTD3D同樣的網格大小，時間大概在18秒左右，性能在700Mcells/s左右。然後我做了簡單的profiling，偽譜法的性能熱點和想像的一樣，集中在三維FFT的部分。因為FFT在最慢的一維上，要頻繁地跨越很大的stride訪問全局存儲，所以慢也是可以預見的吧。

所以結論是在同樣的網格規模，同樣的時間步數下面，偽譜法的GPU實現大概要比FDTD3D要慢4-5倍。所以偽譜法必須在更大的網格步長下計算，才能體現出他的優勢。

最後呢就是這些日子，利用空閑的時間，研究了一下提升CUDA kernel性能的一些技巧吧，並把這些技巧應用到了一些數值演算法上(比如：矩陣乘法和stencil computation)，算是對GPU的架構和並行計算模式有了更進一步的認識，希望對以後的工作能有一些幫助。

不過對CUDA高性能計算優化技術的學習以後應該不會投入太多精力了吧，感覺太偏、太艱深，畢竟對硬體架構還是欠缺一些吧。