郝 琛，韓立會(huì)，程有瑩

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

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

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

1.1 并行GMRES算法特性分析

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

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

13)end for

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

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

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

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

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

2.1 集中點(diǎn)積的計(jì)算

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

fori>0,

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

7) forj=0,…,i-2

9)end for

A2vi-1=Azi

(1)

(2)

Hessenberg矩陣元素的計(jì)算可表示為

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

(3)

(4)

(5)

2.2 非阻塞通信

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

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

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

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

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

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

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

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

3.1 計(jì)算平臺(tái)及測(cè)試算例

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

3.2 數(shù)值分析

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

(a)C5G7 benchmark problem

(b)VENUS2 benchmark problem

p1-GMRES算法和MGS-GMRES算法求解CMFD的用時(shí)對(duì)比，如表1所列。

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

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

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

(6)

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

(a)C5G7 benchmark problem

(b)VENUS2 benchmark problem

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

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

(a)N=224

(b)N=392

(c)N=512

(d)N=768

(e)N=864

(f)N=918

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

(7)

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

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

4 結(jié)論

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