基于CUDA并行技術加速2維矩陣MOC方法的研究

鄭 勇 彭敏俊 安 萍 強勝龍 蘆 韡

（1.中國核動力研究設計院 核反應堆系統設計技術重點實驗室，四川 成都610213；2.哈爾濱工程大學，黑龍江 哈爾濱150001）

0 引言

隨著反應堆物理計算對幾何適應性和計算性能要求的提高，特征線方法在中子輸運計算領域應用越來越廣泛，而且由傳統特征線方法出發，可以推得矩陣形式的特征線方程（MMOC）[1]，該方法僅需在構造系數矩陣階段掃描特征線信息，而在隨后的內迭代和源迭代過程中，均是對矩陣方程組的數值求解。由于減少了繁復的特征線掃描次數，因而，MMOC具有減少輸運計算時間的潛力，但是否具備效率優勢還依賴于矩陣方程的求解過程。以往的研究結果表明[2]，該矩陣方程具有較強的稀疏性和不規則的非零元分布，這些特點使得直接求解和常規迭代算法束手無策，而隸屬于Krylov子空間方法的極小殘余算法（GMRES）具有較好的收斂性能和計算效率[3]。

重啟型的GMRES算法具有良好的收斂特性和數值穩定性，目前國際上已有大量的研究以期在算法層次上進一步減少GMRES算法的迭代次數。針對矩陣特征線方法中系數矩陣的特點，本文基于GPU并行技術和CUDA程序開發規范，從代碼層面上對重啟型GMRES算法進行并行化，以提升該算法的計算效率，從而加速矩陣MOC程序的求解。

1 矩陣特征線方法概述

特征線方法是將偏微分方程簡化為常微分方程，然后使用成熟的數學方法求解常微分方程的一種方法。在反應堆中子輸運計算領域，特征線方法在每個離散方向上生成一族射線覆蓋待求解區域，并跟蹤每條射線來研究中子通量的分布規律。沿著每條射線，中子輸運方程可以寫成如下的特征線形式：

式中，ψ（s）為中子角通量；s為沿射線的當地坐標；∑t和q（s）分別為宏觀總截面和中子源項。假設宏觀截面和源項在計算網格內呈均勻分布時，可以解析求解上式，得到中子出射角通量和入射角通量之間的關系：

由方程（2）不難看出，根據特征線在相鄰網格之間的角通量連續條件，可以從外邊界入射角通量依次求得任意網格的入射和出射角通量，從而獲得中子角通量的傳播方程：

結合中子角通量和通量密度之間的加權求和關系，可以建立起任意計算網格上中子通量密度同外邊界入射角通量和網格源項的方程；同時利用外邊界條件，還可以建立各條特征線外邊界處入射角通量之間的關系。從而，構成了待求解的矩陣方程：

式中，S1、S2、S′1、S′2、S″1、S″2均為掃描特征線得到的系數矩陣，需要在迭代計算開始前計算完畢；Φ和Ψ分別代表由中子通量密度和外邊界入射角通量組成的向量；E為單位矩陣；Q表示無自散射源項的其他源項之和組成的向量。

2 極小殘余算法的分析

由特征線方法得到的矩陣方程通常規模較大，且系數矩陣具有較強的稀疏性，難以應用直接求解算法，同時非零元的非對稱分布也限制了迭代算法的可選擇范圍。以往的研究表明，重啟型極小殘余算法（GMRES）具有比較穩定的收斂特性，若能選用恰當的預條件技術，還能夠進一步加速其收斂過程。下面給出帶左預條件的重啟型GMRES算法的標準過程（其中，M為預條件矩陣）：

在矩陣特征線方法中，當子空間維度m不大于50時，算法即可收斂。因此，該算法中最耗時的部分是用于構造正交基的Arnold循環，而最小二乘問題的求解時間幾乎可以忽略不計，表1分析了Arnold過程中涉及向量、矩陣的各種基本操作的計算量。

表1 Arnold過程基本操作計算量

假設矩陣方程的階數為n，而矩陣A及M-1的非零元個數分別為nza、nzm，那么向量-向量與向量-標量運算的總浮點計算量為（2m2+7m+4）n，矩陣-向量乘操作需要的總計算量為4（m+1）（nza+nzm）。不妨定義稀疏度sp=（nza+nzm）n，則可以得到兩種運算量的比值為：

如前所述，當子空間維度m∈（0，50）時，f∈（1，26），因此當稀疏度sp小于26時，向量運算的浮點計算量可能會更大，反之則必然是SpMV的浮點計算量更大。通常情況下，特征線所穿過的計算網格數量以及覆蓋網格的特征線條數都大于26，因此矩陣的稀疏度sp通常會遠大于26。

3 基于CUDA技術的GMRES算法并行優化

SpMV運算可以分成如算法1所示的兩步，首先將矩陣非零元素與向量對應元素相乘，并把乘積保存在臨時矩陣temp中；然后對每個行向量進行累加求和，得到目標向量對應的元素。

表2 CSR格式的稀疏矩陣-向量乘算法

在第1步中，由于對向量x的訪問模式是隨機的，因此，可以將它綁定到訪存效率較高且適合隨機訪問的紋理內存中，其訪存延遲可以通過激活海量線程而有效隱藏。該算法的并行優化關鍵在于第2步，由于稀疏矩陣每一行的非零元個數都是不相等的（如圖1所示），如何同時高效地求得臨時矩陣temp每個行向量的累加和是優化的難點。

圖1 系數矩陣中各行非零元個數

針對稀疏矩陣不同部分稀疏程度的不同，本文結合文獻[4]提出合并訪存的方法，采取分部并行優化策略：第2部分由于每行平均非零元個數較少，因此，可以直接采用合并訪存和共享內存技術實現CUDA并行；對第1部分提出了兩種CUDA并行計算優化方案，在合并訪存的基礎上將每行的求和操作分別布局到線程和線程塊上。以上三種并行方案，在后文中分別簡寫為：gpu_cv/sm、gpu_thread及gpu_block。

4 數值驗證與討論

為了評估并行程序的加速性能，本節對2維C5G7問題進行了計算。計算環境參數配置如表3所示，雖然兩塊GPU均支持CUDA編程，但本文僅采用Tesla專業計算顯卡作分析比較。

表3 計算環境參數配置表

國際上，C5G7基準題被廣泛用于輸運程序計算精度的驗證，本文針對C5G7基準題的1/4堆芯進行計算，以比較上述幾種方法在處理較大規模問題時的加速效果。計算區域共包含4個17×17的燃料組件和外圍的5個反射層組件，采用文獻[5]給出的7群截面參數。計算條件為16個方位角，2個優化的Leonard極角[6]，特征線間距為0.05 cm，采用粗網有限差分求解器加速；求解器參數設置為子空間維度為10，允許的最大迭代次數為10。

計算結果表明，無論采取串行還是并行計算，2維C5G7基準題的計算結果都是完全相同的，且并行方案對計算精度幾乎沒有影響，部分結果如表3所示，可以發現與國際上其他輸運計算結果相比[5]，本文給出的計算精度完全在可以接受的范圍內。

表5給出了采取不同并行方案的耗時，同時CMFD加速方法和Krylov子空間維度對計算效率的影響也在表中給出。從表中可以看出，無論是cpu串行計算還是gpu并行計算，CMFD都能帶來超過6倍的加速比。當子空間維度越小，GMRES（m）求解器需要的重啟動次數越多，自然也就需要更多的SpMV操作；而子空間維度增大后，對向量的操作也會呈平方關系增加，這種平方增加關系不允許選擇太大的子空間維度。與gpu并行程序比較，在cpu串行計算中，子空間維度對計算效率的影響更為顯著，這是由于gpu并行程序對向量的操作耗時幾乎可以忽略不計。

表5 2維C5G7問題運行時間比較/單位：秒

圖2給出了以cpu串行計算為參考，不同計算條件下獲得的加速比，圖中的Case編號對應于表4中的第2列。由圖中可以看到，不同的并行策略對加速比有顯著影響，優化后的并行方案獲得的加速比成倍增加，且gpu_block具有最高的計算效率，能夠達到6倍以上的加速比。在每種并行策略中，CMFD加速方案對加速比均無明顯影響。隨子空間維度增加，加速比有輕微的變化，這種變化主要是因為cpu串行計算效率對子空間維度更為敏感。

圖2 2維C5G7問題不同計算條件下的加速比

表4 C5G7問題計算精度比較

5 結論

本文首先簡要介紹了矩陣特征線方法的基本模型，引出了求解矩陣方程的極小殘余算法，繼而對該算法中涉及矩陣和向量的基本運算量進行了分析。理論分析結果表明，稀疏矩陣-向量乘操作是算法中的熱點，開發并行求解器時時需要重點優化。基于CUDA并行技術，采用訪存合并和共享內存并行方案，編制了GMRES算法的GPU并行版本，并在此基礎上提出了兩種并行優化策略。為了驗證并行策略的有效性，評估并行求解器的加速效果，對2維C5G7基準題進行了對比計算。計算結果表明，并行方案的不同對計算精度沒有顯著影響，而優化后的并行程序較原始并行程序有更高的加速比。