武 錚,許 樂(lè),安 虹,金 旭,文 可

(中國(guó)科學(xué)技術(shù)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,合肥 230027) E-mail： zhengwu@mail.ustc.edu.cn

1 引 言

矩陣乘是許多關(guān)鍵科學(xué)計(jì)算算法的核心,包括線(xiàn)性方程組的求解、最小二乘問(wèn)題、奇異和特征值計(jì)算[1],這一本質(zhì)使得它成為許多科學(xué)應(yīng)用領(lǐng)域的關(guān)鍵高性能kernel之一,包括深度學(xué)習(xí)、信號(hào)處理、天體物理學(xué)和先進(jìn)物理分析等諸多領(lǐng)域[2,3].除此之外,矩陣乘還經(jīng)常被用于評(píng)估針對(duì)科學(xué)計(jì)算的新型處理器架構(gòu)的性能和效率[4],以及探究如何在新架構(gòu)上進(jìn)行高性能優(yōu)化[5-7].矩陣乘的廣泛使用以及計(jì)算量稠密且巨大的本質(zhì),使得關(guān)于它的高性能計(jì)算研究一直是學(xué)術(shù)界和工業(yè)界的熱門(mén)話(huà)題之一.

當(dāng)前矩陣乘性算法的研究主要可以分為兩個(gè)方面[8],一個(gè)是宏觀層次的分塊方案,基于存儲(chǔ)器層次結(jié)構(gòu),通過(guò)合適的矩陣分塊策略,探索數(shù)據(jù)局部性,協(xié)調(diào)計(jì)算和訪(fǎng)存的平衡;一個(gè)是微觀層次的計(jì)算kernel,基于高級(jí)語(yǔ)言(如：C,C++等)或者底層匯編,通過(guò)SIMD運(yùn)算、指令多發(fā)射和亂序執(zhí)行等微體系結(jié)構(gòu)技術(shù),力求最大化硬件執(zhí)行單元的計(jì)算性能.

在學(xué)術(shù)界,K?gstr?等人[9]最早提出所有的三級(jí)BLAS(basic linear algebra subprograms)操作都可以在一個(gè)高性能矩陣乘的基礎(chǔ)上實(shí)現(xiàn),而現(xiàn)實(shí)中許多科學(xué)應(yīng)用都與頻繁的稠密線(xiàn)性操作緊密相關(guān),大算力的需求使得眾核處理器上的矩陣乘優(yōu)化一直備受?chē)?guó)內(nèi)外研究學(xué)者的關(guān)注.Aberdeen等人[10]考慮到Intel Pentium等商用處理器在矩陣乘設(shè)計(jì)時(shí)多級(jí)緩存有效使用的困難,針對(duì)Intel Pentium III系列處理器詳細(xì)探究了矩陣乘實(shí)現(xiàn)時(shí)的寄存器分配以及L1、L2緩存的分塊策略,最終平均性能是先進(jìn)的公開(kāi)矩陣乘例程的2.09倍.Nakasato[11]在Cypress GPU上優(yōu)化并實(shí)現(xiàn)了單精度和雙精度矩陣乘的kernel,分別可以達(dá)到理論峰值性能的73%和87%.同時(shí)證明了相比于共享內(nèi)存,紋理緩存對(duì)于Cypress架構(gòu)更有效.Kurzak等人[12]針對(duì)Fermi GPU提出了一種通過(guò)啟發(fā)式自動(dòng)調(diào)優(yōu)生成高效矩陣乘kernel的方法,其搜索空間的修剪機(jī)制,可以根據(jù)不同的需求很容易地將搜索空間調(diào)小或者調(diào)大.實(shí)驗(yàn)證明其自動(dòng)生成的kernel可以等同甚至超過(guò)cuBLAS和MAGMA,同時(shí)在雙精度復(fù)數(shù)的數(shù)據(jù)類(lèi)型下有顯著的性能提升.Lim等人[13]在矩陣乘分塊優(yōu)化數(shù)據(jù)重用的基礎(chǔ)上,通過(guò)數(shù)據(jù)預(yù)取、循環(huán)展開(kāi)和Intel AVX-512優(yōu)化分塊后的矩陣乘,實(shí)現(xiàn)了Intel KNL上的高效矩陣乘算法,性能可以與官方最優(yōu)庫(kù)Intel MKL相媲美.Lim等人[14]針對(duì)Intel KNL和Intel Skylake-SP處理器提出了一種啟發(fā)式自動(dòng)調(diào)校生成矩陣乘kernel的方法,調(diào)校參數(shù)包括寄存器級(jí)/緩存級(jí)的矩陣塊大小、預(yù)期距離和循環(huán)展開(kāi)深度.結(jié)果顯示kernel性能優(yōu)于開(kāi)源的BLAS庫(kù),矩陣乘性能同Intel MKL相當(dāng).Chen等人[15]針對(duì)GPU架構(gòu)提出了一種TSM2的矩陣乘算法,主要關(guān)注于不規(guī)則形狀的輸入矩陣.該算法在Kepler、Maxwell和Pascal 3種GPU架構(gòu)上,計(jì)算性能加速1.1x～3x,訪(fǎng)存帶寬提升8%～47.6%,算力利用率提升7%～37.3%.總體來(lái)說(shuō),學(xué)術(shù)界當(dāng)前對(duì)于矩陣乘的研究主要集中于通用眾核處理器Intel Xeon、Intel Xeon Phi和NVIDIA GPU展開(kāi).

在工業(yè)界,許多眾核處理器廠商提供了便于用戶(hù)使用的高效BLAS庫(kù),比如：Intel的oneMKL[16],NVIDIA的cuBLAS[17],AMD的AOCL[18]等,這些數(shù)學(xué)庫(kù)對(duì)于使得各自的處理器平臺(tái)廣泛應(yīng)用于現(xiàn)實(shí)世界起到了巨大的推動(dòng)作用.

申威眾核處理器SW26010是世界上第一臺(tái)峰值運(yùn)算性能超過(guò)100PFlops的超級(jí)計(jì)算機(jī)神威·太湖之光算力的主要貢獻(xiàn)者,也是當(dāng)前眾核處理器領(lǐng)域一顆冉冉升起的新星[19].不同于發(fā)展歷史悠久的通用眾核處理器,申威處理器上矩陣乘的研究仍處在初步階段,現(xiàn)有的研究工作[20-22]也僅僅是對(duì)雙精度矩陣乘算法設(shè)計(jì)進(jìn)行了討論.然而僅僅針對(duì)雙精度矩陣乘的研究工作是難以覆蓋現(xiàn)實(shí)中全部需求的,尤其是近些年迅速發(fā)展的人工智能領(lǐng)域,其使用的深度神經(jīng)網(wǎng)絡(luò)(deep neural networks,DNNs)技術(shù)中單精度矩陣乘是性能的核心[23].另一方面,申威處理器雙精度數(shù)據(jù)運(yùn)算和單精度數(shù)據(jù)運(yùn)算共用雙精度浮點(diǎn)運(yùn)算單元的微體系結(jié)構(gòu)特征,使得單精度矩陣乘算法設(shè)計(jì)無(wú)法直接依賴(lài)于雙精度矩陣乘的研究工作.

綜上所述,為了滿(mǎn)足現(xiàn)實(shí)生活中人工智能等應(yīng)用領(lǐng)域的需求,彌補(bǔ)申威處理器SW26010上單精度矩陣乘性能不足的問(wèn)題,本文主要貢獻(xiàn)如下所示：

1) 面向國(guó)產(chǎn)申威眾核處理器的架構(gòu)特征,在先前研究工作[22]的基礎(chǔ)上,提出了3種基于不同存儲(chǔ)層次的并行單精度矩陣乘算法;

2)對(duì)于國(guó)產(chǎn)申威處理器上的單精度計(jì)算過(guò)程中數(shù)據(jù)類(lèi)型轉(zhuǎn)換的可能情況,從計(jì)算層面和訪(fǎng)存層面進(jìn)行了綜合討論,能夠?yàn)槲磥?lái)該處理器上涉及到單精度計(jì)算的研究工作提供有意義的參考借鑒;

3)從峰值性能和通用性?xún)蓚€(gè)方面對(duì)本文工作進(jìn)行了綜合評(píng)估.實(shí)驗(yàn)證明,本文設(shè)計(jì)的算法能夠有效提升申威處理器上單精度矩陣乘的性能,其運(yùn)行時(shí)峰值性是理論峰值性能的86.17%,對(duì)比最先進(jìn)的官方庫(kù)xMath中的實(shí)現(xiàn),提升了6.8%.同時(shí)且重要的是,在通用性的實(shí)驗(yàn)分析中,相比xMath的實(shí)現(xiàn),在95.33%的場(chǎng)景中能夠獲得性能提升,最高性能加速比達(dá)到247.9%,平均性能加速比為61.66%.

可見(jiàn)本文提出并實(shí)現(xiàn)的單精度矩陣乘算法,相比于當(dāng)前SW26010處理器上的最優(yōu)實(shí)現(xiàn),不僅能夠?qū)Σ煌仃嚦藞?chǎng)景起到更好的適應(yīng)性,而且性能提升顯著.

2 背景介紹

2.1 矩陣乘

BLAS[24,25]有三級(jí)子程序,分別是一級(jí)的向量-向量操作、二級(jí)的矩陣-向量操作和三級(jí)的矩陣-矩陣操作.矩陣乘屬于三級(jí)BLAS,是一個(gè)二維數(shù)據(jù)乘累加的例行程序,如公式(1)所示：

C=αop(A)op(B)+βC,op(X)=XorXT

(1)

其中A∈M×K和B∈K×N為輸入矩陣,有轉(zhuǎn)置和不轉(zhuǎn)置兩種情況,C=M×N為輸出矩陣;α和β是標(biāo)量,為矩陣乘的運(yùn)算系數(shù);M、N和K表示矩陣的維度大小.

基礎(chǔ)的矩陣乘算法是一個(gè)關(guān)于M、N和K的三層循環(huán),如算法1所示.其是所有BLAS操作中最關(guān)鍵的部分[26],因?yàn)槠渌鄶?shù)操作可以依據(jù)矩陣乘實(shí)現(xiàn),同時(shí)也是最具優(yōu)化空間的部分,因?yàn)樗兄罡叩挠?jì)算復(fù)雜度O(n3).本文主要面向單精度矩陣乘在申威處理器上的并行算法進(jìn)行研究,主要目的在于彌補(bǔ)該處理器上單精度研究工作的不足,以及推動(dòng)其在人工智能領(lǐng)域的發(fā)展.

算法1.基礎(chǔ)的矩陣乘算法

輸入：矩陣A、B和C,矩陣乘系數(shù)α和β;

輸出：矩陣C.

1. form=0∶1∶Mdo

2. forn=0∶1∶Ndo

3.C[m,n]×=β

4. fork=0∶1∶Kdo

5.C[m,n]+=α(A[m,k]×B[k,n])

6. end for

7. end for

8. end for

2.2 SW26010眾核處理器

SW26010處理器[27]是由我國(guó)自主研發(fā)的一款異構(gòu)眾核處理器,支持64位自主申威指令級(jí)系統(tǒng),采用分布式共享存儲(chǔ)和片上計(jì)算陣列的異構(gòu)眾核體系架構(gòu).

如圖1所示,整個(gè)處理器包含4個(gè)核組(core-group,CG),每個(gè)核組由一個(gè)主要負(fù)責(zé)任務(wù)管理的主核(management processing element,MPE)和64個(gè)主要負(fù)責(zé)計(jì)算任務(wù)的從核(computing processing element,CPE)構(gòu)成,共計(jì)260個(gè)運(yùn)算核心.四個(gè)核組通過(guò)片上網(wǎng)絡(luò)(network on chip,NoC)互連,每個(gè)核組直連一塊8GB的DDR3主存,主核和從核通過(guò)主存控制器(memory controller,MC)進(jìn)行相應(yīng)的訪(fǎng)存操作.

圖1 SW26010眾核處理器架構(gòu)Fig.1 Architecture of SW26010 many-core processor

主核擁有L1、L2兩級(jí)緩存,L1分為32KB的指令緩存和32KB的數(shù)據(jù)緩存,L2緩存大小為256KB.每個(gè)從核配置有一塊16KB的L1指令緩存和一塊64KB用戶(hù)可控的局部設(shè)備內(nèi)存(local device memory,LDM).整個(gè)從核陣列共享一塊大小為64KB的L2指令緩存.

該處理器為了能夠盡量簡(jiǎn)化微架構(gòu)設(shè)計(jì),從而使得有限的芯片面積能夠集成更多的計(jì)算核心,所以沒(méi)有采用通用處理器的傳統(tǒng)的多級(jí)緩存設(shè)計(jì)方案,而是為每個(gè)從核配置了一塊64KB的軟件控制的片上存儲(chǔ)空間LDM.基于此,該處理器支持兩種主從核間的數(shù)據(jù)訪(fǎng)問(wèn),一種是全局內(nèi)存訪(fǎng)問(wèn),也稱(chēng)為gld/gst離散訪(fǎng)主存,即是從核直接對(duì)主存進(jìn)行讀寫(xiě)操作,這種方式的好處是簡(jiǎn)單易用,缺陷是訪(fǎng)存延遲高達(dá)278個(gè)時(shí)鐘周期;一種是以L(fǎng)DM為橋梁進(jìn)行主存的訪(fǎng)問(wèn),也稱(chēng)為DMA(direct memory access)批量式訪(fǎng)主存,整個(gè)過(guò)程的訪(fǎng)存延遲較低,大約29個(gè)時(shí)鐘周期.但由于LDM大小有限、DMA操作復(fù)雜多變等原因?qū)е滤惴ㄔO(shè)計(jì)難度大,要想獲取高性能對(duì)開(kāi)發(fā)人員自身體系結(jié)構(gòu)基礎(chǔ)要求較高.盡管如此,后者仍是SW26010處理器上研究開(kāi)發(fā)高性能工作的首選.

為盡可能降低訪(fǎng)存頻率,SW26010處理器為用戶(hù)提供了單核組從核陣列上不同從核間的數(shù)據(jù)共享機(jī)制——寄存器通信.為了有效支持這一機(jī)制,每個(gè)從核上配備了發(fā)送緩沖區(qū)、行接收緩沖區(qū)和列接收緩沖區(qū).寄存器通信機(jī)制需要注意3點(diǎn)：1) 通信數(shù)據(jù)的大小固定為256-bits;2) 只能行通信或者列通行;3) 通信的過(guò)程是匿名的,基于先到先得的策略接收消息.

最后,SW26010處理器每個(gè)從核支持兩條硬件流水線(xiàn)P0和P1,其中P0流水支持浮點(diǎn)和整數(shù)的標(biāo)量/向量操作,P1流水支持?jǐn)?shù)據(jù)遷移、比較、跳轉(zhuǎn)和整數(shù)標(biāo)量操作.兩條流水共享一個(gè)指令解碼器(Instruction Decoder,ID)和一個(gè)指令隊(duì)列,每個(gè)時(shí)鐘周期ID對(duì)隊(duì)列中前兩條指令進(jìn)行檢測(cè),當(dāng)滿(mǎn)足這些情況時(shí),兩條指令可以同時(shí)被加載到兩條流水中：1) 兩條指令同已發(fā)射未完成的指令不存在沖突;2) 兩條指令間不存在(read after write)和(write after write)沖突;3) 兩條指令可以分別被兩條流水處理.如何混合交差程序中兩種類(lèi)型的指令,使P0和P1兩條流水充分并行工作,對(duì)于SW26010處理器性能的發(fā)揮起著重要的作用.

3 算法設(shè)計(jì)

不同于通用眾核處理器,既支持雙精度浮點(diǎn)運(yùn)算單元,又支持單精度浮點(diǎn)運(yùn)算單元.申威處理器只具備雙精度浮點(diǎn)運(yùn)算單元,因此在進(jìn)行單精度運(yùn)算時(shí),需要顯式地或者隱式地進(jìn)行數(shù)據(jù)類(lèi)型轉(zhuǎn)換,以間接完成雙精度浮點(diǎn)運(yùn)算單元計(jì)算單精度數(shù)據(jù)的任務(wù).

基于上述硬件特征,如圖2所示,在申威處理器上進(jìn)行矩陣乘算法設(shè)計(jì)時(shí),單精度矩陣乘算法有必要在雙精度矩陣乘算法實(shí)現(xiàn)的基礎(chǔ)上,額外地增加一個(gè)數(shù)據(jù)類(lèi)型轉(zhuǎn)換模塊的設(shè)計(jì).

圖2 單精度矩陣乘算法設(shè)計(jì)的基本思想Fig.2 Base idea of single-precision matrix multiplication algorithm

本文將在先前研究工作[22]的基礎(chǔ)上,對(duì)單精度矩陣乘算法進(jìn)行深入研究.依據(jù)申威處理器的存儲(chǔ)器層次結(jié)構(gòu)可以將雙精度矩陣乘劃分為3個(gè)層次,分別是寄存器級(jí)的矩陣乘、LDM級(jí)的矩陣乘和主存級(jí)的矩陣乘,具體實(shí)現(xiàn)細(xì)節(jié)可以參考[22].同理,數(shù)據(jù)類(lèi)型轉(zhuǎn)換模塊也可以插入到上述3個(gè)不同的層次中,分別對(duì)應(yīng)寄存器級(jí)的數(shù)據(jù)類(lèi)型轉(zhuǎn)換、LDM級(jí)的數(shù)據(jù)類(lèi)型轉(zhuǎn)換和主存級(jí)的數(shù)據(jù)類(lèi)型轉(zhuǎn)換,三者的核心區(qū)別在于轉(zhuǎn)換后數(shù)據(jù)的暫存存儲(chǔ)空間不同,也由此提出了本文的3種不同的單精度矩陣乘算法——基于寄存器級(jí)數(shù)據(jù)類(lèi)型轉(zhuǎn)換的單精度矩陣乘,標(biāo)記為RSMM;基于LDM級(jí)數(shù)據(jù)類(lèi)型轉(zhuǎn)換的單精度矩陣乘,標(biāo)記為L(zhǎng)SMM;基于主存級(jí)數(shù)據(jù)類(lèi)型轉(zhuǎn)換的單精度矩陣乘,標(biāo)記為MSMM.

表1 統(tǒng)一符號(hào)聲明Table 1 Descriptor of symbols

為了方便后續(xù)敘述,本文對(duì)相關(guān)符號(hào)統(tǒng)一聲明,如表1所示.

3.1 RSMM

單精度矩陣乘算法最直接的設(shè)計(jì)方式就是將[22]中所有的雙精度操作全部替換為單精度操作,如算法2所示.

算法2.基于寄存器級(jí)數(shù)據(jù)類(lèi)型轉(zhuǎn)換的單精度矩陣乘算法

輸入：矩陣ASM、BSM和CSM,矩陣乘系數(shù)α和β;

輸出：矩陣CSM.

1.loadAASM[0,0]toASL[cmpt]

2.loadBBSM[0,0]toBSL[cmpt]

3.end loadA

4.end loadB

5.form=0∶1∶mMdo

6. forn=0∶1∶nMdo

7. loadCCSM[m,n] toCSL

8.CSL×=β

9. fork=0∶1∶kMdo

10. computem′,n′,k′of next memory access aboutASMandBSM

11. loadAASM[m′,k′]toASL[ldst]

12. loadBBSM[k′,n′]toBSL[ldst]

13.ASL[cmpt]×=α

14. kernel(ASL[cmpt],BSL[cmpt],CSL)

15. end loadA

16. end loadB

17. exchange the value ofcmptandldst

18. end for

19. storeCCSLtoCSM[m,n]

20. end storeC

21. end for

22.end for

首先對(duì)主存上的單精度矩陣數(shù)據(jù)ASM、BSM和CSM進(jìn)行分塊,分別劃分mM×kM、kM×nM和mM×nM個(gè)矩陣塊.然后將劃分后的矩陣塊按計(jì)算需要依次加載到LDM的單精度矩陣空間ASL、BSL和CSL中,其中ASL又被分為ASL[cmpt]和ASL[ldst],BSL又被分為BSL[cmpt]和BSL[ldst],cmpt和ldst分別用于存放當(dāng)前計(jì)算需要的矩陣塊以及下一次計(jì)算需要的矩陣塊.通過(guò)對(duì)ASM和BSM矩陣塊的加載和核心計(jì)算kernel進(jìn)行一次循環(huán)的錯(cuò)開(kāi),實(shí)現(xiàn)了申威處理器上異步DMA訪(fǎng)存操作下的計(jì)算和訪(fǎng)存的并行.最后更新后的CSL會(huì)被寫(xiě)回CSM中對(duì)應(yīng)的主存空間.

因?yàn)樯晖幚砥鞑淮嬖趩尉雀↑c(diǎn)運(yùn)算單元的緣故,算法2中行14的初始實(shí)現(xiàn)如圖3左側(cè)所示,ldse和vlds指令包含單精度數(shù)據(jù)到雙精度數(shù)據(jù)類(lèi)型的隱式轉(zhuǎn)換過(guò)程,前者加載轉(zhuǎn)換單精度ASL并向量擴(kuò)展到雙精度ADR,后者向量加載并轉(zhuǎn)換單精度BSL到雙精度BDR,整個(gè)計(jì)算過(guò)程由vmad指令運(yùn)行于雙精度浮點(diǎn)運(yùn)算單元.此處展示的是RSMM中kernel的最內(nèi)層循環(huán)部分的指令序列,也是指令重排的部分,該循環(huán)前的CSL到CDR的轉(zhuǎn)換和循環(huán)后的CDR到CSL的轉(zhuǎn)換不在指令重排的考慮范圍.

初始情況下整條指令序列的執(zhí)行時(shí)間是60個(gè)時(shí)鐘周期,IPC(instruction per cycle)為0.65,遠(yuǎn)低于理論值2,造成這一糟糕結(jié)果的原因主要有兩點(diǎn)——1)P0和P1兩條流水線(xiàn)的運(yùn)行幾乎是串行的;2)訪(fǎng)存指令(ldse和vlds)和寄存器通信指令(putr和putc)間的數(shù)據(jù)依賴(lài)造成了大量的流水線(xiàn)氣泡.為此本文參考[28,29],通過(guò)手動(dòng)重排指令序列,解決上述問(wèn)題,提升整個(gè)指令序列的并行度.整體的重排思想有兩點(diǎn)——1)循環(huán)交叉適用于P0和P1流水線(xiàn)的指令,盡可能地掩藏向量計(jì)算外的其他指令的開(kāi)銷(xiāo);2)ldse和vlds指令延遲是4個(gè)時(shí)鐘周期[27],通過(guò)使具有相關(guān)數(shù)據(jù)依賴(lài)的putr和putc指令的理論開(kāi)始時(shí)間滯后至少4個(gè)時(shí)鐘周期,從而消除不必要的流水線(xiàn)氣泡.結(jié)果如圖3右側(cè)所示,重排后指令序列的執(zhí)行時(shí)間是23個(gè)時(shí)鐘周期,IPC約為1.7,性能是初始時(shí)的2.61倍.

圖3 計(jì)算kernel的指令重排Fig.3 Reorder the instructions of the computational kernel

3.2 LSMM

RSMM雖然是最直觀的一種雙精度矩陣乘映射到單精度矩陣乘的算法,但其計(jì)算kernel最內(nèi)層循環(huán)指令重排后的開(kāi)銷(xiāo)高達(dá)23個(gè)時(shí)鐘周期,相比于[22]的16個(gè)時(shí)鐘周期,計(jì)算性能下降了30.4%.針對(duì)這一問(wèn)題,本文提出了LSMM算法,通過(guò)將數(shù)據(jù)類(lèi)型轉(zhuǎn)換操作從寄存器轉(zhuǎn)移到片上存儲(chǔ)LDM上,從而實(shí)現(xiàn)雙精度矩陣乘計(jì)算kernel的無(wú)縫鏈接.

圖4 優(yōu)化LDM空間使用Fig.4 Optimize the usage of the LDM

如圖4左側(cè)所示,為每一個(gè)ASL、BSL和CSL分別匹配一個(gè)雙倍大小的ADL、BDL和CDL,用于暫存轉(zhuǎn)換后的雙精度數(shù)據(jù).此時(shí)LSMM相比于RSMM,雖然計(jì)算kernel的性能得到了極大的改善,但是卻需要額外兩倍的LDM開(kāi)銷(xiāo),這對(duì)于申威處理器64KB的有限片上存儲(chǔ)來(lái)說(shuō)是難以接受的.為了提高LDM的使用效率,如圖4右側(cè)所示,分配一個(gè)CDL空間,通過(guò)嵌套CSL和CDL,實(shí)現(xiàn)兩者LDM空間物理上的共用和邏輯上的分離;分配一個(gè)ASL/BSL和一個(gè)ADL/BDL,固定前者被用于DMA讀主存數(shù)據(jù),固定后者被用于暫存類(lèi)型轉(zhuǎn)換后的數(shù)據(jù).此時(shí),LSMM的額外LDM開(kāi)銷(xiāo)下降至約60%,整體設(shè)計(jì)思路如算法3所示.

算法3.基于LDM級(jí)數(shù)據(jù)類(lèi)型轉(zhuǎn)換的單精度矩陣乘算法

輸入：矩陣ASM、BSM和CSM,矩陣乘系數(shù)α和β;

輸出：矩陣CSM.

1.loadAASM[0,0]toASL

2.loadBBSM[0,0]toBSL

3.end loadA

4.end loadB

5.ADL=(double)(ASL×α)

6.BDL=(double)(BSL)

7.form=0∶1∶mMdo

8. forn=0∶1∶nMdo

9. loadCCSM[m,n]toCSL

10.CDL=(double)(CSL×β)

11. fork=0∶1∶kMdo

12. computem′,n′,k′of next memory access aboutASMandBSM

13. loadAASM[m′,k′]toASL

14. loadBBSM[k′,n′]toCDL

15. kernel(ADL,BDL,BDL)

16. end loadA

17. end loadB

18.ADL=(double)(ASL×α)

19.BDL=(double)(BSL)

20. end for

21.CSL=(float)(CDL)

22. storeCCSLtoCSM[m,n]

23. end storeC

24. end for

25.end for

對(duì)于算法3可以概述如下：1)DMA讀主存數(shù)據(jù)和DMA寫(xiě)主存數(shù)據(jù)都以單精度數(shù)據(jù)類(lèi)型進(jìn)行;2)數(shù)據(jù)類(lèi)型轉(zhuǎn)換操作緊鄰DMA實(shí)時(shí)進(jìn)行,轉(zhuǎn)換結(jié)果暫存于LDM;3)計(jì)算kernel直接使用轉(zhuǎn)換后的雙精度數(shù)據(jù).對(duì)于CSL和CDL的嵌套,物理上CDL的后半部分也是CSL,從而使得算法3的行10順序遍歷時(shí),不會(huì)造成轉(zhuǎn)換后的雙精度數(shù)據(jù)覆蓋后續(xù)轉(zhuǎn)換所需要的單精度數(shù)據(jù)的錯(cuò)誤.相應(yīng)的,算法3中行21則要逆序遍歷,以保證轉(zhuǎn)換數(shù)據(jù)的準(zhǔn)確性.同時(shí),對(duì)于ASL/BSL和ADL/BDL的轉(zhuǎn)換操作需要在DMA讀主存數(shù)據(jù)完成后立即進(jìn)行,從而閑置ASL/BSL,為下一次的DMA讀主存數(shù)據(jù)做準(zhǔn)備.

3.3 MSMM

LSMM雖然保證了計(jì)算kernel的高效運(yùn)行,但數(shù)據(jù)類(lèi)型轉(zhuǎn)換的開(kāi)銷(xiāo)仍是不可忽視的.其中,CSM的數(shù)據(jù)類(lèi)型轉(zhuǎn)換僅進(jìn)行了一次,ASM和BSM卻分別重復(fù)轉(zhuǎn)換了nM次和mM次.由此可見(jiàn),隨著矩陣規(guī)模逐漸增大,ASM和BSM的數(shù)據(jù)類(lèi)型轉(zhuǎn)換對(duì)LSMM性能的影響將會(huì)越來(lái)越明顯.為了解決這一問(wèn)題,本文提出了MSMM算法,通過(guò)消除不必要的數(shù)據(jù)類(lèi)型轉(zhuǎn)化操作,進(jìn)一步提升大規(guī)模矩陣乘運(yùn)算時(shí)的性能.

圖5 消除不必要的數(shù)據(jù)類(lèi)型轉(zhuǎn)換開(kāi)銷(xiāo)Fig.5 Eliminate the unnecessary cost of data type conversion

如圖5所示,在主存上為ASM和BSM分別分配一個(gè)雙倍大小的ADM和BDM,用于暫存轉(zhuǎn)換后的雙精度數(shù)據(jù).此時(shí)MSMM執(zhí)行包含兩個(gè)階段：1)DMA讀主存數(shù)據(jù)ASM和BSM,轉(zhuǎn)換單精度數(shù)據(jù)到雙精度數(shù)據(jù),再DMA寫(xiě)回主存ADM和BDM;2)基于ADM、BDM和CSM執(zhí)行矩陣乘運(yùn)算.

通過(guò)用主存暫存轉(zhuǎn)換后的雙精度數(shù)據(jù),消除了LSMM中不必要的轉(zhuǎn)換開(kāi)銷(xiāo).但是此時(shí)也在轉(zhuǎn)換階段引入了新的DMA讀主存和DMA寫(xiě)主存操作.為此,本文對(duì)MSMM算法進(jìn)一步優(yōu)化,如算法4所示.

算法4.基于主存級(jí)數(shù)據(jù)類(lèi)型轉(zhuǎn)換的單精度矩陣乘算法

輸入：矩陣ASM、BSM和CSM,矩陣乘系數(shù)α和β;

輸出：矩陣CSM.

1.loadAASM[0,0]toASL[cmpt]

2.loadBBSM[0,0]toBSL[cmpt]

3.end loadA

4.end loadB

5.ADL[cmpt]=(double)(ASL[cmpt]×α)

6.BDL[cmpt]=(double)(BSL[cmpt])

7.form=0∶1∶mMdo

8. forn=0∶1∶nMdo

9. loadCCSM[m,n]toCSL

10.CDL=(double)(CSL×β)

11. fork=0∶1∶kMdo

12. computem′,n′,k′of next memory access aboutASMandBSM(ADMandBDM)

13. ifADM[m′,k′]未暫存轉(zhuǎn)換后的結(jié)果

14. loadAASM[m′,k′]toASL[ldst]

15. else

16. loadAADM[m′,k′]toADL[ldst]

17. end if

18. ifASM[m,k]轉(zhuǎn)換后的數(shù)據(jù)未被寫(xiě)回過(guò)主存

19. storeAADL[cmpt]toADM[m,k]

20. end if

21. ifBDM[k′,n′]未暫存轉(zhuǎn)換后的結(jié)果

22. loadABSM[k′,n′]toBSL[ldst]

23. else

24. loadABDM[k′,n′]toBDL[ldst]

25. end if

26. ifBSM[k,n]轉(zhuǎn)換后的數(shù)據(jù)未被寫(xiě)回過(guò)主存

27. storeBBDL[cmpt]toBDM[k,n]

28. end if

29. kernel(ADL[cmpt],BDL[cmpt],CDL)

30. end loadA

31. end loadB

32. ifADM[m′,k′]未暫存轉(zhuǎn)換后的結(jié)果

33.ADL[ldst]=(double)(ASL[ldst]×α)

34. end if

35. ifBDM[k′,n′]未暫存轉(zhuǎn)換后的結(jié)果

36.BDL[ldst]=(double)(BSL[ldst])

37. end if

38. ifASM[m,k]轉(zhuǎn)換后的數(shù)據(jù)未被寫(xiě)回過(guò)主存

39. end storeA

40. end if

41. ifBSM[k,n]轉(zhuǎn)換后的數(shù)據(jù)未被寫(xiě)回過(guò)主存

42. end storeB

43. end if

44. exchange the value ofcmptandldst

45. end for

46.CSL=(float)(CDL)

47. storeCCSLtoCSM[m,n]

48. end storeC

49. end for

50.end for

通過(guò)將MSMM的數(shù)據(jù)類(lèi)型轉(zhuǎn)換階段和矩陣乘運(yùn)算階段融合到一起,雙緩沖消除了額外的DMA訪(fǎng)存操作.在算法4中,本文在LDM上分配雙精度矩陣空間ADL、BDL和CDL,同時(shí)將ASL、BSL和CSL分別嵌套在它們的后半部分,從而實(shí)現(xiàn)物理上的共用和邏輯上的分離.其中,ADL/BDL又被分為ADL[cmpt]/BDL[cmpt]和ADL[ldst]/BDL[ldst],分別用于存放當(dāng)前計(jì)算需要的矩陣塊以及下一次計(jì)算需要的矩陣塊.類(lèi)似的,對(duì)應(yīng)的ASL/BSL也被分為ASL[cmpt]/BSL[cmpt]和ASL[ldst]/BSL[ldst].在每次DMA讀下一次kernel計(jì)算所需的矩陣塊ADM[m′,k′]和BDM[k′,n′]時(shí),判斷其主存空間中是否已經(jīng)存入轉(zhuǎn)換好的雙精度數(shù)據(jù).如果沒(méi)有,則DMA讀ASM[m′,k′]和BSM[k′,n′],并在加載之后轉(zhuǎn)換為雙精度數(shù)據(jù);如果有,則DMA直接讀ADM[m′,k′]和BDM[k′,n′].類(lèi)似的,在每次計(jì)算kernel前,判斷當(dāng)前計(jì)算所對(duì)應(yīng)的矩陣塊ASM[m,k]和BSM[k,n]轉(zhuǎn)換后的數(shù)據(jù)是否已經(jīng)被寫(xiě)回主存.如果沒(méi)有,則將ADL[cmpt]和BDL[cmpt]寫(xiě)回對(duì)應(yīng)的主存空間ADM[m,k]和BDM[k,n].通過(guò)上述操作,本文利用DMA異步技術(shù)下計(jì)算和訪(fǎng)存的雙緩沖,盡可能地減少了不必要的訪(fǎng)存開(kāi)銷(xiāo).雖然該算法能夠最大化運(yùn)行時(shí)峰值性能,但是卻需要額外的雙倍ASM+BSM字節(jié)大小的主存空間開(kāi)銷(xiāo).因此,在使用該算法時(shí)應(yīng)考慮到SW26010單芯片主存容量32GB的架構(gòu)特性.

4 實(shí)驗(yàn)結(jié)果及分析

4.1 實(shí)驗(yàn)設(shè)置

本文使用的申威眾核處理器SW26010,配置如表2所示.

本文主要在SW26010處理器單核組實(shí)現(xiàn)了上述3種不同的單精度矩陣乘算法,因?yàn)榭绮煌私M的并行通常由更高編程級(jí)別的用戶(hù)自己處理.對(duì)比對(duì)象為該平臺(tái)上最先進(jìn)的官方數(shù)學(xué)庫(kù)xMath中的單精度矩陣乘,標(biāo)記為SWSMM.后續(xù)將針對(duì)兩個(gè)方面進(jìn)行實(shí)驗(yàn)設(shè)計(jì)和結(jié)果分析,一個(gè)是峰值性能,一個(gè)是通用性,前者反映了本文工作在硬件資源上發(fā)揮出的最大利用率,后者驗(yàn)證了本文工作能否更好地應(yīng)用于現(xiàn)實(shí)世界中多變的矩陣乘配置.

表2 SW26010配置參數(shù)Table 2 Configure parameters for SW26010

本文所有測(cè)試均是迭代10次,去掉一個(gè)最優(yōu)值和最差值,取剩余測(cè)試結(jié)果的平均值.

4.2 峰值性能測(cè)試

圖6 申威處理器上不同單精度矩陣乘實(shí)現(xiàn)的運(yùn)行時(shí)峰值性能比較Fig.6 Peak performance comparison of different single-precision matrix multiplication implementations on SW26010

方陣(M=N=K)是矩陣乘運(yùn)行時(shí)峰值性能測(cè)試的最佳選擇,因此本文選取M=N=K=128X,其中X=1,2,…,100,共100組不同的矩陣乘配置測(cè)試.

如圖6所示,用當(dāng)前最先進(jìn)的單精度矩陣乘實(shí)現(xiàn)SWSMM對(duì)比本文提出的3種算法RSMM、LSMM和MSMM.其中,SWSMM的最高性能可以達(dá)到598.85GFlops,而RSMM、LSMM和MSMM分別可以達(dá)到489.03GFlops、609.11GFlops和639.66GFlops.RSMM遠(yuǎn)低于SWSMM,主要原因是計(jì)算kernel指令執(zhí)行周期更長(zhǎng),限制了運(yùn)行時(shí)的峰值性能;LSMM與SWSMM較為接近,可以猜想兩者的實(shí)現(xiàn)原理大致相同;MSMM則明顯高于SWSMM,可見(jiàn)其極大地降低了大規(guī)模單精度矩陣乘中數(shù)據(jù)類(lèi)型轉(zhuǎn)換的開(kāi)銷(xiāo).

單獨(dú)對(duì)比本文提出的3種算法RSMM、LSMM和SWSMM,從圖6可以看出,當(dāng)方陣規(guī)模在896以?xún)?nèi)時(shí),RSMM性能較高;隨著規(guī)模的增加,LSMM性能變得更為突出;最后當(dāng)規(guī)模達(dá)到4608之后,MSMM的性能基本最優(yōu).這是因?yàn)殡S著方陣規(guī)模的增加,訪(fǎng)存的影響逐漸減小,計(jì)算kernel的影響逐漸增大,這使得小規(guī)模時(shí)RSMM較優(yōu);反之,LSMM的性能更好,這是因?yàn)長(zhǎng)SMM是以額外的LDM占用來(lái)?yè)Q取更高的性能.LSMM和MSMM計(jì)算kernel相同,但是隨著方陣規(guī)模的再度增加,數(shù)據(jù)類(lèi)型轉(zhuǎn)換的影響越發(fā)突出,這也使得MSMM的性能優(yōu)勢(shì)開(kāi)始凸顯.

圖7 不同矩陣乘配置的性能加速比Fig.7 Speedup of different matrix multiplication configurations

最終本文工作最高可以發(fā)揮申威處理器單精度理論峰值性能的86.17%,相比于SWSMM的運(yùn)行時(shí)峰值性能提升了6.8%.同時(shí),從圖6性能的變化趨勢(shì)中,可以看出本文工作實(shí)現(xiàn)的單精度矩陣乘性能波動(dòng)較小,穩(wěn)定性更好.

4.3 通用性測(cè)試

矩陣形狀的不同使得矩陣乘在現(xiàn)實(shí)世界中也是多種多樣的,依據(jù)M、N和K相對(duì)大小的不同,可以分為8種不同的類(lèi)型,標(biāo)記為BMBNBK、SMSNSK、SMBNBK、BMSNBK、BMBNSK、BMSNSK、SMBNSK和SMSNBK.‘B’表示相對(duì)值小,‘S’表示相對(duì)值大,舉例說(shuō)明：SMBNBK就表示小M大N大K的矩陣乘.在3.2節(jié)中方陣的測(cè)試,代表著B(niǎo)MBNBK和SMSNSK兩種場(chǎng)景,但是比起方陣這種矩陣乘,剩余的6種長(zhǎng)條形矩陣乘場(chǎng)景更為普遍.本節(jié)以SWSMM為基準(zhǔn),SMBNBK等6種場(chǎng)景各取25例,測(cè)試并分析本文工作在各個(gè)矩陣乘配置下的最優(yōu)加速比.

如圖7所示,從單類(lèi)場(chǎng)景來(lái)看,BMSNBK、BMBNSK、BMSNSK和SMSNBK中,本文工作相比SWSMM在所有的矩陣乘測(cè)試中都具有明顯性能提升,最低提升分別為49.49%、154.53%、12.93%和12.1%,最高提升分別為67.73%、179.02%、247.9%和54.43.

SMBNBK和SMBNSK是存在性能下降的兩類(lèi)場(chǎng)景,而造成這一現(xiàn)象的主要原因很可能是由于在某些情況下本文算法的矩陣乘分塊方式不夠完善,沒(méi)有發(fā)揮算法的最大潛力,使得其性能略低于SWSMM.這兩種場(chǎng)景雖然都有性能下降的情況,但是無(wú)論在程度還是數(shù)量上都遠(yuǎn)遠(yuǎn)低于性能提升的情況,兩者最差時(shí)性能分別下降0.7%和23.5%,而性能最優(yōu)時(shí)分別提升18.6%和116%,且性能提升的情況在數(shù)量上分別占據(jù)92%和80%.所有場(chǎng)景來(lái)看,平均性能提升達(dá)到了61.66%,并在95.33%的情況下性能更高.由此可見(jiàn),本文工作因既考慮了計(jì)算層面的優(yōu)化技術(shù)又兼顧了訪(fǎng)存層面的優(yōu)化技術(shù),從而相比SWSMM不僅性能提升顯著,而且通用性更強(qiáng).

5 總 結(jié)

隨著我國(guó)自主研發(fā)的申威眾核處理器SW26010在科學(xué)計(jì)算和人工智能領(lǐng)域的快速發(fā)展,對(duì)于該眾核處理器上高性能矩陣乘算法實(shí)現(xiàn)提出了更高的要求.而現(xiàn)有的關(guān)于矩陣乘的研究工作寥寥無(wú)幾,僅有的工作也僅僅是針對(duì)雙精度矩陣乘算法的討論.考慮到SW26010特有的架構(gòu)特征——雙精度數(shù)據(jù)運(yùn)算和單精度數(shù)據(jù)運(yùn)算共用雙精度浮點(diǎn)運(yùn)算單元,使得單精度矩陣乘算法直接套用雙精度矩陣乘算法難以充分發(fā)揮該處理器的硬件性能.本文針對(duì)上述問(wèn)題,結(jié)合SW26010的體系結(jié)構(gòu)特征,提出了3種不同存儲(chǔ)層次的單精度矩陣乘并行算法,分別是RSMM、LSMM和MSMM.在進(jìn)行算法設(shè)計(jì)時(shí),綜合考慮了計(jì)算和訪(fǎng)存兩個(gè)核心因素,在計(jì)算方面,結(jié)合該處理器的從核雙流水,從匯編層面手動(dòng)控制核心計(jì)算任務(wù)的指令序列,保證了高效的指令級(jí)并行;在訪(fǎng)存方面,綜合考慮了有限片上存儲(chǔ)資源的有效使用,以及訪(fǎng)存任務(wù)和計(jì)算任務(wù)的交叉并行.在實(shí)驗(yàn)分析中,以當(dāng)前最先進(jìn)的官方庫(kù)xMath中的實(shí)現(xiàn)SWSMM作為基準(zhǔn),從峰值性能和通用性?xún)蓚€(gè)方面進(jìn)行測(cè)試分析,最終證明本文工作不僅能夠不同程度上提升單精度矩陣乘運(yùn)算的性能,而且對(duì)于不同的矩陣乘場(chǎng)景也有著更好的適應(yīng)性.

下一步將重點(diǎn)研究本文優(yōu)化算法同人工智能應(yīng)用的具體結(jié)合,從而探索申威眾核處理上的深度學(xué)習(xí)并行優(yōu)化技術(shù).