蔡 雨,孫成國,杜朝暉,劉子行,康夢博,李雙雙

(信息技術(shù)有限公司,江蘇 蘇州 215000)

BLAS(basic linear algebra subprograms)是基本線性代數(shù)子程序的縮寫,是目前應(yīng)用廣泛的核心線性代數(shù)數(shù)學(xué)庫.隨著不斷的發(fā)展完善,BLAS 在高性能計(jì)算、科學(xué)和工程領(lǐng)域都得到了廣泛的應(yīng)用.由于不同廠商、不同體系結(jié)構(gòu)下的CPU 差異較大,通用BLAS 函數(shù)庫不能達(dá)到很好的性能.因此,在大規(guī)模矩陣運(yùn)算時,針對不同CPU 的架構(gòu)特點(diǎn)來優(yōu)化BLAS,可以充分發(fā)揮處理器的計(jì)算性能,從而使計(jì)算效率大大提升.目前主流的BLAS分為專用和通用兩種.專用型BLAS 是由特定的CPU 公司開發(fā),針對特定平臺芯片進(jìn)行代碼優(yōu)化,極大地提高了性能,Intel 的MKL 和AMD 的ACML 就是此種類型.通用BLAS 由非盈利機(jī)構(gòu)開發(fā),不同平臺都可以對開源代碼進(jìn)行優(yōu)化,最具代表性的是ATLAS[1]、GotoBLAS[2,3]和OpenBLAS[4,5].

BLIS(BLAS-like library instantiation software)數(shù)學(xué)庫是美國Texas 大學(xué)超算中心高性能組開發(fā)的開源架構(gòu),是在GotoBLAS 的基礎(chǔ)上進(jìn)行了改寫和優(yōu)化,具有可移植性、易用性和模塊化設(shè)計(jì)易于開發(fā)的特點(diǎn),支持混合類型矩陣存儲和多線程[6].文獻(xiàn)[7]提出了BLIS 架構(gòu),并做了全面的介紹;文獻(xiàn)[8]研究了如何使用BLIS 在通用、低功耗和多核架構(gòu)中實(shí)現(xiàn)高效的level-3 級BLAS 函數(shù);文獻(xiàn)[9]系統(tǒng)地探討了矩陣乘法算法在BLIS 五層循環(huán)中并行化的機(jī)會.BLIS 具有的這些特性使其更容易針對具體平臺進(jìn)行BLAS 優(yōu)化開發(fā).BLAS 包含三級函數(shù),分別為第1 級向量與向量計(jì)算,第2 級向量與矩陣計(jì)算,第3 級矩陣與矩陣的計(jì)算,其中第3 級函數(shù)DGEMM 計(jì)算量最大,可以達(dá)到理論計(jì)算量的95%以上[10],因此三級函數(shù)的優(yōu)化是性能提升的關(guān)鍵.

高性能計(jì)算系統(tǒng)可以分為同構(gòu)或異構(gòu)兩類,同構(gòu)系統(tǒng)是由大量CPU 構(gòu)建,計(jì)算和數(shù)據(jù)處理都在CPU 上進(jìn)行;異構(gòu)系統(tǒng)是CPU 和加速部件組成,加速部件一般由Cell(cell broadband engine architecture)、FPGA(Field programmable gate array)和GPU(graphic processing unit)等為代表,集成大量的浮點(diǎn)計(jì)算單元的同時舍棄了一些CPU 上復(fù)雜的控制單元.異構(gòu)系統(tǒng)中CPU 負(fù)責(zé)任務(wù)調(diào)度和簡單計(jì)算,主要計(jì)算在加速部件上進(jìn)行,計(jì)算能力更強(qiáng).國際上對高性能系統(tǒng)的評測標(biāo)準(zhǔn)程序是HPL(high performance Linpack),由美國田納西大學(xué)教授Dongarra 設(shè)計(jì)提出[11],是在高性能計(jì)算系統(tǒng)中使用高斯消元法求解一個隨機(jī)的N 元一次線性方程組問題,用以評價高性能計(jì)算機(jī)的浮點(diǎn)性能.

異構(gòu)HPL 浮點(diǎn)計(jì)算在CPU端是通過BLAS 函數(shù)完成的,BLAS庫的性能會影響整體系統(tǒng)性能.本文異構(gòu)系統(tǒng)計(jì)算單元由CPU+協(xié)處理器(co-processor)組成(如圖1 所示).CPU 計(jì)算核心之間和協(xié)處理器通過高速互聯(lián)進(jìn)行數(shù)據(jù)交換和通信.協(xié)處理器的計(jì)算能力一般都比通用CPU 強(qiáng)大很多,往往造成計(jì)算負(fù)載嚴(yán)重不均衡,損失性能.因此在CPU端利用更高效的數(shù)學(xué)庫,能夠充分地發(fā)揮CPU 性能,使負(fù)載更均衡,是提高計(jì)算能力的可行方式之一.文獻(xiàn)[12?15]分別針對ARM 架構(gòu)處理器以及異構(gòu)平臺的GEMM 優(yōu)化做了深入研究.文獻(xiàn)[16,17]研究了在申威眾核處理器上,1 級、2 級和3 級BLAS 的優(yōu)化方法,文獻(xiàn)[18]對異構(gòu)多核平臺上的數(shù)學(xué)庫寄存器分配優(yōu)化方法進(jìn)行了研究.

Fig.1 Heterogeneous computing node architecture圖1 異構(gòu)計(jì)算節(jié)點(diǎn)架構(gòu)

本文在BLIS算法庫基礎(chǔ)上,分析異構(gòu)HPL算法中調(diào)用的各級BLAS 函數(shù),針對體系結(jié)構(gòu)特點(diǎn),利用訪存優(yōu)化、指令集優(yōu)化和多線程并行等技術(shù)手段形成了優(yōu)化的HBLIS庫,提高了異構(gòu)系統(tǒng)中CPU端的計(jì)算效率,后文詳細(xì)分析了優(yōu)化策略和方法.本文第1 節(jié)首先對異構(gòu)HPL 的BLIS庫進(jìn)行分析,特別是異構(gòu)環(huán)境的HPL 特點(diǎn)以及各個數(shù)學(xué)庫函數(shù)調(diào)用關(guān)系.第2 節(jié)～第4 節(jié)對高性能BLIS庫的優(yōu)化策略做詳細(xì)的分析和介紹,優(yōu)化HPL 中涉及的各級BLAS 函數(shù),實(shí)現(xiàn)高性能的HBLIS算法庫.第5 節(jié)在異構(gòu)平臺上做詳細(xì)的性能測試與分析,通過性能對比,驗(yàn)證優(yōu)化的HBLIS庫性能有大幅度的提高.第6 節(jié)對全文進(jìn)行總結(jié)和展望.

1 算法分析和CPU 微體系結(jié)構(gòu)

1.1 HPL算法與BLAS函數(shù)

HPL算法核心計(jì)算公式如式(1)所示,其求解過程實(shí)際是對系數(shù)矩陣[A,b]使用高斯消元法進(jìn)行LU 分解,變成[A,b]=[[L,U],Y],其中L是下三角矩陣,U為上三角矩陣,隨著因式分解的進(jìn)行,方程等價于求解Ux=Y.LU 分解過程中每一次迭代計(jì)算都要經(jīng)過panel 分解、panel 廣播、行交換和尾矩陣更新.HPL 程序在進(jìn)程通信和矩陣計(jì)算上需要MPI[19]庫和BLAS 或VSIPL[20]庫的支持,除基本算法不可以改變外,可以通過配置文件調(diào)節(jié)問題規(guī)模N(矩陣大小)、進(jìn)程數(shù)等測試參數(shù),以及使用各種優(yōu)化方法來執(zhí)行HPL 程序,以獲取最佳的性能.當(dāng)求解問題規(guī)模為N時,浮點(diǎn)運(yùn)算次數(shù)為2/3×N3–2×N2.因此,給出問題規(guī)模N,測得系統(tǒng)計(jì)算時間T,計(jì)算系統(tǒng)的浮點(diǎn)計(jì)算能力=計(jì)算量(2/3×N3–2×N2)/計(jì)算時間(T),測試結(jié)果以每秒浮點(diǎn)運(yùn)算次數(shù)(FLOPS)為單位.

HPL算法過程分為LU 分解和回代過程.異構(gòu)HPL 環(huán)境下,LU 分解過程中大量計(jì)算要在CPU 和協(xié)處理器之間進(jìn)行任務(wù)分配,CPU 負(fù)責(zé)調(diào)度和適當(dāng)?shù)妮o助計(jì)算,而主要計(jì)算要在加速部件上進(jìn)行.當(dāng)問題規(guī)模N非常大時,HPL 的最大計(jì)算量是更新尾矩陣的計(jì)算,包括雙精度的三角矩陣求解和稠密矩陣乘 DGEMM,因此將DGEMM 放在協(xié)處理器上計(jì)算,CPU 負(fù)責(zé)計(jì)算量不大的panel 分解、panel 廣播和行交換.異構(gòu)高性能計(jì)算系統(tǒng)優(yōu)化的難點(diǎn)主要是在主從芯片對計(jì)算任務(wù)的劃分和通信開銷這兩個方面,一些學(xué)者從任務(wù)靜態(tài)劃分和動態(tài)劃分實(shí)現(xiàn)負(fù)載均衡[12,21?23],還有學(xué)者在數(shù)據(jù)重用、存儲優(yōu)化、軟件流水和應(yīng)用雙緩沖等方面深入研究[24?26],達(dá)到實(shí)現(xiàn)更高效的通信隱藏的目的.異構(gòu)HPL 主程序算法中任務(wù)分配調(diào)用加速部件專用編程接口進(jìn)行計(jì)算,CPU端則調(diào)用BLAS庫進(jìn)行選主元、panel 分解計(jì)算.

本文針對的異構(gòu)HPL算法由中國科學(xué)院軟件研究所開發(fā)[27],其LU 分解過程會將panel 的一部分矩陣(N×NBMIN)分解為單位下三角L1、上三角U和L2,這3 個部分的運(yùn)算在CPU 上進(jìn)行,主要涉及一級BLAS 函數(shù)中的DCOPY、IDMAX、DSCAL;二級BLAS 函數(shù)中的DTRSM、DTRSV、DGEMV;三級BLAS 函數(shù)中的DGEMM.各個BLAS 函數(shù)作用如下:

(1) DCOPY 執(zhí)行數(shù)據(jù)拷貝操作,將一個向量拷貝到另一個向量.

(2) IDMAX 用于選主元操作,在一列元素中選擇最大的元素.

(3) DSCAL 用于向量更新,用一個標(biāo)量乘以一個向量.

(4) DTRSM 用于三角矩陣求解操作.

(5) DTRSV 用于求解三角方程操作.

(6) DGEMV 用于計(jì)算向量矩陣乘法.

(7) DGEMM 用于每個進(jìn)程上剩余矩陣數(shù)據(jù)更新.

同時,panel 分解的遞歸部分也會調(diào)用三級BLAS 函數(shù)DGEMM 加速panel 分解過程.因此DGEMM 的效率極大地影響了panel 分解的效率.在異構(gòu)HPL算法中,參數(shù)NBMIN表示遞歸分解算法可以分解到的最小方陣大小,隨著它的增大,一級BLAS 在整體過程中的操作次數(shù)與規(guī)模不變,二級BLAS 操作次數(shù)保持不變,而運(yùn)算規(guī)模將增加,導(dǎo)致panel 分解非遞歸部分的整體時間延長.當(dāng)NBMIN減小時,HPL_pdpanllN[28]函數(shù)的耗時將減少,但是也會增加外層非遞歸部分對三級BLAS 的調(diào)用次數(shù).因此需要通過權(quán)衡三級與二級BLAS 的調(diào)用次數(shù)與效率來調(diào)整NBMIN的大小,獲得更好的性能.另一方面,通過DGEMV 的優(yōu)化可以增大NBMIN來減少三級BLAS 的調(diào)用次數(shù).針對這些函數(shù),本文應(yīng)用多種優(yōu)化策略對BLIS庫中的各級BLAS 函數(shù)進(jìn)行了優(yōu)化,獲得了性能的大幅度提升,提高了異構(gòu)HPL 的效率,更好地發(fā)揮了異構(gòu)系統(tǒng)的性能.

1.2 異構(gòu)單元處理器概述

異構(gòu)計(jì)算單元CPU 是64 位X86 架構(gòu)的通用服務(wù)器芯片,采用14nm 工藝制造.從圖1 中可以簡單地了解CPU 架構(gòu),共32 個計(jì)算核心.其基本的微體系結(jié)構(gòu)如下:

(1) 支持所有標(biāo)準(zhǔn)X86 指令集,如AVX、AVX2、SMEP、SSE 和SSE2 等.

(2) 每4 個核心組成一組,每個核心獨(dú)享L1 和L2 級緩存,共享L3 級緩存,支持SMT 技術(shù),每個物理核心支持同步2 線程.

(3) 三級緩存結(jié)構(gòu)如圖2 所示,L1 cache 由64K 指令緩存和32K 數(shù)據(jù)緩存組成,其中L2 cache 大小為512K,共享L3 cache 大小為8M.所有Cache line 大小為64 字節(jié).

(4) X86 處理器有三級指令TLB 和兩級數(shù)據(jù)TLB.訪存頁面大小為4KB,映射的物理地址空間最大可以達(dá)到1GB.

Fig.2 CPU cache hierarchy圖2 CPU 緩存層次結(jié)構(gòu)

2 訪存優(yōu)化

HPL算法中DGEMM 運(yùn)算占據(jù)了BLAS 運(yùn)算中最大的計(jì)算量和訪存耗時.優(yōu)化的HBLIS 通過實(shí)現(xiàn)互補(bǔ)的兩類GEMM kernels 優(yōu)化了DGEMM 性能,即大kernel 和小kernel.大kernel 用于處理矩陣A和B的規(guī)模都相對較大的情況,而小kernel 用于處理矩陣A更加“高窄”和矩陣B規(guī)模較小的情況.在DGEMM 函數(shù)接口內(nèi)根據(jù)矩陣規(guī)模切換兩類kernel.同時引入Auto-tuning 自適應(yīng)調(diào)優(yōu)技術(shù)對矩陣分塊參數(shù)進(jìn)行優(yōu)化,優(yōu)化后的分塊參數(shù)可以更充分地利用CPU 緩存,提高cache 命中率,從而提高DGEMM 訪存性能.

2.1 矩陣分塊

GEMM 大kernel 采用經(jīng)典高效的GEMM 分塊算法[8],如圖3 左側(cè)的6 層循環(huán)偽代碼所示;GEMM 小kernel的算法如圖3 右側(cè)的偽代碼所示.

Fig.3 High performance GEMM algorithm of big and little kernels圖3 高性能GEMM 的大kernel算法和小kernel算法

大kernel算法描述如下:

? 最外層Loop1 遍歷n維度(以jc索引),以寬度nc為單位,對矩陣C和矩陣B進(jìn)行分塊,形成列向panels.

? 第2 層Loop2 遍歷k維度(以pc索引),以高度kc為單位,將矩陣A劃分成列向panels,同時將Loop1 中劃分出的矩陣B的列向panel 進(jìn)一步劃分成block;Loop2 中,還會將Bblock(B(pc:pc+kc–1,jc+nc–1)表示)進(jìn)行打包(packing).關(guān)于大kernel 中的packing 操作見第2.2 節(jié).

? 第3 層Loop3 遍歷m維度(以ic索引),以高度mc為單位,將Loop1 中劃分出的矩陣C的列向panel 進(jìn)一步劃分出block,以Cc表示(未進(jìn)行packing 操作),將Loop2 中劃分出的矩陣A的列向panel 也進(jìn)一步劃分成block.同樣地,Loop3 中,還會將Ablock 進(jìn)行打包,以Ac表示打包后的矩陣塊.

? 第4 層Loop4 再次遍歷n維度(以jr索引),以寬度nr為單位,將packed blockBc和blockCc再次進(jìn)一步劃分成micro-panels.

? 第5 層Loop5 再次遍歷m維度(以ir索引),以高度mr為單位,將packed blockAc劃分成micro-panels,將Loop4 中Cc中的micro-panel 進(jìn)一步劃分成micro-tiles.

? 最內(nèi)層Loop6 再次遍歷k維度(以pr索引),以1 為單位,計(jì)算Ac中的micro-panel 與Bc中的micro-panel的點(diǎn)積(dot product),并將結(jié)果累加到Cc中的micro-panel,即所謂的“rank-1 update”.

小kernel算法描述如下:

? 前兩次loops 的作用在于將小矩陣B進(jìn)行打包,其目的是便于inner-loops 中最內(nèi)層kernel 對B矩陣的連續(xù)、對齊訪問,關(guān)于小kernel 中的packingB操作見第2.2 節(jié).

? inner-loops 中的3 層循環(huán)類似于大kernel 中的實(shí)現(xiàn),只是此時循環(huán)操作的上限不再是大kernel 中分塊后的nc、mc和kc,而直接是整個矩陣的大小.

與大kernel算法相比,小kernel算法有如下優(yōu)點(diǎn):

(1) 省去了對矩陣A的packing 操作.

異構(gòu)HPL 中的矩陣A往往是“高窄陣”,其規(guī)模比矩陣B大非常多,對矩陣A進(jìn)行打包操作的耗時占整個GEMM 操作的比例非常高.

(2) 極大地避免了大kernel 中out-loops 的循環(huán)次數(shù).

大kernel算法與小kernel算法間的動態(tài)切換,主要依據(jù)輸入矩陣的規(guī)模,即計(jì)算量的大小.但并不存在完美的切換閾值(threshold),因?yàn)? 個矩陣A、B和C的行、列維度千變?nèi)f化,某一個特定的切換閾值無法做到全面適用于所有的使用場景,需要在具體的應(yīng)用中實(shí)際測試來找到比較合適的閾值設(shè)置.根據(jù)本文系統(tǒng)平臺實(shí)測,切換進(jìn)入小kernel 的閾值條件滿足下面兩者之一即可:

(1)C矩陣的面積,即M×N≤700×700 時.

(2) 維數(shù)M≤160 同時維數(shù)K≤128 時.

2.2 Auto-tuning自適應(yīng)分塊參數(shù)調(diào)優(yōu)

DGEMM 計(jì)算量大,計(jì)算訪存比高,除上述kernel算法層面的優(yōu)化外,對BLIS 框架內(nèi)矩陣分塊參數(shù)進(jìn)行優(yōu)化還可以提高cache 命中率,提高性能.矩陣乘法循環(huán)分塊的大小是通過BLIS庫中DGEMM 的分塊參數(shù)來決定的,分塊參數(shù)的選擇可以影響cache 命中率,不合適的參數(shù)可能導(dǎo)致流水線停頓,從而使DGEMM 性能下降.因此,DGEMM 分塊參數(shù)的選擇至關(guān)重要.BLIS 中DGEMM 分塊是通過開發(fā)人員的經(jīng)驗(yàn)和測試選擇的一組數(shù)值,手動測試不能保證測試覆蓋率,也不能保證得到的是最優(yōu)解.因此本文采用了Auto-tuning 遺傳算法自適應(yīng)調(diào)優(yōu)技術(shù),更高效地優(yōu)化DGEMM 分塊參數(shù).

遺傳算法基本算法思想是從可能具有潛在解集的一個種群開始的,種群是基因編碼后形成的一組個體集合,在這個種群基礎(chǔ)上通過選擇(selection)、交叉(crossover)和變異(mutation)操作,優(yōu)勝劣汰,不斷地進(jìn)化出最優(yōu)的個體.Auto-tuning 自適應(yīng)調(diào)優(yōu)就是應(yīng)用遺傳算法對DGEMM 分塊參數(shù)進(jìn)行智能優(yōu)化選擇,將分塊參數(shù)作為個體,一定范圍內(nèi)的多個個體形成一個種群進(jìn)行自適應(yīng)迭代進(jìn)化,從而針對具體平臺體系結(jié)構(gòu)產(chǎn)生最優(yōu)的DGEMM 矩陣分塊參數(shù).文獻(xiàn)[29]在Intel 平臺上對遺傳算法在DGEMM 的分塊參數(shù)智能優(yōu)化選擇上的實(shí)現(xiàn)做了詳細(xì)介紹.如圖4 所示是遺傳算法的基本流程.

? 首先初始化種群,以分塊參數(shù)[mc,kc,nc]為一個個體,隨機(jī)生成一組個體的集合.

? 然后計(jì)算種群每個個體的適應(yīng)度,適應(yīng)度大小決定遺傳概率.

? 接著判斷迭代次數(shù),如果沒有達(dá)到預(yù)設(shè)的最大值,則進(jìn)行selection、crossover 和mutation 的遺傳算法進(jìn)化步驟.

? 最后,經(jīng)過selection、crossover 和mutation 這3 個步驟后生成一個新的種群,與原始種群大小一致,但個體都經(jīng)過優(yōu)勝劣汰,進(jìn)化到了更好的水平.重復(fù)上述步驟,直到滿足終止條件退出.

Fig.4 Generic algorithm flow chart圖4 遺傳算法流程圖

根據(jù)上述遺傳算法的介紹,自適應(yīng)程序的實(shí)現(xiàn)需要經(jīng)過編碼個體、建立種群、計(jì)算適應(yīng)度、選擇(selection)、交叉(crossover)、變異(mutation)和終止條件判斷這幾個步驟.

(1) 個體編碼

編碼方法有二進(jìn)制編碼、格雷編碼、浮點(diǎn)數(shù)編碼和符號編碼等多種方式,根據(jù)DGEMM 分塊參數(shù)特點(diǎn),采取浮點(diǎn)數(shù)編碼更適合.也就是將分塊參數(shù)[mc,kc,nc]這3 個浮點(diǎn)數(shù)直接作為參數(shù)傳入.

(2) 建立種群(population)

個體集合的建立過程就是要在一個確定范圍內(nèi)去隨機(jī)生成一定數(shù)量的個體.種群范圍的確定可以根據(jù)cache 的大小進(jìn)行理論計(jì)算得出[7],HBLIS 分塊參數(shù)原始值是[1080,120,8600],在原始參數(shù)值的基礎(chǔ)上擴(kuò)大范圍至[512

(3) 計(jì)算適應(yīng)度(fitness)

適應(yīng)度函數(shù)就是DGEMM 測試程序,其性能作為判斷適應(yīng)度高低的標(biāo)準(zhǔn).性能高的分塊參數(shù)會遺傳給下一代種群,性能低的分塊參數(shù)就會被淘汰.

(4) 選擇(selection)

基本遺傳算法僅僅使用selection、crossover 和mutation 這3 種算子是沒有辦法收斂于全局最優(yōu)解的,因?yàn)楹唵蔚剡M(jìn)行雜交,可能反而會把較好的組合破壞掉.本文采用精英保留策略,也就是把目前最優(yōu)個體不進(jìn)行任何操作直接復(fù)制到下一代種群中.選擇方法有很多種,包括比例選擇法、排序選擇等.本文采用比例選擇法和精英保留策略結(jié)合,種群中每一個個體的適應(yīng)度占總適應(yīng)度的比例越大,其被選擇的概率就越大,同時這一代種群最優(yōu)秀個體無條件復(fù)制到下一代.

(5) 交叉(crossover)

單點(diǎn)交叉、雙點(diǎn)交叉、均勻交叉等都適合二進(jìn)制編碼和浮點(diǎn)數(shù)編碼,根據(jù)GEMM 分塊參數(shù)只有3 個,選擇單點(diǎn)交叉,按某一概率將兩個個體的某一位置參數(shù)互換.交叉率設(shè)置越大,新個體產(chǎn)生越快,其全局搜索能力就會越強(qiáng),但是優(yōu)秀個體被破壞的可能性也就越高,綜合考慮,設(shè)置交叉概率為0.8.

(6) 變異(mutation)

變異的主要目的是改善局部搜索能力以及維持種群個體的多樣性,變異方法采用均勻變異,按某一概率在一個范圍內(nèi)隨機(jī)生成一個數(shù)替換原來個體中某一位置的參數(shù).變異概率不宜設(shè)置過大,導(dǎo)致丟失優(yōu)秀個體而變成隨機(jī)搜索.考慮設(shè)置參數(shù)范圍已經(jīng)在合理范圍內(nèi),本文按0.1 的概率對每一個個體的某一參數(shù)進(jìn)行變異.

(7) 終止條件判斷(iteration)

種群中所有的個體適應(yīng)度計(jì)算完成后要判斷設(shè)置的終止條件是否滿足,以判斷程序是否結(jié)束退出.如果不滿足終止條件就進(jìn)行selection、crossover 和mutation 操作,形成下一個同等規(guī)模的個體集合(種群).

Auto-tuning 自適應(yīng)調(diào)優(yōu)算法完成后,需要進(jìn)行DGEMM 測試,BLIS庫中config 目錄是存放不同平臺kernel的配置信息,其中DGEMM 分塊參數(shù)通過bli_kernel.h 頭文件中的BLIS_DEFAULT_MC_D、BLIS_DEFAULT_KC_D 和BLIS_DEFAULT_NC_D 這3 個宏定義配置.

Auto-tuning 的具體測試流程如下:

(1) 在包含初始分塊的更大的范圍內(nèi)隨機(jī)生成一組滿足遺傳算法個體要求的分塊參數(shù)[mc,kc,nc],并用其替換掉blis_kernel.h 中的分塊參數(shù).

(2) 自動編譯BLIS庫,生成此分塊參數(shù)下的動態(tài)鏈接庫libblis.so.

(3) 調(diào)用libblis.so 執(zhí)行測試程序DGEMM,主程序?qū)⑿阅苤涤涗洸⒎祷?

通過Auto-tuning 的不斷迭代,根據(jù)算法模型,我們會得到一個全局最優(yōu)解,這組參數(shù)就是性能更高的DGEMM分塊參數(shù),實(shí)測數(shù)據(jù)分塊參數(shù)由[1080,120,8400]優(yōu)化為[792,822,8628],將這組參數(shù)替換進(jìn)BLIS中,并使用高性能計(jì)算標(biāo)準(zhǔn)測試程序HPL 和單獨(dú)的DGEMM 測試程序進(jìn)行測試,具體測試對比數(shù)據(jù)參考第5 節(jié)性能分析.

2.3 矩陣打包

連續(xù)的內(nèi)存訪問比非連續(xù)的內(nèi)存訪問要快.異構(gòu)HPL 中傳給DGEMM 的矩陣A、B和C的索引是初始化矩陣的某個位置指針,矩陣元素通過leading-dimension 索引,leading-dimension 表示列優(yōu)先存儲方式的矩陣同一行(或行優(yōu)先存儲方式的矩陣同一列)中兩相鄰元素的距離,也就是實(shí)際存儲矩陣的二維數(shù)組的第1 維(或第2維)的大小.矩陣A、B和C的行或列間跨度通過LDA、LDB和LDC表示,這3 個數(shù)值可能會非常大.比如異構(gòu)HPL 求解問題Ax=b,其矩陣規(guī)模N的最大取值可以使矩陣A的存儲空間占用內(nèi)存80%左右,則LDA不能小于N.在算法實(shí)現(xiàn)中,大kernel算法會不斷地對A和B的不同分塊進(jìn)行打包,形成Ac和Bc;而小kernel算法只將原始矩陣B進(jìn)行一次打包,形成行優(yōu)先存儲(row-major)的Bpacked.如圖5 所示,其中假設(shè)nr等于6,矩陣B是以列優(yōu)先存儲(column-major)方式存儲數(shù)據(jù),且行數(shù)k和列數(shù)n都是nr的整數(shù)倍.

Fig.5 Packing matrix B in little GEMM kernel圖5 GEMM 小kernel 中矩陣B 的打包操作

Packing 數(shù)據(jù)打包操作的優(yōu)點(diǎn)有3 個.

(1) 創(chuàng)造連續(xù)的內(nèi)存訪問方式.

BLAS 接口傳入的矩陣存儲方式是不確定的,可能是row-major,也可能是column-major,且矩陣的行或列間的跨度會非常大,最內(nèi)層的kernel 為了適應(yīng)不同的數(shù)據(jù)存儲方式以及不連續(xù)的數(shù)據(jù)分布,需要通過packing 制造統(tǒng)一的、連續(xù)的內(nèi)存布局視圖.

(2) 數(shù)據(jù)對齊.

對于大多數(shù)CPU 微架構(gòu)實(shí)現(xiàn)而言,對齊到寄存器寬度,cache 邊界和頁邊界的數(shù)據(jù)訪問往往比非對齊的數(shù)據(jù)訪問方式更快.

(3) 提前l(fā)oad 數(shù)據(jù)到cache.

通過packing 操作,最內(nèi)層kernel 可以在盡可能接近CPU 寄存器的cache 上獲取到需要進(jìn)行運(yùn)算的操作數(shù),這比直接從內(nèi)存中獲取數(shù)據(jù)要快很多.

2.4 數(shù)據(jù)預(yù)取

X86 CPU 提供了硬件取功能和符合X86 ISA 的軟件預(yù)取指令.其中硬件預(yù)取對開發(fā)者是透明的,軟件預(yù)取指令PREFETCHn 指示處理器將后續(xù)需要訪問的數(shù)據(jù)提前讀入到指定的某級緩存中,當(dāng)程序使用這些數(shù)據(jù)時,可以直接讀取緩存.PREFETCHn 指令格式見表1.

Table 1 X86 software prefetch instructions表1 X86 軟件預(yù)取指令

軟件預(yù)取指令不要求數(shù)據(jù)行對齊,如果預(yù)取數(shù)據(jù)行已經(jīng)存在于比指定的n更低的某層cache 中,那么指令將被當(dāng)作NOP 空指令處理.

在最內(nèi)層的核心匯編代碼中,PREFETCHn 指令的使用主要考慮兩點(diǎn):

(1) 何處插入PREFETCHn 指令.

(2) 預(yù)取的地址跨度設(shè)置.

X86 CPU 的cache line 大小為64B,一個緩存行可以存放8 個雙精度浮點(diǎn)數(shù)據(jù).如圖6 所示的匯編代碼第1行的軟件預(yù)取操作,表示把Ac的第9 個雙精度數(shù)據(jù)開始的一個緩存行預(yù)取到L1 緩存中.此處的預(yù)取操作是為第22 行的load 指令準(zhǔn)備數(shù)據(jù).第5 行的load 操作會在L1 緩存發(fā)生cache miss,但是因?yàn)閷進(jìn)行了packing 操作,理論上能夠在L2 cache 命中.同理,第9 行、第14 行和第18 行處的load 操作會在L1 中命中,第11 行預(yù)取操作將Ac的第17 個數(shù)據(jù)預(yù)取到內(nèi)存.

Fig.6 Code segment of inner-most DGEMM kernel圖6 DGEMM 核心匯編代碼片段

軟件預(yù)取指令的插入位置需要綜合考慮體系結(jié)構(gòu)特點(diǎn),合理安排匯編指令.X86 CPU 的L2 load-to-use 的latency 至少是12 個core cycles,而每時鐘周期只可以執(zhí)行一條FMA 指令,在CPU 流水線能夠連續(xù)retire FMA指令的理想狀態(tài)下,預(yù)取指令和后續(xù)load 指令之間至少插入12 條FMA 指令.

3 指令集優(yōu)化

高效BLAS kernel 的實(shí)現(xiàn)方法具有通用性,不同于其他數(shù)學(xué)庫大部分代碼通過匯編實(shí)現(xiàn)優(yōu)化,BLIS 只在最內(nèi)層的kernel 核心代碼處使用了匯編代碼,框架內(nèi)的其他代碼使用C 實(shí)現(xiàn).這給可移植性帶來了極大的便利,在適配不同架構(gòu)的CPU 時,開發(fā)者可以將主要精力集中在最內(nèi)層kernel 的匯編代碼的優(yōu)化上.

3.1 X86向量指令

不同架構(gòu)的CPU 大多會提供高效的SIMD(single instruction multiple data)指令,X86 CPU 支持AVX2 指令集,可以使用128 位的XMM 寄存器和256 位的YMM 寄存器,但由于浮點(diǎn)運(yùn)算部件的數(shù)據(jù)路徑寬度是128 位,每個時鐘周期最多可以同時處理2 個雙精度浮點(diǎn)數(shù),單核理論雙精度浮點(diǎn)性能是(128/64)×2=4 DP-FLOPs/cycle,實(shí)際測試結(jié)果表明,使用XMM 寄存器的效率要優(yōu)于YMM 寄存器.AVX2 指令集中主要有兩類指令:數(shù)據(jù)傳輸類指令和數(shù)據(jù)操作類指令.其中,VMOVDDUP、VMOVAPD、VMOVUPD、VMULPD 和VFMADD231PD 是在X86 架構(gòu)CPU 平臺上高效BLAS kernels 的實(shí)現(xiàn)中使用頻率非常高的幾條指令.

HPL算法中傳遞給BLAS 函數(shù)的參數(shù)會有固定值的情況,可以通過適當(dāng)?shù)闹噶钐鎿Q來精簡指令,舉例如下.

(1) 對于Level-3 級函數(shù)DGEMM

系數(shù)Alpha 固定為–1.0,Beta 固定為1.0.此時的DGEMM 計(jì)算公式由C:=α·AB+β·C變成了C:=C?AB可以在packingB(見第2.3 節(jié))的過程中使用xorpd 指令進(jìn)行符號位翻轉(zhuǎn),精簡了對系數(shù)Alpha 的乘法操作.矩陣C不變,精簡了對系數(shù)Beta 的乘法操作.精簡指令后,原來需要4 個core cycle 的乘法操作和5 個core cycle 的乘加指令只需3 個core cycle 的加法指令即可完成.在矩陣規(guī)模較大,最內(nèi)層kernel 被循環(huán)多次的情況下,指令時鐘周期的減少可以優(yōu)化整體性能.

(2) 對于Level-2 級函數(shù)DGEMV

轉(zhuǎn)置參數(shù)TransA 固定為NoTrans,可以省略對矩陣A進(jìn)行轉(zhuǎn)置操作的代碼,提高CPU 流水線分支預(yù)測的準(zhǔn)確率,避免分支預(yù)測失敗的懲罰.Alpha 固定為–1.0,Beta 固定為 1.0,則可以在內(nèi)層 kernel 中使用一條VFNMADD231PD 指令取代VMULPD 和VFMADD231PD 兩條指令.

(3) 對于Level-1 級函數(shù)DSCAL

IncX 固定為1,可以省去對數(shù)據(jù)成員間隔較大時所需的軟件預(yù)取操作的指令(見第2.4 節(jié)).

異構(gòu)HPL 傳遞給BLAS 接口的參數(shù)還有其他的固定值情況存在,這些固定值的存在可以起到精簡指令的作用,為BLAS庫提供了優(yōu)化空間.

3.2 循環(huán)展開

循環(huán)展開是編譯器常用的優(yōu)化方式之一,可以減少分支預(yù)測失敗帶來的開銷,提高指令級的并行度.但高級語言通過編譯器循環(huán)展開后的匯編代碼在BLAS 這種訪存密集型的程序特征下,往往性能表現(xiàn)較差.因此BLAS的核心循環(huán)代碼通常采用手寫匯編的方式進(jìn)行展開.

手寫B(tài)LAS 核心匯編代碼,不僅需要考慮底層CPU 架構(gòu)層面的架構(gòu)寄存器的數(shù)量和復(fù)用方式,而且需要考慮底層CPU 的微架構(gòu)特點(diǎn),比如浮點(diǎn)部件的寬度、存儲層級的大小和層級間的延遲等.這些CPU 特性共同決定了手寫匯編循環(huán)展開的程度,也就是“循環(huán)展開因子”.本文第2.1 節(jié)中提到的5 個參數(shù)nc、kc、mc、nr和mr,其中內(nèi)核參數(shù)nr和mr與匯編代碼循環(huán)展開直接相關(guān),而矩陣分塊參數(shù)nc,kc和mc則與CPU 存儲層級特性直接相關(guān).針對具體平臺體系結(jié)構(gòu)可以理論計(jì)算出這5 個參數(shù)的大致取值[30],但實(shí)際代碼測試中發(fā)現(xiàn),理論數(shù)據(jù)往往并不是最優(yōu)設(shè)置,X86 CPU 有16 個128 位的向量寄存器XMM0～XMM15,內(nèi)核參數(shù)nr,mr取值的首要策略是最大化利用架構(gòu)寄存器,因此綜合考慮,nr和mr分別取值為6 和4.同時矩陣參數(shù)[nc,kc,mc]依據(jù)第2.2 節(jié)中的自動調(diào)優(yōu)技術(shù)同樣獲得了比理論值更好的參數(shù)設(shè)置.

如圖6 所示,通用寄存器RAX 用于索引Ac的micro-panel 中的數(shù)據(jù),不斷獲取一個雙精度數(shù)據(jù)然后復(fù)制成兩份填充到XMM0 中;RBX 用于索引Bc的micro-panel 中的數(shù)據(jù),一次獲取并存儲6 個雙精度數(shù)據(jù)到XMM1～XMM3 中.XMM4～XMM15 用于存儲不斷累加的中間結(jié)果.在此過程中,向量寄存器XMM0 實(shí)際共存儲過4 個Ac數(shù)據(jù),XMM1～XMM3 共存儲過6 個Bc數(shù)據(jù).

4 多線程并行

常用的多線程并行手段有OpenMP[31]和POSIX Threads(pthread)[32].BLIS 同時支持兩種并行方式,且當(dāng)前只在Level-3 實(shí)現(xiàn)了并行.與OpenMP 相比,Pthreads 技術(shù)在實(shí)際編程中需要考慮臨界區(qū)、線程同步原語等非常底層的操作,故我們在GEMM 以及TRSM 等小kernel 和Level-1 與Level-2 的并行化實(shí)現(xiàn)中,選擇使用更加便利的OpenMP 進(jìn)行并行化.

4.1 Control-tree優(yōu)化

針對Level-3 級BLAS 的并行實(shí)現(xiàn)方式,BLIS 提出了新穎的“control-tree”結(jié)構(gòu)[7].該結(jié)構(gòu)以統(tǒng)一的方式實(shí)現(xiàn)了如圖3 所示的GEMM 大kernel算法的并行以及HERK、TRMM 和TRSM 運(yùn)算的并行.Control-tree 結(jié)構(gòu)最突出的優(yōu)點(diǎn)是可以在任意層loop 內(nèi)分別使用環(huán)境變量BLIS_JC_NT、BLIS_IC_NT、BLIS_JR_NT和BLIS_IR_NT進(jìn)行多線程并行的控制,甚至是多個維度同時并行.這種統(tǒng)一結(jié)構(gòu)的便利性的代價就是在所有能夠發(fā)生并行的地方設(shè)置同步屏障(barrier),以便在當(dāng)前l(fā)oop 的下層所有操作全部返回后才能繼續(xù)回溯.雖然文獻(xiàn)[7]中的性能測試結(jié)果表明control-tree 結(jié)構(gòu)的開銷是非常小的,但在異構(gòu)HPL 環(huán)境中,control-tree 的開銷變得尤為突出.BLIS原生control-tree 的結(jié)構(gòu)大致如圖7 所示.

Fig.7 Original BLIS GEMM control-tree圖7 BLIS 原生GEMM control-tree

對于異構(gòu)HPL 中的GEMM 運(yùn)算,我們選擇只在loop3 處(Mc維度)實(shí)現(xiàn)并行,當(dāng)各個線程的Acpacking 打包操作和后續(xù)操作都完成后,需要在loop3 的Icbarrier 處進(jìn)行第1 次同步,隨后的Pcbarrier 和global barrier 處再次進(jìn)行同步,才可以保證最終GEMM 運(yùn)算的正確性.可以看到,Pcbarrier 和global barrier 的同步開銷其實(shí)是可以避免的,只需將loop3 放到算法的最前端,然后使用OpenMP 自身的同步操作實(shí)現(xiàn)各個線程的同步即可.基于此分析,對GEMM 運(yùn)算的control-tree 結(jié)構(gòu)做出了如圖8 所示的調(diào)整.

通過異構(gòu)HPL 的實(shí)際測試,未優(yōu)化control-tree 時DGEMM 和DTRSM 運(yùn)算的同步barrier 等待的時間開銷為2.29s,約占CPU端計(jì)算總時間的20%.優(yōu)化control-tree 結(jié)構(gòu)后,GEMM 函數(shù)同步barrier 的等待時間減少為1.49s,TRSM 運(yùn)算應(yīng)用類似于優(yōu)化思路,同步barrier 的等待時間減少到基本可以忽略不計(jì).

Fig.8 Optimized HBLIS GEMM control-tree圖8 優(yōu)化后的HBLIS GEMM control-tree

4.2 Level-1和Level-2級BLAS并行優(yōu)化

原BLIS庫沒有實(shí)現(xiàn)Level-1 級和Level-2 級BLAS 函數(shù)運(yùn)算的并行化.本文使用如下相似的策略實(shí)現(xiàn)了DGEMV、DTRSV、DSCAL、DAXPY 以及IDAMAX 運(yùn)算的并行化:

(1) 通過運(yùn)行腳本中設(shè)置的環(huán)境變量“OMP_NUM_THREADS”獲取當(dāng)前MPI 進(jìn)程中的線程數(shù).

(2) 將矩陣或向量中的某個維度按設(shè)置的線程數(shù)進(jìn)行切分.

(3) 最后使用OpenMP 的Pragma(編譯指示)啟動多線程,在不同的線程空間中復(fù)用同一段匯編kernel 來處理不同的數(shù)據(jù).在Pragma 包裹的并行域(parallel region)中需要合理而準(zhǔn)確地為每個線程劃分任務(wù)空間.可以通過omp_get_thread_num()函數(shù)獲得當(dāng)前線程的線程編號(TID)劃分任務(wù),避免計(jì)算沖突.

多線程并行對于BLAS 運(yùn)算的性能提升是非常顯著的,矩陣或向量計(jì)算達(dá)到一定規(guī)模后,性能往往隨著線程數(shù)的增加而提升,但提升并不會完全是線性的,通常是計(jì)算越密集的運(yùn)算,提升的線性程度越高,因此計(jì)算量最大的Level-3 級BLAS 運(yùn)算通過多線程并行會大大提升性能.如圖9 所示,在N=43776 時的Level-1 運(yùn)算IDMAX 隨著線程數(shù)的增加表現(xiàn)出性能的提升.

Fig.9 Performance improvement by multi-threading of IDAMAX圖9 多線程對IDAMAX 的性能提升

4.3 并行化問題

多線程并行化并不是總會對計(jì)算帶來幫助,為了進(jìn)一步提高性能,需要綜合考慮以下問題.

(1) 動態(tài)多線程

創(chuàng)建和啟動多線程的開銷在計(jì)算量很小的情況下會變得突出,從而導(dǎo)致程序性能下降.所以需要設(shè)置一個或多個特定的閾值來進(jìn)行線程數(shù)的切換.具體是設(shè)置一個閾值將多線程直接切換到單線程,還是設(shè)置多個閾值將線程數(shù)分段遞減,需要具體的BLAS 函數(shù)結(jié)合特定的輸入數(shù)據(jù)進(jìn)行大量的測試歸納(profiling).在DGEMM 的多線程優(yōu)化中,我們設(shè)計(jì)了精確的運(yùn)算量模型.DGEMM 需要計(jì)算一個m×k階矩陣A和k×n階矩陣B的乘積,其主要計(jì)算量和乘積mkn成正比.考慮到矩陣A的規(guī)模m·k,矩陣B的規(guī)模k·n等都對最終運(yùn)算性能造成影響,我們使用最小二乘法擬合矩陣計(jì)算的性能公式,假設(shè)不同線程數(shù)目下運(yùn)算時間根據(jù)式(2)計(jì)算.

于是對于輸入?yún)?shù)為mi、ki和ni,運(yùn)行時間為ti的測試數(shù)據(jù),對應(yīng)擬合數(shù)據(jù)中的一行如式(3)所示.

所有這些擬合數(shù)據(jù)組成N×8 階矩陣P,其中N是采樣數(shù),而所有的采集時間構(gòu)成對應(yīng)的列向量b,于是構(gòu)成式(4)的擬合方程.

解方程得出的8 階變量x的各分量就是a1,a2,…,a8的最小二乘擬合結(jié)果.分別對單線程、4 線程、8 線程和32 線程測試擬合公式,在運(yùn)行中根據(jù)擬合結(jié)果動態(tài)切換線程數(shù)據(jù),得到如圖10 所示的性能結(jié)果.

Fig.10 Performance speedup of DGEMM in 32 threads圖10 DGEMM 在32 個線程中的性能加速比

實(shí)踐中還發(fā)現(xiàn),不同的kernel 實(shí)現(xiàn)方式以及硬件設(shè)置等因素也會影響到閾值的選取.比如,在異構(gòu)HPL 的DGEMV 函數(shù)的優(yōu)化過程中,根據(jù)不同的“OMP_NUM_THREADS”、不同的CPU 主頻以及矩陣A的列數(shù)N(M方向做并行)設(shè)置不同的多個閾值.

(2) 均衡余量(remainder)處理

當(dāng)使用多線程將求解問題切分成多塊后,往往會遇到不能均分的情況.所謂不能“均分”,是指整體的運(yùn)算量不能均勻地分配到每個線程中去.此時有兩個選擇:或者將所有多出的運(yùn)算量全部放到某一個線程中,比如最后一個線程;或者將余量再次均分到所有線程中.Remainder 處理的問題在線程數(shù)較少時(比如6 個線程)不會那么關(guān)鍵,但在線程數(shù)較多時(比如32 個線程)會變得突出.因?yàn)榇藭r的木桶效應(yīng)更加明顯,整體運(yùn)算的同步結(jié)束取決于擁有最大余量的那個線程.根據(jù)本文HPL算法測試分析后,采取了將remainder 分配給每個線程處理,得到的性能結(jié)果最優(yōu).

(3) 循環(huán)展開與多線程的均衡

循環(huán)展開的程度和線程的數(shù)量需要均衡考慮.在使用動態(tài)多線程去加速運(yùn)算的過程中,往往會有多種求解方法去完成同一個運(yùn)算.考慮以下場景:異構(gòu)HPL 中的DGEMM 運(yùn)算在M較小時,比如M=28,此時可以使用4個線程各處理6 個M維度的行,然后使用1 個線程處理剩下的4 個M行;或者使用5 個線程各處理5 個M維度的行,然后使用3 個線程處理各處理一個M行,如式(5)所示.

不同的求解方法涉及同一個函數(shù)的不同kernel 函數(shù)編寫以及動態(tài)多線程threshold 的選取.通過實(shí)際測試發(fā)現(xiàn):DGEMM 運(yùn)算在M較小時,更多的循環(huán)展開比更多的線程性能要好,如圖11 所示.

Fig.11 In the case that B matrix is square matrix (dimension N=K),more unrolling is better performance than more threading in small DGEMM圖11 B 矩陣是方陣(維數(shù)N=K)的情況下,在小規(guī)模的DGEMM 中,更多的循環(huán)展開比更多的線程性能更優(yōu)

5 性能測試與分析

異構(gòu)HPL算法中綜合使用了以上介紹的訪存優(yōu)化、指令集優(yōu)化和多線程并行等優(yōu)化技術(shù),對算法中調(diào)用的7 個BLAS 函數(shù)進(jìn)行了針對性優(yōu)化,并與OpenBLAS 0.3.6 和MKL(Intel Math Kernel Library 2018.0.128)進(jìn)行了對比測試.實(shí)驗(yàn)測試平臺是高性能計(jì)算集群系統(tǒng)中的單個計(jì)算節(jié)點(diǎn),其配置是如圖1 所示的32 核X86 CPU 與4 個協(xié)處理組成的異構(gòu)計(jì)算單元.在測試平臺上分別進(jìn)行了單核DGEMM 和多核心并行DGEMM 測試,最后測試了異構(gòu)HPL 程序.

5.1 Cache性能分析

通過系統(tǒng)性能分析工具讀取CPU 中各種計(jì)數(shù)器的數(shù)值,可以分析應(yīng)用程序的各項(xiàng)指標(biāo).X86 處理器支持非常多的性能監(jiān)控計(jì)數(shù)器,包括浮點(diǎn)、訪存、指令等CPU 各種資源的監(jiān)控.如表2 所示為Auto-tuning 優(yōu)化后的矩陣分塊參數(shù)與初始分塊參數(shù)進(jìn)行DGEMM 計(jì)算時的cache miss 對比情況,分析得出分塊參數(shù)[mc,kc,nc]從[1080,120,8400]優(yōu)化為[792,822,8628],kc和nc的增大導(dǎo)致矩陣打包到L1 緩存時,因L1 容量小而發(fā)生更多的cache miss.但是,L2 和L3 緩存容量比L1 大得多,其命中率的提高可以掩蓋L1 的性能損失并提高整體DGEMM的性能.

Table 2 Cache miss rate of HPL counted by performance monitor unit表2 通過性能監(jiān)控單元統(tǒng)計(jì)HPL 的緩存失效率

5.2 DGEMM效率

如圖12 和圖13 所示分別為單核和32 核情況下不同數(shù)學(xué)庫的DGEMM 性能對比.因?yàn)镮ntel MKL 只有在Intel 平臺上才能發(fā)揮最佳效果,在X86 異構(gòu)平臺上效率不是最好,僅作為參考;而Netlib 為未經(jīng)過任何優(yōu)化的Fortran 語言實(shí)現(xiàn)的通用標(biāo)準(zhǔn)數(shù)學(xué)庫.可以看到,HBLIS 的DGEMM 實(shí)現(xiàn)無論是在單核心還是在32 核心下,其性能均超過了其他BLAS庫,單核DGEMM 最高效率達(dá)到了98%,32 核最高效率達(dá)到了96.9%.與未優(yōu)化前的BLIS和MKL 相比,單核性能分別提升了19.2%和33.3%,32 核性能分別提升了42.1%和33.6%.

Fig.12 DGEMM efficiency comparison of single-core圖12 單核DGEMM 效率對比

Fig.13 DGEMM efficiency comparison of 32 cores圖13 32 核DGEMM 效率對比

5.3 異構(gòu)HPL測試

HPL 實(shí)測性能值是最終評價標(biāo)準(zhǔn),分別統(tǒng)計(jì)每個BLAS 函數(shù)的時間作為BLAS 優(yōu)化的評價標(biāo)準(zhǔn).為避免網(wǎng)絡(luò)通信的干擾,本文基于單節(jié)點(diǎn)進(jìn)行異構(gòu)HPL 測試,對比包括panel 分解中pfact 時間以及HPL算法中使用的7 個BLAS 函數(shù)DGEMM、DGEMV、DTRSM、DTRSV、IDMAX、DSCAL 和DCOPY 的計(jì)算時間.

異構(gòu)HPL算法panel 分解部分主要在CPU端進(jìn)行計(jì)算.如圖14 左側(cè)所示為HPL 程序panel 分解中的pfact總時間統(tǒng)計(jì),調(diào)用HBLIS 的pfact 與調(diào)用MKL 的相比,計(jì)算時間節(jié)省了53%.如圖14 右側(cè)所示為DGEMM 函數(shù)計(jì)算時間統(tǒng)計(jì),可以看到,HBLIS 的DGEMM 性能比MKL 的DGEMM 性能提升了約25%.

Fig.14 Performance comparison of pfact and DGEMM calling different math libraries in heterogeneous HPL圖14 異構(gòu)HPL 中調(diào)用不同數(shù)學(xué)庫的pfact 和DGEMM 的性能對比

如圖15 所示為DGEMV 和DTRSV 函數(shù)性能對比.DGEMV 在panel 分解中時間占比僅次于DGEMM,因此DGEMV 性能對整體性能的影響也僅次于DGEMM.圖15 左側(cè)為DGEMV 性能,優(yōu)化的HBLIS 比MKL 提升了約62%,比原BLIS 提升了約65%.圖15 右側(cè)為DTRSV 函數(shù)時間統(tǒng)計(jì),原BLIS 性能與MKL 差距較大,經(jīng)過優(yōu)化后的BLIS 性能超過了MKL.因?yàn)镈TRSV 計(jì)算量較小,在總時間中占比也相對很小.

Fig.15 Performance comparison of DGEMV and DTRSV calling different math libraries in heterogeneous HPL圖15 異構(gòu)HPL 中調(diào)用不同數(shù)學(xué)庫的DGEMV 和DTRSV 的性能對比

如圖16 所示為DTRSM 和IDMAX 函數(shù)性能對比.原生BLIS算法庫中的三角矩陣求解DTRSM 性能較差,優(yōu)化后的DTRSM 性能與MKL 基本一致.在pfact 第106 步～第141 步的性能是超過MKL 的,而在panel 分解最后段,MKL 運(yùn)算更快,原因是DTRSM 在切換到小kernel 的threshold 和動態(tài)多線程的threshold 時,需要做進(jìn)一步的調(diào)整,以適應(yīng)運(yùn)算量的變化.IDMAX 函數(shù)主要是用于選主元操作,計(jì)算量可以忽略不計(jì).

Fig.16 Performance comparison of DTRSM and IDAMAX calling different math libraries in heterogeneous HPL圖16 異構(gòu)HPL 中調(diào)用不同數(shù)學(xué)庫的DTRSM 和IDAMAX 的性能對比

如圖17 所示,Level-1 級BLAS 函數(shù)DSCAL 和DCOPY 優(yōu)化后的性能基本與Intel MKL 一致.其中,DSCAL優(yōu)化后的性能比原生BLIS 有大幅度提升,但在panel 分解后段的性能還是略低于MKL,原因是MKL 使用的多線程庫IOMP 比開源的GOMP 在多線程啟動、調(diào)度和動態(tài)多線程的切換上更加高效.而DCOPY 代碼實(shí)現(xiàn)簡單,優(yōu)化空間非常有限,硬件平臺性能直接決定了其時間開銷,故可以看到三者DCOPY 的性能基本一致.DAXPY 運(yùn)算只在回代(back substitution)過程中調(diào)用一次,性能影響可以忽略,故未列出詳細(xì)對比.

Fig.17 Performance comparison of DSCAL and DCOPY calling different math libraries in heterogeneous HPL圖17 異構(gòu)HPL 中調(diào)用不同數(shù)學(xué)庫的DSCAL 和DCOPY 的性能對比

6 總結(jié)與展望

本文針對CPU 平臺系統(tǒng)架構(gòu)的特點(diǎn),在開源BLIS算法庫的基礎(chǔ)上開發(fā)了優(yōu)化的HBLIS 數(shù)學(xué)庫,詳細(xì)介紹了HPL 調(diào)用的各級BLAS 函數(shù)優(yōu)化方法.通過訪存優(yōu)化、指令集優(yōu)化和多線程并行優(yōu)化,HPL 調(diào)用的各級BLAS函數(shù)的計(jì)算性能與MKL 的相應(yīng)函數(shù)相比有著普遍和顯著的提升,個別計(jì)算量小的函數(shù)性能基本與MKL 相應(yīng)函數(shù)性能一致,其中最重要的DGEMM 函數(shù)提升了約25%,DGEMV 函數(shù)性能提升了約62%.在異構(gòu)單節(jié)點(diǎn)HPL測試中,與調(diào)用MKL庫的HPL 相比,采用HBLIS 的HPL 整體效率值提升了11.8%,更好地發(fā)揮了異構(gòu)高性能平臺的計(jì)算能力.同時本文也有不足之處,雖然HBLIS庫性能測試發(fā)現(xiàn)DGEMM、DGEMV、DTRSV 和IDMAX有著很好的優(yōu)化效果,其性能與MKL 和開源BLIS庫相應(yīng)函數(shù)相比都有大幅度提升,但是DTRSM 和DSCAL優(yōu)化后性能并沒有優(yōu)于Intel MKL,這可能與計(jì)算量、多線程庫支持和編譯器優(yōu)化支持都有密切關(guān)聯(lián).底層子例程庫優(yōu)化是一個系統(tǒng)工程,不僅需要硬件架構(gòu)的優(yōu)化設(shè)計(jì),而且需要編譯器和線程庫等底層軟件的優(yōu)化.本文未涉及的其他BLAS 函數(shù)的優(yōu)化將是HBLIS庫下一步優(yōu)化工作的方向.