“神威·太湖之光”上Tend_lin并行優化

， ， ，，,

(1.山東科技大學 計算機科學與工程學院，山東 青島 266590；2.山東省計算中心(國家超級計算濟南中心)，山東 濟南 2501013.中國科學院計算技術研究所 計算機體系結構國家重點實驗室，北京 100190)

地球系統模式已成為研究和預測全球變化的有力工具，也是當前大氣科學領域的研究熱點之一。中科院大氣物理所研發的大氣環流模式(atmospheric general circulation models，AGCM)，利用偏微分方程求解物理問題，是典型的高性能應用，是整個中科院地球系統模式CAS-ESM(earth system model)中最復雜的部分。其動力框架基于準一致網格[1]實現，在分辨率提高的同時，計算量增大，進程間通信量也隨之增加。如何提升較高分辨率下的AGCM模式性能，已成為重要的研究課題。

“神威·太湖之光”是由我國自主研發的高性能異構眾核處理器“申威(SW)26010”構建的超級計算機系統，其峰值性能超過100Pflops，曾4次蟬聯世界超級計算機Top500榜首，目前在世界超算平臺排名第3。申威眾核處理器支持眾核加速編程模型(神威OpenACC，簡稱SWACC)和加速線程庫(Athread庫)兩種并行編程模型。徐金秀等[2]在神威平臺上開展耦合算法的多級并行優化，并使用SWACC實現耦合算法的眾核并行；李億淵等[3]基于神威平臺，利用Athread在線程級層面，對稀疏矩陣向量乘法進行并行算法設計和優化實現。眾核異構平臺相比傳統的高性能硬件平臺能提供更強大的計算能力和訪存帶寬[4]，在國內當前主流的異構眾核高性能平臺上對AGCM進行眾核并行化，成為目前的重要研究方向之一。

本研究面向國產眾核異構平臺“神威·太湖之光”，以中科院大氣物理研究所研發的IAP AGCM4.0為基礎，對其熱點程序動力框架的適應過程構造Tend_lin應用，利用神威OpenACC編程模型實現Tend_lin應用的眾核并行化。通過循環分布、循環分塊并行將計算任務分配到計算核心，以高效DMA的方式將計算需要的數據批量傳輸到從核局存，消除計算過程中低效的全局離散訪存，實現了Tend_lin應用在國產異構眾核平臺上的OpenACC并行化。

圖1 申威眾核處理器結構圖

1 神威異構眾核處理器體系結構

“神威·太湖之光”整機由40 960塊SW26010處理器組成，并通過計算插件板、計算超節點和計算機倉等模式進行系統擴展。“申威(SW)26010”為主從核結構的異構眾核處理器，集成4個運算核組共260個運算核心，采用片上計算陣列集群和分布式共享存儲相結合的體系結構(圖1)。

每個CPU包括4個異構核組(core-groups，CGs)，每個異構核組由1個管理核心(management processing element，MPE，主核)、1個運算核心簇(computing processing elements clusters，CPE Cluster)、1個存儲控制器(memory controller，MC)和系統接口(system interface，SI)組成[5]。

主核在大規模計算過程中通常用于完成管理計算核心以及與其他核組進行通信的任務[7]。CPE Cluster包含64個運算核心(computing processing elements，CPE，從核)，以8×8的Mesh布局通過簇通信網絡進行連接。每個從核具備私有的L1指令Cache，一個運算核心簇共享L2指令Cache，從核只能運行于用戶模式并且不支持中斷，適用于處理邏輯簡單而計算相對密集的過程[8]。核組內共享內存。主/從核通過MC完成與主存的數據傳輸。每個從核均包括一個以SPM(scratch pad memory)方式組織的64 kB大小的局部數據存儲器(local data memory，LDM)，該空間需要用戶顯式的設計規劃使用。

SW26010處理器從核LDM空間有限，往往無法滿足使用大數組的科學計算程序存儲要求，大部分數據存儲在主存中。在計算過程中，必須通過一定的方式完成主存到從核LDM的數據傳輸。從核支持兩種主存訪問模式：1)全局讀入/寫出(gld/gst)，直接完成主存到從核寄存器間的零散數據訪問；2)直接內存訪問(direct memory access，DMA)，批量完成從內存到LDM的大量數據拷貝，通常延時較大，但帶寬利用率較高。相比其他異構平臺，從核由于存儲空間和帶寬限制，數據傳輸往往成為程序性能瓶頸。Xu等[9]在主存數據128字節對齊的前提下，在一個核組內64個從核中進行Stream Triad，PE模式下測得DMA的帶寬可以達22.6 GB/s，而gld/gst方式每個核組帶寬利用只有1.5 GB/s左右；延遲方面，從核gld單核延遲需要177個cycles，而從核訪問局存的1d只需要4個cycles。可見，DMA方式比gld/gst能夠更有效的利用帶寬，但啟動一次DMA需要300 cycles左右，開銷較大，若DMA次數過多會大大降低程序性能。因此在SW26010異構眾核處理器上，獲得更好的數據傳輸性能的關鍵在于一次完整數據傳輸的數據量要盡可能大，并減少DMA的啟動次數。

圖2 單核組訪存周期

2 AGCM及其適應過程Tend_lin應用的構成及特點分析

2.1 AGCM動力框架

CAS-ESM是以耦合器為核心的模塊化軟件，集成了大氣、陸面、陸冰、海冰和海洋等分量模式的復雜系統，AGCM是其中最復雜的分量模式。

IAP AGCM主要由物理過程和動力框架兩個過程組成。其中，物理過程計算量占比小，并行可擴展性較高；動力框架主要是用來求解關于時間的偏微分方程組，其計算在空間上是三維的，數值方法復雜，是主要的優化部分。IAP AGCM 4.0[12]的動力框架采用均勻經緯網格，并在高緯度采用靈活跳點格式[13]。格點模式中極區計算難以處理，模式引入濾波模塊，將緯度劃分為3個緯度帶，不同緯度采用不同濾波方法以增加時間步長，避免計算不穩定。本研究提取IAP AGCM 4.0的熱點動力框架中的適應過程Tend_lin作為研究對象。

2.2 Tend_lin構成及特點分析

圖3 Tend_lin應用中fft計算函數調用層次圖

Tend_lin的計算過程主要包括兩個不同的階段，兩個階段的計算數據來自不同的數組，但都主要由stencil和1維fft兩類計算組成，主要包括積分、濾波和平滑等計算類型；通信模式主要是近鄰通信和集合通信。其中stencil計算涉及22組計算循環，占應用串行運行時間的71%，計算中用到臨近若干點，每個計算循環的總訪問足跡較大，其中14個緯度-高度-經度三維空間全遍歷，計算訪問3維數組較多。濾波和平滑功能的計算分別在filt和smoother中，不同緯度選擇了不同的濾波方法，函數調用層次深(圖3)。該階段計算循環面對的是同一個高度上的數組區域，包含緯度-經度上遍歷的兩層循環，每個高度上要進行多階段處理，不同高度間無數據依賴關系。雖然數據訪問量小，但計算循環個數較多。

Tend_lin應用涉及73個計算循環和單入單出的代碼塊，數據訪問多以數組形式在函數內動態分配內存空間。Tend_lin程序具有以下特點：①沒有熱點循環，任務并行代碼散落于73個并行循環和單入單出的代碼塊上；②氣候模式算法現狀導致網格點分辨率不高，每個循環的并行度都不是很大，充分利用眾核資源難度大；③循環間大小不一，并行度差異大；④計算/訪存比值較低，訪存優化為應用性能提升的關鍵；⑤各進程的計算行為不同構，在緯度方向上，不同進程采用不同濾波方法，且南北極需要特別處理。

為方便并行化，除適應過程的主函數Tend_lin之外，本研究增加了緯度和高度方向的區域分解，以及變量初始化部分，并為應用配備3個不同的問題規模A、B、C以匹配低、中、高不同分辨率，建立了包含59個源文件、2萬行左右的benchmark程序。

下面，將AGCM中的Tend_lin串行代碼移植到SW26010的主核上，并選取網格中典型分辨率即經度×緯度×高度為1.4°×1.4°×30L(下稱規模B)闡述并行化過程；將B規模在經緯度的分辨率擴大2倍，得到計算規模是B規模4倍的C規模，用來進行全程序評估。B規模下，單個緊嵌計算循環最大訪問足跡達62.8 M(stencil_11)；2個計算循環為緯度-經度的二維計算，保持了緯度上的并行度，但缺少高度維，計算量較小；4個單入單出的代碼塊，包含經度-高度方向的規約操作；2個經度方向的單層循環計算。

3 Tend_lin應用的神威OpenACC眾核并行化和優化

眾核并行化主要目標就是將最耗時的循環迭代分散到多個從核線程上并行執行。下面主要從循環分布、循環分塊及并行化、數據傳輸的優化以及函數調用的從核化4個方面描述Tend_lin應用在神威平臺上的眾核并行化方案。

3.1 循環分布

循環分布是在保證應用程序正確性的前提下，將源循環分解為多個獨立的目標循環。Tend_lin應用中，并行循環均為典型的DOALL循環，DO循環中所有的循環實例都可以并行的執行，每個循環實例中的所有語句均串行執行。應用中包含了兩類需要特殊處理的并行循環：一類是不同分支的數組讀區域，難以被SWACC編譯器合并，導致編譯錯誤；另一類循環體中包含了多個訪問區間不同的子循環且數據足跡較大，其循環體內訪問了大量數據重用較低的數組，計算訪存比過低。由于LDM尺寸有限，若直接并行該類循環，由于分塊尺寸較小，數據重用低，會出現較大的數據傳輸開銷，需要在并行化之前進行循環分布。

為避免出現激進的循環分布帶來的開銷可能會大于并行收益的情況，在對以上兩種并行循環進行循環分布的同時要進行空間優化。

1) 分支結構的循環分布

stencil計算中存在帶有分支結構的循環，不同分支間數組的J維訪問區域沖突，如GHI(IB:IE, beglev：endlev, JB)和GHI(IB:IE, beglev：endlev, loc_JB:loc_JE)，將其抽象為如圖4左側代碼。由于并不知道JB和loc_JB:loc_JE的大小關系，編譯器無法給出統一規則完成合并，進而不能完成編譯。

圖4 分支結構分布示意

經過對該循環的不同分支分布為圖4右側所示循環代碼，程序能夠正常編譯且數據傳輸總量并沒有變化。

2) 大數據足跡循環的循環分布

Tend_lin應用的02_stencil階段存在非緊嵌循環(如圖5左側)，并且：①最外層J循環內，有多個非緊嵌K-I子循環，這些子循環訪問變量眾多，對從核空間需求較大(約47 360 kB)，即便最外層J維分塊大小為1(約370 kB)，仍遠超從核LDM空間(64 kB)；②不同的K-I子循環間訪問數據足跡大小不一，不同循環間訪問的數組不同，且循環間訪問的數據無任何依賴或者重用；③K維訪問區間各異，難以合并為緊嵌循環，且K維空間過小，不能提供有效的并行度，不適合K維并行。

對內層循環進行循環分布時，必須兼顧合法性和收益性問題。分布合法性指循環分布的原有數據依賴語義必須保持，既要保證源循環套并列循環間沒有跨迭代依賴，又要保證兩個語句組之間沒有反向的數據依賴。而為保證循環分布后的收益，有數據重用的盡量放到一個循環套內[14]，減少因拆分導致的同一數組的多次讀取，避免額外開銷。

本研究在保證J維并行的前提下進行循環分布，將一個執行次數較多的非緊嵌循環分割成若干個執行次數較少、訪問足跡較小的完美緊嵌子循環(如圖5右側)，以求每個小循環并行時，獲得更大的循環分塊尺寸，提高DMA效率。

圖5 大數據足跡循環分布示意

3.2 循環分塊及并行化

神威OpenACC按照OpenACC 2.0[11]標準進行循環并行化，使用parallel和loop語句指示加速區代碼，加載到從核并行執行。但與標準OpenACC不同的是，由于硬件結構的區別，神威上gang、worker、vector的三層循環設計沒有分層的需求，默認gang設置成64，worker為1，vector也并未做SIMD功能支持。將硬件特征與程序特征進行匹配的循環分塊可以充分的開發硬件架構的潛能，將串行程序生成適合粒度的并行分塊程序。

1) 線程并行維度的選擇

Tend_lin應用中進程間保留了IAPAGCM4.0中的緯度-高度二維的進程劃分，本研究在線程間使用神威OpenACC進行眾核并行化，支持一維的并行。應用中的并行DOALL循環除stencil_09外，數組索引與循環索引緊耦合，可直接選取循環索引最外層，即選擇數組的最高維為并行維，以有效避免低效的低維跨步傳輸，提升DMA傳輸性能。

應用中并行循環之間訪問足跡大小差距較大，導致對LDM空間的需求差異大，若使用全局一致的循環分塊策略，則導致大部分的并行區LDM空間利用率不高，且OpenACC作為同步式的并行模型不需要全局一致的分塊，故本文采取局部最優劃分策略，保證每個循環執行時間最優。由于LDM空間有限，對于數據傳輸量較大的循環，即便并行維最小分塊，也無法滿足LDM空間需求。為了避免低效的跨步傳輸，本文選取最高維(緯度)和次高維(高度)使用循環分塊子句tile進行二維劃分，即緯度維循環并行分塊，高度維做DMA流水分塊。對于只劃分并行維就可以滿足空間需求的并行循環，則不再做次高維的流水劃分。對于非耦合循環stencil_09，由于存在最內層K維的規約操作，不宜調整循環順序，僅選取次高維(經度)流水分塊。

2) 循環分塊尺寸的選擇

為充分利用從核LDM空間，并保證最少的DMA通信次數，編譯器在對循環進行分塊時，既要保證內循環塊中訪問的數據總量不超過LDM的大小，又要保證使用盡可能大的LDM空間。根據文獻[14]中的面向異構多核處理器的循環分塊方案，以訪問足跡最大的緊嵌并行循環stencil_11為例，求解循環分塊的可行解。分析過程如下：

step 1：計算并行循環的總空間需求

stencil_11使用的數組數據類型為雙精度浮點類型，每個元素8個字節，計算該循環用到8個三維數組、3個二維數組和4個一維數組，另外還訪問1個緯度方向跨度為1的3維數組和1個緯度方向跨度為2的2維數組作為陰影區，則共需要的數組尺寸為((256×128×30×8+256×128×3+256×4+256×1×30×1+256×2×1)×8)=63 774 720字節=62 280 kB，遠超LDM空間64 kB，須進行循環分塊后才能夠放入局存當中。

step 2：確定最高維緯度j層分塊大小

令該層的分塊大小為1，這時循環訪問的數組所占空間即數據足跡為558 kB，仍超過LDM空間，說明即便使分塊大小為1，僅進行最外層(最高維)的循環分塊，仍不能滿足LDM空間需求的，需要對下一維進行劃分。

step 3：確定次高維高度k層分塊大小

當最外層大小為1時，數據足跡是LDM空間的8.7倍。由于循環訪問的數組形狀差異大，數據足跡和k層分塊之間并不是嚴格成比例。當k層循環分塊大小取3(30/8.7=3)時，數據足跡大小72 kB，超過LDM空間；當k層循環次外層分塊大小取2時，數據足跡為54 kB，為本分塊策略下最大可用空間，并為最大次外層分塊尺寸。

通過以上3步，得到循環分塊維度劃分方案：最外層(并行維)按分塊大小1劃分，次外層按照分塊大小2或1進行劃分，最內層不劃分。

3.3 數據傳輸優化

在第1部分介紹的從核線程的兩種訪存方式中，直接訪存比全局離散訪存能更好的利用帶寬，因此從應用特征出發，最大可能地將計算所需的數據以DMA方式批量拷貝到LDM存儲，消除計算過程中的離散全局訪存是性能優化步驟的關鍵。

本研究采用將計算中所需的數據都存儲在從核LDM中的數據管理策略，Tend_lin應用計算過程中訪問的數組多為動態數組，編譯器無法直接獲取其區域信息；循環間訪問行為各異，需要對數據傳輸進行差異化表達。經過3.1節循環分布，實現了循環中訪問數組索引變量與循環的索引變量緊耦合，這些循環直接使用tile子句由編譯器自動完成數組劃分，搭配copy/copyin/copyout語句，在每個并行循環分塊計算前，以DMA傳輸方式將所需數據傳輸到各個從核的LDM中，實現數據的提前拷貝。具體進行數據傳輸時需要解決以下幾個問題：

1) 動態數組區域信息的表達

應用中的大量動態數組，在編譯時編譯器不能獲得其維度信息，其形狀不能被分析出來，無法在從核LDM中申請對應循環分塊尺寸的空間，需要使用暗示制導語句annotate(dimension(array))來指明數組維度，為編譯器劃分數組提供依據。如圖6所示，并行循環中DGH和GHI1數組為動態數組，利用copy/copyin子句搭配dimension暗示制導語句的形式實現動態數組的DMA傳輸。

圖6 Dimension使用示意

2) 最低維的傳輸區域表達

Tend_lin應用的并行循環中，存在對非并行維非全維度的訪問，如圖7所示，J維為并行維，VT為動態分配內存空間的三維數組，在計算中I維只訪問了區間[IB:IE]，占全維度的98%，并沒有訪問全維度。

系統提供的parallel copy子句僅支持標量和數組，不支持子數組形式，因此copy子句只能傳輸全維度數據，而不是該維度的子區間(IB:IE)。若要實現非全維度的準確低維區間傳輸，只能利用數組聲明進行手工修改。但由于精確傳輸的dma_get為跨步傳輸，單從核帶寬為6.33 GB/s，64核帶寬為30.48 GB/s。而直接使用copy子句進行I維全維度傳輸為非跨步傳輸，雖然數據量增加了2%，但單從核帶寬可達8.14 GB/s，64核帶寬為30.18 GB/s。導致兩個方案在64線程時的DMA差別不到1%，即全維度傳輸雖然多使用了2%的存儲空間，但未明顯影響傳輸時間。因此最終選用copy子句實現I維全維度傳輸。

3) 循環不變數組的拷貝

圖8中，相對于并行維J，只讀數組A是循環不變數組，在并行循環中的每個實例是相同的，存在讀讀重用。循環中數組A的訪問區域只與流水循環索引有關，編譯器依照流水分塊尺寸在流水循環內生成的多次DMA傳輸。這時利用annotate(entire(array))將整個數組在流水循環外以一次DMA的方式傳輸到每個從核LDM中，使每個從核的LDM中都保有該數組的一個完整副本。但由于LDM空間有限，該方式只適用于傳輸較小的數組。

圖7 非并行維非全維度訪問示意

圖8 循環不變數組用法示意

4) 數組的轉置后傳輸

如圖9中的stencil_09循環，并行循環維度從高到低順序為J、I、K，I循環內存在子循環K，對應TT1數組的最高維，循環內對TT1數組的訪問是跨步訪問，需要利用轉置子句swapin對原始數據進行重新布局，改善數據傳輸效率。如圖9所示，在數據傳輸前，主核使用swapin語句數據進行轉置，使循環索引與數組索引的高低維順序對應，以轉置后的數組代替原數組訪問，再按照循環分塊進行傳輸。

圖9 Swapin子句使用示意

對于該計算循環，轉置前由于不能使用copyin將TT1拷入LDM，計算過程只能使用gld訪存，計算時間占總執行時間的近80%；使用swapin轉置后，該數組訪問由低效的gld變為高效的ld操作，計算時間縮短了13倍，這樣雖然轉置時間占到了轉置后計算循環總執行時間的30%，但計算時間的大幅縮短使整個循環執行時間減少2倍以上，該循環在整個程序中占計算比重20%，為tend_lin全程序帶來5%左右的影響。

5) 標量的傳輸

每個并行循環中訪問的標量最多有12個，平均5個左右，按照讀寫類型分為只寫和只讀兩類。只寫標量都為線程私有變量，可直接使用local子句本地化。

只讀標量多帶有parameter屬性，即便不在制導語句中有任何體現，編譯器仍然可以根據其parameter屬性自動添加到copyin和annotate readonly子句中，在從核LDM中生成對等的parameter變量，直接賦值。對于不具備parameter屬性的只讀循環邊界量，編譯器也可自動添加到copyin和annotate readonly子句中，在每個從核LDM中生成一份拷貝，每個變量以一次gld的形式load到LDM中；對于不具備parameter屬性的只讀循環內部變量，編譯器無法分析出該變量是否需要私有化，必須手動添加到copyin子句中才能在LDM中生成拷貝并以gld形式拷入。經過以上處理，將從核計算過程中所需要的數據在計算前存入從核LDM中。

為避免計算過程中不存在的低效gld/gst操作，利用“神威·太湖之光”上的性能分析工具penv_slave2_ gld_count和penv_slave2_gst_count，在編譯器生成的中間代碼中統計計算過程中從核gld/gst請求次數，確保完全消除了計算過程中的全局離散訪存。

3.4 函數調用的從核并行

應用中，加速區存在大量routine函數調用，且關系復雜，若直接在routine函數調用層的循環外使用制導語句，循環內的routine函數調用無法直接生成從核版本，導致生成的從核代碼找不到對應被調routine函數。

對于加速區的函數調用，在被調用函數內使用函數指示routine處理，為加速代碼區中被調用的函數自動創建可在從核上執行的版本。

4 并行化方法的測評

本研究以Tend_lin應用在SW26010處理器1個主核上，使用-O3優化選項的串行版本為測試基準，并以主核串行版本輸出為標準，相對誤差控制在10-10以下。針對第3節提出的并行化方案，分別利用B、C規模進行實驗性能評估。

4.1 不同分塊大小對程序性能影響

根據3.2節，在B規模下，以最高維緯度為并行維，分塊大小為1，次高維高度為流水維，最大訪問足跡的并行循環流水分塊分別使用1(36 kB)或者2(54 kB)對比測試；在最大訪問足跡循環流水分塊為2的基礎上，為每個并行循環選取空間需求不超過LDM空間的最大流水分塊尺寸，記為n。依此共設置三組對照實驗，分別取分塊大小為1×1、1×2和1×n，多線程相比主核單線程獲得的加速比結果如圖10所示。

圖10 神威OpenACC版多線程不同分塊尺寸性能對比

由圖10(a)可看出，隨著線程增加，純計算階段按比例加速，具有強擴展性，應用性能提升的關鍵在于優化訪存性能。由圖10(b)可知，隨著流水分塊的增大，應用性能不斷提升，說明在滿足LDM空間限制的前提下，DMA傳輸總量不變，隨著流水分塊尺寸的增大，DMA傳輸啟動次數減少，降低了DMA開銷；同時隨著流水分塊尺寸的增加，每次dma_get傳輸的向量長度增加，帶寬利用率增大。以上兩方面訪存性能的提升使全程序性能隨之提升。

圖11 Tend_lin應用不同規模多進程下加速比

4.2 全程序評估

在B、C兩規模下使用二維進程布局，每個進程選取最優線程數下表現的性能，得到Tend_lin應用最終獲得的全程序性能(如圖11)，在B規模下最終獲得相比主核串行25倍的加速；隨著分辨率的進一步擴大，在C規模下64進程時最高獲得69倍加速。

5 總結

本研究基于國產異構眾核平臺“神威·太湖之光”，使用神威OpenACC并行編程模型對Tend_lin應用進行了并行化，分別從循環分布、循環分塊、數據傳輸表達和函數調用從核化4個方面說明并行化方案。為保證高效的數據傳輸效率，從核并行選擇數組的最高維為并行維，對于最高維并行不能滿足LDM空間需求的循環，在次高維進行流水劃分，并為每個并行循環選取最大的流水分塊尺寸；從核并行化使用高效的DMA完成主存與LDM之間的數據傳輸，保證計算過程中只訪問從核LDM空間，并重點對典型數據傳輸的表達進行討論。最終應用在國產異構眾核平臺上，B規模單核組多線程獲得6.8倍的全程序加速，多進程獲得25倍加速；隨著分辨率擴大，加速更加明顯，C規模多進程下最多獲得69倍加速。