(山東科技大學(xué) 計算機科學(xué)與工程學(xué)院，山東 青島 266590)

分子動力學(xué)(molecular dynamics，MD)模擬是指使用數(shù)值方法，利用計算機模擬原子核和電子所構(gòu)成的多體系統(tǒng)的運動過程，已被廣泛應(yīng)用于物理、化學(xué)、生物、材料、醫(yī)學(xué)等多個領(lǐng)域，用來研究系統(tǒng)的結(jié)構(gòu)和性質(zhì)[1]。在材料領(lǐng)域虛擬過程工程中,Silicon-Crystal應(yīng)用是研究硅晶體熱傳導(dǎo)性的MD模擬應(yīng)用，有限的計算能力一直是制約模擬效率的瓶頸[2]。近年來，高性能計算技術(shù)的發(fā)展為材料領(lǐng)域的虛擬過程工程提供了可能[3]。

“神威·太湖之光”是世界上首臺運行速度超過十億億次的超級計算機，也是我國第一臺全部采用國產(chǎn)處理器SW26010構(gòu)建的超級計算機[4]。清華大學(xué)付昊桓等[5]在“神威·太湖之光”上，利用OpenACC移植大氣模型CAM應(yīng)用，單核組內(nèi)實現(xiàn)2倍加速，但未進行移植后優(yōu)化；上海交通大學(xué)王一超等[6]利用OpenACC移植并優(yōu)化了磁約束聚變領(lǐng)域GTC-P應(yīng)用，單核組內(nèi)實現(xiàn)2.5倍加速，但缺少對訪存密集型應(yīng)用帶寬訪存優(yōu)化；中國科學(xué)院計算應(yīng)用研究中心張帥等[7]在GPU平臺上對MD模擬進行訪存優(yōu)化，但未提供模擬中跨時間迭代問題的解決方法。CPU、GPU架構(gòu)與SW26010架構(gòu)存在著差異。SW26010采用片上計算陣列群和分布式共享存儲相結(jié)合的異構(gòu)眾核體系架構(gòu)，使得MD模擬應(yīng)用的移植具有更大的靈活性，但也使得移植難度加大，目前對MD模擬移植到“神威·太湖之光”超級計算機上的相關(guān)研究尚未見到。

本文設(shè)計了一種SW26010主從計算并行化方案，實現(xiàn)對Silicon-Crystal應(yīng)用的神威OpenACC移植與優(yōu)化；以數(shù)據(jù)流驅(qū)動的任務(wù)圖并行化方法解決任務(wù)間的峰值訪存、跨時間迭代問題，針對該應(yīng)用訪存密集型特點進行帶寬訪存優(yōu)化。

1 背景介紹

1.1 SW26010處理器架構(gòu)及神威OpenACC執(zhí)行模型

“神威·太湖之光”是中國自主研發(fā)的超級計算機，峰值性能為125.4 PFlops，實測峰值約為93 PFlops。采用新一代的眾核異構(gòu)處理器SW26010(架構(gòu)如圖1)。神威OpenACC程序的執(zhí)行模型是在主核指導(dǎo)下，主從核協(xié)同工作，其加速執(zhí)行模型如圖2所示。

圖1 “SW26010”異構(gòu)眾核架構(gòu)Fig.1 Heterogeneous multi-core processor architecture of “SW26010”

SW26010異構(gòu)眾核架構(gòu)中，各核組之間采用片上網(wǎng)絡(luò)互連，每個核組包含1個主核(management processing element MPE)、1個從核簇(8×8=64個，computing processing element，CPE)、1個協(xié)議處理單元和1個內(nèi)存控制器。核組內(nèi)采用共享存儲架構(gòu)，內(nèi)存與主、從核之間可通過內(nèi)存控制器傳輸數(shù)據(jù)，處理器可通過系統(tǒng)接口與外部設(shè)備相連[8]。申威眾核處理器旨在用少量具備指令級并行能力的管理核心集成眾多面向計算開發(fā)的精簡運算核心高效處理線程級并行，從而大幅提高芯片性能[9]。

程序首先在MPE上啟動，以一個主線程串行執(zhí)行，計算密集區(qū)域則在主線程的控制下作為加速任務(wù)被加載到加速設(shè)備CPE上執(zhí)行[10]。任務(wù)的執(zhí)行過程包括：在CPE設(shè)備內(nèi)存上分配所需的數(shù)據(jù)空間；加載任務(wù)代碼至CPE；任務(wù)將所需的數(shù)據(jù)從MPE傳輸至CPE內(nèi)存；等待數(shù)據(jù)傳輸完成；CPE進行計算并將計算結(jié)果傳送回主存；釋放設(shè)備上的數(shù)據(jù)空間等步驟。

圖2 神威OpenACC執(zhí)行模型Fig.2 Execution model of the Sunway OpenACC

MPE加載一系列任務(wù)到加速設(shè)備上同時執(zhí)行，但這種fork-join模式在訪存帶寬有限的SW26010處理器上易產(chǎn)生峰值訪存問題，使CPE之間相互爭搶帶寬，從而影響計算性能。

1.2 AceMesh編程框架

AceMesh編程框架是面向網(wǎng)格應(yīng)用[11-12]、以數(shù)據(jù)為中心，應(yīng)用于多核、眾核平臺上的數(shù)據(jù)流驅(qū)動并行編程框架。AceMesh并行編程框架通過底層的任務(wù)調(diào)度系統(tǒng)[13](運行時庫)對網(wǎng)格應(yīng)用進行任務(wù)圖并行，其核心設(shè)計思想來源于圖論中的有向無環(huán)圖(directed acyclic graph，DAG)。

任務(wù)調(diào)度系統(tǒng)采用探測-執(zhí)行(inspector-executor)兩階段執(zhí)行的并行模式[14]，該模式對并行區(qū)域進行代碼級調(diào)度。探測階段將代碼區(qū)域的控制流和數(shù)據(jù)流信息提交給運行時系統(tǒng)，由運行時系統(tǒng)根據(jù)任務(wù)間依賴關(guān)系建立任務(wù)依賴圖。執(zhí)行階段以構(gòu)建的任務(wù)圖為基礎(chǔ)，依據(jù)資源配置及利用率搭配不同的任務(wù)調(diào)度策略和算法，動態(tài)的調(diào)度并行任務(wù)。

圖3 AceMesh編程框架任務(wù)調(diào)度系統(tǒng)Fig.3 Task scheduling system ofAceMesh programming framework

AceMesh編程框架任務(wù)調(diào)度系統(tǒng)結(jié)構(gòu)圖如圖3所示，該調(diào)度系統(tǒng)包括四層：

1) 用戶接口層，收集任務(wù)粒度[15]的描述、數(shù)據(jù)流信息、任務(wù)構(gòu)造等信息；

2) 任務(wù)構(gòu)建層，根據(jù)上層用戶提供的信息，在系統(tǒng)內(nèi)部產(chǎn)生任務(wù)、建立依賴和進行任務(wù)圖管理；

3) 任務(wù)調(diào)度層，通過靜態(tài)調(diào)度和動靜結(jié)合調(diào)度兩種方式提供任務(wù)調(diào)度支持。靜態(tài)調(diào)度采用輪詢法按權(quán)值將任務(wù)分配至線程；動靜結(jié)合調(diào)度指靜態(tài)調(diào)度策略與任務(wù)竊取調(diào)度算法[13，16]相結(jié)合，提高任務(wù)數(shù)據(jù)重用率和線程間負(fù)載均衡性；

4) 隊列調(diào)度層，利用線程庫對線程私有并發(fā)任務(wù)調(diào)度隊列進行任務(wù)級調(diào)度。

2 Silicon-Crystal應(yīng)用分析及移植方案設(shè)計

2.1 應(yīng)用算法和模擬過程數(shù)據(jù)特征分析

MD模擬中，通過差分求解牛頓運動方程可得到系統(tǒng)中原子的一系列位形。由于模擬過程中力的計算工作量很大，常用的龍格-庫塔法已不再適用，Silicon-Crystal應(yīng)用中的TP(Tersoff Potent)模塊利用leap-frog算法[17]模擬原子在Tersoff勢能作用下的運動軌跡，在所有的線性微分方程的求解器中都有應(yīng)用。

基于有限差分法leap-frog算法，求解線性常微分方程公式如下：

(1)

(2)

其中，r、V、m、F分別為原子的位置矢量、速度、質(zhì)量、所受勢能力，Δt為計算時間步長。

2.2 TP模塊并行方案

加速線程庫(athread庫)是針對主從加速編程模型所設(shè)計的程序加速庫，旨在使用戶能夠方便、快捷地使用核組內(nèi)的線程進行控制和調(diào)度，從而更好地發(fā)揮組內(nèi)多計算核的性能。本研究使用加速線程庫將TP模塊移植到從核的運算模式如圖4所示。

圖4 TP模塊并行方案設(shè)計Fig.4 Parallel scheme design of TP module

TP模塊的移植主要分為：

1) 計算網(wǎng)絡(luò)劃分。MPE端沿三維空間x、y和z三個方向?qū)?shù)據(jù)區(qū)域按比例分成若干矩形體，每一矩形體計算視為一個任務(wù)。這樣的劃分方式有兩個好處：其一，CPE端得到的數(shù)據(jù)在空間上是連續(xù)的，數(shù)據(jù)塊訪問開銷比較??；其二，分塊內(nèi)中心原子占比相對較高，減少分塊間的原子通信量，提升計算效率。

2) 初始化環(huán)境。CPE端對劃分后任務(wù)內(nèi)的原子信息進行初始化，初始化信息包括原子的位置矢量、加速度、速度等。

3) 計算參數(shù)初始化。初始化MPE端對Tersoff勢能下的離散計算參數(shù)。

4) 力場計算。以任務(wù)為基本單位將原子信息加載至CPE端進行加速計算，首先進行MPE端至CPE端的數(shù)據(jù)拷貝，其次利用CPE端的計算陣列群加速核心計算，最后將計算后的各個任務(wù)原子信息由CPE端傳回MPE端。

5) 同步力場數(shù)據(jù)。MPE端按照鄰居關(guān)系索引表進行任務(wù)間數(shù)據(jù)更新操作，保證數(shù)據(jù)全局一致性。

6) 更新殘量和輸出文件信息。CPE端進行每個時間步計算后的殘量更新，MPE端將計算范數(shù)值輸出至文件系統(tǒng)。

3 神威OpenACC從核移植

神威OpenACC并行編程模型，用編譯指示的方式把應(yīng)用中可并行化的計算循環(huán)移植到申威處理器從核以加速計算。具體到Silicon-Crystal應(yīng)用的從核移植，主要分為以下三個步驟：

1) 循環(huán)并行化。Silicon-Crystal應(yīng)用以任務(wù)分片存儲的數(shù)據(jù)為基本單位進行模擬計算，在分塊級的for循環(huán)上添加相應(yīng)的指導(dǎo)語句#pragma acc parallel loop，將計算部署在64個從核上并行執(zhí)行。

需要注意的是，gang、worker、vecotr是OpenACC2.0中的3層循環(huán)設(shè)計，由于神威眾核架構(gòu)在物理上并沒有分層需求，所以神威OpenACC的實現(xiàn)是把gang設(shè)置成64，worker設(shè)為1。

2) 基于計算數(shù)據(jù)優(yōu)先的數(shù)據(jù)管理。神威眾核架構(gòu)中存在訪存帶寬較小的問題，故從核移植并行化過程最為關(guān)鍵的是將加速計算需要的數(shù)據(jù)提前拷貝到訪問延遲低的SPM(scratch pad memory)。本研究采用計算數(shù)據(jù)優(yōu)先傳輸策略即將所有計算涉及的數(shù)據(jù)優(yōu)先傳至SPM。

計算數(shù)據(jù)優(yōu)先傳輸過程如下：

i) 按循環(huán)索引劃分傳輸。若數(shù)組的索引變量與循環(huán)索引變量緊耦合時，神威OpenACC編譯器將數(shù)組劃分為64份，然后利用DMA的方式將劃分后的數(shù)據(jù)集中傳輸至各從核SPM中，并將任務(wù)內(nèi)的鄰居關(guān)系表、打包后的計算參數(shù)順序傳遞給從核。具體使用copy/copyin/cpoyout等指導(dǎo)語句完成(如圖5)。

①#pragma acc parallel loop?②local(pi) //變量局存私有化③copy(ax,ay,az,pox,poy,poz)annotate(dimension(ax(PN,N),ay(PN,N),az(PN,N))) //按循環(huán)索引傳輸④copyin(maplis) annotate(dimension(maplist(4?N))) //鄰居關(guān)系表傳輸⑤packin(deltaT,soma,Mu,Lamda,Beta,nexp,InteriorCutoff,ExteriorCutoff) //離散計算參數(shù)打包傳輸⑥for (pi = 0; pi < PN ; pi++){…}

圖5 數(shù)據(jù)管理過程的函數(shù)指導(dǎo)語言實現(xiàn)
Fig.5 Implementation of functional instruction language based on data management

ii) 變量局存私有化。對于并行循環(huán)索引變量等線程私有變量，既可使用private子句也可使用local子句將變量私有化，考慮到private是線程私有化變量，變量值仍在主存中，而local是線程私有化的局存變量，存儲在SPM中，數(shù)據(jù)訪問更加高效，故采用local子句進行變量的私有化。

iii) 離散計算參數(shù)打包傳輸。Silicon-Crystal應(yīng)用存在多個離散標(biāo)量的模擬參數(shù)需要傳送至從核，若一一傳輸需要頻繁的使用DMA方式，會大大增加訪存開銷。在此情況下，本文利用pack/packin/packout等神威定制的指導(dǎo)語句將離散數(shù)據(jù)打包后一次傳遞，以更充分有效地利用有限的訪存帶寬。

綜上，得到 Silion-Crystal應(yīng)用移植中的數(shù)據(jù)管理過程的函數(shù)指導(dǎo)語言實現(xiàn)如圖5。

3) 加速代碼區(qū)約束處理。SWACC編譯器進行OpenACC并行化過程中，對并行區(qū)的代碼有一定的要求。如在加速區(qū)代碼中存在函數(shù)調(diào)用時，需在函數(shù)定義處添加routine子句指示，否則生成從核代碼將找不到函數(shù)的位置。但目前routine子句只適用Fortran程序，C代碼暫不支持。Silicon-Crystal應(yīng)用程序是C代碼程序，無法利用routine子句修飾從核函數(shù)。本研究通過利用宏定義實現(xiàn)力場計算的內(nèi)聯(lián)函數(shù)，來解決移植過程中加速區(qū)函數(shù)返回值異常的問題。

圖6 神威OpenACC移植性能數(shù)據(jù)Fig.6 Performance data of the Sunway OpenACC transplant

將運行在1個主核上的Silicon-Crystal作為測試基準(zhǔn)，分別與循環(huán)并行化、基于計算數(shù)據(jù)優(yōu)先的訪存和離散計算參數(shù)打包傳輸3個方面在單核組上進行性能測試(圖6)。測試問題規(guī)模：回環(huán)中存在131 072個粒子，迭代計算次數(shù)為1 000次。

可以看出，對于訪存密集型的應(yīng)用，僅進行循環(huán)并行化將計算過程移至CPE端，性能反而會降低；按照計算數(shù)據(jù)優(yōu)先方式通過DMA方式放入從核SPM中，性能開始超越主核；通過pack子句對離散計算參數(shù)打包后再傳輸，性能進一步提升，整體應(yīng)用較主核版實現(xiàn)了2.26倍的加速。

4 數(shù)據(jù)流驅(qū)動的任務(wù)圖并行化

AceMesh任務(wù)調(diào)度系統(tǒng)的設(shè)計思想來源于數(shù)據(jù)結(jié)構(gòu)中的有向無環(huán)圖，即任務(wù)依賴圖。任務(wù)依賴圖在圖論中是指：如果一個有向圖無法從某個頂點出發(fā)經(jīng)過若干條邊回到該點，則這個圖是一個任務(wù)依賴圖。任務(wù)依賴圖中的頂點代表任務(wù)，圖中的邊代表任務(wù)間的依賴關(guān)系。根據(jù)任務(wù)依賴圖的特點，將并行計算中的大規(guī)模計算問題劃分為N(N≥1)個任務(wù)，并根據(jù)各個任務(wù)的依賴關(guān)系建立任務(wù)依賴圖，圖中所有沒有后繼的頂點都執(zhí)行完后，任務(wù)依賴圖的執(zhí)行完成。

在神威眾核處理器上任務(wù)圖并行化過程分為任務(wù)構(gòu)圖期和任務(wù)執(zhí)行期。構(gòu)圖期是任務(wù)構(gòu)建的探測過程，旨在根據(jù)注冊的數(shù)據(jù)地址去建立任務(wù)間的依賴關(guān)系，在不改變串行序結(jié)果的情況下以數(shù)據(jù)流調(diào)整執(zhí)行序列；執(zhí)行期是指按照構(gòu)圖期間構(gòu)建的DAG圖，搭配任務(wù)調(diào)度系統(tǒng)的不同調(diào)度策略執(zhí)行任務(wù)的過程。故任務(wù)圖并行化總時間等于構(gòu)圖時間(graph time)加執(zhí)行時間(execution time)。

Silicon-Crystal應(yīng)用在太湖之光上使用任務(wù)圖并行化主要分為以下3個步驟：

1)主核構(gòu)建任務(wù)依賴圖。根據(jù)不同并行區(qū)內(nèi)劃分的任務(wù)按照對內(nèi)部訪問的數(shù)據(jù)依賴關(guān)系進行地址注冊，構(gòu)建出任務(wù)執(zhí)行序DAG圖。

圖7為Silicon-Crystal應(yīng)用2線程4任務(wù)依賴圖。其中，每個橢圓代表包裝后的一個任務(wù)，橢圓中第一個數(shù)字為并行區(qū)編號，第二個數(shù)字為任務(wù)編號；箭頭代表任務(wù)間的數(shù)據(jù)訪問先后的依賴關(guān)系；陰影、非陰影圓圈代表執(zhí)行時不同的線程；實線邊是任務(wù)垂直后繼依賴邊，虛線邊是普通后繼依賴邊。在任務(wù)執(zhí)行期間，采用的調(diào)度策略使垂直后繼任務(wù)優(yōu)先于普通后繼任務(wù)執(zhí)行，旨在使任務(wù)間的數(shù)據(jù)重用得到最大化。此外，截斷并行區(qū)間任務(wù)執(zhí)行的依賴關(guān)系，按照與神威OpenACC相同的控制流驅(qū)動的fork-join執(zhí)行模式，稱為任務(wù)圖單步執(zhí)行。

圖7 2線程4任務(wù)時TP模塊任務(wù)依賴圖Fig.7 Task dependency graph of TP module under 2 threads 4 tasks

圖8 神威OpenACC與任務(wù)圖并行化性能Fig.8 Performance of the Sunway Open ACC andtask graph parallelization

2)從核包裝任務(wù)函數(shù)，將應(yīng)用主要的計算代碼包裝成任務(wù)函數(shù)。從核任務(wù)函數(shù)根據(jù)構(gòu)圖期分配的函數(shù)參數(shù)、循環(huán)劃分尺寸、數(shù)據(jù)區(qū)劃分尺寸等信息，包裝任務(wù)函數(shù)，放入從核陣列并行計算。

3)從核數(shù)據(jù)管理。SW26010主從核間的數(shù)據(jù)傳輸通過DMA實現(xiàn)，DMA只能由從核發(fā)起，主核被動進行數(shù)據(jù)傳輸。從核制定傳輸?shù)哪Ｊ綍r，數(shù)據(jù)傳輸依據(jù)數(shù)據(jù)在主存數(shù)據(jù)區(qū)內(nèi)存儲地址的連續(xù)性和從核計算實際需要的數(shù)據(jù)尺寸進行傳輸。DMA數(shù)據(jù)傳輸方式分為跨步式數(shù)據(jù)傳輸和非跨步式數(shù)據(jù)傳輸。兩種傳輸模式下，軟件開銷主要體現(xiàn)在傳輸?shù)膯雍蛯MA傳輸回答字的處理。本研究采用數(shù)據(jù)分片存儲的數(shù)據(jù)結(jié)構(gòu)，將任務(wù)訪問的數(shù)據(jù)進行分片存儲，并以數(shù)據(jù)塊編號為索引劃定數(shù)據(jù)區(qū)，通過athread_get/athread_put接口進行非跨步傳輸；對于離散數(shù)據(jù)訪問，在主核代碼中對離散數(shù)據(jù)打包后，使用DMA方式進行數(shù)據(jù)傳輸，最后在從核代碼中對數(shù)據(jù)進行解包，提高離散數(shù)據(jù)的訪問效率。

本節(jié)將第3節(jié)中OpenACC優(yōu)化后的版本(ACC)作為基礎(chǔ)版，同等的優(yōu)化條件下，與任務(wù)圖單步版(single step of DAG)、任務(wù)圖亂序版(Unordered DAG)進行實驗對比。迭代時間步長為1時，Silicon-Crystal應(yīng)用性能如圖8所示。

可以看出，由于任務(wù)圖單步版和OpenACC采用相同的fork-join模式，即并行區(qū)開始spawn線程，并行區(qū)結(jié)束wait所有線程，故二者的執(zhí)行時間一致，說明兩者具有計算一致性。任務(wù)圖亂序版比任務(wù)單步版提升27%，驗證了任務(wù)間的亂序執(zhí)行，可以錯開峰值帶寬競爭，充分的利用從核訪存帶寬。但是，由于神威OpenACC采用fork-join模式，其執(zhí)行時間即為總時間；數(shù)據(jù)流驅(qū)動的任務(wù)圖并行需要在構(gòu)圖期構(gòu)建任務(wù)執(zhí)行的依賴關(guān)系，故其總時間為構(gòu)圖時間與執(zhí)行時間之和。實驗結(jié)果表明，Silicon-Crystal應(yīng)用的任務(wù)圖并行化存在相對總時間8%的構(gòu)圖時間，加上此部分構(gòu)圖開銷，總時間上任務(wù)圖并行比ACC性能提升11.5%。

表1 多時間步擴展下任務(wù)圖并行化構(gòu)時間Tab.1 Times of task graph parallelization based on multiple time steps

圖9 多時間步下任務(wù)圖并行化加速比Fig.9 Acceleration ratio of task graph parallelizationunder multiple time steps

傳統(tǒng)的fork-join模式無法擴展多時間步的迭代計算，任務(wù)圖卻可打通迭代時間步間的并行區(qū)域，即在構(gòu)圖期依據(jù)多個時間步下任務(wù)的數(shù)據(jù)流構(gòu)建出任務(wù)依賴圖，執(zhí)行期按多時間步下任務(wù)亂序調(diào)度方式執(zhí)行任務(wù)，從而將迭代時間步由單時間步擴展至多時間步，如表1所示，隨著任務(wù)圖并行在多時間步的擴展，執(zhí)行時間進一步降低，構(gòu)圖時間逐漸降低，使得任務(wù)圖并行性能進一步提升。

多時間步下任務(wù)圖并行化加速比如圖9所示，以主核版作為基準(zhǔn)版，使用神威OpenACC移植利用從核加速，實現(xiàn)2.26倍加速比；時間步為1時任務(wù)圖并行加速比為2.52；隨著時間步的擴展，任務(wù)圖規(guī)模隨之增加，任務(wù)的亂序使錯峰訪存的優(yōu)勢進一步擴大，時間步擴展至20時趨于平穩(wěn)，加速比達到3.2。

5 總結(jié)與未來工作

本研究為Silicon-Crystal應(yīng)用設(shè)計了一套在SW26010上實現(xiàn)主從計算的并行化方案，利用OpenACC完成了向目標(biāo)平臺“神威·太湖之光”上的移植，在單核組內(nèi)實現(xiàn)了2.52倍加速；針對該應(yīng)用訪存密集的行為特點，以數(shù)據(jù)流驅(qū)動的任務(wù)圖并行化方法解決任務(wù)間的峰值訪存和跨時間迭代問題，結(jié)果表明，Silicon-Crystal應(yīng)用在數(shù)據(jù)流驅(qū)動的任務(wù)圖并行在單時間步下性能提升11.5%，多時間步下性能提升42%，總體較主核實現(xiàn)3.2倍加速。

數(shù)據(jù)流驅(qū)動的任務(wù)圖并行編程模型采用AceMesh任務(wù)調(diào)度系統(tǒng)中的低級接口對程序源碼進行優(yōu)化，隨著本課題組數(shù)據(jù)驅(qū)動的并行調(diào)度系統(tǒng)自動轉(zhuǎn)譯器的完善，未來將使用指導(dǎo)語言的高級形式對代碼進行自動源源變換，從而實現(xiàn)通過指導(dǎo)語言方式對應(yīng)用任務(wù)圖并行的自動化過程。