黎雷生 ,楊文浩 ,馬文靜 ,張 婭 ,趙 慧 ,趙海濤 ,李會元,孫家昶

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

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

HPL(high performance Linpack)是評測計算系統(tǒng)性能的程序,是早期Linpack 評測程序的并行版本,支持大規(guī)模并行超級計算系統(tǒng)[1],其報告的每秒浮點運算次數(shù)(floating-point operations per second,簡稱FLOPS)是世界超級計算機Top500 列表排名的依據(jù).

基于Linpack 的Top500 排名開始于1993 年,每年發(fā)表兩次.2007 年11 月及之前的列表,排名前10 位的超級計算機的計算能力全部由同構(gòu)CPU提供.2008年6月Top500首臺性能超過1 PFLOPS 的超計算機Roadrunner使用異構(gòu)CPU 結(jié)構(gòu),通用CPU 執(zhí)行調(diào)度、管理等任務(wù),加速CPU 主要負責計算密集任務(wù).2009 年11 月排名第5 位的Tianhe-1 使用了CPU+GPU 的異構(gòu)架構(gòu).此后榜單上排名前10 的系統(tǒng)CPU+加速器的架構(gòu)成為趨勢.2019年6 月最新的排名顯示,前10 位中有7 臺系統(tǒng)使用CPU+加速器架構(gòu),其中使用GPU 加速器的5 臺,使用XEON Phi 的1 臺,使用Matrix-2000 的1 臺[2].

HPL 浮點計算集中在BLAS 函數(shù),特別是DGEMM.對于同構(gòu)CPU 架構(gòu),優(yōu)化BLAS 函數(shù)特別是DGEMM性能是提高HPL 的浮點效率的關(guān)鍵.Dongarra 等人[1]總結(jié)了2002 年之前Linpack 發(fā)展歷史和BLAS 函數(shù)的優(yōu)化方法.

對于CPU+加速器架構(gòu),優(yōu)化方向集中于CPU端BLAS 函數(shù)、加速器端BLAS 函數(shù)、CPU 與加速器之間負載分配和數(shù)據(jù)傳輸?shù)?Kurzak 等人[3]優(yōu)化多核CPU 多GPU 的HPL,其4 個GPU 的DGEMM 浮點性能達到1 200 GFLOPS.Bach 等人[4]面向AMD CPU 與Cypress GPU 架構(gòu)優(yōu)化了DGEMM 和HPL,HPL 效率達到70%.Wang 等人[5]采用自適應(yīng)負載分配方法優(yōu)化CPU 與GPU 混合系統(tǒng)HPL 性能,并調(diào)優(yōu)HPL 參數(shù),目標系統(tǒng)Tianhe1A 浮點性能達到0.563 PFLOPS.Heinecke 等人[6]基于Intel Xeon Phi Knights Corner 優(yōu)化了DGEMM,使用動態(tài)調(diào)度和改進的look-ahead 優(yōu)化HPL,100 節(jié)點集群HPL 效率達到76%.Gan 等人[7]優(yōu)化基于加速器的DGEMM,并應(yīng)用于Tianhe-2 HPL 評測.

隨著以GPU 為代表的加速器技術(shù)的發(fā)展,加速器浮點性能越來越高,CPU 與加速器的浮點性能差距越來越大.2019 年6 月,Top500 排名第一的Summit 系統(tǒng)[8]的1 個節(jié)點裝備2 個CPU,理論浮點性能1 TFLOPS,裝備6個GPU,理論浮點性能42 TFLOPS.本文研究的目標系統(tǒng)使用CPU+GPU 異構(gòu)架構(gòu),每個節(jié)點裝備1 個32 核CPU,4 個GPU,CPU 浮點計算性能約是GPU 的1/61.同時加速器本身的結(jié)構(gòu)也變得越來越復(fù)雜,通過增加特定的硬件滿足特定領(lǐng)域的需求,如Nvidia GPU 的Tensor Core 等.已有研究使用Tensor Core 的強大的半精度運算能力混合雙精度計算開發(fā)了HPL-AI[9],報告Summit 的HPL-AI 性能是全雙精度HPL 的2.9 倍.并且已有應(yīng)用采用混合精度算法加速計算,從HPL 和應(yīng)用角度來看,混合精度都是值得研究的方向.面對這種新的計算架構(gòu),內(nèi)存、總線、網(wǎng)絡(luò)以及系統(tǒng)設(shè)計都要與之相適應(yīng),形成復(fù)雜的異構(gòu)計算系統(tǒng),這為HPL 評測帶來很多機遇與挑戰(zhàn).

本文在基礎(chǔ)HPL 代碼之上,針對目標系統(tǒng)實現(xiàn)HPL 并研究其優(yōu)化方法.第1 節(jié)簡要介紹基礎(chǔ)算法.第2 節(jié)介紹目標系統(tǒng)基本情況.第3 節(jié)描述復(fù)雜異構(gòu)系統(tǒng)HPL 平衡點理論.第4 節(jié)介紹復(fù)雜異構(gòu)系統(tǒng)HPL 高效并行算法.第5 節(jié)介紹基礎(chǔ)模塊的性能優(yōu)化,包括panel 分解優(yōu)化和行交換long算法的優(yōu)化.第6 節(jié)介紹目標系統(tǒng)HPL實驗和結(jié)果分析.第7 節(jié)總結(jié)并展望未來的工作.

1 基礎(chǔ)算法

HPL算法使用64 位浮點精度矩陣行偏主元LU 分解加回代求解線性系統(tǒng).矩陣是稠密實矩陣,矩陣單元由偽隨機數(shù)生成器產(chǎn)生,符合正態(tài)分布.

線性系統(tǒng)定義為:

行偏主元LU 分解N×(N+1)系數(shù)矩陣[A,b]:

其中,Pr表示行主元交換矩陣,分解過程中下三角矩陣因子L已經(jīng)作用于b,解x通過求解上三角矩陣系統(tǒng)得到:

HPL 采用分塊LU算法,每個分塊是一個NB列的細長矩陣,稱為panel.LU 分解主循環(huán)采用right-looking算法,單步循環(huán)計算panel 的LU 分解和更新剩余矩陣.基本算法如圖1 所示,其中A1,1和A2,1表示panel 數(shù)據(jù).需要特別說明的是,圖示矩陣是行主順序,HPL 代碼中矩陣是列主存儲的.

Fig.1 Block LU algorithm圖1 分塊LU算法

計算公式如下:

第1 個公式表示panel 的LU 分解,第2 個公式表示求解U,一般使用DTRSM 函數(shù),第3 個公式表示矩陣更新,一般使用DGEMM 函數(shù).

對于分布式內(nèi)存計算系統(tǒng),HPL 并行計算模式基于MPI,每個進程是基本計算單元.進程組織成二維網(wǎng)格.矩陣A被劃分為NB×NB的邏輯塊,以Block-Cycle 方式分配到二維進程網(wǎng)格,數(shù)據(jù)布局示例如圖2 所示.

Fig.2 Block-Cycle algorithm圖2 Block-Cycle算法

對于具有多列的進程網(wǎng)格,單步循環(huán)只有一列進程執(zhí)行panel 分解計算,panel 分解過程中每一列執(zhí)行一次panel 的行交換算法選擇并通信最大主元行.Panel 分解計算完成后,把已分解數(shù)據(jù)廣播到其他進程列.HPL 基礎(chǔ)代碼包含6 類廣播算法,可以通過測試選擇較優(yōu)的算法.

HPL 采用行主元算法,單步矩陣更新之前,要把panel 分解時選出的最大主元行交換到U矩陣中,需要執(zhí)行未更新矩陣的主元行交換和廣播.主元行交換和廣播后,每個進程獲得完整的主元行數(shù)據(jù).

矩陣更新包括兩部分計算,一是使用DTRSM 求解U,二是使用DGEMM 更新矩陣數(shù)據(jù).

LU 分解完成后,HPL 使用回代求解x,并驗證解的正確性.

2 系統(tǒng)介紹

2.1 系統(tǒng)硬件

優(yōu)化HPL 的目標系統(tǒng)是典型的復(fù)雜異構(gòu)系統(tǒng),計算部件由通用CPU 和計算加速器GPU 組成,網(wǎng)絡(luò)接口和互聯(lián)系統(tǒng)是高速Infiniband 網(wǎng)絡(luò).

2.2 基礎(chǔ)軟件

目標CPU 兼容x86 指令集.編譯器使用gcc,MPI 使用OpenMPI.

HPL 基礎(chǔ)代碼使用2.3 版本.HPL 會調(diào)用多個BLAS 函數(shù).CPU端BLAS 使用針對目標CPU 優(yōu)化的OpenBLIS,支持OpenMP 實現(xiàn)多核并行.GPU端BLAS 使用針對目標GPU 優(yōu)化的BLAS.

3 復(fù)雜異構(gòu)系統(tǒng)HPL 平衡點理論

從硬件的角度來看,CPU 的每個Die 與對應(yīng)的內(nèi)存系統(tǒng)、GPU 和網(wǎng)絡(luò)接口組成相對獨立的系統(tǒng).因此,在選擇并行方案時,采用1 個進程使用1 個Die,并且與相應(yīng)的內(nèi)存、GPU 和網(wǎng)絡(luò)接口綁定的架構(gòu).HPL 基礎(chǔ)代碼中只有1 列或1 行進程時執(zhí)行流程不同,較大規(guī)模并行一般是多行多列進程,因此使用1 個節(jié)點4 進程作為研究對象,覆蓋多行多列進程執(zhí)行流程.

對于復(fù)雜異構(gòu)系統(tǒng)HPL 計算,需要協(xié)同調(diào)度CPU、GPU、PCIe 和通信網(wǎng)絡(luò)等資源.首先研究CPU 與GPU計算任務(wù)分配方式,進程內(nèi)CPU 與GPU 以及進程間的數(shù)據(jù)傳輸,并且建立線程模型控制CPU、GPU 和網(wǎng)絡(luò)接口協(xié)同計算.然后分析HPL 各部分的性能,提出復(fù)雜異構(gòu)系統(tǒng)的平衡點理論指導(dǎo)性能優(yōu)化.

3.1 計算任務(wù)分配

文獻[4,6]等研究的異構(gòu)系統(tǒng)中,CPU 執(zhí)行panel 分解計算,加速器執(zhí)行矩陣更新的DTRSM 和DGEMM 計算.CPU 計算能力可以達到加速器計算能力的10%或更高,除了panel 分解,CPU 還可以計算部分矩陣更新.并且由于加速器內(nèi)存限制,矩陣不能完全存儲在加速器,矩陣更新一般采用流水的方式.

經(jīng)過分析和實驗對比,由于panel 分解控制流程復(fù)雜,并且大量調(diào)用GPU 計算效率較低的DGEMV、DSCAL等函數(shù),不適合在GPU端實現(xiàn),因此,仍然采用CPU 計算panel 分解的方案.

對于矩陣更新來說,目標系統(tǒng)CPU 和GPU 之間的數(shù)據(jù)傳輸成為瓶頸,導(dǎo)致矩陣更新流水線方法不再適用.針對單個GPU 分析,假設(shè)矩陣規(guī)模N=58368(占用約80%的系統(tǒng)內(nèi)存),NB=384,矩陣更新的DGEMM 效率為85%.可以計算出矩陣更新計算時間0.52s.PCIe 單向傳輸矩陣數(shù)據(jù)的理論時間1.70s.數(shù)據(jù)傳輸時間超過GPU 的計算時間,導(dǎo)致GPU 超過2/3 處于空閑狀態(tài),GPU 計算能力不能充分利用.

因此,采用矩陣數(shù)據(jù)常駐GPU 內(nèi)存的方式,矩陣更新時避免傳輸整個剩余矩陣.這種模式下,矩陣的規(guī)模受限于GPU 設(shè)備內(nèi)存的容量.為了運算方便,CPU端內(nèi)存保留GPU端數(shù)據(jù)存儲的鏡像.

3.2 數(shù)據(jù)傳輸

數(shù)據(jù)傳輸包括進程內(nèi)主機內(nèi)存與GPU 設(shè)備內(nèi)存之間數(shù)據(jù)傳輸和進程間數(shù)據(jù)傳輸.主機內(nèi)存與GPU 設(shè)備內(nèi)存數(shù)據(jù)傳輸使用PCIe 總線,一般采用DMA 異步模式,數(shù)據(jù)傳輸?shù)耐瑫rCPU 和GPU 可以執(zhí)行計算任務(wù).進程間數(shù)據(jù)傳輸使用MPI.如果數(shù)據(jù)傳輸?shù)倪M程位于同一個節(jié)點,MPI 使用共享內(nèi)存等技術(shù)提高數(shù)據(jù)傳輸性能.節(jié)點之間的數(shù)據(jù)傳輸通過網(wǎng)絡(luò)接口.

由于矩陣數(shù)據(jù)常駐GPU 內(nèi)存,panel 分解前把當前的panel 數(shù)據(jù)從GPU 傳輸?shù)紺PU.Panel 分解在CPU端執(zhí)行panel 的主元選取與廣播,使用基礎(chǔ)代碼的算法.已完成panel 分解的數(shù)據(jù),進程內(nèi)部通過PCIe 傳輸?shù)紾PU.進程之間通過MPI 傳輸,非當前panel 分解進程接收到panel 數(shù)據(jù)后,再通過PCIe 傳輸?shù)紾PU.

相比于基礎(chǔ)代碼,矩陣更新行交換和廣播增加了CPU 與GPU 數(shù)據(jù)傳輸過程.矩陣的行在GPU 內(nèi)存不連續(xù),數(shù)據(jù)傳輸前后需要執(zhí)行封裝和重排.數(shù)據(jù)傳輸關(guān)系如圖3 所示.

Fig.3 Data transmisson圖3 數(shù)據(jù)傳輸

3.3 多線程模型

在1 個節(jié)點4 個進程的架構(gòu)下,每個進程使用8 個CPU 核心,使用多線程管理和調(diào)度這些核心.Panel 分解和panel 廣播共享數(shù)據(jù),執(zhí)行順序存在依賴關(guān)系,放在同一個線程.矩陣行交換和矩陣更新共享數(shù)據(jù),執(zhí)行順序存在依賴關(guān)系,放在同一個線程.兩組操作之間相對獨立,分別綁定到特定核心.兩個線程之間通過信號量同步.

Panel 分解線程調(diào)用和管理OpenBLIS 函數(shù),OpenBLIS 使用OpenMP 線程模型.Panel 分解調(diào)用OpenBLIS函數(shù)時,與panel 分解的其他過程沒有資源共享和沖突,OpenBLIS 的OpenMP 線程綁定的核心可以包括panel分解主線程的核心,也就是最多可綁定7 個核心.

3.4 復(fù)雜異構(gòu)系統(tǒng)協(xié)同計算平衡點理論

目標系統(tǒng)是復(fù)雜異構(gòu)系統(tǒng),涉及CPU、GPU、PCIe 和網(wǎng)絡(luò)接口等部件的協(xié)同計算.為了更好分析和評測各部分對HPL 性能的影響,指導(dǎo)性能優(yōu)化,提出平衡點理論.

HPL算法中,矩陣更新的浮點計算占據(jù)絕大部分,這部分計算由GPU 完成.性能分析過程中,矩陣更新被認為是有效計算,其他功能模塊包括panel 分解和矩陣行交換的目標是為矩陣更新提供數(shù)據(jù)準備.為了獲得較高的效率,要盡可能保證GPU 處于工作狀態(tài),發(fā)揮GPU 浮點計算能力.對于panel 分解和矩陣行交換,盡可能讓它們與GPU 計算并行,使用GPU 計算隱藏panel 分解和矩陣行交換.

研究沒有使用特殊優(yōu)化、單步循環(huán)各部分時間以及串行執(zhí)行的總時間,如圖4 所示.采用2×2 進程布局,圖中是進程0 的時間曲線.進程0 在計算過程中處于不同的位置,包括當前列-當前行和非當前行-非當前列兩種情況,圖中”A”代表當前列-當前行,”B”代表非當前行-非當前列.兩種情況下執(zhí)行的計算有差異,每種運算有兩條時間曲線.panel 分解時間在不同規(guī)模下稍有波動.可以看出,單步循環(huán)串行執(zhí)行時間遠大于GPU 執(zhí)行有效計算的時間.那么首先必須做好矩陣更新與panel 分解和行交換的并行,GPU 計算時間隱藏panel 分解與和行交換時間.

Fig.4 Time of panel factorization,row swap,matrix update and each loop圖4 Panel 分解、行交換、矩陣更新和三部分串行執(zhí)行累加時間

另一方面,HPL 的LU 分解過程分為兩個階段,矩陣更新時間大于panel 分解和行交換的時間的階段和矩陣更新時間小于panel 分解直至小于行交換時間的階段.在panel 分解與矩陣行交換并行的情況下,兩個階段的分界點即矩陣更新與panel 分解時間曲線的交點就是CPU 與GPU 協(xié)同計算的平衡點.

針對1 節(jié)點4 進程進行分析,N=88320 情況下,從計算規(guī)模來看,平衡點在37.82%,LU 計算量占76.00%,之后GPU 不能充分發(fā)揮性能.有必要繼續(xù)優(yōu)化panel 分解和行交換,減少執(zhí)行時間,平衡點后移,延長GPU 高效率計算時間,進而提高整體的計算效率.平衡點理論可以很好地指導(dǎo)HPL 的優(yōu)化工作.

4 復(fù)雜異構(gòu)系統(tǒng)HPL 高效并行算法

根據(jù)第3.4 節(jié)的分析,首先做好GPU 計算與panel 分解和行交換的并行算法,實現(xiàn)GPU 計算重疊其他部分的時間.

4.1 GPU參與計算的look-ahead算法

基礎(chǔ)代碼用look-ahead算法實現(xiàn)panel 分解與panel 廣播重疊執(zhí)行.計算過程中保存兩個panel 的數(shù)據(jù)結(jié)構(gòu).執(zhí)行當前panel 分解的進程列,優(yōu)先執(zhí)行當前panel 列數(shù)據(jù)的行交換和矩陣更新,然后執(zhí)行panel 分解計算并發(fā)起panel 廣播,廣播數(shù)據(jù)傳輸?shù)耐瑫r更新剩余矩陣.對于當前不計算panel 分解的進程列,執(zhí)行上一次循環(huán)的行交換和矩陣更新,行交換和矩陣更新計算過程中監(jiān)測本次panel 廣播數(shù)據(jù)通信,數(shù)據(jù)到達后,執(zhí)行廣播通信流程.

由于加入了GPU 計算資源,look-ahead算法相比于基礎(chǔ)版本有所改變.按照第3.3 節(jié)的線程模型,在執(zhí)行當前panel 分解計算的進程中,panel 分解線程CPU 優(yōu)先執(zhí)行panel 列的行交換,然后GPU 優(yōu)先更新panel 分解的列.CPU 等待panel 數(shù)據(jù)更新后,發(fā)起panel 數(shù)據(jù)從GPU 內(nèi)存到CPU 內(nèi)存的傳輸,傳輸結(jié)束后執(zhí)行panel 分解計算,最后發(fā)起panel 廣播.矩陣更新線程在panel 分解列行交換完成后發(fā)起剩余矩陣的行交換,然后調(diào)用GPU 矩陣更新計算.對于非當前panel 分解計算進程,矩陣更新線程執(zhí)行行交換和調(diào)用GPU 矩陣更新.Panel 分解線程只負責處理panel 廣播的操作.

4.2 行交換連續(xù)流水線算法

矩陣行交換的前提是已接收panel 分解計算的行交換索引數(shù)據(jù)(包含在panel 廣播數(shù)據(jù)中),并且上一次迭代矩陣更新已完成.HPL 基礎(chǔ)代碼中,行交換完成后啟動矩陣更新計算,如果采用這種算法,GPU 在行交換執(zhí)行的過程是空閑的.文獻[6]提出行交換流水線算法優(yōu)化行交換流程.在使用GPU 的環(huán)境,也可以采用這一算法.使用行交換流水線算法,對矩陣更新的行分段,首先執(zhí)行第1 段的行交換,然后GPU 執(zhí)行第1 段的矩陣更新,同時執(zhí)行第2 段的行交換,依次類推.采用這種算法,GPU 計算與行交換數(shù)據(jù)傳輸重疊執(zhí)行,減少了GPU 的空閑時間.但是,單步循環(huán)執(zhí)行第1 段行交換的這段時間內(nèi)GPU 是空閑的.

為了避免GPU 等待,提出了連續(xù)流水線算法.算法中,第1 個分段完成上一次循環(huán)矩陣更新和當前進程接收到下一次循環(huán)的行交換信息之后,不需要等待上一次循環(huán)矩陣更新全部完成,就執(zhí)行第1 分段的下一次循環(huán)的行交換.采用這種算法,下一次循環(huán)第1 分段的行交換被上一次循環(huán)矩陣更新隱藏,下一次循環(huán)第1 分段的矩陣更新隱藏第2 分段的行交換,后續(xù)分段繼續(xù)流水.

使用連續(xù)流水算法避免單步循環(huán)的流水線啟動過程,如果接收到行交換信息早于矩陣更新,則GPU 可以一直處于工作狀態(tài).從第3.4 節(jié)可以看出,行交換信息在平衡點之前是早于矩陣更新的,而這一階段正是浮點運算量較大的情況,在這一階段充分發(fā)揮GPU 計算能力,可以提高整體的計算效率.

行交換分段方法選擇倍數(shù)遞增,也就是第2 段列數(shù)是第1 段列數(shù)的倍數(shù),依次類推,增長因子可調(diào).連續(xù)行交換流水算法下實現(xiàn)look-ahead算法,對于當前panel 分解,增加panel 分解數(shù)據(jù)提前行交換和矩陣更新的分段.

5 基礎(chǔ)模塊性能優(yōu)化

從平衡點理論來看,平衡點之前的部分GPU 計算時間可以隱藏panel 分解和行交換時間.隨著計算的推進,panel 分解和行交換成為系統(tǒng)性能的瓶頸.為了進一步提高效率,需要針對panel 分解和行交換進行優(yōu)化.

5.1 Panel分解優(yōu)化

5.1.1 基本參數(shù)調(diào)優(yōu)

Panel 分解浮點計算集中在BLAS 函數(shù),使用針對目標系統(tǒng)優(yōu)化的Hygon OpenBLIS庫.Panel 分解使用遞歸算法,中間遞歸層次的浮點運算集中在BLAS 的DTRSM 和DGEMM 兩個函數(shù).當遞歸層次包含的列數(shù)小于等于閾值時,使用非遞歸算法,浮點運算集中在BLAS 的DGEMV、DTRSV、DSCAL、IDAMAX 等函數(shù).對于每一列的LU 分解,需要通過交換和廣播通信主元所在的panel 行,并記錄主元行交換信息.交換和廣播通信使用binary-exchange算法.

Panel 分解計算有關(guān)參數(shù)有NB、NBMIN、PFACT、RFACT、DIV等.NB取決于GPU 執(zhí)行矩陣更新的效率,同時考慮CPU 與GPU 計算的平衡.當節(jié)點規(guī)模較小時,NB＝384;當節(jié)點規(guī)模較大時,NB＝256.通過參數(shù)調(diào)優(yōu),選擇優(yōu)化的參數(shù)組合,見表1.

Table 1 Parameters of panel factorization表1 Panel 分解參數(shù)

進一步分析panel 分解各部分時間,0 號進程主要計算函數(shù)時間見表2.

Table 2 Time of BLAS functions in panel factorization表2 Panel 分解BLAS 函數(shù)時間

5.1.2 GPU 加速panel 分解DGEMM

從第5.1.1 節(jié)可以看出,panel 分解中,DGEMM 時間占有最大比例,需要進一步優(yōu)化.DGEMM 在panel 分解遞歸層次調(diào)用.根據(jù)panel 分解的left-looking算法,從左到右執(zhí)行LU 分解,左側(cè)的subpanel 完成分解之后,執(zhí)行DGEMM,更新相同層次的右側(cè)的subpanel,GPU 可以加速這部分DGEMM.

GPU 加速DGEMM 的算法過程,首先把左側(cè)subpanel 和相應(yīng)的U數(shù)據(jù)傳輸?shù)紾PU 設(shè)備內(nèi)存,然后執(zhí)行DGEMM 更新GPU 內(nèi)存的右側(cè)subpanel,更新后數(shù)據(jù)傳輸?shù)紺PU端內(nèi)存,繼續(xù)執(zhí)行后續(xù)的panel 分解.算法增加了CPU 和GPU 之間的數(shù)據(jù)傳輸,但由于GPU 執(zhí)行DGEMM 的速度遠高于CPU,在一定規(guī)模DGEMM 的情況下,這一過程總的執(zhí)行時間少于CPU端執(zhí)行DGEMM 時間,整體性能提高.

并不是所有的panel 分解DGEMM 都可以使用GPU 加速,如果傳輸?shù)臅r間開銷大于DGEMM 加速的效果,就不能采用這種方式.實際計算中,使用參數(shù)控制采用GPU 加速DGEMM 的閾值,通過調(diào)優(yōu)獲得最佳的性能.

5.1.3 Panel 廣播優(yōu)化

5.1.3.1 避免數(shù)據(jù)封裝

基礎(chǔ)代碼中,panel 分解數(shù)據(jù)使用CPU端矩陣A的內(nèi)存,A是列存儲,panel 分解的列是分段連續(xù)的.使用MPI傳輸panel 數(shù)據(jù)之前,需要執(zhí)行數(shù)據(jù)封裝,把數(shù)據(jù)復(fù)制到發(fā)送緩沖區(qū),而復(fù)制操作會帶來一定的時間開銷.

為此提出一種避免數(shù)據(jù)封裝的方法.待分解的panel 數(shù)據(jù)是從GPU 內(nèi)存復(fù)制的,通過使用二維復(fù)制接口把panel 數(shù)據(jù)復(fù)制到連續(xù)存儲區(qū)域.Panel 分解計算完成后,MPI 接口直接使用緩沖區(qū)數(shù)據(jù),避免了數(shù)據(jù)封裝.

5.1.3.2 Panel 廣播流水

基礎(chǔ)代碼中,panel 分解完成后才會執(zhí)行panel 廣播,計算與數(shù)據(jù)傳輸串行,并且一次性數(shù)據(jù)傳輸較大,傳輸時間較長.通過觀察可以看出,left-looking算法中,左側(cè)subpanel計算完成后,主體部分不再變化,后續(xù)只有發(fā)生行交換的數(shù)據(jù)變化.

因此,提出一種廣播流水算法,對已分解的subpanel 數(shù)據(jù)提前發(fā)起廣播,這種情況下,panel 數(shù)據(jù)廣播與后續(xù)panel 分解計算并行.采用廣播流水算法,panel 計算完成后,只需要傳輸最后的subpanel 和行交換發(fā)生變化的數(shù)據(jù),相比于傳輸整個panel 數(shù)據(jù),縮短了傳輸時間.

廣播流水算法與GPU 加速panel 分解DGEMM 協(xié)同使用,subpanel 的數(shù)據(jù)廣播與panel 加速DGEMM 過程中的CPU 與GPU 之間數(shù)據(jù)傳輸、GPU 執(zhí)行DGEMM 計算并行,充分利用系統(tǒng)的CPU、GPU、PCIe 和網(wǎng)絡(luò)接口資源.

HPL 包含的6 種廣播算法都可以使用廣播流水優(yōu)化.

第5.1.2 節(jié)描述的GPU 加速panel 分解DGEMM 和panel 廣播流程緊密耦合,表3 給出了優(yōu)化前后偽代碼對比.

Table 3 Pseudocode of accelerated panel DGEMM and broadcast pipeline表3 Panel DGEMM 加速和廣播流水偽代碼

5.2 行交換long算法優(yōu)化

矩陣更新行交換有兩種算法可選,bianary-exchange算法和long算法,經(jīng)過分析和評測選擇long算法[10].第4.2 節(jié)從行交換與GPU 并行計算的角度針對行交換流程做了優(yōu)化,本節(jié)優(yōu)化具體的行交換算法.

Long算法包括兩個步驟,spread 和roll.基礎(chǔ)代碼中,首先是spread 過程,當前行進程將需要換出的行分發(fā)給其他所有進程.然后,非當前行進程把接收的數(shù)據(jù)與本地需要換出的數(shù)據(jù)交換.接著是roll 過程,總共需要進行P–1 次傳輸.每次傳輸所有的進程都要參與,每個進程只與自己的鄰居交換.第1 步,先把自己所擁有的一部分與左(右)鄰居交換,然后,把自己剛剛得到的數(shù)據(jù)與右(左)鄰居交換.這樣左右交替,直到所有的數(shù)據(jù)都交換完畢,每個進程都擁有了全部的U.最后所有進程還需要把所有的數(shù)據(jù)進行重排,放置到矩陣A相應(yīng)的位置.

基礎(chǔ)代碼中,spread 和roll 使用同一塊數(shù)據(jù)緩沖區(qū),為了避免數(shù)據(jù)沖突,spread 和roll 是嚴格串行的.通過觀察發(fā)現(xiàn),spread 的接收數(shù)據(jù)與roll 的發(fā)送數(shù)據(jù)并沒有依賴關(guān)系,spread 和roll算法上是可并行性.基于以上觀察與分析,將spread 的接收緩沖與roll 的發(fā)送緩沖分離.這樣,非當前進程行在等待接收spread 數(shù)據(jù)時,就可以執(zhí)行roll數(shù)據(jù)封裝,并把數(shù)據(jù)傳輸?shù)桨l(fā)送緩沖區(qū),這部分時間就可以被隱藏.

在使用上述的緩沖區(qū)分離算法之后,spread 和roll 已經(jīng)沒有依賴關(guān)系,因此提出了一種新的行交換方法,即交換spread 和roll 的執(zhí)行順序,如圖5(b)所示.對當前行進程來說,在roll 操作之前,只需將本地需要交換的數(shù)據(jù)拷貝到發(fā)送緩沖區(qū)即可,也就是減少了行交換的啟動時間.在roll 執(zhí)行網(wǎng)絡(luò)傳輸?shù)耐瑫r,當前行進程可以將spread 所需的數(shù)據(jù)封裝,并異步傳輸?shù)桨l(fā)送緩沖區(qū).roll 執(zhí)行完成后,開啟執(zhí)行spread.同時,把roll 過程接收的數(shù)據(jù)異步傳輸?shù)紾PU,并執(zhí)行數(shù)據(jù)交換.在spread 過程結(jié)束后,非當前行進程再將換入的數(shù)據(jù)傳輸?shù)紾PU,并交換到相應(yīng)位置,當前行進程對需要進行內(nèi)部交換的數(shù)據(jù)執(zhí)行GPU 上的本地交換.這種算法可以實現(xiàn)網(wǎng)絡(luò)通信傳輸、GPU 數(shù)據(jù)封裝和數(shù)據(jù)交換計算以及CPU 與GPU 數(shù)據(jù)傳輸相互重疊,減少行交換執(zhí)行時間.

Fig.5 Original row swap algorithm (a),and modified row swap procedure (b)圖5 原始行交換算法(a)及新的行交換算法(b)

6 實驗與分析

6.1 單節(jié)點性能分析

首先針對單節(jié)點實驗,分析評估各種優(yōu)化方法的效果.根據(jù)平衡點理論,優(yōu)化的目標是GPU 盡可能地處于高效工作狀態(tài),執(zhí)行有效的計算,也即GPU 盡可能地連續(xù)執(zhí)行矩陣更新.圖6 對比了優(yōu)化后單步循環(huán)的時間和GPU 矩陣更新時間.圖6 給出了1 個節(jié)點內(nèi)0 號進程的時間曲線.”A”表示當前行-當前列,”B”表示非當前行-非當前列.可以看出,N=88320 情況下,LU 分解矩陣規(guī)模56.09%之前單步循環(huán)時間接近GPU 矩陣更新時間.56.09%之后,當前進程執(zhí)行panel 分解,單步循環(huán)時間仍然接近GPU 矩陣更新時間,GPU 大部分時間處于工作狀態(tài),說明算法優(yōu)化已經(jīng)達到較好的效果.對于非當前panel 分解進程,增加了等待panel 廣播數(shù)據(jù)的時間開銷,GPU 有一段時間處于空閑狀態(tài).

Fig.6 GPU time vs.total time of each loop圖6 GPU 時間與單步循環(huán)總時間

可以看出,圖6 優(yōu)化版本的平衡點是56.09%,圖4 未優(yōu)化版本版本是37.82%.從LU 分解的計算量來看,前56.09%占總的LU 分解計算量的91.53%,大部分時間充分發(fā)揮了GPU 的計算能力.優(yōu)化后單節(jié)點HPL 效率達到了79.51%.

和優(yōu)化有關(guān)的參數(shù)主要是GPU 加速DGEMM 和panel 廣播流水的層次,這里都是3.行交換非當前列第1段列數(shù)是待更新列數(shù)的1/8,增長因子是2.Panel 分解參數(shù)和表1 一致.

6.2 可擴展分析

多節(jié)點HPL 較大規(guī)?？蓴U展,目前最大規(guī)模4 096 個節(jié)點.HPL 較大規(guī)模并行有關(guān)的參數(shù)和效率見表4,表中列出了針對相應(yīng)節(jié)點數(shù)執(zhí)行的實驗中效率最高的配置.

Table 4 Parameters and efficiency of multi-nodes HPL表4 多節(jié)點HPL 參數(shù)和效率

參數(shù)N是較大規(guī)模節(jié)點HPL 評測中比較重要的參數(shù),每個節(jié)點計算的矩陣列數(shù)一般是NB的倍數(shù),多節(jié)點的N正比于節(jié)點數(shù)平方根.N值越大,計算效率較高的矩陣更新占比越大,一般來說HPL 效率越高,但N受限于GPU 的內(nèi)存容量,同時還要考慮系統(tǒng)的負載.NB使用384 或256.P和Q確定二維進程網(wǎng)格的組織方式,經(jīng)過驗證,P=Q時目標系統(tǒng)性能最優(yōu).Panel 廣播和矩陣行交換分別使用優(yōu)化的long算法和廣播流水算法.

從HPL 效率來看,4 節(jié)點效率高于單節(jié)點效率,原因在于4 節(jié)點矩陣規(guī)模大于單節(jié)點,效率較高的計算占比較大,并且可以抵消4 節(jié)點規(guī)模下網(wǎng)絡(luò)傳輸路徑增加帶來的性能損失.隨著規(guī)模的擴大,網(wǎng)絡(luò)傳輸對性能的影響增大,HPL 效率逐步降低,但降低的趨勢比較平緩.

7 結(jié)論與未來的工作

針對復(fù)雜異構(gòu)系統(tǒng),分析HPL算法的特點,提出CPU 與GPU 協(xié)同計算方法,提出平衡點理論指導(dǎo)各方面優(yōu)化.實現(xiàn)CPU 與GPU 協(xié)同計算的look-ahead、行交換連續(xù)流水線算法隱藏panel 分解和行交換時間,并優(yōu)化panel 分解和行交換算法,延長GPU 計算占優(yōu)的時間,最終提高整個系統(tǒng)的HPL 效率,單節(jié)點效率79.51%.

當前的超級計算機正向百億億級邁進,目標系統(tǒng)采用的節(jié)點內(nèi)通用CPU 裝備加速器,配合高速的節(jié)點內(nèi)總線被認為是較有潛力的一種架構(gòu).未來工作主要是針對百億億級計算架構(gòu),開展加速器DGEMM 優(yōu)化、CPU 與加速器協(xié)同計算和混合精度等相關(guān)研究.同時,百億億級超級計算機節(jié)點規(guī)模繼續(xù)增加,節(jié)點數(shù)量有可能達到數(shù)萬甚至更多,節(jié)點間的互連也更加重要.需要進一步研究HPL 大規(guī)模可擴展能力,評測行交換算法和廣播算法的性能,并進行相應(yīng)的優(yōu)化.