張 淼，周 宇，陳建海+，何欽銘，徐 順，宮 明

1.浙江大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，杭州310012

2.中國(guó)科學(xué)院 計(jì)算機(jī)網(wǎng)絡(luò)信息中心，北京100190

3.中國(guó)科學(xué)院 高能物理研究所，北京100049

+通訊作者E-mail:chenjh919@zju.edu.cn

1 引言

格點(diǎn)量子色動(dòng)力學(xué)（lattice quantum chromodynamics，LQCD）是粒子物理最前沿的研究課題之一。自然界存在四種基本相互作用——強(qiáng)相互作用、弱相互作用、電磁相互作用及引力相互作用。物理學(xué)家將強(qiáng)相互作用、弱相互作用及電磁相互作用統(tǒng)一起來(lái)稱之為“標(biāo)準(zhǔn)模型”。量子色動(dòng)力學(xué)（quantum chromodynamics，QCD）規(guī)范理論是標(biāo)準(zhǔn)模型中描述夸克、膠子組成的粒子強(qiáng)相互作用的基本理論，其重要特征是漸進(jìn)自由，這使得高能強(qiáng)相互作用過(guò)程可利用微擾論進(jìn)行處理，但此方法針對(duì)低能強(qiáng)相互作用過(guò)程就無(wú)能為力了。1974 年，Wilson 創(chuàng)建了格點(diǎn)規(guī)范理論[1]，其將連續(xù)時(shí)空用離散晶格來(lái)替代，夸克放置于格點(diǎn)上，它們之間用規(guī)范場(chǎng)鏈建立聯(lián)系。此理論不借助任何QCD以外的參數(shù)或假設(shè)從第一原理出發(fā)用非微擾方法來(lái)解決上述問(wèn)題。

在用計(jì)算機(jī)對(duì)LQCD 理論進(jìn)行數(shù)值模擬時(shí)需使用蒙特卡羅算法，而這一過(guò)程是高度計(jì)算密集型的，故現(xiàn)代LQCD 數(shù)值模擬過(guò)程均使用超級(jí)計(jì)算機(jī)來(lái)進(jìn)行。目前，已經(jīng)有很多工作在大規(guī)模超級(jí)計(jì)算機(jī)上實(shí)現(xiàn)了LQCD數(shù)值模擬。例如：2006年，Vranas等人[2]基于IBM 的BlueGene/L 超級(jí)計(jì)算機(jī)，在32 000 個(gè)CPU核心下，模擬QCD實(shí)現(xiàn)了12.2 TFLOPS的性能，并以此獲得了2006年的戈登貝爾獎(jiǎng)。2010年，Babich等人[3]在32 個(gè)由InfiniBand 互連的GPU 上，采用MPI（message passing interface）并行化QUDA 庫(kù)，實(shí)現(xiàn)了4 TFLOPS的性能。2011年，Ammendola等人[4]用多核CPU（central process unit）、GPU（graphics process unit）及FPGA（field-programmable gate array）構(gòu)建了一個(gè)由42U racks 組成的3D torus 網(wǎng)絡(luò)原型系統(tǒng)——QUonG 用于LQCD 計(jì)算；同年，Babich 等人[5]通過(guò)使用多維度的并行化策略和域分解的預(yù)調(diào)節(jié)器將由post-Monte Carlo 算法實(shí)現(xiàn)的LQCD 擴(kuò)展到256 個(gè)GPU上。2013年，Joó等人[6]將Dslash內(nèi)核移植到Intel Xeon Phi 加速卡上，單節(jié)點(diǎn)取得了280 GFlops 的性能，整個(gè)求解器獲得了大約215 GFlops 的性能，當(dāng)將節(jié)點(diǎn)擴(kuò)展到64 個(gè)Intel KNC 時(shí)，整個(gè)求解器獲得了3.6 TFlops 的性能。2017 年，Bonati 等 人[7]使用OpenACC 實(shí)現(xiàn)了可移植的LQCD 蒙特卡羅算法代碼。2018 年，Kanamori 等人[8]將Wilson 內(nèi)核在Oakforest-PACS超級(jí)計(jì)算機(jī)上實(shí)現(xiàn)了并行化。

從“天河一號(hào)”初登TOP500榜首，到“天河二號(hào)”連續(xù)3 年登頂，世界首臺(tái)10 億億次超級(jí)計(jì)算機(jī)“神威·太湖之光”連續(xù)4次奪冠，我國(guó)超算硬件架構(gòu)已達(dá)世界先進(jìn)水平，軟件水平也有了一定的提升，如在神威平臺(tái)上實(shí)現(xiàn)的千萬(wàn)核可擴(kuò)展全球大氣動(dòng)力學(xué)全隱式模擬[9]和高可擴(kuò)展非線性地震模擬[10]分獲了2016年與2017 年的戈登貝爾獎(jiǎng)，這是歷史性的突破。但目前尚無(wú)工作針對(duì)LQCD 在這些國(guó)內(nèi)超算平臺(tái)上進(jìn)行移植優(yōu)化。2018年最新一期TOP500榜單，“神威”已被美國(guó)最新的“Summit”超級(jí)計(jì)算機(jī)超越。未來(lái)，我國(guó)不僅要重奪超算硬件架構(gòu)水平的制高點(diǎn)，還需加快自身軟件生態(tài)的建設(shè)。故為提升我國(guó)自主超算平臺(tái)的軟件水平，針對(duì)神威平臺(tái)的特有架構(gòu)，本文工作主要是將LQCD中的熱點(diǎn)代碼Dslash進(jìn)行移植，通過(guò)代碼向量化、循環(huán)展開(kāi)、并行文件讀寫(xiě)、主-從核DMA（direct memory access）通訊、主-主核MPI 通訊等優(yōu)化手段，實(shí)現(xiàn)了Dslash代碼在申威26010芯片上的異構(gòu)眾核并行,并實(shí)現(xiàn)了一定的優(yōu)化效果，單核組從核陣列優(yōu)化取得了165倍的加速比，使用MPI連并16 個(gè)核組優(yōu)化取得了25 倍的加速比，同時(shí)發(fā)現(xiàn)使用MPI 引發(fā)了性能瓶頸，通過(guò)實(shí)驗(yàn)定位了瓶頸所在，為后續(xù)優(yōu)化工作打下基礎(chǔ)。

本文組織結(jié)構(gòu)如下：第2章將介紹神威平臺(tái)的系統(tǒng)架構(gòu)與通信模式；第3章詳細(xì)闡述了LQCD在神威平臺(tái)上的優(yōu)化過(guò)程；第4章對(duì)優(yōu)化效果進(jìn)行了數(shù)據(jù)測(cè)試并加以分析；第5章對(duì)全文內(nèi)容進(jìn)行了總結(jié)并對(duì)未來(lái)工作給予展望。

2 神威架構(gòu)及其通信模式

2.1 系統(tǒng)架構(gòu)

神威·太湖之光超級(jí)計(jì)算機(jī)采用的是基于申威26010 處理器的異構(gòu)眾核架構(gòu)，其處理器架構(gòu)如圖1所示。每個(gè)處理器包含有4個(gè)核組（core group，CG），而每個(gè)核組又包括一個(gè)主核（management processing element，MPE）、64 個(gè) 從 核（computing processing element，CPE）組成的8×8計(jì)算處理元件陣列，以及一個(gè)協(xié)議處理單元（protocol processing unit，PPU）和一個(gè)DDR3 存儲(chǔ)控制器（memory controller，MC）。4 個(gè)核組通過(guò)片上網(wǎng)絡(luò)（network on chip，NOC）連接，且每個(gè)核組有自己的內(nèi)存空間，主核和從核均通過(guò)存儲(chǔ)控制器來(lái)訪問(wèn)。主核是具有256 位向量單元的全功能64 位精簡(jiǎn)指令集（reduced instruction set computing，RISC）核心，而從核是精簡(jiǎn)過(guò)的具有相同長(zhǎng)度向量單元的64位RISC核心。每個(gè)從核都配置了一個(gè)64 KB大小的本地設(shè)備存儲(chǔ)器（local device memory，LDM），其可以配置為用戶可編程的高速緩沖區(qū)或用于數(shù)據(jù)自動(dòng)緩存的軟件模擬緩存。

Fig.1 Architecture of SW26010圖1 申威26010芯片架構(gòu)

神威·太湖之光超級(jí)計(jì)算機(jī)總共有40 960 個(gè)處理器，它們之間通過(guò)中央交換機(jī)網(wǎng)絡(luò)互連，且每256 個(gè)處理器組成一個(gè)完全連接的超節(jié)點(diǎn)。每個(gè)申威處理器具有260 個(gè)核心，其雙精度峰值性能為3.06 TFlops，整機(jī)核心數(shù)目超過(guò)一千萬(wàn)個(gè)，總體性能達(dá)125 PFlops[11]。

2.2 通信模式

神威平臺(tái)目前支持三種主要的編程模型，根據(jù)應(yīng)用的特點(diǎn)，本文工作將研究重點(diǎn)放在了MPI+Athread混合編程模型上，其中Athread是由設(shè)備供應(yīng)商提供的面向細(xì)粒度眾核編程的輕量級(jí)線程庫(kù)。利用Athread在從核陣列上進(jìn)行眾核高效率計(jì)算及從核間通訊，MPI用于不同核組的主核間數(shù)據(jù)傳遞。

基于此，本文工作所涉及的算法和優(yōu)化策略將主要充分利用申威處理器的三種高級(jí)架構(gòu)特性來(lái)獲得高性能，即用戶控制的LDM、直接存儲(chǔ)器訪問(wèn)（DMA）和寄存器數(shù)據(jù)通信機(jī)制。每個(gè)從核上用戶可控制的LDM 可視為私有的用戶可編程存儲(chǔ)空間。DMA 主要用于主存儲(chǔ)器和LDM 之間的高效數(shù)據(jù)移動(dòng)，以便用戶可快速?gòu)腖DM 中訪問(wèn)數(shù)據(jù)而不是從外面的主存儲(chǔ)器中。低延遲的片上寄存器數(shù)據(jù)通信機(jī)制使得在同一核組中的從核之間可以進(jìn)行快速有效的數(shù)據(jù)交換，而無(wú)需從核陣列和主存儲(chǔ)器之間的額外數(shù)據(jù)移動(dòng)。利用MPI的消息傳遞機(jī)制進(jìn)行不同核組的邊界數(shù)據(jù)的傳輸。

3 數(shù)學(xué)描述與優(yōu)化方案設(shè)計(jì)及實(shí)現(xiàn)

3.1 數(shù)學(xué)描述

LQCD 描繪了一個(gè)四維(x,y,z,t)的格點(diǎn)陣列，并把每一維的首尾連接起來(lái)，如圖2 所示，其為二維的示意圖，類似一個(gè)“甜甜圈”。在每個(gè)格點(diǎn)上放置表示夸克的費(fèi)米子場(chǎng)量，其是含有3個(gè)色分量和4個(gè)旋量分量的格拉斯曼數(shù)，這樣每個(gè)格點(diǎn)就成為一個(gè)旋子，其在計(jì)算機(jī)中具體操作時(shí)以12 個(gè)復(fù)數(shù)構(gòu)成的向量表示。在相鄰旋子之間的連接上放置表示膠子的規(guī)范場(chǎng)量，其可以表示為色空間的3×3復(fù)數(shù)矩陣。由于是四維空間，且旋子間的連接又分為正、負(fù)兩個(gè)方向，故每個(gè)格點(diǎn)上將會(huì)有8 個(gè)3×3 的復(fù)矩陣，以連接到8個(gè)鄰居格點(diǎn)。

Fig.2 2-dimensional schematic diagram of LQCD圖2 LQCD二維示意圖

通過(guò)對(duì)每個(gè)格點(diǎn)附加費(fèi)米子常量，每個(gè)連接附加規(guī)范場(chǎng)量來(lái)對(duì)格點(diǎn)的狀態(tài)進(jìn)行描述。不同的狀態(tài)出現(xiàn)的實(shí)際概率不盡相同，此概率正比于e-βSg(U)+ψˉM(U)ψ，其中Sg(U)是全部規(guī)范場(chǎng)量Uμ(x)的一個(gè)泛函，M(U)是與規(guī)范場(chǎng)量相關(guān)的一個(gè)矩陣，其物理含義為全空間中任意兩個(gè)局部夸克場(chǎng)之間的相互作用，這二者定義了LQCD在格點(diǎn)上的物理模型。

對(duì)四維空間中的所有狀態(tài)進(jìn)行積分，通過(guò)歐式空間的路徑積分得到配分函數(shù)Z：

任意物理量O(U,Gq(U))的期望值可寫(xiě)為：

其中，Gq(U)=M-1(U)代表夸克傳播子。

使用蒙特卡羅算法進(jìn)行式（2）積分的計(jì)算僅需兩個(gè)步驟：第一步，隨機(jī)地生成規(guī)范場(chǎng)組態(tài)Uμ(x)，并使其出現(xiàn)的概率為P(U)∝e-βSg(U)detM(U)；第二步，計(jì)算可觀測(cè)量值。其中，由于難以直接計(jì)算detM(U)與M-1(U)，因此通過(guò)采用隨機(jī)偽費(fèi)米子場(chǎng)的方法來(lái)估算detM(U)，從而將其轉(zhuǎn)化為求解M-1(U)，即稀疏矩陣的求逆操作。M-1(U)的物理含義是為了表示任意兩個(gè)局部夸克場(chǎng)之間的傳播振幅，對(duì)于求解M-1(U)，每次只需計(jì)算它的若干列而非整個(gè)逆矩陣，問(wèn)題便轉(zhuǎn)化為求解線性方程組：

其中，b表示四維格點(diǎn)空間，通過(guò)x的求解，方便對(duì)M-1(U)的計(jì)算。

對(duì)于求解上述方程組，只需在{b,Mb,M2b,M3b,…}張開(kāi)的Krylov子空間中迭代逼近即可達(dá)到滿足要求的精度。在如共軛梯度法等Krylov 子空間算法中，關(guān)鍵的操作是矩陣乘以向量和向量?jī)?nèi)積的計(jì)算，其中矩陣向量乘是主要的計(jì)算熱點(diǎn)，而向量?jī)?nèi)積涉及全局規(guī)約操作。

對(duì)于矩陣向量乘主要是計(jì)算以下的公式，即Wilson費(fèi)米子矩陣：

其中，μ表示4 個(gè)時(shí)空方向；x和y是兩個(gè)4-矢量時(shí)空坐標(biāo)；Uμ(x)表示從x點(diǎn)向μ方向的規(guī)范場(chǎng)連接，其為3×3的復(fù)矩陣，稱為“組態(tài)”，U?是它的轉(zhuǎn)置復(fù)共軛；γμ是4×4的稀疏復(fù)矩陣，形式為：

此費(fèi)米子矩陣對(duì)角元部分的計(jì)算相對(duì)簡(jiǎn)單，非對(duì)角元上的計(jì)算才是關(guān)鍵部分，定義：

Wilson費(fèi)米子矩陣每經(jīng)過(guò)一次迭代計(jì)算，空間中的所有旋子就會(huì)更新自身的狀態(tài)，上式定義的Dslash即為重點(diǎn)優(yōu)化的熱點(diǎn)部分。

3.2 代碼結(jié)構(gòu)

本文重點(diǎn)優(yōu)化的包括Dslash 熱點(diǎn)部分在內(nèi)的Wilson費(fèi)米子矩陣主要由3個(gè)函數(shù)所構(gòu)成，它們的依賴關(guān)系如圖3 所示（箭頭指向的是被調(diào)用者）。這3個(gè)函數(shù)分 別 為mr_solver(…)、m_wilson(…)、dslash(…)，它們的代碼結(jié)構(gòu)如下：

Fig.3 Dependences of 3 important functions圖3 3個(gè)主要函數(shù)的依賴關(guān)系

LQCD的整個(gè)計(jì)算流程可簡(jiǎn)化為：在求解線性方程組的迭代條件下，將表示費(fèi)米子場(chǎng)量和規(guī)范場(chǎng)量的復(fù)矩陣通過(guò)Wilson費(fèi)米子矩陣調(diào)用Dslash運(yùn)算以遍歷更新所有格點(diǎn)的數(shù)據(jù)，然后進(jìn)行規(guī)約運(yùn)算，最后計(jì)算殘差，輸出結(jié)果解。如圖4展示了上述基本流程。

Fig.4 Flow of LQCD圖4 LQCD計(jì)算流程

3.3 單核組優(yōu)化方案

從式（6）可以看出，Dslash 運(yùn)算可以歸為一種四維的9-points stencil運(yùn)算，如圖5所示為二維示意，即每計(jì)算一個(gè)格點(diǎn)的數(shù)據(jù)均需要訪問(wèn)4個(gè)維度上相鄰8個(gè)格點(diǎn)的數(shù)據(jù)及8個(gè)方向上的連接的數(shù)據(jù)。

Fig.5 2-dimensional schematic diagram of 4-dimensional 9-points stencil圖5 四維9-points stencil的二維示意圖

關(guān)于stencil 計(jì)算，目前也已有很多的研究成果。如：Datta等人[12]針對(duì)stencil計(jì)算開(kāi)發(fā)了一套行之有效的優(yōu)化策略，并以此構(gòu)建了一個(gè)自動(dòng)調(diào)優(yōu)的機(jī)制，其可通過(guò)搜索優(yōu)化策略組及其參數(shù)來(lái)選出最優(yōu)的策略，從而最大限度地減少運(yùn)行時(shí)間，同時(shí)提高性能。此機(jī)制在Intel Clovertown、IBM QS22 PowerXCell 8等多核架構(gòu)上進(jìn)行了驗(yàn)證，并取得了很好的效果。Ho等人[13]通過(guò)利用本地存儲(chǔ)器和異步軟件預(yù)取技術(shù)來(lái)優(yōu)化三維Lattice Boltzmann（LBM）stencil 模型，并在Kalray MPPA-256 多核處理器上獲得了33%的性能提升。Bondhugula 等人[14]提出了一種全新的平鋪技術(shù)——Diamod Tiling，其解決了現(xiàn)有自動(dòng)平鋪框架流水啟動(dòng)和負(fù)載不均衡的問(wèn)題，使得可以并發(fā)啟動(dòng)且保證負(fù)載均衡，從而最大化stencil 計(jì)算的并行化。Ao 等人[15]在算法層面提出了一種計(jì)算-通信重疊技術(shù)，以減少進(jìn)程間的通信開(kāi)銷；一種局部感知的阻塞方法，以增強(qiáng)數(shù)據(jù)局存性，從而充分利用片上并行性；以及一種針對(duì)不同線程中共享數(shù)據(jù)的聯(lián)合訪問(wèn)模式。以此，他們將三維非靜力大氣模擬模型在神威·太湖之光上進(jìn)行了實(shí)現(xiàn)，并獲得了26 PFlops 的stencil計(jì)算性能。

根據(jù)Dslash運(yùn)算的特點(diǎn)，目前實(shí)現(xiàn)的Dslash代碼是根據(jù)從核LDM 的大小專門設(shè)計(jì)的，希望盡可能地把每個(gè)從核需要處理的數(shù)據(jù)存放在它上面的64 KB大小的LDM 上，從而最小化主存儲(chǔ)DMA 開(kāi)銷。在此設(shè)計(jì)下，每個(gè)核組負(fù)責(zé)84大小的子格子空間，每個(gè)從核負(fù)責(zé)64 個(gè)格點(diǎn)，且這64 個(gè)點(diǎn)采用指定z和t的一個(gè)xy平面，如圖6所示。

3.3.1 代碼向量化及循環(huán)展開(kāi)

Fig.6 Schematic diagram of 4-demensional lattice data spreading on CPE cluster圖6 四維格點(diǎn)數(shù)據(jù)在從核陣列上的分布示意圖

在實(shí)際開(kāi)發(fā)過(guò)程中，發(fā)現(xiàn)申威處理器針對(duì)簡(jiǎn)單循環(huán)以外的其他結(jié)構(gòu)并沒(méi)有提供自動(dòng)向量化（如SLP（successive linear programming）算法）功能，因此在無(wú)法利用編譯器進(jìn)行向量化的情況下只能手動(dòng)向量化代碼。申威處理器每次能夠處理4個(gè)浮點(diǎn)數(shù)值，且旋量的4個(gè)分量剛好可以一并處理，于是修改了數(shù)據(jù)的存儲(chǔ)格式，將旋量指標(biāo)放到最內(nèi)層，這樣矩陣元素在計(jì)算時(shí)被逐個(gè)處理，轉(zhuǎn)置也不會(huì)帶來(lái)額外的開(kāi)銷。Dslash中有8項(xiàng)不同的操作，不同的項(xiàng)需要考慮轉(zhuǎn)置或不轉(zhuǎn)置、1±γμ的符號(hào)等情況進(jìn)行分別處理。另外，SU(3)矩陣有18個(gè)浮點(diǎn)數(shù)，而向量的存儲(chǔ)是4個(gè)浮點(diǎn)數(shù)對(duì)齊的，這就意味著如果某個(gè)方向上的矩陣是對(duì)齊的，那么隨后的一個(gè)就剛好不是對(duì)齊的，因此對(duì)齊與不對(duì)齊的情況也要分別處理。通過(guò)使用宏嵌套定義來(lái)分別對(duì)8個(gè)項(xiàng)生成對(duì)應(yīng)的函數(shù)定義，從而對(duì)數(shù)據(jù)對(duì)齊或不對(duì)齊等各種不同的情況分別進(jìn)行處理。

除了Dslash以外，其他涉及旋量場(chǎng)操作的函數(shù)也進(jìn)行了向量化，主要是把循環(huán)內(nèi)的操作改寫(xiě)成向量的形式。

指令流水線是否一直通暢并不斷發(fā)射指令是影響從核計(jì)算效率的重要因素之一。循環(huán)是計(jì)算中最常見(jiàn)的結(jié)構(gòu)，每次循環(huán)都需要進(jìn)行分支判斷，而條件判斷會(huì)采用分支預(yù)測(cè)技術(shù)，但申威編譯器的分支預(yù)測(cè)功能尚不健全，若直接采用其默認(rèn)的分支預(yù)測(cè)功能就會(huì)造成流水線的不通暢，因此直接將Dslash中的64個(gè)時(shí)空點(diǎn)相關(guān)的循環(huán)全部展開(kāi)。代碼徹底展開(kāi)后增加了代碼的靈活性，其中有關(guān)和主存進(jìn)行DMA傳輸?shù)拇a部分是異步的，這樣就可以通過(guò)巧妙地安排代碼次序來(lái)實(shí)現(xiàn)隱藏通訊時(shí)間。

3.3.2 從核寄存器通訊機(jī)制

從核寄存器通訊機(jī)制是申威處理器的一個(gè)非常特殊的機(jī)制。在實(shí)踐過(guò)程中，感覺(jué)此機(jī)制使得從核上的程序并不像線程，而更接近于進(jìn)程，這就直接改變了很多編程的基本思路。寄存器通訊只能發(fā)生在同行或同列的從核上，因此安排每個(gè)從核上的計(jì)算任務(wù)時(shí)，這種通訊關(guān)系必須滿足。也就是說(shuō)，對(duì)于任何一個(gè)時(shí)空點(diǎn)而言，相鄰的8個(gè)時(shí)空點(diǎn)要么在同一個(gè)從核上，要么在同行或同列的其他從核上。若不滿足這個(gè)條件，數(shù)據(jù)傳輸就需要中轉(zhuǎn)，代碼因此會(huì)變得很復(fù)雜，效率也會(huì)大大降低。如前所述，在目前的版本中，采用簡(jiǎn)單的分配方案，將從核的行和列映射到四維模型的z和t維度，然后給每個(gè)從核分配一個(gè)xy平面的8×8的薄片。此方案的好處是代碼清晰直觀，但效率不是最好的，一方面由于薄片的面積較大，另一方面是與主存通信的DMA 數(shù)量不平衡。在后續(xù)版本中，將采用參數(shù)化代碼自動(dòng)生成方案，為每個(gè)從核生成其特定的任務(wù)代碼，使得從核間協(xié)調(diào)工作，并使用遺傳算法、模擬退火算法等隨機(jī)算法來(lái)優(yōu)化從核的代碼設(shè)計(jì)，使得從核上的任務(wù)操作順序最優(yōu)?；诒∑桨福陂_(kāi)發(fā)過(guò)程中了解到寄存器通信涉及到3 套隊(duì)列，分別是“發(fā)送”“路上”“接收”，且它們有各自不同的大小和共用規(guī)則，非常復(fù)雜。如果接收從核的接收隊(duì)列滿了，路上隊(duì)列也滿了，那么有些數(shù)據(jù)就會(huì)停留在發(fā)送從核的發(fā)送隊(duì)列中，在同步過(guò)后，不同從核的隊(duì)列轉(zhuǎn)移進(jìn)路上隊(duì)列的次序并不能保證是最初的發(fā)送次序，從而導(dǎo)致次序錯(cuò)亂。這個(gè)原理弄清楚后，讓后一個(gè)發(fā)送從核與當(dāng)前接收從核同步，保證發(fā)送的時(shí)候接收隊(duì)列是空的，從而確保發(fā)送次序的穩(wěn)定。

3.3.3 主-從核DMA異步通信

從核LDM與主存之間的數(shù)據(jù)傳輸采用的是異步DMA。神威架構(gòu)定制的Athread 庫(kù)包含了athread_get()等API（application programming interface）接口，這些接口一般都是一些宏，包裝了編譯器實(shí)現(xiàn)的“內(nèi)置函數(shù)”，如dma()等。這些內(nèi)置函數(shù)一般會(huì)被發(fā)射為若干行匯編代碼，大致相當(dāng)于嵌入?yún)R編。這三層接口的效率由慢到快，在程序的不同部分會(huì)使用不同層級(jí)的接口，其中對(duì)計(jì)算要求較高的地方采用直接嵌入?yún)R編代碼或內(nèi)置宏接口的方式調(diào)用所需功能。在DMA 操作方面，經(jīng)測(cè)試Athread 的效率明顯較低，因此直接采用dma_set_*()系列宏接口來(lái)實(shí)現(xiàn)DMA 傳輸。在主核和從核相互傳輸數(shù)據(jù)時(shí)，雙方需要保證對(duì)“進(jìn)度”的統(tǒng)一認(rèn)知，這可以通過(guò)在主存上設(shè)置信號(hào)量來(lái)實(shí)現(xiàn)。具體編寫(xiě)代碼時(shí)，無(wú)論是普通變量還是指針變量，都需要用volatile 關(guān)鍵字來(lái)修飾（其中對(duì)于指針變量的修飾，放置位置不同含義也不同），以保證緩存一致性。對(duì)于信號(hào)量的更新操作，申威處理器提供了3個(gè)指令，它們各自有相應(yīng)的包裝宏來(lái)調(diào)用，在實(shí)際使用時(shí)，重寫(xiě)了這些宏，以保證代碼正確性。

3.4 基于狀態(tài)機(jī)模型的多核組MPI通訊

在整個(gè)LQCD程序優(yōu)化過(guò)程中，設(shè)計(jì)的整體思路是從核負(fù)責(zé)計(jì)算，主核負(fù)責(zé)通信。故在從核優(yōu)化的基礎(chǔ)上，使用MPI消息傳遞機(jī)制來(lái)進(jìn)行核組間通信，以擴(kuò)展并行規(guī)模。

核與核之間的通信主要出現(xiàn)在熱點(diǎn)計(jì)算dslash函數(shù)中，每個(gè)核組負(fù)責(zé)存儲(chǔ)84數(shù)量的格點(diǎn)，并負(fù)責(zé)本區(qū)域內(nèi)格點(diǎn)的計(jì)算。邊界從核在對(duì)四維空間邊界上的格點(diǎn)進(jìn)行計(jì)算時(shí)，需要從存儲(chǔ)相鄰格點(diǎn)的相鄰核組中的從核中獲取數(shù)據(jù)，同時(shí)也需要將自身的邊界數(shù)據(jù)發(fā)送給對(duì)方。由于是四維空間，每個(gè)維度上都需要交換數(shù)據(jù)，且分前、后兩個(gè)方向，那么就需要8次發(fā)送數(shù)據(jù)與8次接收數(shù)據(jù)，每次發(fā)送和接收的數(shù)據(jù)量為四維空間某個(gè)邊界上的數(shù)據(jù)，即83數(shù)量的格點(diǎn)。整個(gè)計(jì)算需要多次執(zhí)行dslash函數(shù)，因此相應(yīng)的主核之間也需要進(jìn)行多次MPI通信。

3.4.1 建立笛卡爾通信域

由于每個(gè)核組對(duì)應(yīng)著四維空間的數(shù)據(jù)，那么核組之間就有了空間關(guān)系，為了更好地表示這種空間位置關(guān)系，建立了對(duì)應(yīng)維度和維度長(zhǎng)度的MPI 笛卡爾四維torus通信域，即此通信域的維度為4(x、y、z、t)，且每個(gè)維度有其對(duì)應(yīng)的進(jìn)程數(shù)。這樣在通信中就可以很方便地通過(guò)進(jìn)程自身在笛卡爾通信域中的坐標(biāo)（MPI_Cart_coords）找到其對(duì)應(yīng)邊界相鄰進(jìn)程的坐標(biāo)（MPI_Cart_shift）。

3.4.2 非阻塞輪詢監(jiān)聽(tīng)及持久通信

在一次dslash 執(zhí)行過(guò)程中，需要進(jìn)行8 次數(shù)據(jù)發(fā)送和接收操作，如果使用阻塞函數(shù)（MPI_Send/MPI_Recv）來(lái)進(jìn)行發(fā)送和接收，那么就需要很繁瑣的設(shè)計(jì)來(lái)避免進(jìn)程間的死鎖。同時(shí)存放在從核陣列中的格點(diǎn)數(shù)據(jù)的狀態(tài)在同一時(shí)刻不盡相同，因此主核并不能隨時(shí)對(duì)邊界數(shù)據(jù)進(jìn)行發(fā)送，亦不能隨時(shí)提供從核所需要的邊界數(shù)據(jù)?；诖耍瑢?duì)發(fā)送和接收采用非阻塞函數(shù)（MPI_Isend以及MPI_Irecv）[16]，并考慮到需要多次進(jìn)行dslash計(jì)算，對(duì)非阻塞發(fā)送和接收采用持久通信方式（MPI_Send_init 和MPI_Recv_init）以減小通信開(kāi)銷。但是采用非阻塞方式，就需要等待或者查看非阻塞函數(shù)是否完成，還涉及到主從核如何彼此相互通知數(shù)據(jù)已交予對(duì)方，因此采用模擬網(wǎng)絡(luò)socket 通信中的I/O 復(fù)用機(jī)制select 方式來(lái)輪詢監(jiān)聽(tīng)發(fā)送和接收數(shù)據(jù)事件。算法如下：

將8 個(gè)方向的數(shù)據(jù)的發(fā)送和接收看成是16 個(gè)事件（8 個(gè)發(fā)送事件、8 個(gè)接收事件）并進(jìn)行監(jiān)聽(tīng)。主核層面設(shè)立8種狀態(tài)來(lái)進(jìn)行事件描述，分別為發(fā)送過(guò)程的4種狀態(tài)NOT_READY（從核并沒(méi)有準(zhǔn)備好主核需要發(fā)送的數(shù)據(jù)）、NOT_SEND（從核數(shù)據(jù)已準(zhǔn)備好，但主核尚未發(fā)送）、SENDING（數(shù)據(jù)發(fā)送中）、SEND_SUCCESS（數(shù)據(jù)發(fā)送成功）和接收過(guò)程的4 種狀態(tài)NOT_RECV（還未接收數(shù)據(jù)）、RECVING（接收中）、RECV_READY（主核接收完成，等待從核接收）、RECV_SUCCESS（從核接收完成）。

每次輪詢前對(duì)發(fā)送和接收事件的狀態(tài)進(jìn)行初始化。輪詢時(shí)，對(duì)每個(gè)事件查看其狀態(tài)，并執(zhí)行相應(yīng)的操作。若某個(gè)事件狀態(tài)對(duì)應(yīng)的動(dòng)作未執(zhí)行完，那么就立即查看下一個(gè)事件的狀態(tài)；若某個(gè)事件的狀態(tài)對(duì)應(yīng)的函數(shù)或者變量發(fā)生變化，那么就改變此事件的狀態(tài)并查看下一個(gè)事件的狀態(tài)。發(fā)送和接收事件的狀態(tài)都顯示為成功時(shí)即標(biāo)識(shí)著一次輪詢的結(jié)束。

另外主從核之間需要特定標(biāo)志來(lái)判斷從核已準(zhǔn)備好主核需要發(fā)送的數(shù)據(jù)，主核已接收到從核所需的數(shù)據(jù)。在設(shè)計(jì)方案中為這些事件各設(shè)立了一個(gè)變量，用以表示主核需發(fā)送的數(shù)據(jù)有多少?gòu)暮艘呀?jīng)放置到主核上，以及主核已經(jīng)接收到數(shù)據(jù)并通知從核進(jìn)行獲取。

主核間除上述的通過(guò)MPI 通訊來(lái)傳輸格點(diǎn)數(shù)據(jù)，還需在求解殘差過(guò)程中對(duì)每個(gè)核組中的殘差數(shù)據(jù)進(jìn)行全局規(guī)約。

3.5 并行文件讀寫(xiě)

MPI-2 中提供了一套并行I/O 接口，實(shí)現(xiàn)了進(jìn)程間并行讀寫(xiě)文件的方法。LQCD 中每個(gè)進(jìn)程對(duì)應(yīng)的數(shù)據(jù)在數(shù)據(jù)文件中并不是連續(xù)存放的。進(jìn)程數(shù)據(jù)分布方式為循環(huán)塊分布，通過(guò)建立新的數(shù)據(jù)類型構(gòu)造符（MPI_Type_create_darray）來(lái)定義分布式數(shù)組類型。打開(kāi)文件之后（MPI_File_open），設(shè)置進(jìn)程的文件視圖（MPI_File_set_view），使得每個(gè)進(jìn)程通過(guò)各自的文件視圖訪問(wèn)數(shù)據(jù)時(shí)，文件中只包含視圖對(duì)應(yīng)的數(shù)據(jù)，且數(shù)據(jù)之間是沒(méi)有空隙的。

4 性能測(cè)試與分析

為測(cè)試上述兩部分的并行優(yōu)化效果，現(xiàn)分別設(shè)置幾組實(shí)驗(yàn)進(jìn)行驗(yàn)證。首先，測(cè)試從核優(yōu)化版本的效果；其次，測(cè)試MPI版本的擴(kuò)展性效果。測(cè)試標(biāo)準(zhǔn)為熱點(diǎn)計(jì)算部分（即涉及Dslash 運(yùn)算的上述3 個(gè)函數(shù)）迭代150次的運(yùn)行時(shí)間。

4.1 單主核版本與單核組從核優(yōu)化版本

基于84 大小的格點(diǎn)數(shù)量，首先將程序在單個(gè)主核上編譯成功并運(yùn)行，這相當(dāng)于純CPU 版本；然后，將使用從核優(yōu)化策略改進(jìn)過(guò)的異構(gòu)并行版本在單個(gè)核組，即64 個(gè)從核上運(yùn)行。它們的運(yùn)行時(shí)間如表1所示。

Table 1 Time comparison of single MPE version and single CG version表1 單主核版本與單核組從核版本運(yùn)行時(shí)間對(duì)比

從以上實(shí)驗(yàn)數(shù)據(jù)可以看出，目前的從核優(yōu)化策略取得了不錯(cuò)的效果，從核優(yōu)化版本相較單主核版本的加速比達(dá)到了165倍。

4.2 單主核版本與從核優(yōu)化MPE版本

在從核優(yōu)化取得不錯(cuò)加速效果的基礎(chǔ)上，為測(cè)試使用MPI 擴(kuò)展核組規(guī)模的性能情況，將程序規(guī)模擴(kuò)展到16 個(gè)核組上，每個(gè)核組負(fù)責(zé)的格點(diǎn)數(shù)據(jù)量仍為84，則總數(shù)據(jù)量為16×84。與之對(duì)應(yīng)的單主核版本的數(shù)據(jù)量也擴(kuò)大到16×84。它們的運(yùn)行時(shí)間如表2所示。

Table 2 Time comparison of single MPE version and MPI version表2 單主核版本與MPI版本運(yùn)行時(shí)間對(duì)比

從以上實(shí)驗(yàn)數(shù)據(jù)可以看出，MPI版本相較單主核版本依然有不錯(cuò)的運(yùn)行速度提升，但相比上一節(jié)得出的加速比卻有了較大的下降。由于單主核版本的數(shù)據(jù)量增大了16 倍，理論上運(yùn)行時(shí)間也應(yīng)增大16倍，實(shí)際的數(shù)據(jù)也符合這一情況（57.73/3.31=17.44,接近16）。但MPI版本相較單核組從核版本而言，整體數(shù)據(jù)量雖增大了16 倍，但每個(gè)核組負(fù)責(zé)的數(shù)據(jù)量沒(méi)有變，且處理器的數(shù)量也相應(yīng)擴(kuò)大了16倍，那么理想情況下這兩個(gè)版本的運(yùn)行時(shí)間應(yīng)相差不大，然而實(shí)際的數(shù)據(jù)卻相差較大，MPI 版本慢了99.1%。這說(shuō)明，MPI 版本在實(shí)際運(yùn)行過(guò)程中出現(xiàn)了性能瓶頸，具體原因需要進(jìn)一步探究。

進(jìn)一步地，測(cè)試了隨著核組規(guī)模的擴(kuò)大，運(yùn)行時(shí)間的變化趨勢(shì)，如表3和圖7所示。

Table 3 Running time of MPI version of different CG numbers表3 不同核組數(shù)MPI版本的運(yùn)行時(shí)間

Fig.7 Trend of CGs scale and running time圖7 核組規(guī)模與運(yùn)行時(shí)間的變化趨勢(shì)

從圖7 可分析得出，使用MPI 擴(kuò)展核組規(guī)模后，程序運(yùn)行時(shí)間與單核組版本相比都有較大的差距，且隨著核組規(guī)模的擴(kuò)展，運(yùn)行時(shí)間有增大的趨勢(shì)。這進(jìn)一步說(shuō)明了使用MPI 誘發(fā)了性能瓶頸，且隨著核組規(guī)模的擴(kuò)大，程序的整體性能有下降的趨勢(shì)。具體原因?qū)⒃谙乱还?jié)中詳細(xì)討論。

4.3 從核代碼拍數(shù)與MPI代碼拍數(shù)對(duì)比

為了詳細(xì)說(shuō)明使用MPI 在主核間傳輸數(shù)據(jù)引發(fā)性能下降的原因，需要更準(zhǔn)確地分析各部分程序運(yùn)行的時(shí)間。使用神威平臺(tái)提供的程序運(yùn)行拍數(shù)計(jì)數(shù)接口rpcc進(jìn)行測(cè)量，將拍數(shù)乘以處理器的主頻即可準(zhǔn)確得到該段程序運(yùn)行的墻鐘時(shí)間（wall clock time）。

利用神威平臺(tái)提供的計(jì)數(shù)器接口，分別在從核關(guān)鍵代碼和MPI代碼中進(jìn)行插樁，分別得到這些代碼部分的拍數(shù)統(tǒng)計(jì)，測(cè)試數(shù)據(jù)如表4和表5所示。

Table 4 rpcc results of CPE codes表4 從核代碼拍數(shù)測(cè)試結(jié)果

Table 5 rpcc results of MPI codes表5 MPI代碼拍數(shù)測(cè)試結(jié)果

表4中，m_wilson_dslash 表示從核Dslash熱點(diǎn)部分計(jì)算過(guò)程的拍數(shù)；slave_doit 表示整個(gè)從核代碼運(yùn)算過(guò)程的拍數(shù)。

表5中，在3 408次輪詢下，共進(jìn)行了27 264次的MPI 發(fā)送（send）和接收（recv）測(cè)試，它們之間的關(guān)系為TPC=TSC(orTRC)/8。在此基礎(chǔ)上，測(cè)得單次輪詢過(guò)程需17 981 513拍數(shù)，其中一次調(diào)用MPI發(fā)送指令需134 716 拍數(shù)，測(cè)試數(shù)據(jù)發(fā)送完成的過(guò)程需8 131拍數(shù)，且需大約9 次測(cè)試完成此過(guò)程；接收過(guò)程同理。從核發(fā)送到主核上的數(shù)據(jù)進(jìn)行轉(zhuǎn)化需1 791 068拍數(shù)。它們之間的關(guān)系為：

從以上兩個(gè)表格數(shù)據(jù)中可以看出，MPI代碼運(yùn)行的平均拍數(shù)比從核代碼運(yùn)行的平均拍數(shù)大了3 個(gè)數(shù)量級(jí)，這就充分說(shuō)明了MPI 代碼占用了大量的運(yùn)行時(shí)間，且從核代碼的運(yùn)算不能將其隱藏，這是目前LQCD 神威版本最大的性能瓶頸。再具體地分析，MPI 運(yùn)行拍數(shù)中的Data transfering 部分又占據(jù)了大量的運(yùn)行時(shí)間，其占比達(dá)到了895 534×16/17 981 513×100%=79.7%，這部分的功能主要是為了主核間傳輸數(shù)據(jù)，將從核上準(zhǔn)備好的數(shù)據(jù)通過(guò)DMA傳輸給對(duì)應(yīng)的主核，然后主核根據(jù)接收的次序?qū)?shù)據(jù)進(jìn)行重排。

經(jīng)過(guò)實(shí)驗(yàn)數(shù)據(jù)測(cè)試與分析，單核組從核優(yōu)化取得了很好的效果，使用MPI 進(jìn)行規(guī)模擴(kuò)展后亦取得了一定的效果，但同時(shí)也造成了性能下降，通過(guò)利用拍數(shù)測(cè)量工具對(duì)代碼的精確分析，將性能瓶頸定位到了主核對(duì)從核數(shù)據(jù)進(jìn)行重排轉(zhuǎn)化的過(guò)程，這促使下一步的工作需要改進(jìn)這部分的算法設(shè)計(jì)。另一方面，減去這部分的性能障礙，MPI 傳輸過(guò)程的拍數(shù)也要比從核代碼的拍數(shù)高出兩個(gè)數(shù)量級(jí)，因此優(yōu)化MPI傳輸也是下一步工作的重點(diǎn)。

5 總結(jié)與展望

LQCD 是當(dāng)今高能物理領(lǐng)域最前沿的研究方向之一，在大規(guī)模超級(jí)計(jì)算機(jī)上實(shí)現(xiàn)其數(shù)值模擬，對(duì)于一個(gè)國(guó)家占領(lǐng)科技制高點(diǎn)具有重要意義。在本文中，總結(jié)了目前在神威·太湖之光超級(jí)計(jì)算機(jī)上對(duì)LQCD 中的Dslash 熱點(diǎn)部分所做的移植和并行優(yōu)化工作。主要貢獻(xiàn)包含如下幾點(diǎn)：通過(guò)分析LQCD的應(yīng)用特點(diǎn)及數(shù)值特征，首次將其在神威平臺(tái)上實(shí)現(xiàn)了成功移植及運(yùn)行；通過(guò)使用向量化、指令流水線、寄存器通訊機(jī)制等手段在申威26010 處理器上實(shí)現(xiàn)了異構(gòu)眾核并行，并實(shí)現(xiàn)了不錯(cuò)的加速比；在實(shí)現(xiàn)從核陣列并行化的基礎(chǔ)上，進(jìn)一步使用MPI 實(shí)現(xiàn)了多核組連并運(yùn)行，以此實(shí)現(xiàn)一定的并行規(guī)模，并著重分析了在目前版本中其成為性能瓶頸的原因。將LQCD在全自主國(guó)產(chǎn)超算平臺(tái)上實(shí)現(xiàn)并行化進(jìn)行了有力嘗試。測(cè)試數(shù)據(jù)表明，神威平臺(tái)對(duì)實(shí)現(xiàn)大規(guī)模的LQCD并行化具有巨大的潛力，對(duì)我國(guó)未來(lái)在高能物理強(qiáng)相互作用領(lǐng)域的數(shù)值模擬能力建設(shè)與提升具有重大意義。

接下來(lái)將會(huì)在以下幾方面繼續(xù)進(jìn)行研究和開(kāi)發(fā)：

（1）在神威平臺(tái)上對(duì)LQCD 進(jìn)行全面的眾核優(yōu)化。目前版本所使用的格點(diǎn)數(shù)量尚且較小，后面將優(yōu)化所采用的stencil 模型，增大數(shù)據(jù)量，進(jìn)一步發(fā)掘LDM 與寄存器通訊機(jī)制的功能，以更加充分地利用從核陣列的并行計(jì)算能力，提高運(yùn)行效率。

（2）通過(guò)進(jìn)一步的分析，極大地消除MPI 通訊瓶頸阻礙，或?qū)⑹褂肕PI-2 等其他版本進(jìn)行開(kāi)發(fā)，以期進(jìn)一步擴(kuò)大并行規(guī)模，充分挖掘神威平臺(tái)的整體計(jì)算能力。

（3）在并行規(guī)模和效率達(dá)到滿意的程度后，將程序與Chroma框架進(jìn)行整合，便于科研人員使用，為高能物理領(lǐng)域做出些許貢獻(xiàn)。