王一超 林新華, 蔡林金 Tang William Ethier Stephane 王 蓓 施忠偉,4 松崗聰

1(上海交通大學(xué)高性能計(jì)算中心 上海 200240) 2(東京工業(yè)大學(xué) 日本東京 1528550) 3(普林斯頓大學(xué)等離子體物理實(shí)驗(yàn)室 美國新澤西州普林斯頓 08540) 4 (英偉達(dá)公司 新加坡 138522) (wangyichao@sjtu.edu.cn)

國家超級計(jì)算無錫中心的神威“太湖之光”超級計(jì)算機(jī)在2016年6月發(fā)布的世界Top500超算榜單中排名第一，其理論峰值性能達(dá)到125.4PFLOPS，Linpack實(shí)測效率74.15%，約為93PFLOPS[1].值得一提的是，“太湖之光”搭載的處理器芯片是完全由中國自主研發(fā)的神威異構(gòu)眾核處理器SW26010.該芯片的架構(gòu)面向高性能計(jì)算，通過統(tǒng)一指令系統(tǒng)、統(tǒng)一執(zhí)行模型和支持一致性的主存共享，實(shí)現(xiàn)異構(gòu)核心的深度融合[2].

在軟件層面，“太湖之光”的操作系統(tǒng)和編譯環(huán)境等也全部由中國自主研發(fā).該系統(tǒng)實(shí)現(xiàn)的并行編程模型神威OpenACC，兼容OpenACC 2.0標(biāo)準(zhǔn)，還添加了部分定制功能.基于指導(dǎo)語句的并行編程模型可以幫助用戶在應(yīng)用移植過程中，避免硬件架構(gòu)不同帶來的差異，提供高可移植性.本文將研究在神威眾核處理器上利用OpenACC對于實(shí)際應(yīng)用進(jìn)行移植及優(yōu)化.

由美國普林斯頓大學(xué)開發(fā)的高性能計(jì)算應(yīng)用GTC-P(the Princeton gyrokinetic toroidal code)利用 PIC(particle-in-cell)算法求解Vlasov-Poisson方程，模擬粒子在托卡馬克裝置內(nèi)的運(yùn)動(dòng)，具有重要的物理意義.又由于其極佳的可擴(kuò)展性，GTC-P已在Top500排名前10的6臺超級計(jì)算機(jī)上進(jìn)行了性能測試[3]，并入選了美國NERSC國家超算中心的基準(zhǔn)測試集[4]，是最具代表意義的PIC應(yīng)用.

綜上所述，本文使用OpenACC在“太湖之光”超級計(jì)算機(jī)上成功移植了GTC-P應(yīng)用.針對神威OpenACC編譯器尚無法解決的性能瓶頸，本文提出了3種基于中間代碼二次開發(fā)的優(yōu)化方法.實(shí)驗(yàn)證明優(yōu)化方法有效，同時(shí)我們還對神威OpenACC結(jié)合MPI的多節(jié)點(diǎn)版本進(jìn)行了強(qiáng)可擴(kuò)展性的評估.

簡而言之，本文有以下2點(diǎn)貢獻(xiàn):

1) 首次在“太湖之光”上利用OpenACC成功并行移植了基于PIC算法的實(shí)際應(yīng)用GTC-P.

2) 提出了3種基于中間代碼的手動(dòng)優(yōu)化方法：①消除原子操作；②避免低效的全局訪存操作；③手動(dòng)添加SIMD intrinsics指令.實(shí)驗(yàn)證明優(yōu)化方法對于PIC算法在“太湖之光”上的性能提升效果顯著.

1 相關(guān)工作

本文的相關(guān)工作主要有3類:

1) “太湖之光”上的應(yīng)用移植.清華大學(xué)付昊桓等人[5]利用OpenACC在“太湖之光”上成功移植了大氣模型CAM應(yīng)用.實(shí)驗(yàn)結(jié)果顯示，利用OpenACC移植后的應(yīng)用，相比1個(gè)主核，在64個(gè)從核上整體實(shí)現(xiàn)了2倍的加速比.但該文尚未提出在“太湖之光”上針對OpenACC移植后的性能優(yōu)化方法.

2) GTC-P的移植與調(diào)優(yōu).普林斯頓大學(xué)Wang等人[6]對GTC-P應(yīng)用在Mira,Sequoia,Hopper等超級計(jì)算機(jī)上的移植與調(diào)優(yōu)做了介紹，并評估了可擴(kuò)展性.在異構(gòu)系統(tǒng)上的移植方面，美國勞倫斯伯克利國家實(shí)驗(yàn)室Ibrahim等人[7]介紹了利用CUDA在基于GPU的Titan超級計(jì)算機(jī)上移植GTC-P的方法.但這些工作都是在基于CPU,GPU,MIC加速卡構(gòu)建的超算系統(tǒng)上完成的，與神威眾核的架構(gòu)存在差異.

3) 利用OpenACC移植應(yīng)用. Calore等人[8]和Hariri等人[9]分別利用OpenACC對基于Lattice Boltzmann和particle-in-cell算法的科學(xué)應(yīng)用進(jìn)行了并行移植，并介紹了常用的優(yōu)化手段，移植平臺為GPU和CPU.東京工業(yè)大學(xué)Hoshino等人[10]對于利用OpenACC在GPU上進(jìn)行并行移植存在的性能瓶頸問題做了深入的研究，并指出缺乏對于GPU shared memory的控制是導(dǎo)致性能不佳的原因.但本文中所提及的神威OpenACC實(shí)現(xiàn)上與OpenACC標(biāo)準(zhǔn)存在差異.

2 背景介紹

本文的研究背景介紹將從神威眾核處理器架構(gòu)、神威OpenACC以及GTC-P應(yīng)用展開.

2.1 神威眾核處理器架構(gòu)

神威眾核處理器旨在用少量具備指令級并行能力的管理核心，集成眾多面向計(jì)算開發(fā)的精簡運(yùn)算核心高效處理線程級并行，從而大幅提高芯片性能[2].

“太湖之光”超級計(jì)算機(jī)上所使用的神威SW26010芯片架構(gòu)如圖1所示，1塊芯片由4個(gè)核組組成，每個(gè)核組內(nèi)64個(gè)從核按照8×8的mesh拓?fù)浣Y(jié)構(gòu)，由片上內(nèi)部網(wǎng)絡(luò)互聯(lián).

Fig. 1 The architecture of Sunway many-core processor圖1 神威眾核處理器架構(gòu)圖

此外，主從核均支持256b向量化指令，以及乘加融合(fused multiply add, FMA)指令.每個(gè)從核支持一條浮點(diǎn)數(shù)DP(double precision)流水線，主核支持雙浮點(diǎn)數(shù)流水線.由此可得，主從核雙精度理論峰值性能的計(jì)算為

1×1.45×8×2=23.2GFLOPS;

64×1.45×4×2=742.4GFLOPS.

不難發(fā)現(xiàn)，SW26010上逾98%的理論性能均源自于從核，所以應(yīng)用必須在從核上得到并行化移植才能充分發(fā)揮“太湖之光”超級計(jì)算機(jī)的性能.

存儲架構(gòu)方面，主核擁有2級緩存：第1級緩存指令數(shù)據(jù)分離，即指令和數(shù)據(jù)各自擁有32 KB緩存；第2級則為指令數(shù)據(jù)共享的256 KB緩存.而每個(gè)從核擁有16 KB的一級指令緩存.從核數(shù)據(jù)存儲的設(shè)計(jì)是一種類CELL處理器的用戶控制的草稿本方式存儲，即SPM(scratch pad memory)，每個(gè)從核的SPM容量為64 KB.這種存儲方式不僅省去了緩存實(shí)現(xiàn)上的控制開銷，還避免了眾多運(yùn)算核心間一致性處理帶來的設(shè)計(jì)復(fù)雜性和性能降級[2].但也為應(yīng)用開發(fā)中從核的數(shù)據(jù)管理帶來了新挑戰(zhàn)，即需要人為對于應(yīng)用在SPM上的訪存進(jìn)行設(shè)計(jì)與規(guī)劃.

SW26010片上集成了4個(gè)內(nèi)存控制器，每個(gè)核組擁有8 GB內(nèi)存.但在具體實(shí)現(xiàn)上，內(nèi)存管理器仍是基于DDR3的，芯片的實(shí)測訪存帶寬為136.51 GBps.

基于上述設(shè)計(jì)，相比NVIDIA GPU K20X或者Intel Xeon Phi 7110P在Stream (Triad)測試中實(shí)測180 GBps左右的訪存帶寬，SW26010的帶寬差了近25%.由此，在SW26010上移植應(yīng)用，性能更容易受限于訪存帶寬.充分利用SPM局存訪問，緩解主存訪問壓力，就成了應(yīng)用移植和優(yōu)化的關(guān)鍵之一.

2.2 神威OpenACC

神威OpenACC是“太湖之光”上基于定制OpenACC 2.0標(biāo)準(zhǔn)實(shí)現(xiàn)的并行編程模型，面向從核的并行化移植.其指導(dǎo)語句的設(shè)計(jì)與OpenACC 2.0標(biāo)準(zhǔn)存在一定區(qū)別，具體情況如表1所示：

Table 1 Sunway OpenACC vs OpenACC 2.0 Standard表1 神威OpenACC與OpenACC 2.0標(biāo)準(zhǔn)的對比

這樣的設(shè)計(jì)主要出于3方面的考慮：

1) 由于神威架構(gòu)與GPU和MIC等主流異構(gòu)加速器存在區(qū)別，導(dǎo)致如tile,cache等語句在實(shí)現(xiàn)上與標(biāo)準(zhǔn)OpenACC的實(shí)現(xiàn)產(chǎn)生了差異.

3) 自主研發(fā)的神威OpenACC編譯器，與PGI，Cray等成熟的商業(yè)編譯器相比，還存在不完善之處.尤其在數(shù)據(jù)管理方面，需要指導(dǎo)語句給予編譯器更多的提示，這反映在了定制的annotate語句中.

此外，神威OpenACC編譯器是基于ROSE開源框架開發(fā)的代碼轉(zhuǎn)譯器，主要負(fù)責(zé)中間代碼生成.其后端仍使用的是“太湖之光”的本地編譯器.該編譯器向用戶開放中間代碼，允許用戶基于編譯生成的中間代碼進(jìn)行二次開發(fā).第4節(jié)的性能優(yōu)化中，本文就會(huì)用到該功能.

2.3 GTC-P應(yīng)用

GTC-P是由美國普林斯頓大學(xué)等離子體物理實(shí)驗(yàn)室開發(fā)的應(yīng)用于磁約束聚變領(lǐng)域中全局性大規(guī)模并行模擬的高性能計(jì)算程序.該應(yīng)用在近10年來一直都是研究托卡馬克裝置內(nèi)等離子體湍流運(yùn)動(dòng)等重要非線性物理問題的主要工具.

粒子在1束激光中的運(yùn)動(dòng)可以通過六維的Vlasov-Maxwell方程來模擬[11].而對于托卡馬克這種帶有零階穩(wěn)態(tài)平衡導(dǎo)向場的等離子體系統(tǒng)，回旋動(dòng)理學(xué)的簡化降維又是非常有效的.在這種情況下，磁場中的帶電粒子的回旋運(yùn)動(dòng)將近似于一個(gè)帶電環(huán).而PIC算法則是現(xiàn)如今最為普及的求解該系統(tǒng)方程的數(shù)值計(jì)算方法.

基于回旋動(dòng)理學(xué)的PIC方法的每個(gè)時(shí)間步主要分為5個(gè)計(jì)算步驟，對應(yīng)到代碼中的主要函數(shù).1)charge.使用4點(diǎn)陀螺平均法求解粒子到網(wǎng)格上的電荷沉積.2)poisson.在網(wǎng)格上求解回旋動(dòng)理學(xué)的Poisson方程.3)field.計(jì)算網(wǎng)格上的電場.4)push.用矢量場和其他衍生物助推粒子.5)smooth.使電荷密度和潛在矢量平滑.其中，charge和push的計(jì)算強(qiáng)度均小于2.結(jié)合SW26010的峰值性能和訪存帶寬可知，計(jì)算強(qiáng)度低于20的應(yīng)用在神威平臺上的性能就會(huì)受限于訪存性能，因此這2個(gè)函數(shù)的訪存實(shí)現(xiàn)將會(huì)成為移植過程中的主要難點(diǎn).

3 利用OpenACC移植GTC-P

在使用神威OpenACC進(jìn)行從核并行化移植之前，本文首先將GTC-P在x86平臺上實(shí)現(xiàn)的串行代碼移植到SW26010的1個(gè)主核上.并記錄了各部分函數(shù)的執(zhí)行時(shí)間，分析了程序熱點(diǎn)，結(jié)果如圖2所示.

Fig. 2 Profiling the hotspots of GTC-P code圖2 GTC-P代碼熱點(diǎn)分析

可見，主核上的熱點(diǎn)測試結(jié)果與2.3節(jié)的分析相印證，即訪存密集型的函數(shù)charge和push是程序的主要熱點(diǎn).本文將對于GTC-P應(yīng)用中的這2個(gè)函數(shù)進(jìn)行從核并行化移植與優(yōu)化，移植工作分為3個(gè)步驟：數(shù)據(jù)管理、循環(huán)并行化、原子操作.

3.1 數(shù)據(jù)管理

如2.1節(jié)中提到的，在神威眾核架構(gòu)中，存在著訪存帶寬性能有限的潛在瓶頸.所以從應(yīng)用性能角度出發(fā)，將計(jì)算涉及的數(shù)據(jù)提前拷貝至訪問延遲低的SPM局存中，是從核并行化移植的關(guān)鍵步驟.OpenACC 2.0標(biāo)準(zhǔn)中的data指導(dǎo)語句提供了如copy,present,cache等多種數(shù)據(jù)管理模式，在神威OpenACC中為提升DMA帶寬添加了swap,pack等定制語句.本文采用了將所有在計(jì)算中涉及的數(shù)據(jù)優(yōu)先傳輸至SPM的實(shí)現(xiàn)策略，具體實(shí)現(xiàn)方法為：

1) 按循環(huán)索引劃分傳輸.若數(shù)組的索引變量與循環(huán)的索引變量緊耦合時(shí)，神威OpenACC編譯器會(huì)自行將數(shù)組劃分為64份，然后利用DMA的方式分別將64份數(shù)據(jù)傳輸至各從核的SPM中.具體使用如圖3所示，其中 annotate(dimension(array(size))) 指導(dǎo)語句是為了指明數(shù)組長度，方便編譯器劃分?jǐn)?shù)據(jù).

2) 完整傳輸.當(dāng)訪問的數(shù)組與循環(huán)索引無關(guān)聯(lián)時(shí)，由于SW26010缺少共享的局部存儲，編譯器將無法判斷如何劃分?jǐn)?shù)據(jù)，從而默認(rèn)不會(huì)把該數(shù)組傳輸至SPM中.在此情況下，可以通過annotate(entire(array))指導(dǎo)語句提示編譯器將數(shù)組的完整副本分別傳輸至64個(gè)從核的SPM上.

4) 使用緩存.OpenACC 2.0中的cache指導(dǎo)語句在SW26010上實(shí)際對應(yīng)的是SPM.又因?yàn)镾PM沒有實(shí)現(xiàn)硬件層的緩存管理機(jī)制，所以神威OpenACC的cache指導(dǎo)語句背后是一個(gè)軟件實(shí)現(xiàn)緩存機(jī)制.

① #pragmaaccparallelloop② copyin(z0,z1,z2,z4,z5)annotate(dimension(z0(mi),z1(mi),z2(mi),z4(mi),z5(mi)))③ copyin(delt,mtheta,igrid,qtinv)annotate(dimension(delt(90+1),mtheta(90+1),igrid(90+1),qtinv(90+1)),entire(delt,mtheta,igrid,qtinv))④ copyout(wpion,wtion0,wtion1,jtion0,jtion1)annotate(dimension(wpion(4*mi),wtion0(4*mi),wtion1(4*mi),jtion0(4*mi),jtion1(4*mi)))⑤ cache(pgyro,tgyro)⑥ for(m=0;m

Fig. 3 charge function implementation based on 4 data management strategies of tile
圖3 tile基于4種數(shù)據(jù)管理策略的函數(shù)charge實(shí)現(xiàn)

基于以上4種數(shù)據(jù)管理方法的函數(shù)charge實(shí)現(xiàn)如圖3所示.其中，copyin語句中的z0，z1，z2等數(shù)組編譯后會(huì)基于循環(huán)索引m均等劃分，分別傳輸至64個(gè)從核的SPM中.而delt,mtheta,igrid,qtiny這4個(gè)數(shù)組則將被完整地傳輸至64個(gè)從核的SPM中，等同于在每個(gè)從核的SPM中保存了一份這4個(gè)數(shù)組的副本，方便從核隨時(shí)進(jìn)行訪問.copyout指導(dǎo)語句中的數(shù)組傳輸方法與z0等相同，按照循環(huán)索引進(jìn)行劃分.值得一提的是cache指導(dǎo)語句中的pgyro和tgyro，這2個(gè)數(shù)組的訪存是不規(guī)則的，無法根據(jù)循環(huán)索引進(jìn)行均等劃分；并且由于數(shù)組所占的空間大于500 KB，遠(yuǎn)大于SPM的64 KB空間，所以也無法完整傳輸.

3.2 循環(huán)并行化

在循環(huán)并行化方面，神威OpenACC與OpenACC 2.0標(biāo)準(zhǔn)保持一致，使用parallel和loop語句使代碼在64個(gè)從核上并行執(zhí)行.

需要注意的是，gang,worker,vector是OpenACC 2.0標(biāo)準(zhǔn)中的3層循環(huán)設(shè)計(jì)，由于神威眾核架構(gòu)在物理上并沒有分層需求，所以神威OpenACC的實(shí)現(xiàn)是默認(rèn)把gang設(shè)置為64，worker設(shè)為1，而vector也沒有對SIMD功能做支持.

3.3 原子操作

為保證數(shù)據(jù)競爭情況下的結(jié)果正確性，神威OpenACC與OpenACC 2.0一樣，都支持了原子操作.atomic指導(dǎo)語句可以保證數(shù)據(jù)在同一時(shí)間只會(huì)被某一個(gè)進(jìn)程讀寫.在GTC-P的函數(shù)charge中，數(shù)組densityi_part由于不規(guī)則訪問會(huì)產(chǎn)生數(shù)據(jù)競爭.為保證結(jié)果的正確性，需要添加atomic指導(dǎo)語句.

針對charge中的原子操作，本文利用“太湖之光”上提供的性能分析工具[12]測得執(zhí)行一次原子操作的開銷約為750個(gè)周期.同時(shí)，由于數(shù)據(jù)競爭問題的存在，導(dǎo)致原子操作的數(shù)組無法通過copy指導(dǎo)語句傳輸至64個(gè)從核的SPM中，使得計(jì)算部分必須進(jìn)行主存訪問，增加了額外的訪存開銷.

Table 2 Performance Comparison of Implementation on

表2為利用本節(jié)所提方法移植GTC-P后，主核與從核版本之間的執(zhí)行時(shí)間對比，從核版本比主核版本慢了將近92倍.造成性能瓶頸的原因有以下2點(diǎn)：

1) 神威眾核處理器上的原子操作由鎖機(jī)制實(shí)現(xiàn)，實(shí)際的運(yùn)行開銷大.

2) 從核之間缺乏共享局部存儲，為配合原子操作，無法將數(shù)據(jù)預(yù)存入各從核的SPM中.

以上原因與神威眾核的硬件架構(gòu)相關(guān)，目前僅依靠神威OpenACC指導(dǎo)語句是無法解決該性能瓶頸的.

4 性能優(yōu)化

針對3.3節(jié)提出的性能瓶頸原因，本文提出不依賴于神威OpenACC編譯器直接編譯得到的可執(zhí)行程序，而是基于中間代碼的二次開發(fā)來優(yōu)化性能.

4.1 消除原子操作

為了解決多從核之間的數(shù)據(jù)讀寫競爭問題，同時(shí)又避免原子操作，我們先對函數(shù)charge中相關(guān)的代碼段進(jìn)行分析，如圖4所示:

① #pragmaaccatomic② densityi_part[ij1]+=d1;③ #pragmaaccatomic④ densityi_part[ij1+1]+=d2;⑤ #pragmaaccatomic⑥ densityi_part[ij1+mzeta+1]+=d3;

Fig. 4 charge function atomic operation code
圖4 函數(shù)charge原子操作代碼段

這部分的數(shù)組操作為累加，在從核并行時(shí)并不涉及共享同步機(jī)制，因此本文采用的優(yōu)化策略是在每個(gè)從核的SPM中保存一個(gè)數(shù)組副本，在代碼段執(zhí)行完之后對數(shù)組副本進(jìn)行1次統(tǒng)一的歸約求和計(jì)算.

由于需要分別對64個(gè)從核SPM中的數(shù)組副本進(jìn)行維護(hù)，這部分優(yōu)化需要利用神威OpenACC編譯器在中間代碼中提供的acc_myid接口才能實(shí)現(xiàn).通過用線程號維護(hù)densityi_part[acc_myid][]來實(shí)現(xiàn)數(shù)組副本的管理，代碼轉(zhuǎn)換如圖5所示：

① densityi_part[acc_myid][ij1]=d1;② densityi_part[acc_myid][ij1+1]=d2;③ densityi_part[acc_myid][ij1+mzeta]=d3;④ densityi_part[acc_myid][ij1+mzeta+1]=d4;

Fig. 5 Array copies maintained by thread id on CPEs
圖5 數(shù)組根據(jù)從核進(jìn)程號維護(hù)副本

因?yàn)閿?shù)組densityi_part大于500 KB無法完整存放到SPM中，為了在不規(guī)則訪問的情況下盡可能避免主存訪問，本文在該代碼段前添加copyin語句，并結(jié)合變量索引進(jìn)行DMA數(shù)據(jù)傳輸，即在每次對于數(shù)組副本densityi_part進(jìn)行累加前，先將涉及到的數(shù)據(jù)傳輸?shù)絊PM中.優(yōu)化后的性能如表3所示：

Table 3 Results after Eliminating Atomic Operations表3 消除原子操作的優(yōu)化效果 s

在消除了原子操作和主存訪問帶來的額外開銷之后，函數(shù)charge實(shí)現(xiàn)了約180倍的加速，性能接近于主核版本，但仍慢了將近13，需要進(jìn)一步優(yōu)化.

4.2 提升DMA傳輸帶寬

如2.1節(jié)所述，在神威眾核架構(gòu)下，64個(gè)從核同時(shí)發(fā)起DMA訪問時(shí)，訪存帶寬可以達(dá)到30 GBps.然而神威OpenACC中的copycopyincopyout指導(dǎo)語句默認(rèn)卻是按照單次循環(huán)涉及的數(shù)組元素進(jìn)行傳輸?shù)?，若不進(jìn)行優(yōu)化，則無法提高DMA傳輸帶寬.

指導(dǎo)語句tile是OpenACC 2.0標(biāo)準(zhǔn)中的語句，在CPU和GPU上的實(shí)現(xiàn)是為了通過循環(huán)合并達(dá)到更好的數(shù)據(jù)局部性.神威OpenACC實(shí)現(xiàn)的tile指導(dǎo)語句則是通過設(shè)定循環(huán)合并的層數(shù)來控制每次DMA傳輸?shù)臄?shù)據(jù)塊大小，從而達(dá)到提升傳輸帶寬的目的.根據(jù)“pragma acc data copy(z0) tile(250)”的指導(dǎo)語句編譯出的中間代碼示例如圖6所示.其中acc_sync_pe_m2l_nostride_copy為實(shí)際的DMA數(shù)據(jù)傳輸接口，而acc_cp_4_size則為傳輸數(shù)據(jù)的尺寸參數(shù)，則此次傳輸?shù)臄?shù)據(jù)塊為250×8 B，即一次DMA傳輸約為2 KB的數(shù)據(jù).

① for(acc_blockindex_m0_1=acc_myid*250;acc_blockindex_m0_1<=_ldm_mi-1;acc_blockindex_m0_1+=250*acc_corenum)② {acc_loop_ub_tmp_var0=_ldm_mi-1-acc_blockindex_m0_1;③ acc_loop_ub_m=min(249,acc_loop_ub_tmp_var0);④ acc_cp_4_size=(acc_loop_ub_m+1)*8;⑤ _swacc_ret=acc_sync_pe_m2l_nostride_copy(_ldm_z0,&z0[acc_blockindex_m0_1],acc_cp_4_size). ?}

Fig. 6 Intermediate code generated by tile directive
圖6 tile指導(dǎo)語句生成的中間代碼

根據(jù)上述原理，本文針對charge和push中的tile尺寸設(shè)定做了測試實(shí)驗(yàn)，結(jié)果如圖7所示.伴隨tile尺寸增大、執(zhí)行時(shí)間變小，性能得到優(yōu)化，但受限于實(shí)際可用的SPM大約為58 KB，在這2個(gè)函數(shù)中可用的tile尺寸最大分別為251和336.

與消除原子操作的版本相比，在提升DMA傳輸帶寬后，函數(shù)charge和push分別實(shí)現(xiàn)了2.7及4倍的加速比.

Fig. 7 Test for tile size in Sunway OpenACC圖7 神威OpenACC的tile尺寸測試

4.3 避免低效的訪存操作

“太湖之光”上提供的性能分析工具中的penv2_slave_gld_count接口可以收集中間代碼里的gldgst指令數(shù)目，幫助用戶發(fā)現(xiàn)冗余指令.該工具提供了一系列的函數(shù)調(diào)用接口，可以幫助用戶收集硬件相關(guān)的統(tǒng)計(jì)數(shù)據(jù)，包括總指令數(shù)、loadstore指令數(shù)、TLB miss數(shù)等[12]，使用代碼如下所示：

penv2_slave_gld_init();

?

?

penv2_slave_gld_count(&ic1);

?

penv2_slave_gld_count(&ic2).

在函數(shù)push中發(fā)現(xiàn)的冗余gldgst指令，主要分2種情況：

1) 直接訪問保存在主存上的指針.對于這類情況的優(yōu)化方法是:在SPM上重新聲明一個(gè)局存指針，然后在代碼段初始部分將指針地址賦值給局存指針，使其后的指針操作都是直接訪SPM，而非通過gld指令的全局離散訪主存.

2) “太湖之光”編譯環(huán)境中默認(rèn)鏈接的數(shù)學(xué)庫需要進(jìn)行離散訪主存操作，而函數(shù)push需要使用sin,cos,exp這3個(gè)數(shù)學(xué)庫計(jì)算.本文調(diào)用了加載在SPM上的優(yōu)化從核數(shù)學(xué)庫[13].該數(shù)學(xué)庫完全在SPM上實(shí)現(xiàn)，所以不會(huì)產(chǎn)生gldgst指令，但會(huì)占用更多的SPM空間，需要提前規(guī)劃好SPM的使用.最新版本中雙精度的函數(shù)cos和sin 各需要4 KB空間，而函數(shù)exp 占用9 KB空間.

在對函數(shù)push進(jìn)行消除低效的訪存操作后，對比上一步的優(yōu)化結(jié)果，進(jìn)一步取得了6.7倍的加速比.

4.4 手動(dòng)添加SIMD intrinsics

正如2.1節(jié)中的理論峰值性能公式所示，256b的SIMD向量化(相當(dāng)于4個(gè)雙精度數(shù))操作對于性能起了重要作用.若所有計(jì)算都僅使用64b的雙精度標(biāo)量計(jì)算，實(shí)際上即是損失了75%的計(jì)算性能.神威OpenACC中尚未提供針對SIMD的支持，且自動(dòng)SIMD向量化一直以來是編譯器的實(shí)現(xiàn)難點(diǎn).所以本文采用在中間代碼中手動(dòng)添加SIMD intrinsics指令的方法，不依賴編譯器，手動(dòng)向量化計(jì)算部分的代碼.

本文針對函數(shù)charge中的部分計(jì)算進(jìn)行了向量化優(yōu)化，實(shí)現(xiàn)示例如圖8所示：

①for(_ldm_i=0;_ldm_i<=acc_loop_ub_i;_ldm_i+=1){② vtmp=simd_set_doublev4(0,0,0,0);③ for(_ldm_j=0;_ldm_j<=63;_ldm_j+=4){④simd_load(vd,&ldm_swap_mydensityi_local0[_ldm_i][_ldm_j]);⑤ vtmp=simd_vaddd(vd,vtmp);⑥ _ldm_denloc_sum+=_ldm_swap_mydensityi_local0[_ldm_i][_ldm_j];⑦}⑧simd_store(vtmp,&tmp[0]);⑨_ldm_densityi[_ldm_i]=tmp[0]+tmp[1]+tmp[2]+tmp[3];⑩}

Fig. 8 SIMD optimization code in charge function
圖8 函數(shù)charge向量化優(yōu)化代碼

通過在中間代碼中添加SIMD intrinsics指令、對比執(zhí)行周期數(shù)，該部分代碼實(shí)現(xiàn)了5.6倍的加速比，結(jié)果證明該優(yōu)化方法效果明顯.但由于該代碼段在全局中所占的時(shí)間尚不足1%，所以對于整體的GTC-P性能并未產(chǎn)生影響.

5 實(shí)驗(yàn)驗(yàn)證

本節(jié)在“太湖之光”上針對第4節(jié)提出的優(yōu)化方法進(jìn)行了實(shí)驗(yàn)驗(yàn)證.其中，手動(dòng)添加SIMD intrinsics指令的方法對整體性能影響不大，所以未加入本實(shí)驗(yàn).

5.1 性能對比分析

本文將運(yùn)行在1個(gè)主核上的串行版本GTC-P作為測試基準(zhǔn)，分別與消除原子操作、提升DMA傳輸帶寬、消除冗余gldgst指令這3個(gè)優(yōu)化方法在單核組的從核上進(jìn)行了性能對比，結(jié)果如圖9所示.該測試的問題規(guī)模為A，即回環(huán)中有3 235 896個(gè)粒子.

Fig. 9 Performance comparison of GTC-P step by step圖9 GTC-P性能對比

除了4.1節(jié)中提及的消除原子操作優(yōu)化效果顯著外，本節(jié)的實(shí)驗(yàn)結(jié)果還驗(yàn)證了2點(diǎn)：

1) 通過tile指導(dǎo)語句優(yōu)化DMA傳輸效率，對于神威眾核上訪存密集型的應(yīng)用，是有效且必要的.該優(yōu)化對比消除原子操作的版本，在函數(shù)charge和push上分別實(shí)現(xiàn)了2.7及4倍的加速比，并使得GTC-P在從核上的整體性能要高于主核上的.

在單核組實(shí)驗(yàn)中，函數(shù)charge和push在從核上的最終版本對比主核版本分別實(shí)現(xiàn)了約1.6倍和8.6倍的加速比，對于整個(gè)GTC-P應(yīng)用而言，實(shí)現(xiàn)了2.5倍的加速比.

5.2 強(qiáng)可擴(kuò)展性分析

本文在利用OpenACC實(shí)現(xiàn)的單核組版本基礎(chǔ)上，結(jié)合MPI[6]在“太湖之光”上測試了OpenACC+MPI版本，分析了強(qiáng)可擴(kuò)展性.問題規(guī)模較之前的性能對比，使用的輸入規(guī)模為B(1 235 644 416個(gè)粒子)，步長仍為100.其中，圖10(a)為輸入規(guī)模A的結(jié)果，圖10(b)為規(guī)模B的結(jié)果.

圖10中的橫軸代表MPI的線程數(shù)，1個(gè)MPI線程控制1塊SW26010芯片上的1個(gè)核組，即4個(gè)MPI線程控制整塊芯片(即單節(jié)點(diǎn)).實(shí)驗(yàn)的線程數(shù)從4一直擴(kuò)展到了128，即從1個(gè)節(jié)點(diǎn)擴(kuò)展到了32個(gè)節(jié)點(diǎn).

伴隨并行計(jì)算節(jié)點(diǎn)數(shù)的上升，函數(shù)push的強(qiáng)可擴(kuò)展性表現(xiàn)要優(yōu)于函數(shù)charge的，執(zhí)行時(shí)間基本呈理想的線性遞減趨勢.而GTC-P應(yīng)用的整體強(qiáng)可擴(kuò)展性則未與理想的線性趨勢相吻合.主要原因是由于GTC-P的整體執(zhí)行時(shí)間由函數(shù)charge主導(dǎo)(占80%左右)，而charge函數(shù)則受限于不規(guī)則訪存導(dǎo)致的性能瓶頸，尚需進(jìn)一步優(yōu)化.

6 總 結(jié)

本文利用OpenACC在“太湖之光”上并行移植了GTC-P應(yīng)用.在實(shí)現(xiàn)了完全基于指導(dǎo)語句的移植版本后，我們發(fā)現(xiàn)神威OpenACC編譯器尚無法解決由于原子操作開銷大以及離散訪存操作代價(jià)高所導(dǎo)致的性能瓶頸.本文提出，不依賴于編譯器，而是基于中間代碼進(jìn)行二次開發(fā)，通過消除原子操作和冗余gldgst指令、在中間代碼里手動(dòng)添加SIMD intrinsics指令的優(yōu)化方法，提升應(yīng)用性能.實(shí)驗(yàn)結(jié)果證明所提出的優(yōu)化方法有效.

此外，本文還通過與OpenACC 2.0標(biāo)準(zhǔn)進(jìn)行對比，提出了利用OpenACC移植應(yīng)用到“太湖之光”上需要注意的地方，總結(jié)為以下3點(diǎn)：

1) 在data指導(dǎo)語句使用過程中，需要注意數(shù)組索引與循環(huán)索引的關(guān)聯(lián)性，盡可能使數(shù)組可以按索引劃分并傳輸至從核的SPM中.若無法劃分且所占空間小于64 KB，則建議可利用annotate(entire(array))指導(dǎo)語句將數(shù)組完整傳輸至SPM.

3) tile指導(dǎo)語句對于應(yīng)用在“太湖之光”上的優(yōu)化效果顯著，GTC-P應(yīng)用中部分函數(shù)通過該指導(dǎo)語句優(yōu)化實(shí)現(xiàn)了3～4倍的加速比.

我們下一步將探索適應(yīng)神威眾核架構(gòu)，尤其是針對64 KB SPM局存設(shè)計(jì)的新算法，并利用底層編程方法對GTC-P進(jìn)行深度優(yōu)化.

致謝感謝國家并行計(jì)算機(jī)工程技術(shù)研究中心的劉鑫、尤洪濤、何香老師在本研究中提供的大力支持!

[1]Dongarra J. Report on the Sunway TaihuLight system[EBOL]. [2016-06-24]. http:www.netlib.orgutkpeopleJackDongarraPAPERSsunway-report-2016.pdf

[2]Zheng Fang, Zhang Kun, Wu Guiming, et al. Architecture techniques of many-core processor for energy-efficient in high performance computing[J]. Chinese Journal of Computers, 2014, 37(10): 2176-2186 (in Chinese)(鄭方, 張昆, 鄔貴明, 等. 面向高性能計(jì)算的眾核處理器結(jié)構(gòu)級高能效技術(shù)[J]. 計(jì)算機(jī)學(xué)報(bào)，2014，37(10): 2176-2186)

[3]Tang William, Wang Bei, Ethier S. Scientific discovery in fusion plasma turbulence simulations at extreme scale[J]. Computing in Science & Engineering, 2014, 16(5): 44-52

[4]Ethier S, Chang Choon-Seock, Ku Seung-Hoe, et al. NERSC’s impact on advances of global gyrokinetic PIC codes for fusion energy research[J]. Computing in Science & Engineering, 2015, 17(3): 10-21

[5]Fu Haohuan, Liao Junfeng, Yang Jinzhe, et al. The Sunway TaihuLight supercomputer: System and applications[J]. Science China Information Sciences, 2016, 59(7): 072001

[6]Wang Bei, Ethier S, Tang William, et al. Kinetic turbulence simulations at extreme scale on leadership-class systems[C]Proc of the Int Conf on High Performance Computing, Networking, Storage and Analysis. New York: ACM, 2013: 82:1-82:12

[7]Ibrahim K, Madduri K, Williams S, et al. Analysis and optimization of gyrokinetic toroidal simulations on homogenous and heterogenous platforms[J]. International Journal of High Performance Computing Applications, 2013, 27(4): 454-473

[8]Calore E, Kraus J, Schifano S, et al. Accelerating lattice boltzmann applications with OpenACC[C]Proc of the European Conf on Parallel Processing. Berlin: Springer, 2015: 613-624

[9]Hariri F, Tran T, Jocksch A, et al. A portable platform for accelerated PIC codes and its application to GPUs using OpenACC[J]. Computer Physics Communications, 2016, 207: 69-82

[10]Hoshino T, Maruyama N, Matsuoka S, et al. CUDA vs OpenACC: Performance case studies with kernel benchmarks and a memory-bound CFD application[C]Proc of the 13th Int Symp on Cluster, Cloud, and Grid Computing. Piscataway, NJ: IEEE, 2013: 136-143

[11]Madduri K, Ibrahim K, Williams S, et al. Gyrokinetic toroidal simulations on leading multi-and manycore HPC systems[C]Proc of Int Conf for High Performance Computing, Networking, Storage and Analysis. New York: ACM, 2011: 23:1-23:12

[12]He Xiang, Zhou Mingzhong, Liu Xin. Design and implementation of multi-level heterogeneous parallel algorithm of 3D acoustic wave equation forwarded[J]. Computer Applications and Software, 2014, 31(1): 264-267 (in Chinese)(何香, 周明忠, 劉鑫. 三維聲波方程正演多級異構(gòu)并行算法設(shè)計(jì)與實(shí)現(xiàn)[J]. 計(jì)算機(jī)應(yīng)用與軟件，2014，31(1): 264-267)

[13]Xu Jinchen, Guo Shaozhong, Huang Yongzhong, et al. Access optimization technique for mathematical library of slave processors on heterogeneous many-core architectures[J]. Computer Science, 2014, 41(6): 12-17 (in Chinese)(許瑾晨, 郭紹忠, 黃永忠, 等. 面向異構(gòu)眾核從核的數(shù)學(xué)函數(shù)庫訪存優(yōu)化方法[J]. 計(jì)算機(jī)科學(xué), 2014, 41(6): 12-17)WangYichao, born in 1990. Master. Member of CCF. His main research interests include high performance computing and GPU computing.

LinXinhua, born in 1979. PhD. Member of CCF. His main research interests include computer architecture and high performance computing.

CaiLinjin, born in 1994. Master candidate. His main research interests include high performance computing (caijinjin4@sjtu.edu.cn).

TangWilliam, born in 1944. Professor, APS fellow. His main research interests include plasma physics and scientific computing (tang@pppl.gov).

EthierStephane, born in 1966. PhD. His main research interests include computa-tional physics and plasma physics (ethier@pppl.gov).

WangBei, born in 1980. PhD. Her main research interests include plasma physics and high performance computing (beiwang@princeton.edu).

SeeSimon, born in 1965. PhD. Member of CCF. His main research interests include high performance computing and GPU architecture (ssee@nvidia.com).

SatoshiMatsuoka, born in 1963. Professor and ACM fellow. His main research interests include high performance computing and computer architecture.