鄭 勇 彭敏俊 安 萍 強(qiáng)勝龍 蘆 韡

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

0 引言

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

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

1 矩陣特征線方法概述

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

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

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

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

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

2 極小殘余算法的分析

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

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

表1 Arnold過程基本操作計(jì)算量

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表5 2維C5G7問題運(yùn)行時(shí)間比較/單位：秒

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

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

表4 C5G7問題計(jì)算精度比較

5 結(jié)論

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