求解CMFD的改進(jìn)流水線并行GMRES方法

郝 琛，韓立會，程有瑩

(1.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院，哈爾濱150001；2.國家核應(yīng)急響應(yīng)技術(shù)支持中心，北京100080)

基于特征線方法的全堆芯3維精細(xì)化物理計算方法是國際研究熱點之一。目前，全堆芯精細(xì)化中子輸運計算中多應(yīng)用2維或1維聚合算法，并應(yīng)用粗網(wǎng)有限差分方法(CMFD)加速3維輸運計算[1]。在求解CMFD大型非對稱線性稀疏系統(tǒng)時，首爾大學(xué)開發(fā)的全堆芯輸運計算程序nTRACER應(yīng)用BICGSTAB算法求解CMFD線性系統(tǒng)，可降低并行矩陣向量乘積計算的耗時[2]。先進(jìn)輕水反應(yīng)堆模擬仿真聯(lián)盟開發(fā)的中子輸運程序MPACT中采用無Jacobian矩陣的Newton-Krylov方法，將非線性牛頓法和有限差分方法相結(jié)合代替Jacobian矩陣向量乘積[3]。我國自主開發(fā)的堆芯高保真中子輸運計算程序HNET中應(yīng)用GMRES方法求解CMFD線性系統(tǒng)，但在多核并行環(huán)境下，由于傳統(tǒng)GMRES算法中點積的計算比矩陣向量乘積計算更加耗時，所以CMFD的加速效果受到限制。例如，用標(biāo)準(zhǔn)的并行GMRES方法求解C5G7-3D多群CMFD線性系統(tǒng)，當(dāng)計算核數(shù)為768時，全局通信用時高達(dá)GMRES求解總時間的95%[4]。因此，減少點積計算引起的全局通信用時是GMRES方法并行算法設(shè)計與實現(xiàn)的研究重點。s步GMRES方法通過減少內(nèi)積次數(shù)可實現(xiàn)算法優(yōu)化[5-6]，但該方法的數(shù)值穩(wěn)定性較差，且計算和全局通信連續(xù)執(zhí)行，處理器在通信階段處于空閑狀態(tài)，會造成處理器計算資源浪費。本文采用流水線思想[7]，重構(gòu)GMRES算法，實現(xiàn)每次迭代僅一次全局通信，并以最大限度的計算任務(wù)隱藏全局通信造成的延遲，充分利用處理器的資源，開發(fā)了流水線式并行GMRES求解算法及線性求解器。數(shù)值結(jié)果表明，流水線式并行GMRES求解器的并行計算效率大大提高，可有效加速CMFD線性系統(tǒng)的求解。

1 求解CMFD的標(biāo)準(zhǔn)并行GMRES算法

1.1 并行GMRES算法特性分析

GMRES算法通常用于求解非對稱線性方程組，可在一定迭代步數(shù)后將殘差降到最小。并行GMRES算法涉及的主要計算包括矩陣向量乘積(SpMV)、向量點積(dot product)、向量校正(AXPY)及預(yù)處理。

CMFD方法的并行實現(xiàn)基于區(qū)域分解，矩陣向量乘積(通常包括預(yù)處理部分)僅涉及相鄰子區(qū)域之間的數(shù)據(jù)交換，在相鄰的處理器之間進(jìn)行通信，因此，矩陣向量乘積不會造成使性能嚴(yán)重降低的通訊問題。向量校正不涉及處理器之間的通信，而向量點積的并行計算需要來自每個結(jié)點、每個計算核的信息，是一個全局通信的過程，需要所有處理器保持同步。隨著并行核數(shù)的增加，全局通信的用時越來越大，對算法的并行可擴(kuò)展性產(chǎn)生嚴(yán)重的負(fù)面影響，點積計算引起的全局通信成為限制并行GMRES算法的主要因素。

1.2 Gram-Schmidt正交化方法

GMRES算法利用Gram-Schmidt方法構(gòu)造標(biāo)準(zhǔn)正交基，為提高數(shù)值穩(wěn)定性，常用修正的Gram-Schmidt方法(MGS)取代經(jīng)典的Gram-Schmidt方法(CGS)?；贛GS的GMRES算法(MGS-GMRES)為

1)r0=b-Ax0

3)fori=0,…，m-1,do

4)w=Avi

5)forj=0,…，i,do

6)hj,i=(w,vj)

7)w=w-hj,ivj

8)end for

12)對Hi+1,i應(yīng)用Givens變換

13)end for

14)xm=x0+Vmym，其中，ym是使‖Hm+1,mym-‖r0‖2e1‖2最小的解。

該過程中，進(jìn)行點積運算的w向量在內(nèi)迭代過程中不斷更新，且只有得到更新后的w向量才能進(jìn)行下一次點積運算，因此，對于第m次迭代正交化過程，需要逐次進(jìn)行m+1次點積運算。在并行環(huán)境下，這意味著要進(jìn)行m+1次全局通信，造成巨大的通信用時，因此，MGS-GMRES算法不適用于多核并行計算。MGS-GMRES算法在4個處理機(jī)上的迭代過程，如圖1所示。其中，綠色表示本地計算，在典型的并行實現(xiàn)中不需要處理機(jī)之間進(jìn)行通信；藍(lán)色表示矩陣向量乘積計算，只需要相鄰處理機(jī)之間進(jìn)行通信；紅色實線表示全局通信。

圖1 MGS-GMRES算法在4個處理器上的迭代過程Fig.1 Schematic representation of a single iteration of the MGS-GMRES algorithm on four processors

CGS方法在正交化過程中的點積運算彼此獨立，不存在數(shù)據(jù)相關(guān)性。通過合理的算法設(shè)計，首先進(jìn)行局部點積運算，然后將每次點積運算需要傳遞的數(shù)據(jù)合并，提高通訊粒度，進(jìn)而通過一次全局通信得到正交化部分所有點積的計算結(jié)果，極大地降低了全局通信的次數(shù)。針對CGS方法的不穩(wěn)定性，可以采取迭代CGS方法進(jìn)行改善[8]。本研究中采用重啟動GMRES算法提高計算精度和穩(wěn)定性，保證基向量的正交性，測試表明，效果良好。

2 并行GMRES算法的優(yōu)化

2.1 集中點積的計算

基于CGS的GMRES算法(CGS-GMRES)進(jìn)行程序設(shè)計，實現(xiàn)正交化部分點積合并計算。CGS-GMRES算法正交化與標(biāo)準(zhǔn)化過程點積計算之間存在數(shù)據(jù)依賴性，本質(zhì)上是串行過程，因而正交化部分與標(biāo)準(zhǔn)化部分的點積計算不能進(jìn)行合并。分析每次迭代過程發(fā)現(xiàn)，通過算法重構(gòu)可以將一次迭代正交化過程包含的點積與上一次迭代標(biāo)準(zhǔn)化過程中的點積集中進(jìn)行計算，從而實現(xiàn)1次全局通信得到所有點積計算結(jié)果。該算法的主要過程為

fori>0,

1)forj=0,…,i-2，

7) forj=0,…,i-2

9)end for

A2vi-1=Azi

(1)

(2)

Hessenberg矩陣元素的計算可表示為

hj,i-1=(Avi-1,vj)=(zi,vj)=

(3)

(4)

(5)

2.2 非阻塞通信

大規(guī)模并行計算環(huán)境下，全局通信耗時巨大，在阻塞通信結(jié)束之前，各個計算核心必須等待同步通信結(jié)束才能進(jìn)行下一項工作，因此，浪費了計算資源。流水線技術(shù)是高效率利用硬件資源的經(jīng)典方法，當(dāng)不同階段需要不同資源時，多任務(wù)重疊可以同時利用不同資源[9]。本文采用流水線思想，應(yīng)用非阻塞通信，將通信任務(wù)交給特定的通信硬件去執(zhí)行，在該通信硬件進(jìn)行通信操作期間，充分應(yīng)用各個計算核心的計算資源，處理與該通信無關(guān)的計算任務(wù)，提高程序執(zhí)行效率。

非阻塞通信，一般需要調(diào)用2個MPI函數(shù)。首先要啟動非阻塞通信，但啟動并不代表該通信過程的完成。為了確保通信的完成，還必須調(diào)用與該通信相關(guān)聯(lián)的通信完成接口，如MPI_Wait或MPI_Test，以等待或查詢非阻塞通信是否完成，完成了檢測調(diào)用才真正將非阻塞通信完成。

2.3 并行GMRES算法優(yōu)化

采用非阻塞通信方式，在確保程序正確性的基礎(chǔ)上，調(diào)整程序各部分執(zhí)行的先后順序，確保盡早啟動和盡晚結(jié)束非阻塞通信[10]。用最大限度的計算任務(wù)屏蔽通信延遲，充分利用處理機(jī)通信期間的空閑等待時間。

本文基于已有的GMRES求解器，采用MPI并行編程模型，對重構(gòu)算法部分程序執(zhí)行順序進(jìn)行調(diào)整，將第m+1次迭代的矩陣向量積計算調(diào)整至第m次迭代點積計算后進(jìn)行。先發(fā)出點積計算非阻塞通信組歸約命令，即調(diào)用MPI_Iallreduce函數(shù)，后進(jìn)行下一次迭代矩陣向量乘積的計算，完成計算后再調(diào)用MPI_Wait函數(shù)，等待非阻塞通信的完成。這樣上一次迭代還未結(jié)束，下一次迭代已經(jīng)開始，2次迭代在時間上重疊了一部分。改進(jìn)的流水線式并行GMRES算法簡稱為p1-GMRES算法。其中，p代表流水線，1表示在執(zhí)行非阻塞組歸約的同時進(jìn)行一個矩陣向量乘積的計算。

應(yīng)用非阻塞組歸約操作，將2次迭代融合。圖2為流水線式并行GMRES算法迭代示意圖。雖然每次迭代構(gòu)造標(biāo)準(zhǔn)正交基的用時未發(fā)生變化，但以流水線的方式工作，減少了總的求解時間。

圖2 流水線式并行GMRES算法迭代示意圖Fig.2 Schematic representation of pipelineparallel GMRES algorithm iteration

3 數(shù)值結(jié)果

3.1 計算平臺及測試算例

數(shù)值試驗均在“天河一號”超算平臺上進(jìn)行。選取C5G7[11]及VENUS2[12]2個基準(zhǔn)例題為研究對象，分別構(gòu)建相應(yīng)的3維多群CMFD線性系統(tǒng)。

3.2 數(shù)值分析

圖3給出了GMRES求解計算時間t隨核數(shù)N的變化關(guān)系。由圖3可見，p1-GMRES算法具有較高的計算性能，GMRES求解計算時間明顯減小。

(a)C5G7 benchmark problem

(b)VENUS2 benchmark problem

p1-GMRES算法和MGS-GMRES算法求解CMFD的用時對比，如表1所列。

表1 p1-GMRES算法和MGS-GMRES算法求解CMFD的用時對比Tab.1Comparison between p1-GMRES algorithmand MGS-GMRES algorithm for CMFD

由表1可見，在N<392時，即MGS-GMRES加速效果達(dá)到峰值前，p1-GMRES算法求解CMFD的用時比MGS-GMRES算法的用時減小約35%，可以有效加速CMFD線性系統(tǒng)的求解。

定義加速比SN為同一算法在單個計算核心上的運行時間與在N個計算核心構(gòu)成的并行系統(tǒng)上的運行時間之比[13]，即

(6)

圖4給出了并行計算的加速比。由圖4可見，p1-GMRES算法具有較高的加速比。對于C5G7基準(zhǔn)題的CMFD線性系統(tǒng)求解，MGS-GMRES算法在核數(shù)為128附近加速比達(dá)到峰值，此后隨著核數(shù)增加，加速比呈現(xiàn)下降趨勢，運行時間逐漸增加；而p1-GMRES算法在核數(shù)大于128后加速比仍在增大，直到核數(shù)為416附近加速比達(dá)到峰值。對于規(guī)模更大的VENUS2基準(zhǔn)題的求解，MGS-GMRES算法在核數(shù)為512附近加速比達(dá)到峰值；而p1-GMRES算法在核數(shù)為864附近加速比達(dá)到峰值，極大地提高了加速比。

(a)C5G7 benchmark problem

(b)VENUS2 benchmark problem

以稍微增加計算時間為代價，降低全局通信次數(shù)，并將全局通信與矩陣向量乘積的計算重疊，從而將全局通信用時降為最低。以求解C5G7基準(zhǔn)題的CMFD線性系統(tǒng)為例，圖5給出了并行GMRES算法各部分的執(zhí)行時間。

由圖5可見，p1-GMRES算法重疊計算后的全局通信用時比MGS-GMRES算法大大降低；而且，隨著核數(shù)的增加，計算時間逐漸降低，全局通信用時逐漸增加，當(dāng)加速計算的時間小于全局通信增加的時間時，運行時間開始增加。結(jié)合圖3可知，隨著核數(shù)增加，計算用時呈反比例關(guān)系迅速下降。與MGS-GMRES算法相比，在核數(shù)為768時p1-GMRES算法增加的計算量很小，且遠(yuǎn)遠(yuǎn)小于降低通信時所獲得的時間收益。

(a)N=224

(b)N=392

(c)N=512

(d)N=768

(e)N=864

(f)N=918

并行效率反映一個計算核心的計算能力被有效利用的比率。并行效率的表達(dá)式為

(7)

表2給出了MGS-GMRES算法和p1-GMRES算法并行效率的對比。由表2可見，在對 C5G7基準(zhǔn)題計算中，核數(shù)為128時，p1-GMRES算法的并行效率高于50%，而MGS-GMRES算法的并行效率小于30%；對規(guī)模更大的VENUS2基準(zhǔn)題計算中，核數(shù)為224時，p1-GMRES算法的并行效率高于60%，而MGS-GMRES算法的并行效率小于35%。表2數(shù)據(jù)表明，p1-GMRES算法的并行效率顯著提高，改進(jìn)的流水線式并行GMRES求解器可以高效地利用處理器的計算能力，具有良好的可擴(kuò)展性。

表2 MGS-GMRES算法和p1-GMRES算法并行效率的對比Tab.2Comparison of parallel efficiency betweenMGS-GMRES and p1-GMRES

4 結(jié)論

本文對并行GMRES方法中通信用時最長的向量點積并行計算進(jìn)行了深入研究，采用MPI非阻塞通信方式，優(yōu)化和重構(gòu)了GMRES算法，提出了一種改進(jìn)的流水線式并行GMRES算法，提高了通信粒度，開發(fā)了線性求解器，并將其成功耦合于精細(xì)化中子物理計算程序HNET中。以C5G7和VENUS2的3維CMFD線性系統(tǒng)求解為例，對MGS-GMRES算法與p1-GMRES算法的運行時間、加速比及并行效率進(jìn)行了對比分析。數(shù)值結(jié)果表明，改進(jìn)的流水線式并行GMRES算法極大地降低了向量內(nèi)積并行計算對并行性能的影響，具有更高的加速比和并行效率，計算時間更短，可擴(kuò)展性良好，實現(xiàn)了加速求解CMFD線性系統(tǒng)的目的。