2 GPU 并行加速算法

　　GPU 具有高度并行的多流水線架構(gòu)，使其非常適宜于FDTD 加速運算。與CPU 運算的多次循環(huán)逐網(wǎng)格迭代更新方式不同，GPU 可以實現(xiàn)多網(wǎng)格的同時迭代更新，配合GPU 的線程集指令執(zhí)行機制，可以高效地利用流水線資源，隱藏流處理器與設(shè)備內(nèi)存間的場量讀寫延時，從而實現(xiàn)FDTD 運算加速。

　　2.1 GPU 函數(shù)的并行

　　2.1.1 數(shù)據(jù)并行優(yōu)化

　　OpenCL 支持按數(shù)據(jù)并行的編程模型和按任務(wù)并行的編程模型。數(shù)據(jù)并行是一種普遍意義上的并行方式，在程序中，常有一些循環(huán)操作，如果在循環(huán)的內(nèi)部各個計算過程之間不存在相關(guān)性，則可以轉(zhuǎn)換成并行循環(huán)，即數(shù)據(jù)并行。

　　由 FDTD 中電磁場的時間推進(jìn)計算公式可知，F(xiàn)DTD方法中每一網(wǎng)格點的電磁場值只與緊鄰的網(wǎng)格點上電磁場值及本網(wǎng)格點上一時間步的電磁場值有關(guān)，而與遠(yuǎn)端網(wǎng)格點的 場值無直接關(guān)系。也就是說，時間步長上的迭代過程是相互關(guān)聯(lián)、相互影響的，迭代需要上迭代的計算結(jié)果，故迭代之間不宜實現(xiàn)并行化。而在迭代內(nèi) 部，電場和磁場的計算僅需要前一時刻的計算結(jié)果，各個計算過程之間沒有影響，因此FDTD 并行化適合于用空間分區(qū)、時間分段的方式來分配運算任務(wù)，可以使用數(shù)據(jù)并行編程模型實現(xiàn)并行計算。

　　在 GPU 上實現(xiàn)FDTD 加速算法思路為：將每個數(shù)據(jù)元素映射到每個線程，就是將FDTD 三維網(wǎng)格劃分成為一個二維工作組，然后每個工作組中包含一維的工作項，通過工作組的id（groupID）和工作項的id（localID）來索引工作項 在緩沖對象中的位置，用一個工作項來計算一個場值，如圖3 所示。這樣，在一個工作組中的多個工作項可以相互協(xié)作、共享計算單元的內(nèi)存空間、同步、共享結(jié)果以節(jié)省計算，同時每個工作項當(dāng)中的工作組和周圍的兩個 Yee 元胞進(jìn)行數(shù)據(jù)交換，因此可以提供很快的數(shù)據(jù)訪問，大大節(jié)省設(shè)備內(nèi)存帶寬的消耗。

　　此外，由于 OpenCL 的緩沖對象只支持存儲一維元素集合，所以要把CPU 算法中相應(yīng)的三維數(shù)組和二維數(shù)組轉(zhuǎn)換為一維數(shù)組（包括電場參數(shù)和磁場參數(shù)）。

　　2.1.2 局部存儲器和內(nèi)建矢量使用

　　SAR_EH_CPX_Kernel 對一個一維緩沖對象進(jìn)行了更新操作。CPU 代碼中每次循環(huán)更新12 個變量，意味著每個計算單元需要讀/寫12 次全局內(nèi)存。我們設(shè)計的OpenCL并行方案訪存通過OpenCL 內(nèi)建矢量，每次訪存循環(huán)讀入3個cl_float4 變量到速度更快的局部存儲器中。將逐單元的訪存改為逐數(shù)據(jù)塊訪存，如圖4 所示，同時利用了速度與寄存器相當(dāng)?shù)木植看鎯ζ?。分析代碼發(fā)現(xiàn)，對規(guī)模為M*N*K 的三維數(shù)據(jù)塊減少了K*8 次全局內(nèi)存的讀/寫操作。

　　此外，對CPU 上耗時長、計算邏輯也比較簡單的對E 和H 進(jìn)行修正計算的函數(shù)，線程和數(shù)據(jù)之間是一一映射的對應(yīng)關(guān)系，為了提高訪存效率，在kernel 內(nèi)部也使用了內(nèi)建矢量數(shù)據(jù)類型對訪存進(jìn)行了優(yōu)化。

　　3.1.3 并行歸約

　　snefar 函數(shù)中對四個變量進(jìn)行求和操作，在OpenCL實現(xiàn)了兩個kernel 對該函數(shù)進(jìn)行了并行歸約的優(yōu)化。由于實際求和的個數(shù)不是2 的冪，所以設(shè)置工作項為2 的冪次方后，在Kernel 中首先初始化局部內(nèi)存為0,然后將數(shù)據(jù)從全局內(nèi)存順序讀入局部存儲器中，在局部存儲器進(jìn)行運算。實現(xiàn)歸約求和的主要代碼如下：

　　for（unsigned int s = get_local_size（0）/ 2; s > 0; s 》= 1）

　　{

　　if（tx < s）

　　sdata[tx] += sdata[tx + s];

　　barrier（CLK_LOCAL_MEM_FENCE）；

　　}

　　其中 tx 為get_local_id（0），為工作項在工作組中的一維下標(biāo)。歸約如圖5 所示。此并行歸約可以進(jìn)行循環(huán)展開，實現(xiàn)進(jìn)一步的優(yōu)化。

　　2.2 多GPU 系統(tǒng)的并行算法

　　2.2.1 區(qū)域分割與負(fù)載均衡

　　我們使用 MPI 進(jìn)行數(shù)據(jù)通信和調(diào)度，實現(xiàn)多GPU 核并行算法。步就是要對數(shù)據(jù)塊進(jìn)行區(qū)域分割，建立拓?fù)浣Y(jié)構(gòu)。

　　測試所用的 FDTD 數(shù)據(jù)塊為三維數(shù)據(jù)，故區(qū)域分割方式主要有三種：一維分割、二維分割和三維分割。觀察代碼發(fā)現(xiàn)，如果進(jìn)行一維分割，無論沿XYZ 哪個方向進(jìn)行，數(shù)據(jù)塊與數(shù)據(jù)塊之間的通信量極少，只涉及邊界交換，故為代碼移植簡單起見，我們實現(xiàn)的MPI 是基于一維沿X 方向進(jìn)行數(shù)據(jù)分割。

　　負(fù)載平衡是區(qū)域分割的一個重要原則，我們實現(xiàn)的代碼中首先獲得節(jié)點中所有具備計算能力的GPU 核，給每個核分配固定的id,根據(jù)id 分配相同大小數(shù)據(jù)塊，使得具有相同計算能力的處理器獲得相同的工作量。對于update_E_PML 和update_H_PML 這兩個計算邊界而言，由于我們沿X 方向進(jìn)行分割，在計算個節(jié)點和一個節(jié)點需要更新整個X 平面，因此這兩個函數(shù)的計算負(fù)載沒有實現(xiàn)均衡。

　　2.2.2 數(shù)據(jù)通信

　　多 GPU 核FDTD 算法另一個重要問題就是數(shù)據(jù)通信，也是MPI 存在的。在FDTD 計算區(qū)域分割成多個子區(qū)域后，各個子區(qū)域內(nèi)部的場分量按照原有的迭代公式進(jìn)行計算，而子區(qū)域界面上的場分量的迭代計算就需要相鄰子區(qū)域上相應(yīng)的場分量，此 時便需要相鄰子區(qū)域之間的數(shù)據(jù)通信。

　　簡單來說，對于我們實現(xiàn)的三維FDTD多GPU核并行，每個非邊界上的子區(qū)域都需要與左、右方向上的相鄰子區(qū)域進(jìn)行數(shù)據(jù)交換，如圖6 所示。在每一輪迭代中，這種邊界交換一共需要四次，兩次在更新E 和H 分量前，另外兩次在更新E 和H 的邊界前。

相關(guān)資料：

基于CPU-GPU 異構(gòu)機群的FDTD 并行算法加速研究（一）

基于CPU-GPU 異構(gòu)機群的FDTD 并行算法加速研究（三）