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

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

（1.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點實驗室，四川 成都610213；2.哈爾濱工程大學(xué)，黑龍江 哈爾濱150001）

0 引言

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

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

1 矩陣特征線方法概述

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

式中，ψ（s）為中子角通量；s為沿射線的當(dāng)?shù)刈鴺?biāo)；∑t和q（s）分別為宏觀總截面和中子源項。假設(shè)宏觀截面和源項在計算網(wǎng)格內(nèi)呈均勻分布時，可以解析求解上式，得到中子出射角通量和入射角通量之間的關(guān)系：

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

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

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

2 極小殘余算法的分析

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

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

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

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

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

3 基于CUDA技術(shù)的GMRES算法并行優(yōu)化

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

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

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

圖1 系數(shù)矩陣中各行非零元個數(shù)

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

4 數(shù)值驗證與討論

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

表3 計算環(huán)境參數(shù)配置表

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

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

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

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

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

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

表4 C5G7問題計算精度比較

5 結(jié)論

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