國(guó)產(chǎn)異構(gòu)系統(tǒng)上的HPCG 并行算法及高效實(shí)現(xiàn)?

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

1(中國(guó)科學(xué)院 軟件研究所 并行軟件與計(jì)算科學(xué)實(shí)驗(yàn)室,北京 100190)

2(中國(guó)科學(xué)院大學(xué),北京 100049)

3(計(jì)算機(jī)科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室(中國(guó)科學(xué)院 軟件研究所),北京 100190)

4(北京大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,北京 100871)

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

目前超級(jí)計(jì)算機(jī)的體系結(jié)構(gòu)日趨復(fù)雜,早已從同構(gòu)轉(zhuǎn)向異構(gòu),且同節(jié)點(diǎn)異構(gòu)部件的運(yùn)算性能差異越來(lái)越大,這對(duì)異構(gòu)并行算法的設(shè)計(jì)將有很大的影響.本文主要面向一類國(guó)產(chǎn)復(fù)雜異構(gòu)系統(tǒng)研究HPCG 并行算法和高效實(shí)現(xiàn)方法.目前Dongarra 教授課題組研究了HPL-AI 混合精度算法[5,6],我們也正在關(guān)注該方向的研究.

本文工作是在國(guó)產(chǎn)超級(jí)計(jì)算機(jī)研制方大力支持下完成的,在此表示感謝.

1 HPCG 簡(jiǎn)介

HPCG 來(lái)源于半結(jié)構(gòu)網(wǎng)格上的三維熱傳導(dǎo)應(yīng)用,其核心是在三維規(guī)則區(qū)域上對(duì)Poisson 方程進(jìn)行有限差分離散化后,轉(zhuǎn)換成一個(gè)稀疏線性方程組的求解問(wèn)題.其方程如式(1)所示.采用如圖1 所示的27 點(diǎn)stencil 格式進(jìn)行離散化.

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

在大規(guī)模并行環(huán)境中,HPCG 使用三維區(qū)域分解策略,也就是按照3 個(gè)維度將整個(gè)計(jì)算區(qū)域劃分成若干個(gè)子區(qū)域,每個(gè)子區(qū)域被分配一個(gè)MPI 進(jìn)程,其規(guī)模由輸入?yún)?shù)指定.每個(gè)MPI 進(jìn)程處理的計(jì)算區(qū)域中的內(nèi)部點(diǎn)所依賴的26 個(gè)鄰點(diǎn),全部來(lái)自于該計(jì)算區(qū)域,而對(duì)于計(jì)算區(qū)域邊界上的點(diǎn),其依賴的鄰點(diǎn)一部分來(lái)自于該計(jì)算區(qū)域內(nèi)部,一部分則來(lái)自于鄰居進(jìn)程所處理的子區(qū)域.不妨稱計(jì)算區(qū)域?yàn)閮?nèi)區(qū),其對(duì)應(yīng)的點(diǎn)為內(nèi)點(diǎn),計(jì)算區(qū)域所依賴的鄰居進(jìn)程的點(diǎn)構(gòu)成的區(qū)域?yàn)橥鈪^(qū),其對(duì)應(yīng)的點(diǎn)為外點(diǎn).由于我們使用的是27 點(diǎn)stencil,所以每個(gè)進(jìn)程最多有26 個(gè)鄰居進(jìn)程.

Fig.1 3D 27-point stencil in HPCG圖1 HPCG 三維27 點(diǎn)stencil 格式

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

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

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

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

Fig.3 Multigrid pre-conditioner圖3 多重網(wǎng)格預(yù)條件子

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

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

2.鄰居通信.主要發(fā)生在與稀疏矩陣相關(guān)的兩個(gè)核心函數(shù)中:SpMV 和SymGS,使用MPI 點(diǎn)到點(diǎn)的通信實(shí)現(xiàn)halo 區(qū)數(shù)據(jù)交換.HPCG 采用27 點(diǎn)stencil 格式,每個(gè)網(wǎng)格點(diǎn)的計(jì)算依賴周圍緊鄰的26 個(gè)鄰居.在SpMV 和SymGS 計(jì)算開(kāi)始前,需要與周圍的鄰居進(jìn)程交換邊界數(shù)據(jù)(halo),以便將計(jì)算過(guò)程中所需要訪問(wèn)到的鄰居數(shù)據(jù)傳輸?shù)奖镜剡M(jìn)程中.

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

2 國(guó)內(nèi)外進(jìn)展介紹

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

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

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

3 異構(gòu)眾核并行算法設(shè)計(jì)

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

3.1 分塊著色+level scheduling算法

3.1.1 基本算法

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

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

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

3.1.2 數(shù)據(jù)及任務(wù)映射

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

Fig.5 Data dependence and levels inside a block圖5 塊內(nèi)數(shù)據(jù)依賴及分層

3.1.4 SymGS 內(nèi)核優(yōu)化

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

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

在進(jìn)行歸約操作時(shí),我們通過(guò)實(shí)驗(yàn)不同的歸約方法,選出了一種盡量充分利用處理器單元的方式.

3.2 分塊圖著色算法

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

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

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

4 圖著色算法研究

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

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

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

5 異構(gòu)部件任務(wù)劃分

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

1.將所有計(jì)算任務(wù)全部放到GPU 上,CPU 不參與計(jì)算,只負(fù)責(zé)通信.

2.將大部分計(jì)算任務(wù)放到GPU 上,CPU 并行執(zhí)行小部分計(jì)算任務(wù)同時(shí)負(fù)責(zé)通信.

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

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

6 HPCG 優(yōu)化方法

6.1 稀疏矩陣存儲(chǔ)格式

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

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

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

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

6.2 稀疏矩陣重排

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

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

6.3 其他優(yōu)化

分塊后,塊內(nèi)的數(shù)據(jù)存在可能的復(fù)用,本文可以通過(guò)GPU 卡上的程序員可控高速緩存LDS 來(lái)完成這部分的復(fù)用,從而進(jìn)一步提高訪存的效率.

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

7 多進(jìn)程實(shí)現(xiàn)與優(yōu)化

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

由于矩陣、向量等常駐GPU 設(shè)備內(nèi)存,鄰居通信前需要將halo 區(qū)的數(shù)據(jù)先拷貝回CPU 內(nèi)存,待通信結(jié)束后再拷貝到GPU 內(nèi)存.第1 次拷貝是無(wú)法重疊的,因?yàn)橹挥械韧耆截惤Y(jié)束后才能進(jìn)行鄰居通信.

Fig.12 Overlap of computation and communication圖12 計(jì)算通信重疊

8 HPCG 測(cè)試結(jié)果

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

Fig.13 HPCG performance using single process圖13 HPCG 單進(jìn)程性能

本文基于某國(guó)產(chǎn)復(fù)雜異構(gòu)超級(jí)計(jì)算機(jī)系統(tǒng)進(jìn)行測(cè)試,分別測(cè)試了單節(jié)點(diǎn)計(jì)算效率和整機(jī)計(jì)算效率,并測(cè)試了整機(jī)的弱可擴(kuò)展性,最終單節(jié)點(diǎn)計(jì)算效率達(dá)到了1.82%,整機(jī)計(jì)算效率達(dá)到了1.67%,整機(jī)弱可擴(kuò)展性并行效率達(dá)到了92%.

9 結(jié)論及下一步工作

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