劉芳芳 ,王志軍 ,汪 荃 ,吳麗鑫 ,馬文靜 ,楊 超 ,孫家昶

1(中國科學院 軟件研究所 并行軟件與計算科學實驗室,北京 100190)

2(中國科學院大學,北京 100049)

3(計算機科學國家重點實驗室(中國科學院 軟件研究所),北京 100190)

4(北京大學 數學科學學院,北京 100871)

HPCG(high performance conjugate gradients)基準測試程序是一種新的超級計算機排名度量標準[1].該測試基準主要用于衡量超級計算機解決大規(guī)模稀疏線性系統(tǒng)的能力,相比于HPL(high performance Linpack)[2],其計算、訪存與通信模式與當前主流高性能計算機應用程序更為契合,能夠更全面地反映系統(tǒng)的訪存帶寬和延遲的性能.HPL 和HPCG 排名每年發(fā)布兩次,一次在德國的ISC[3]大會上,另一次在美國的SC[4]大會上.位于世界前列的超級計算機均提交HPL 和HPCG 兩個基準測試程序的結果,2019 年德國舉辦的ISC 大會同時發(fā)布了HPL和HPCG 結果.越來越多的應用研究人員開始關注HPCG 的測試結果,期望從中可以了解超級計算機訪存、通信等的特點,一些HPCG 并行和優(yōu)化方法也可能對某些應用有所啟示.

目前超級計算機的體系結構日趨復雜,早已從同構轉向異構,且同節(jié)點異構部件的運算性能差異越來越大,這對異構并行算法的設計將有很大的影響.本文主要面向一類國產復雜異構系統(tǒng)研究HPCG 并行算法和高效實現方法.目前Dongarra 教授課題組研究了HPL-AI 混合精度算法[5,6],我們也正在關注該方向的研究.

本文工作是在國產超級計算機研制方大力支持下完成的,在此表示感謝.

1 HPCG 簡介

HPCG 來源于半結構網格上的三維熱傳導應用,其核心是在三維規(guī)則區(qū)域上對Poisson 方程進行有限差分離散化后,轉換成一個稀疏線性方程組的求解問題.其方程如式(1)所示.采用如圖1 所示的27 點stencil 格式進行離散化.

其中,每個網格點的更新都依賴其周圍緊鄰的點.由于參與計算的點只能在三維網格區(qū)域范圍內,所以每個網格點最多依賴26 個鄰點,具體的鄰點的個數,分別為26(內部點,鄰點都在域內)、17(邊界面上的點,某一鄰面上的9 個點都在域外)、11(邊界線上的點,比邊界面上的點又少了2/3 個面的點,即又少了6 個點)和7(邊界頂點,僅有一個小立方體上的8個點在域內,即圖1的1/8,包括其自身).最終得到的稀疏矩陣A是由27 條對角線組成,并且是一個對稱正定矩陣.

在大規(guī)模并行環(huán)境中,HPCG 使用三維區(qū)域分解策略,也就是按照3 個維度將整個計算區(qū)域劃分成若干個子區(qū)域,每個子區(qū)域被分配一個MPI 進程,其規(guī)模由輸入參數指定.每個MPI 進程處理的計算區(qū)域中的內部點所依賴的26 個鄰點,全部來自于該計算區(qū)域,而對于計算區(qū)域邊界上的點,其依賴的鄰點一部分來自于該計算區(qū)域內部,一部分則來自于鄰居進程所處理的子區(qū)域.不妨稱計算區(qū)域為內區(qū),其對應的點為內點,計算區(qū)域所依賴的鄰居進程的點構成的區(qū)域為外區(qū),其對應的點為外點.由于我們使用的是27 點stencil,所以每個進程最多有26 個鄰居進程.

Fig.1 3D 27-point stencil in HPCG圖1 HPCG 三維27 點stencil 格式

在HPCG 中,該線性系統(tǒng)使用預處理共軛梯度(preconditioned conjugate gradient,簡稱PCG)法來求解,計算流程如圖2 所示.

CG算法中一次迭代計算包含3 個向量點積操作(第4、7、10 行)和3 個向量更新操作(WAXPBY,第6、8、9 行)以及1 個稀疏矩陣向量乘操作(SpMV,第7 行).

HPCG 中使用的預條件子M–1是基于V-cycle 幾何多重網格,其計算流程如圖3 所示.多重網格通過將細網格上頻率變化較緩慢的低頻分量映射到粗網格上來處理,以達到加快收斂速度的目的.作為預條件子,多重網格的計算同樣包含4 個主要操作:前磨光操作(第4 行)和后磨光操作(第8 行)、限制算子(第5 行)、插值算子(第7行)、底部解法器(第2 行).在HPCG 中,多重網格的層數被固定為4 層,前后磨光操作以及底部解法器都是使用對稱Gauss-Seidel 迭代法的一次迭代來實現,并且磨光操作只執(zhí)行一次,不可修改.使用SymGS(x,r)表示一次對稱Gauss-Seidel 迭代,其中r是右端項,x是解向量的初始猜測值.

Fig.2 Conjugate gradient algorithm圖2 共軛梯度法

Fig.3 Multigrid pre-conditioner圖3 多重網格預條件子

在HPCG 中主要存在兩種通信類型:全局通信和鄰居通信.

1.全局通信.主要為MPI_Allreduce 通信,發(fā)生在向量點積操作中.在一次CG 迭代中總共需要計算3 次向量的點積,因此需要調用3 次MPI_Allreduce.

2.鄰居通信.主要發(fā)生在與稀疏矩陣相關的兩個核心函數中:SpMV 和SymGS,使用MPI 點到點的通信實現halo 區(qū)數據交換.HPCG 采用27 點stencil 格式,每個網格點的計算依賴周圍緊鄰的26 個鄰居.在SpMV 和SymGS 計算開始前,需要與周圍的鄰居進程交換邊界數據(halo),以便將計算過程中所需要訪問到的鄰居數據傳輸到本地進程中.

3.目前超級計算機的體系結構日趨復雜,早已從同構轉向異構,且同節(jié)點異構部件的運算性能比差異越來越大,這對異構并行算法的設計將有很大的影響.本文主要面向一類國產復雜異構超級計算機研究HPCG 異構眾核并行算法.首先介紹了國內外進展,然后提出兩種適用于該平臺的HPCG 異構眾核并行算法.對于實際使用的第2 種算法中,還存在圖著色的問題,本文也開展了相關工作,并單獨使用一節(jié)進行介紹.接下來介紹在該超級計算機上的實現問題,包括CPU 和GPU 之間的任務劃分、為了隱藏鄰居通信實現的內外區(qū)劃分方法,以及一些稀疏矩陣重排等優(yōu)化方法,最后介紹了實驗結果.

2 國內外進展介紹

HPCG 的優(yōu)化工作主要是針對熱點函數SpMV 和SymGS 來開展.SpMV 的性能與存儲格式以及計算方式有關.Bell 等人[7]在GPU 平臺上實現了常用的CSR、COO、ELL 等基本格式,提出了一種新的HYB 格式來提高SpMV 的性能,并分析了不同矩陣適合使用的格式以及并行方法.另外,還有一些其他格式包括ELLPACKR[8]、BRC[9]、BCCOO[10]等被提出.還有一些學者研究稀疏矩陣的重排技術[11]及壓縮格式[12],以減少訪存開銷.在計算方式方面,Williams 等人[13]在多種不同架構的多核平臺上使用了線程分塊、緩存分塊、向量化等優(yōu)化手段,并取得了很好的結果.另外還有學者研究自動調優(yōu)技術[14],通過分析稀疏矩陣的特征來選擇參數以提升性能.

SymGS 的求解過程類似于稀疏三角解法器,其關鍵是如何在這種強依賴的限制下獲得并行度.研究人員已經提出了很多不同的并行方法,如Park 等人[15]提出了level scheduling 方法,這種方法使用層次遍歷稀疏矩陣對應的有向無環(huán)圖,保證每一層的節(jié)點間沒有依賴關系,從而可以并行更新.但是這種方法獲得的并行度不高,且各層之間的負載不均衡.Iwashita 等人[16]提出了適合結構化稀疏三角求解器問題的分塊圖著色技術,這種方法以塊為單位對圖進行著色,保證一個塊與其相鄰塊顏色不同,相同顏色的塊一起執(zhí)行.這種方法雖然會破壞部分依賴,但是可以獲得較大的并行度.

HPCG 的整體優(yōu)化工作也受到了很多研究者們的關注.Kumahata 等人[17]在基于CPU 的日本超級計算機K上的HPCG 優(yōu)化工作使用了分塊圖著色排序來挖掘SymGS 中的并行度,并通過改變反向更新時的計算順序來提高cache 命中率,在82 944 個CPU 上取得了0.461 PFLOPS 的成績,達到整機峰值性能的4.34%.Phillips 等人[18]在NVIDIA GPU 平臺上使用圖著色方法對SymGS 進行并行.劉益群等人[19,20]基于天河2 號超級計算機的特征提出了混合的CPU-MIC 協作的異構計算方法,通過區(qū)域劃分充分利用CPU 與MIC 的計算資源,在整機16 000 個節(jié)點上取得了0.623 PFLOPS 的性能,位列2014 年11 月排行榜第一名.敖玉龍等人[21]設計了基于申威架構的并行算法,該算法使用分塊著色的方式挖掘SymGS 中的并行度,并使用神威特色的片上寄存器通信機制優(yōu)化了核間的數據交換,在163 840 個節(jié)點上取得了0.481 PFLOPS 的性能.Ruiz 等人[22]基于ARM 平臺研究了多色重排和分塊多色重排兩個算法的性能,并從cache miss 等角度進行了分析.目前世界排名第一的是IBM 研制的超級計算機Summit,其HPCG 性能達到了2.92 PFLOPS,是整機峰值性能的1.5%,但目前還沒有相關文獻詳細介紹其并行算法和優(yōu)化技術.面向一類國產復雜異構系統(tǒng)研制眾核并行版HPCG 軟件,且峰值性能和整機效率達到或超過Summit 的水平,是一項具有巨大挑戰(zhàn)的工作.

3 異構眾核并行算法設計

GPU 是最近十幾年被廣泛使用的高性能計算機加速部件.一個GPU 由多個計算單元(computing unit,簡稱CU)組成,每個CU 里有多個SIMD,可以同時進行多路向量運算.CU 內部具有高速的軟件可控緩存LDS(local data share).所有CU 共享L2 cache 和設備內存.程序運行時組織成多個workgroup,每個workgroup 可包含一個或若干個wavefront.每個workgroup 只能在一個CU 上運行,在我們使用的平臺上,每個wavefront 包含64 個線程,一個workgroup 內所包含的線程數不能超過1 024 個.在這種高并行度的復雜架構上并行和優(yōu)化HPCG,是一個非常大的挑戰(zhàn).在并行算法層面,不僅要求高并行度,還要求整體迭代次數不能有較大增長.在優(yōu)化方面,不僅要從算法層面盡量減少訪存和通信,并將通信與計算重疊,還需要充分利用硬件的特性減少訪存開銷.

3.1 分塊著色+level scheduling算法

3.1.1 基本算法

如第2 節(jié)所述,HPCG 中SymGS 模塊的經典優(yōu)化方法包括點著色方法、按層計算方法(level scheduling)以及分塊著色算法等.考慮計算平臺的特點及收斂速度的需求,我們選擇了分塊著色算法作為程序的基本框架.我們將SymGS 和SpMV 等基本算子在分塊著色算法的基礎上進行并行化,并提出了針對GPU 異構平臺的優(yōu)化方法.

分塊著色算法的基礎是將三維網格分塊,然后對每一個數據塊進行著色.在預條件子中使用的四重網格里,每一層網格都需要分塊及著色.參考文獻[21]的工作,我們選取了類似的分塊著色算法,既方便數據在GPU 上的分配,又保證了較高的并行度,且不會大幅度增加迭代次數.如圖4 所示,三維空間中相鄰的8 個數據塊使用8 種不同的顏色,每次SymGS 計算時,同樣顏色的數據塊被同時處理.

Fig.4 Block coloring algorithm圖4 分塊著色算法

3.1.2 數據及任務映射

塊內SymGS 需要嚴格保持依賴關系,否則迭代次數將大幅度上升.因此,本文采用level scheduling算法進行并行,方程4x+2y+z=k所定義的平面上的點都可以算作一層(每個點的坐標為(x,y,z)),因為它們所依賴的數據都已經計算完畢.如圖5 所示,從上到下是時間軸,第1 個點的計算不需要本次迭代的任何數據,可以在此數據塊處理過程的第1 次迭代進行更新;然后,下一個點僅依賴于第1 個點,而第1 個點已經有最新數據,因此第2 個點可以進行計算,此為第2 層;再然后,我們可以分析出第3 層,第4 層,...,直到整塊數據被處理完.依照這種排序,每層數據的更新可以并行操作.通過實驗不同的塊大小,我們選出能夠適應高速緩存LDS 大小且保證一定的塊間及塊內并行度的分塊方案.

Fig.5 Data dependence and levels inside a block圖5 塊內數據依賴及分層

3.1.4 SymGS 內核優(yōu)化

在分塊著色+塊內level scheduling算法的基礎上,我們對SpMV 和SymGS 的內核函數進行了優(yōu)化.SpMV函數的優(yōu)化參考文獻[23]的方法,每線程計算一個雙精度乘法,將中間結果保存在LDS 中,然后將LDS 中的乘積進行相加.這種方法保證對矩陣A的讀取都是連續(xù)訪存(coalesced access),因而能夠充分利用內存帶寬.

對SymGS函數,我們對每一level 采用類似的方式在一個workgroup 內運行.這就需要對整個矩陣進行重排,使之在計算開始之前,就按分塊之后所需要的順序排列好.在SymGS 函數被調用時,矩陣A也能實現連續(xù)訪存.由于level 內并行度較小,我們將workgroup 大小定為一個wavefront 的大小,即64 線程.

在進行歸約操作時,我們通過實驗不同的歸約方法,選出了一種盡量充分利用處理器單元的方式.

3.2 分塊圖著色算法

由于HPCG 采用27 點stencil 格式進行離散,第3.1 節(jié)分塊方式使得每個子塊內部網格點依賴關系比較強,只有嚴格保持其依賴關系整個問題才能收斂.而上一節(jié)提到,塊內采用level scheduling 的并行方式不能有效利用硬件資源,所以本節(jié)提出一種新的分塊圖著色算法,如圖6 所示.將網格y方向一條作為一個塊進行處理,這樣塊內依賴較少,可以對其進行并行計算,同時塊內向量x可以進行重用,以提升訪存帶寬利用率.將所有子塊進行著色,每一個顏色的子塊并行執(zhí)行.著色方案需要精細設計,以確保整個算法可以以較少的迭代次數收斂,從而提升整體性能.

Fig.6 Blocked graph coloring algorithm圖6 分塊圖著色算法

該方案可以實現兩級并行:同色子塊并行執(zhí)行和子塊內部網格點并行執(zhí)行.這樣可以更好的利用GPU 的硬件資源,充分發(fā)揮硬件的性能.同時塊內可以利用LDS 對向量x進行重用,從而減少向量x從設備內存訪存的次數,提升訪存帶寬利用率.

4 圖著色算法研究

第3.2 節(jié)的方案中需要對所有子塊進行著色,且著色算法對整個HPCG 的迭代次數有很大的影響.本文首先嘗試了國際上著名的并行圖著色算法JPL[24]、CC[25]等,基于國產處理器進行了并行實現并將其運用于HPCG中.對JPL算法,本文嘗試了貪心著色策略和隨機著色策略.對CC算法,本文也進行了多種嘗試,如每輪迭代的遍歷次數選為1,2,3,13,27 這5 種.對于256×256×256 的測試規(guī)模,迭代次數介于58～67 之間,且最低迭代次數時顏色數為70.這些算法的測試結果無論顏色數還是迭代次數都偏高.

經過分析可知,傳統(tǒng)的JPL 和CC算法只考慮一層依賴,即著色時每個網格點只考慮其周圍最臨近網格點的依賴關系,而實際上依賴關系是會傳遞的,如圖7 中,1 號點和2 號點有依賴,2 號點和3 號點有依賴,則1 號點和3 號點也有依賴.考慮這種依賴傳遞將有可能在并行算法中保持更多的依賴關系,從而降低迭代次數.為了降低迭代次數的同時減少顏色數(增加并行度),本文還考慮了x,y,z這3 個方向的差異,設計了一些實驗,驗證了這3 個方向有一個為強依賴方向,其他兩個為弱依賴方向.對強依賴方向使用較多的顏色,對弱依賴方向使用較少的顏色,期望整體顏色數盡量減少.最終本文使用了33 色對整個網格進行著色.經測試,對于256×256×256 規(guī)模,單進程運行時迭代次數為55 次,相比于JPL 和CC算法的最少迭代次數降低了3 次,即整體性能可提升6%.

Fig.7 Propagation of dependence圖7 依賴傳遞

5 異構部件任務劃分

國產超級計算機每個節(jié)點包括CPU、GPU 等部件,CPU、GPU 均擁有獨立的內存.HPCG 的數據生成部分耗時較多,本文使用GPU 進行矩陣、向量和其他信息生成,并拷貝回CPU 進行參考版的計算.所有優(yōu)化版需要的矩陣、向量等常駐GPU 設備內存.整個任務劃分時需考慮CPU 和GPU 的任務分工.由于GPU擁有比CPU 更高的訪存帶寬,更多的計算資源,所以應該盡量把任務放到GPU 上計算.這樣就有兩種可能的方案:

1.將所有計算任務全部放到GPU 上,CPU 不參與計算,只負責通信.

2.將大部分計算任務放到GPU 上,CPU 并行執(zhí)行小部分計算任務同時負責通信.

第2 種方案實際執(zhí)行時需要CPU 和GPU 不斷地傳輸新的計算數據,而這種傳輸開銷較大,傳輸256×256×256 規(guī)模的向量x約需0.035s,而計算一次SpMV 的時間僅為0.008 78s,傳輸開銷遠大于計算開銷,得不償失.本文采用第1 種方案,并且通過內外區(qū)劃分將CPU 通信與GPU 內區(qū)的計算進行重疊,如圖8 所示,具體介紹見第7 節(jié).

Fig.8 Partitioning of inner areas and outer areas圖8 內外區(qū)劃分

6 HPCG 優(yōu)化方法

6.1 稀疏矩陣存儲格式

HPCG 參考版本里矩陣存儲采用的是CSR 格式,如圖9 所示,矩陣的非零元素和其對應的列索引分別存放在各自的數組里,并用一個指針數組記錄每一行的起始位置.這樣數據的內存分配比較離散,局部性不高.為了解決這個問題,我們采用了如圖10 所示的ELL 格式,按照矩陣中非零元最多的行對其他行進行填充,使每行的長度一樣,這樣數據在內存中占據了一塊連續(xù)的空間,在知道這個長度和矩陣起始地址的情況下可以算出每一行的起始位置,節(jié)省了CSR 格式中每行起始位置的索引數組.

Fig.9 Storage format used by the reference version of HPCG圖9 HPCG 參考版使用的存儲格式

Fig.10 Optimized ELL format圖10 優(yōu)化ELL 格式

對規(guī)模為256×256×256 的矩陣而言,總非零元個數為450 629 374,CSR 格式總存儲量為5.20 GB,ELL 格式總存儲量為5.06 GB,減少了3%.

6.2 稀疏矩陣重排

由于使用了分塊圖著色算法,相同顏色塊一起計算,而在原存儲空間中,相同顏色的塊都是不相鄰的,這樣無法充分利用L2 Cache,因此,我們根據分塊著色的結果對矩陣進行了重排,如圖11 所示,我們將相同顏色的塊放在一起,進一步提高數據局部性.相應的,我們對解向量x和右端項y也做了同樣的重排.

Fig.11 Blocked graph coloring and reordering圖11 分塊圖著色以及重排

6.3 其他優(yōu)化

分塊后,塊內的數據存在可能的復用,本文可以通過GPU 卡上的程序員可控高速緩存LDS 來完成這部分的復用,從而進一步提高訪存的效率.

除此之外,本文還運用了多種優(yōu)化手段,如分支消除、循環(huán)展開、數據預取、數據管理優(yōu)化等.在分支消除方面,考慮到GPU 的分支執(zhí)行特性,本文通過邏輯計算來消除了一些分支指令.在循環(huán)展開方面,綜合考慮寄存器的數量與性能,本文展開了一些熱點函數中的循環(huán).在數據預取方面,本文在計算時先統(tǒng)一預取不規(guī)則的內存區(qū)域中的數據.數據管理方面,本文采用了內存池的方式對內存進行管理,避免了很多不必要的內存申請與釋放,并且在主存分配空間時都聲明為pinned memory 以加速數據拷貝.

7 多進程實現與優(yōu)化

在多進程的實現中,本文采用了內外區(qū)劃分的方式,將整個計算網格分成內區(qū)和外區(qū)兩個部分,外區(qū)寬度為1,如圖8 所示.將內區(qū)按照第4 節(jié)中算法進行著色,將外區(qū)看成最后一個顏色.當內區(qū)部分所有的計算和外區(qū)所需的halo 區(qū)數據通信都完成后開始進行外區(qū)的計算.具體流程如圖12 所示.這樣整個HPCG 中核心函數SpMV、SymGS 中的鄰居通信均可被內區(qū)計算掩蓋,從而減少了整體運行時間.

由于矩陣、向量等常駐GPU 設備內存,鄰居通信前需要將halo 區(qū)的數據先拷貝回CPU 內存,待通信結束后再拷貝到GPU 內存.第1 次拷貝是無法重疊的,因為只有等完全拷貝結束后才能進行鄰居通信.

Fig.12 Overlap of computation and communication圖12 計算通信重疊

8 HPCG 測試結果

本文對第3.1 節(jié)和第3.2 節(jié)的算法均進行了實現,并基于某國產復雜異構超級計算機單節(jié)點單進程進行了測試,測試結果如圖13 所示.圖13 中給出了HPCG參考版、兩個并行方案的性能結果,對第2 個方案,還給出了多個優(yōu)化技術所帶來的性能提升.從圖中可以看出,第3.2 節(jié)算法性能較高(見圖中para2),約是第3.1 節(jié)算法(見圖中para1)的2 倍多,因為第3.2 節(jié)算法能更好的發(fā)揮GPU 的計算能力,并行度更高.圖中opt1 是將setup 部分和OptimizeProblem部分進行優(yōu)化的結果,同時該版本還減少了不必要的同步操作,opt2 是對多重網格中每一層的著色方案進行自適應調整,以增加粗層的并行度,另外該版本中還包括對零元訪存的優(yōu)化.opt2 版本的性能是本文工作單節(jié)點可達到的最高性能,在做完內外區(qū)劃分后,該性能有所下降,見圖中inner-outer 部分.

Fig.13 HPCG performance using single process圖13 HPCG 單進程性能

本文基于某國產復雜異構超級計算機系統(tǒng)進行測試,分別測試了單節(jié)點計算效率和整機計算效率,并測試了整機的弱可擴展性,最終單節(jié)點計算效率達到了1.82%,整機計算效率達到了1.67%,整機弱可擴展性并行效率達到了92%.

9 結論及下一步工作

本文面向國產復雜異構超級計算機研制了異構眾核并行HPCG 軟件,從著色算法、異構任務劃分、稀疏矩陣存儲格式等角度展開研究,提出了一套適用于結構化網格的著色算法,用于HPCG 后,著色質量與現有算法相比有明顯的提升.通過分析異構部件的傳輸開銷,提出了一套異構任務劃分方法,并采用ELL 格式存儲稀疏矩陣,以減少訪存開銷.同時采用分支消除、循環(huán)展開、數據預取等多種優(yōu)化方式進行了優(yōu)化.在多進程實現時,為了盡可能地提升整機性能,本文還采用了內外區(qū)劃分的算法,將鄰居通信與內區(qū)計算進行重疊,以隱藏通信開銷,提高并行效率.最終整機計算效率達到了峰值性能的1.67%,弱可擴展性并行效率更是提高到了92%.下一步將面向其他國產超級計算機開展HPCG 的并行與優(yōu)化工作,研究HPCG 混合精度算法,并與相關應用進行合作,提升應用整體性能.