李 彪，劉 杰，2

(1.國防科技大學(xué)并行與分布處理國家重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410073,2.復(fù)雜系統(tǒng)軟件工程湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410073)

1 引言

1.1 粒子輸運(yùn)簡(jiǎn)介

輸運(yùn)理論是研究微觀粒子在介質(zhì)中遷移統(tǒng)計(jì)規(guī)律的數(shù)學(xué)理論。這里“微觀粒子”指的是中子、光子、分子、電子和離子等。輸運(yùn)理論研究大量粒子在空間或者某種介質(zhì)中運(yùn)行時(shí)，由于各粒子位置、動(dòng)量和其他特征量的變化而引起的各種有關(guān)物理量隨時(shí)空變化所表現(xiàn)出來的非平衡統(tǒng)計(jì)運(yùn)動(dòng)規(guī)律[1]。

19世紀(jì)中后葉的分子運(yùn)動(dòng)論導(dǎo)致粒子輸運(yùn)理論發(fā)展的開端。1872 年，Boltzmann 推導(dǎo)出了Maxwell分子速度分布函數(shù)隨時(shí)空變化的非線性積分-微分方程，該方程也常被稱為粒子輸運(yùn)方程或者Boltzman方程[1]。粒子輸運(yùn)方程描述的分布函數(shù)具有7個(gè)獨(dú)立變量，包括3個(gè)空間變量、2個(gè)方向角度變量、1個(gè)能量變量和1個(gè)時(shí)間變量。由于實(shí)際情況的復(fù)雜性，通常不能采用解析方式進(jìn)行直接求解，只能通過數(shù)值模擬求解。

目前常用的求解粒子輸運(yùn)問題的數(shù)值方法分成2類[2]：一類是“確定性方法”，這類方法通過對(duì)變量進(jìn)行離散，將粒子輸運(yùn)方程轉(zhuǎn)化為一組線性代數(shù)方程，然后通過求解此代數(shù)方程組來獲得精確解或近似解；另一類方法是“蒙特卡羅MC(Monte Carlo)方法”，也稱非確定性方法或概率論方法。它是基于統(tǒng)計(jì)理論的數(shù)值方法，對(duì)所要研究的問題構(gòu)造一個(gè)隨機(jī)概率模型來模擬大量粒子的運(yùn)動(dòng)，以獲得感興趣的物理量分布。

本文研究的是關(guān)于粒子輸運(yùn)中的非確定性方法。

1.2 蒙特卡羅方法在粒子輸運(yùn)上的應(yīng)用

蒙特卡羅MC方法亦稱隨機(jī)模擬法，是20世紀(jì)40年代中期電子計(jì)算機(jī)誕生后發(fā)展起來的一門計(jì)算科學(xué)[3]。它是通過計(jì)算機(jī)來對(duì)中子行為進(jìn)行隨機(jī)模擬的數(shù)值方法，是一種以概率統(tǒng)計(jì)理論為指導(dǎo)的非常重要的數(shù)值計(jì)算方法。

通常蒙特卡羅方法可以粗略地分成2類[3]：一類是所求解的問題本身具有內(nèi)在的隨機(jī)性，利用計(jì)算機(jī)直接模擬這種隨機(jī)的過程。例如，在中子輸運(yùn)理論中，要求建立單個(gè)中子在給定幾何系統(tǒng)中的真實(shí)運(yùn)動(dòng)歷史，通過對(duì)大量中子歷史的跟蹤，得到充分的隨機(jī)試驗(yàn)值(或稱抽樣值)，然后用統(tǒng)計(jì)的方法得出隨機(jī)變量某個(gè)數(shù)值特征的估計(jì)量，用此估計(jì)量作為問題的解。另一類是所求解問題可以轉(zhuǎn)化為某種隨機(jī)分布的特征數(shù)，比如隨機(jī)事件出現(xiàn)的概率，或者隨機(jī)變量的期望值。通過隨機(jī)抽樣的方法，以隨機(jī)事件出現(xiàn)的頻率估計(jì)其概率，或者以抽樣的數(shù)字特征估算隨機(jī)變量的數(shù)字特征，并將其作為問題的解。這種方法多用于求解復(fù)雜的多維積分問題。

1.3 基于MC方法的異構(gòu)協(xié)同計(jì)算

由于功耗和散熱問題日益突出，單塊芯片的性能已經(jīng)不能單純地從提高頻率來獲取提升，摩爾定律將會(huì)消亡[4]，于是工業(yè)界逐漸轉(zhuǎn)向多核處理器。但是，隨著技術(shù)的發(fā)展，多核處理器中每個(gè)內(nèi)核的設(shè)計(jì)越來越復(fù)雜，并且頻率較高，功耗較大，嚴(yán)重阻礙了多核處理器中核數(shù)的擴(kuò)展。為了進(jìn)一步提升性能，同時(shí)滿足功耗和散熱的要求，出現(xiàn)了由眾多相對(duì)簡(jiǎn)單的內(nèi)核構(gòu)建的眾核協(xié)處理器，如 NVIDIA公司推出的 GPU[5]、Intel公司推出的MIC協(xié)處理器[6]和本文所探討的國產(chǎn)加速器Matrix2000等。與多核處理器相比，這些眾核加速器的內(nèi)核相對(duì)簡(jiǎn)單，且頻率較低，但數(shù)量眾多，能夠以合理的功耗實(shí)現(xiàn)更高的性能。

近十年來，異構(gòu)加速器快速發(fā)展，TOP500排行中排名靠前的超算系統(tǒng)的體系結(jié)構(gòu)均采用了異構(gòu)協(xié)同模式，許多研究人員對(duì)已有的求解粒子輸運(yùn)問題的方法在各種加速平臺(tái)上做了許多相關(guān)的研究工作。

楊博[7]提出了一種針對(duì)CPU/GPU混合異構(gòu)系統(tǒng)的深穿透粒子輸運(yùn)MC模擬異構(gòu)協(xié)同算法，針對(duì)GPU的計(jì)算和訪存特點(diǎn)，給出了一種基于粒子數(shù)的任務(wù)劃分方法和高效并行數(shù)據(jù)結(jié)構(gòu)以及線程間的并行規(guī)約方法。相比運(yùn)行在Intel Xeon X5670上的MC程序，該算法獲得了3.4的加速比，并在TianHe-1A的64個(gè)結(jié)點(diǎn)進(jìn)行了測(cè)試，結(jié)果表明該算法具有良好的性能和可擴(kuò)展性。

崔顯濤[8]基于MC方法，提出了一種面向CPU/MIC混合異構(gòu)系統(tǒng)的粒子輸運(yùn)并行算法，針對(duì)MIC訪存特點(diǎn)提出了適于程序并行的高速數(shù)據(jù)結(jié)構(gòu)和基于靜態(tài)分配的任務(wù)劃分方式。相比運(yùn)行在CPU端的程序串行程序，改進(jìn)的MC程序獲得了8.6倍的加速比。

在多能群MC粒子輸運(yùn)領(lǐng)域，近些年來有許多研究和探索,Bergmann 等人[9]基于GPU開發(fā)了一個(gè)MC粒子輸運(yùn)框架WARP(Weaving All the Random Particles)。WARP的目的是使以前的基于事件的輸運(yùn)算法適應(yīng)新的基于GPU硬件而開發(fā)的連續(xù)能量Monte Carlo中子輸運(yùn)代碼。Hamilton等人[10]基于Shift程序包開發(fā)了一個(gè)在GPU上運(yùn)行的連續(xù)能量Monte Carlo中子輸運(yùn)求解器，并把程序移植到了超算系統(tǒng)Summit上，在耗盡燃料的SMR(Small Modular Reactor)核心配置上，進(jìn)行了1 024個(gè)結(jié)點(diǎn)的弱擴(kuò)展測(cè)試，結(jié)果顯示其并行效率接近93%，這說明GPU對(duì)性能提升顯著。

本文的工作與文獻(xiàn)[8,9]的接近，都是只考慮單群的MC程序，不過本文的異構(gòu)算法設(shè)計(jì)是基于不同于GPU和MIC體系結(jié)構(gòu)的國產(chǎn)加速器Matrix2000，且研制的程序可以移植到大規(guī)模結(jié)點(diǎn)上，另外，還針對(duì)MC程序的數(shù)據(jù)通信進(jìn)行了優(yōu)化，取得了顯著的性能提升。

2 粒子輸運(yùn)MC模擬

MC模擬是通過對(duì)大量中子歷史的跟蹤來獲取大量的隨機(jī)試驗(yàn)結(jié)果，并對(duì)這些結(jié)果進(jìn)行統(tǒng)計(jì)方法處理，最終以得到的統(tǒng)計(jì)量作為問題的解。 所謂一個(gè)中子的歷史，是指該中子從源出發(fā)，在介質(zhì)中隨機(jī)游走，經(jīng)過各種核反應(yīng)作用，直到中子歷史結(jié)束或稱中子“死亡”。所謂“死亡”是指中子被吸收、穿出系統(tǒng)、被熱化，或達(dá)到能量權(quán)下限或時(shí)間上限。其中時(shí)間、能量的載斷是無條件的，而權(quán)截?cái)嗍怯袟l件的，由俄羅斯輪盤賭決定[11]。

圖1是一個(gè)中子歷史的循環(huán)過程。初始時(shí)，從粒子源出發(fā)的一個(gè)粒子的位置、方向、能量等參數(shù)是確定的，隨后粒子的運(yùn)動(dòng)方向和與原子核的碰撞類型服從特定的概率分布，在經(jīng)歷了一系列碰撞過程之后，粒子歷史趨向于結(jié)束。有幾種不同類型的結(jié)束，比如達(dá)到時(shí)間與能量界限、逃脫出系統(tǒng)或達(dá)到權(quán)下限。在這一個(gè)過程中，下一次碰撞僅與當(dāng)前碰撞后粒子的位置、能量以及方向相關(guān)，完全獨(dú)立于先前的碰撞。 因此，這一過程是一個(gè)典型的馬爾可夫過程。所以，只要知道粒子與原子核碰撞的規(guī)律，那么粒子的軌跡就可以用蒙特卡羅方法正確地模擬出來，從而得到所關(guān)心的解[12]。

Figure 1 Simulation of particle transport process圖1 模擬粒子輸運(yùn)過程

3 通信模式優(yōu)化

3.1 串行阻塞通信模式

通過測(cè)試MC程序的強(qiáng)擴(kuò)展性發(fā)現(xiàn)，當(dāng)進(jìn)程數(shù)大于256時(shí)，效率急速下降，達(dá)到1 024進(jìn)程時(shí)基本上已經(jīng)沒有加速效果了，這對(duì)移植MC程序到大規(guī)模系統(tǒng)上產(chǎn)生了巨大障礙。經(jīng)過仔細(xì)分析程序的通信模式發(fā)現(xiàn)，由于計(jì)算結(jié)果數(shù)據(jù)收集的通信模式是主從模式，導(dǎo)致等待時(shí)間過長。

首先0號(hào)進(jìn)程初始化任務(wù)，接著把數(shù)據(jù)發(fā)送給其他所有從進(jìn)程，從進(jìn)程接收到數(shù)據(jù)后開始計(jì)算，并在從進(jìn)程的計(jì)算結(jié)束后設(shè)置同步柵欄，直到所有從進(jìn)程的計(jì)算都結(jié)束后，0號(hào)進(jìn)程才開始收集從進(jìn)程發(fā)送過來的結(jié)果數(shù)據(jù)。數(shù)據(jù)的收集過程中所有從進(jìn)程都向主進(jìn)程阻塞發(fā)送數(shù)據(jù)，主進(jìn)程進(jìn)行阻塞接收，先到的信息先被處理。對(duì)于主進(jìn)程來說這是一個(gè)串行處理過程，此通信方式會(huì)造成通信瓶頸，核數(shù)越多，對(duì)程序性能的影響越大。圖2展示了P個(gè)進(jìn)程的通信過程。

Figure 2 Serial data collection mode圖2 串行數(shù)據(jù)收集模式

3.2 優(yōu)化后的二叉樹通信模式

考慮到MC程序0號(hào)進(jìn)程對(duì)收集到的數(shù)據(jù)只進(jìn)行簡(jiǎn)單的統(tǒng)計(jì)處理，所以可以通過逐層進(jìn)行兩兩進(jìn)程通信，實(shí)現(xiàn)局部統(tǒng)計(jì)處理，最終把統(tǒng)計(jì)結(jié)果匯總到0號(hào)進(jìn)程。本文采用二叉樹通信模式來實(shí)現(xiàn)這種通信過程，假設(shè)總進(jìn)程數(shù)為2K+res，其中,K,res都是大于或等于0的整數(shù)，且滿足條件0≤res<2K，為了滿足二叉樹的通信模式，排在前面的2*res個(gè)進(jìn)程先就近奇偶進(jìn)程號(hào)兩兩通信，并把通信過程中已發(fā)送過信息的進(jìn)程剔除，不參與下次通信，從而下一階段通信的進(jìn)程數(shù)變?yōu)?K。接著對(duì)2K個(gè)進(jìn)程號(hào)重排，前2*res個(gè)進(jìn)程號(hào)除以2，其余進(jìn)程依次減去res，得到新的2K個(gè)進(jìn)程的標(biāo)識(shí)號(hào)，下一階段就是進(jìn)行二叉樹通信，算法1描述了具體的通信過程。為了方便描述，圖3展示了7個(gè)進(jìn)程的通信過程，括號(hào)里的數(shù)字表示重排后的進(jìn)程標(biāo)識(shí)號(hào)。

Figure 3 Binary tree data collection mode圖3 二叉樹數(shù)據(jù)收集模式

改進(jìn)后的數(shù)據(jù)收集模式，通信復(fù)雜度由2K+res減少為log (2K+res)，極大地減少了通信時(shí)間，也避免了大規(guī)模通信時(shí)從進(jìn)程同時(shí)向主進(jìn)程發(fā)送數(shù)據(jù)導(dǎo)致的程序阻塞。

4 基于CPU/Matrix2000的異構(gòu)協(xié)同算法

4.1 Matrix2000加速器

天河2A系統(tǒng)中，每個(gè)結(jié)點(diǎn)由2顆Intel Xeon 微處理器和2顆Matrix2000加速器組成，如圖4所示。每個(gè)Intel Xeon微處理器包含12核，工作頻率為2.2 GHz，采用英特爾Ivy Bridge微架構(gòu)，峰值性能0.211 2 TFLPOS。每個(gè)Matrix2000加速器包含128核，由4個(gè)超結(jié)點(diǎn)組成，每個(gè)超結(jié)點(diǎn)包含32個(gè)計(jì)算核，其中超結(jié)點(diǎn)支持64核超線程技術(shù)，峰值性能2.457 6 TFLOPS@1.2 GHz，有8個(gè)DDR4內(nèi)存通道，支持×16 PCIE 3.0 EP工作模式。

算法1二叉樹通信偽代碼

//假設(shè)總進(jìn)程數(shù)為2K+res，0≤res<2K

1mask=1;

2if(myid< 2*res)then/*前2*res個(gè)進(jìn)程先通信，使下一階段通信進(jìn)程個(gè)數(shù)為2K*/

3if(myid%2!=0)then//奇數(shù)號(hào)進(jìn)程

4dest=myid-1;

5new_id=-1/*進(jìn)程號(hào)置為-1，不參與下次通信*/

6myid進(jìn)程SENDmsgtodest;

7else//偶數(shù)號(hào)進(jìn)程

8src=myid+1;

9new_id=myid/2;//重排后新的進(jìn)程號(hào)

10myid進(jìn)程RECVmsgfromsrc

11endif

12else

13new_id=myid-res;/*其余的進(jìn)程重排后的進(jìn)程號(hào)*/

14endif

15if(new_id!=-1)then/*重排后的2K個(gè)進(jìn)程進(jìn)行通信*/

16while(mask<2K)do//二叉樹通信模式

17if(mask&new_id)then/*獲取重排號(hào)之后的偶數(shù)進(jìn)程號(hào)*/

18newsrc=new_id|mask;

19if(newsrc

20src=newsrc*2

21else

22src=newsrc+res

23endif

24 RECVmsgfromsrc

25else

26newdest=new_id& ～mask

27if(newdest

28dest=newdest* 2

29else

30dest=newdest+res

31endif

32 SENDmsgtodest

33 exit//發(fā)送過消息后，此進(jìn)程就退出通信

34endif

35mask= 2*mask

36endwhile

37endif

Figure 4 Structure of Tianhe-2A system node圖4 天河2A系統(tǒng)結(jié)點(diǎn)結(jié)構(gòu)

4.2 異構(gòu)編程模型

天河2A系統(tǒng)支持3種異構(gòu)通信模式OpenMP4.5、ACL和BCL，其中BCL是一種簡(jiǎn)單高效的對(duì)稱傳輸接口，CPU和協(xié)處理器之間數(shù)據(jù)通信利用PCIE總線進(jìn)行傳輸，底層通過SCIF來實(shí)現(xiàn)。BCL相比OpenMP4.5更底層，程序移植也更復(fù)雜，但是傳輸速率更快，移植后的程序靈活性更好。基于BCL接口的異構(gòu)程序需要編譯2套程序，CPU和加速器分別編譯一套程序，2套程序分別在CPU端和加速器端同時(shí)運(yùn)行，其中運(yùn)行時(shí)加速器端的程序需要通過ACL傳送到加速器端并啟動(dòng)。

圖5給出了基于CPU/Matrix2000的MC程序的異構(gòu)模式流程。首先，CPU端0號(hào)進(jìn)程啟動(dòng)MPI進(jìn)行進(jìn)程初始化，0號(hào)進(jìn)程負(fù)責(zé)讀取文件數(shù)據(jù)，并把總計(jì)算任務(wù)按照粒子數(shù)均等的方式分配給從進(jìn)程對(duì)應(yīng)的子任務(wù)，再將子任務(wù)的初始化數(shù)據(jù)傳輸給對(duì)應(yīng)的從進(jìn)程。每個(gè)從進(jìn)程控制一個(gè)Matrix2000超結(jié)點(diǎn)，通過ACL把加速器端的程序傳輸?shù)組atrix2000端并啟動(dòng)程序，利用BCL使CPU端與 Matrix2000超結(jié)點(diǎn)建立連接，CPU端把初始化數(shù)據(jù)傳輸?shù)組atrix2000端，加速器端完成CPU端數(shù)據(jù)接收后就利用多線程技術(shù)并行跟蹤每個(gè)粒子歷史，直到所有粒子歷史全部完成計(jì)算，加速器統(tǒng)計(jì)結(jié)果數(shù)據(jù)并把數(shù)據(jù)通過BCL傳輸給對(duì)應(yīng)的CPU進(jìn)程，再由MPI實(shí)現(xiàn)進(jìn)程間的傳輸。整個(gè)異構(gòu)控制的過程中，CPU主要負(fù)責(zé)傳輸數(shù)據(jù)，Matrix2000主要負(fù)責(zé)計(jì)算。

Figure 5 Flowchart of heterogeneous logic 圖5 異構(gòu)模式邏輯流程圖

4.3 OpenMP線程級(jí)并行

Matrix2000中一個(gè)超結(jié)點(diǎn)包含32核，支持64核超線程技術(shù)，為了充分發(fā)揮Matrix2000的性能，通過OpenMP指令實(shí)現(xiàn)線程級(jí)細(xì)粒度并行，算法2描述了具體算法偽代碼。加速器端接收CPU端發(fā)送過來的數(shù)據(jù)，接著做初始的工作，針對(duì)粒子輸運(yùn)模塊進(jìn)行多線程并行。每個(gè)線程跟蹤一個(gè)粒子，完成當(dāng)前粒子跟蹤后，進(jìn)入time_to_stop判斷是否所有粒子數(shù)都已完成，如果未完成，當(dāng)前線程臨界更新粒子數(shù)變量nps++，更新后的第nps個(gè)粒子分配給當(dāng)前線程進(jìn)行跟蹤模擬。這種調(diào)度方式類似迭代塊大小等于1的dynamic調(diào)度。

算法2Matrix2000加速器上的OpenMP線程級(jí)并行

1 假設(shè)線程數(shù)為T，加速器分配到M個(gè)粒子，粒子編號(hào)區(qū)間為(nps,nps+M-1);

2 #pragma omp parallel;/*開啟多線程，每個(gè)線程跟蹤一個(gè)粒子，完成后繼續(xù)跟蹤下一個(gè)粒子*/

3forall threadsdo

4 各個(gè)線程初始化自己的私有變量;

5omp_set_lock(&lock);/*設(shè)置鎖，當(dāng)前只有獲得鎖的線程才能更新共享數(shù)據(jù)*/

6 更新共享變量;

7while(!time_to_stop())do/*判斷當(dāng)前是否完成M個(gè)粒子數(shù)追蹤*/

8nps=nps+1;/*未完成，更新nps，追蹤第nps個(gè)粒子*/

9 更新共享變量;

10omp_unset_lock(&lock);/*釋放鎖，使得多線程可以調(diào)用hstory并行追蹤粒子*/

11 callhstory();/*線程各自調(diào)用hstory()函數(shù)，互不干涉，完成輸運(yùn)過程*/

12omp_set_lock(&lock);/*設(shè)置鎖，進(jìn)入time_to_stop()中訪問共享數(shù)據(jù)*/

13endwhile

14 更新共享變量;/*完成M個(gè)粒子數(shù)追蹤，進(jìn)行數(shù)據(jù)更新*/

15omp_unset_lock(&lock);

16endfor

相比Intel MIC協(xié)處理器Offload模式多線程技術(shù)以及OpenMP4.5提供的異構(gòu)模式，本文的異構(gòu)模式中的CPU與Matrix2000加速器并不是主從關(guān)系，而是對(duì)等的關(guān)系。Matrix2000加速器端啟動(dòng)的是單獨(dú)的一套程序，所有的數(shù)據(jù)傳輸和初始化工作在啟動(dòng)多線程之前已經(jīng)完成，所以Matrix2000上的OpenMP不需要通過編譯指導(dǎo)語句來傳輸數(shù)據(jù)和設(shè)置各種變量屬性，使得程序結(jié)構(gòu)設(shè)計(jì)更簡(jiǎn)單、更靈活。

5 實(shí)驗(yàn)及結(jié)果分析

5.1 測(cè)試環(huán)境及參數(shù)

測(cè)試平臺(tái)為天河2A系統(tǒng)，由于ACL和BCL指令只提供C/C++的接口，本文需要對(duì)MC程序插入C語言的通信控制接口，CPU端的程序采用Intel的編譯器和基于高速網(wǎng)的mpich3.2進(jìn)行編譯；加速器端采用自定義交叉編譯器進(jìn)行編譯，支持OpenMP指令。

單結(jié)點(diǎn)硬件配置為2顆Intel E5-2692 v2 @2.20 GHz(12 核/顆)+8個(gè)超結(jié)點(diǎn)Matrix2000，每個(gè)超結(jié)點(diǎn)通過超線程技術(shù)實(shí)現(xiàn)45核多線程并行。 MC程序的輸入文件參數(shù)除問題規(guī)模隨具體實(shí)驗(yàn)設(shè)置外，均保持默認(rèn)值。

5.2 通信優(yōu)化測(cè)試

考慮到數(shù)據(jù)的通信僅在CPU端進(jìn)程之間進(jìn)行，而并不涉及到加速器，所以本節(jié)僅在CPU上進(jìn)行對(duì)比測(cè)試，表1給出了通信優(yōu)化前后的測(cè)試結(jié)果。

Table 1 Test results before and after communication optimization

分析圖6可以發(fā)現(xiàn),原始程序測(cè)試結(jié)果中，512核是拐點(diǎn)，小于512核時(shí)程序還有加速效果，大于512核時(shí)測(cè)試時(shí)間反而急速上升，已經(jīng)沒有加速效果，表明此時(shí)計(jì)算能力不是瓶頸，數(shù)據(jù)通信占主要耗時(shí)。相反，優(yōu)化了通信模式的程序測(cè)試時(shí)，運(yùn)行時(shí)間隨進(jìn)程數(shù)遞增逐漸下降，沒有出現(xiàn)進(jìn)程數(shù)越多測(cè)試時(shí)間反而上升的現(xiàn)象，當(dāng)大于2 000核時(shí)耗時(shí)趨于穩(wěn)定，這是由于粒子數(shù)計(jì)算規(guī)模受限，增加進(jìn)程數(shù)的優(yōu)勢(shì)已經(jīng)不明顯導(dǎo)致的?？梢栽O(shè)想如果計(jì)算規(guī)模增加進(jìn)程數(shù)越多，優(yōu)化通信后的程序會(huì)比原始程序效果更好。

Figure 6 Comparison of two communication modes圖6 通信模式對(duì)比圖

5.3 異構(gòu)算法測(cè)試

本節(jié)對(duì)基于CPU/Matrix2000的異構(gòu)MC程序進(jìn)行弱擴(kuò)展測(cè)試，每個(gè)結(jié)點(diǎn)啟動(dòng)10個(gè)CPU進(jìn)程，其中8個(gè)進(jìn)程分別和結(jié)點(diǎn)內(nèi)的8個(gè)Matrix2000超結(jié)點(diǎn)進(jìn)行異構(gòu)計(jì)算，剩余2個(gè)進(jìn)程協(xié)同計(jì)算，每個(gè)超結(jié)點(diǎn)使用48核多線程，表2給出了弱擴(kuò)展測(cè)試結(jié)果，圖7給出了相應(yīng)測(cè)試結(jié)果的柱狀圖。

Table 2 Weakly scalable test

Figure 7 Bar and curve chart of weakly scalabletest result on Tianhe-2A圖7 天河2A弱擴(kuò)展測(cè)試結(jié)果柱狀圖

表2說明弱擴(kuò)展到45萬核時(shí)，相對(duì)5萬核并行效率保持在 22.54%。 分析圖7可以得到，雖然研制的程序可以運(yùn)行到45萬核，但是可以看到，隨著核數(shù)增大并行效率的下降趨勢(shì)很明顯。結(jié)合程序與加速器體系結(jié)構(gòu)分析，導(dǎo)致并行效率下降迅速的原因主要有2點(diǎn)：首先，MC程序的粒子輸運(yùn)部分的計(jì)算不是數(shù)據(jù)密集型，在中子的歷史過程中充滿了大量粒子條件狀態(tài)的判斷，程序的邏輯控制語句居多，這與Matrix2000加速器適合程序邏輯控制簡(jiǎn)單、數(shù)據(jù)計(jì)算密集的條件相沖突；其次，在程序的OpenMP多線程部分，雖然開啟了多線程技術(shù)，但是粒子歷史的模擬追蹤過程中各個(gè)線程都涉及更新共享數(shù)據(jù)，需要設(shè)置臨界區(qū)，在臨界區(qū)內(nèi)同時(shí)只能有一個(gè)線程更新共享數(shù)據(jù)，其他線程只能等待，這會(huì)影響多線程的并行效率。

6 結(jié)束語

本文在現(xiàn)有粒子輸運(yùn)蒙特卡羅模擬程序的基礎(chǔ)上，提出了一種面向CPU/Matrix2000異構(gòu)系統(tǒng)的粒子輸運(yùn)異構(gòu)協(xié)同算法，基于天河2A系統(tǒng)的異構(gòu)通信模式BCL和ACL，針對(duì)CPU和加速器各自設(shè)計(jì)了一套不同的代碼，在CPU端通過ACL實(shí)現(xiàn)加速器端代碼的啟動(dòng)，利用BCL進(jìn)行CPU與Matrix2000之間的通信，進(jìn)而提出了一種CPU與加速器Matrix2000之間的簡(jiǎn)單高效的對(duì)稱通信模式。優(yōu)化了原MC程序的串行數(shù)據(jù)收集通信模式，提出了新的二叉樹通信模式，極大地減少了通信時(shí)間，加速比可達(dá)17.7。通過優(yōu)化通信模式，以及基于MPI-ACL-OpenMP編程框架，本文實(shí)現(xiàn)的基于CPU-Matrix2000異構(gòu)協(xié)同計(jì)算的并行程序，可以弱擴(kuò)展到45萬核，相對(duì)5萬核并行效率保持在22.54%。作為未來工作，將探索如何繼續(xù)增強(qiáng)該異構(gòu)程序的擴(kuò)展性，以更好地發(fā)揮出HPC平臺(tái)的算力優(yōu)勢(shì)。