FitenBLAS：面向FT1000微處理器的高性能線(xiàn)性代數(shù)庫(kù)*

遲利華,劉 杰,晏益慧,謝林川,甘新標(biāo),胡慶豐,蔣 杰,李勝?lài)?guó)

(國(guó)防科技大學(xué) 并行與分布處理重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410073)

FitenBLAS：面向FT1000微處理器的高性能線(xiàn)性代數(shù)庫(kù)*

遲利華?,劉 杰,晏益慧,謝林川,甘新標(biāo),胡慶豐,蔣 杰,李勝?lài)?guó)

(國(guó)防科技大學(xué) 并行與分布處理重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410073)

BLAS庫(kù)是基本線(xiàn)性代數(shù)子程序庫(kù)，是許多大型科學(xué)與工程計(jì)算的核心計(jì)算程序，F(xiàn)itenBLAS庫(kù)是在多核多線(xiàn)FT1000微處理器上開(kāi)發(fā)的基本線(xiàn)性代數(shù)庫(kù)，其研制對(duì)FT1000微處理器在科學(xué)與工程計(jì)算中的應(yīng)用具有重要意義.根據(jù)多級(jí)存儲(chǔ)結(jié)構(gòu)和寄存器的數(shù)目，設(shè)計(jì)了向量與向量、矩陣與向量和矩陣與矩陣運(yùn)算的多級(jí)循環(huán)展開(kāi)方法，采用指令調(diào)度、數(shù)據(jù)預(yù)取等通用優(yōu)化技術(shù)，優(yōu)化BLAS庫(kù)串行程序.對(duì)于BLAS3子程序，設(shè)計(jì)了矩陣乘無(wú)冗余數(shù)據(jù)拷貝分塊算法，采用指令重排、訪(fǎng)存與計(jì)算的重疊、分塊等技術(shù)優(yōu)化矩陣乘子程序，基于矩陣乘子程序?qū)崿F(xiàn)了其他BLAS3子程序.研制了匯編線(xiàn)性代數(shù)程庫(kù)FitenBLAS，其核心子程序矩陣乘的雙精度計(jì)算性能達(dá)到6.91Gflops，是峰值性能的86.4%.

FT1000微處理器；BLAS庫(kù)；性能優(yōu)化

基本線(xiàn)性代數(shù)子程序BLAS(Basic Linear Algebra Subroutines)庫(kù)，提供最基本的線(xiàn)性代數(shù)函數(shù)接口[1]，分為三級(jí)：BLAS 1(Level 1)包括向量與向量操作子程序,如點(diǎn)積、向量相加等.BLAS 2(Level 2)包括矩陣與向量操作子程序，如矩陣向量相乘等.BLAS 3(Level 3)包括矩陣與矩陣操作子程序，如矩陣與矩陣相乘等.

BLAS庫(kù)是每款微處理器要移植和優(yōu)化的數(shù)學(xué)庫(kù)，是許多大型科學(xué)與工程計(jì)算的核心計(jì)算模塊，同時(shí)BLAS庫(kù)子程序可以反映許多應(yīng)用程序的計(jì)算特點(diǎn)，如果BLAS庫(kù)可以在微處理器上獲得高性能，同樣的應(yīng)用程序也可以獲得好的性能，BLAS庫(kù)程序可以驗(yàn)證微處理器的功能和計(jì)算性能.因此各個(gè)廠(chǎng)家在新型號(hào)微處理器推出時(shí)，都會(huì)配套針對(duì)微處理器特點(diǎn)研制、優(yōu)化和推出高性能BLAS庫(kù)，BLAS庫(kù)已經(jīng)成為微處理器的必備數(shù)學(xué)庫(kù)之一.

Intel公司針對(duì)通用CPU開(kāi)發(fā)了MKL基本數(shù)學(xué)運(yùn)算庫(kù)[2]，包含采用多線(xiàn)程進(jìn)行并行計(jì)算的函數(shù)庫(kù)，可以在結(jié)點(diǎn)內(nèi)實(shí)現(xiàn)高性能；AMD公司針對(duì)通用CPU開(kāi)發(fā)了ACML基本數(shù)學(xué)運(yùn)算庫(kù)[3]，具有MKL類(lèi)似的特點(diǎn)；IBM公司開(kāi)發(fā)了ESSL基本數(shù)學(xué)運(yùn)算庫(kù)[4].上述廠(chǎng)商開(kāi)發(fā)的基本數(shù)學(xué)運(yùn)算庫(kù)，均包含了使用匯編語(yǔ)言手工優(yōu)化的高性能BLAS庫(kù).GotoBLAS[5-6]是開(kāi)源的針對(duì)不同類(lèi)型的微處理器開(kāi)發(fā)的BLAS庫(kù)，提供了OpenMP多線(xiàn)程并行版本，是性能最好的BLAS庫(kù)之一，2008后不再更新，不支持最新推出的微處理器；ATLAS[7-8]針對(duì)不同的微處理器提供可以自動(dòng)優(yōu)化的BLAS庫(kù).OpenBLAS是針對(duì)龍芯等微處理器開(kāi)發(fā)的高性能BLAS庫(kù)[9-10].

FT1000是由國(guó)防科學(xué)技術(shù)大學(xué)研制的單芯片多線(xiàn)程(CMT)處理器，是天河-1A計(jì)算結(jié)點(diǎn)采用的微處理器之一[11].BLAS庫(kù)在FT1000上獲得高性能是需要研究的重要問(wèn)題,目前沒(méi)有可以在FT1000上運(yùn)行的高性能BLAS庫(kù)，本文結(jié)合FT1000多核多線(xiàn)微處理器特點(diǎn)，設(shè)計(jì)了并行計(jì)算方法和數(shù)據(jù)結(jié)構(gòu)，研制了手工匯編子程序，進(jìn)行了針對(duì)性的性能優(yōu)化，研制了高性能線(xiàn)性代數(shù)庫(kù)FitenBLAS.

1 FT1000微處理器

FT1000微處理器包含8個(gè)處理器核，每核包含8套硬件現(xiàn)場(chǎng)，支持8個(gè)線(xiàn)程.每個(gè)線(xiàn)程有1個(gè)完整的寄存器文件，大部分ASI，ASR和特權(quán)寄存器都是每個(gè)線(xiàn)程1份.

每個(gè)核包含2條整數(shù)流水線(xiàn)、1條浮點(diǎn)流水線(xiàn)和1條存儲(chǔ)流水線(xiàn).8個(gè)線(xiàn)程共享浮點(diǎn)流水線(xiàn)和存儲(chǔ)流水線(xiàn).8個(gè)線(xiàn)程分成2組，每組4個(gè)線(xiàn)程，共享1條整數(shù)流水線(xiàn).雖然8個(gè)線(xiàn)程同時(shí)運(yùn)行，但是在任意時(shí)刻，最多兩個(gè)線(xiàn)程是活躍的，這兩個(gè)線(xiàn)程發(fā)射的指令只可能是下面的組合：1對(duì)整數(shù)操作、1個(gè)整數(shù)和1個(gè)浮點(diǎn)、1個(gè)整數(shù)和1個(gè)存儲(chǔ)、1個(gè)浮點(diǎn)和1個(gè)存儲(chǔ).每個(gè)組內(nèi)的線(xiàn)程按照LRU算法每個(gè)周期進(jìn)行切換.

每個(gè)核內(nèi)有獨(dú)立的一級(jí)數(shù)據(jù)cache和一級(jí)指令cache.L1I Cache大小為 16 kB，8路組相聯(lián)，塊大小32字節(jié).L1D Cache大小為8 kB，4路組相聯(lián)，塊大小32字節(jié).指令TLB為64項(xiàng)全相聯(lián)，數(shù)據(jù)TLB為128項(xiàng)全相聯(lián).8個(gè)線(xiàn)程共享L1I,L1D和TLB，通過(guò)TLB中的自動(dòng)釋放機(jī)制使多個(gè)線(xiàn)程更新TLB.

FT1000微處理器采用SPARC V9指令集，現(xiàn)有的BLAS線(xiàn)性代數(shù)庫(kù)不能直接運(yùn)行，需要重新設(shè)計(jì).

2 程序優(yōu)化方法

2.1 循環(huán)展開(kāi)方法

循環(huán)展開(kāi)可以減少循環(huán)體的循環(huán)次數(shù)，減少分支執(zhí)行的時(shí)間，為流水線(xiàn)提供更多的并行機(jī)會(huì)，是一種通用的程序優(yōu)化方法，是手工匯編優(yōu)化BLAS子程序的主要方法.在循環(huán)展開(kāi)中，展開(kāi)因子的選擇是核心問(wèn)題，目前的編譯器一般只對(duì)循環(huán)體很小的循環(huán)進(jìn)行循環(huán)展開(kāi)，且使用的展開(kāi)因子是很小的常數(shù)，這可能損失一些計(jì)算性能，特別是對(duì)于多層嵌套循環(huán)，編譯器的優(yōu)化效果不明顯，只能采用手工的循環(huán)展開(kāi)優(yōu)化.

BLAS 1(Level 1)包括向量與向量操作子程序，包含SUM,DOT,AXPY,COPY,SCAL等子程序．當(dāng)跨步為1時(shí)，BLAS1程序中的計(jì)算操作針對(duì)一維數(shù)據(jù)進(jìn)行連續(xù)訪(fǎng)問(wèn)，具有良好的空間數(shù)據(jù)局部性，可以采用手工循環(huán)展開(kāi)來(lái)優(yōu)化性能.采用匯編編程時(shí)寄存器的數(shù)目限制了展開(kāi)次數(shù)，為此，本文根據(jù)寄存器的數(shù)目進(jìn)行多級(jí)展開(kāi)，展開(kāi)因子可以隨意調(diào)節(jié)，在展開(kāi)循環(huán)體內(nèi)，循環(huán)使用寄存器.以AXPY為例來(lái)說(shuō)明多級(jí)展開(kāi)方法，圖1給出了多級(jí)展開(kāi)循環(huán)體，雙精度浮點(diǎn)寄存器的數(shù)目為n，涉及2個(gè)向量運(yùn)算，平均使用向量寄存器，每個(gè)向量使用的寄存器最大數(shù)目為n/2,循環(huán)展開(kāi)因子可以取為m“*”n/2.ldd,fmuld,faddd和std是匯編指令，分別表示取數(shù)、相乘、相加和存數(shù).為了表示方便，圖1中使用了循環(huán)控制變量，在實(shí)際展開(kāi)的循環(huán)體中，沒(méi)有循環(huán)控制變量，而將圖1中的循環(huán)體重復(fù)m“*”n/2次.

AXPY使用寄存器的多級(jí)展開(kāi)循環(huán)體:

浮點(diǎn)寄存器數(shù)目為n，令k=n/2

循環(huán)展開(kāi)因子為m*k

數(shù)組a(*)和b(*)表示不同的寄存器

alpha為寄存器，保存標(biāo)量值

for i ← 0 to m*k-1 step k do

for j ← 0 to k -1 do

ldd x[m*i+j], a[j];

fmuld a[j], alpha, a[j];

ldd y[m*i+j], b[j];

faddd a[j], b[j], b[j];

std b[j], y[m*i+j];

endfor

endfor

圖1 BLAS1子程序AXPY多級(jí)展開(kāi)循環(huán)體

Fig.1 Multi-level loop unrolling for subroutine AXPY of BLAS1

根據(jù)寄存器數(shù)目，圖1中給出的是AXPY多級(jí)循環(huán)展開(kāi)的一般方法，BLAS1其它子程序中可以根據(jù)圖1中給出的展開(kāi)方法，進(jìn)行多級(jí)循換展開(kāi).針對(duì)FT1000微處理器，雙精度寄存器數(shù)目為32，對(duì)于AXPY子程序，每個(gè)向量至多使用16個(gè)寄存器．使用16個(gè)寄存器時(shí)，循環(huán)展開(kāi)因子就是16的倍數(shù)，同樣地循環(huán)展開(kāi)因子可以是12,13,14和15等數(shù)的倍數(shù)，可以根據(jù)實(shí)際測(cè)試情況，進(jìn)行調(diào)整，以找到最佳的循環(huán)展開(kāi)因子.

BLAS2包括矩陣與向量操作子程序，如矩陣-向量乘、rank-1和rank-2矩陣校正等，涉及二維數(shù)組A和一維數(shù)組x和y間的計(jì)算操作，計(jì)算結(jié)果為一個(gè)一維數(shù)組.以矩陣向量乘GEMV子程序?yàn)槔M(jìn)行說(shuō)明，為了復(fù)用一維數(shù)組x，對(duì)數(shù)組x進(jìn)行分段，對(duì)數(shù)組A進(jìn)行分塊.在分配寄存器時(shí)，二維數(shù)組A,一維數(shù)組x和y以及保存臨時(shí)變量的一維數(shù)組z給分配寄存器總數(shù)的1/4.針對(duì)FT1000微處理器，雙精度寄存器數(shù)目為32，對(duì)于GEMV子程序，數(shù)組A至多使用8個(gè)寄存器，一維數(shù)組x和y分別使用8個(gè)寄存器.使用8個(gè)寄存器時(shí)，循環(huán)展開(kāi)因子就是8的倍數(shù)，同樣地循環(huán)展開(kāi)因子可以是4,5,6和7等數(shù)的倍數(shù)，可以根據(jù)實(shí)際測(cè)試情況，進(jìn)行調(diào)整，以找到最佳的循環(huán)展開(kāi)因子.圖2給出了BLAS2子程序GEMV多級(jí)展開(kāi)循環(huán)體.

GEMV使用寄存器的多級(jí)展開(kāi)循環(huán)體:

浮點(diǎn)寄存器數(shù)目為q，令k=q/4

循環(huán)展開(kāi)因子為m*k和n*k

數(shù)組M(*,*)、a(*)和b(*)表示不同的寄存器

alpha為寄存器，保存標(biāo)量值

for i ← 0 to m*k-1 step k do

for j ← 0 to n*k -1step k do

for v ←0 to k-1 do

ldd x[n*j+v], a[v];

fmuld a[v], alpha, a[v];

for u ←0 to k-1 do

ldd A[m*i+u, n*j+v], M[u,v];

ldd y[m*i+u], b[u];

fmuld M[u,v], a[v], c[u];

faddd b[u], c[u], b[u];

endfor

endfor

endfor

for u ←0 to k-1 do

std b[u], y[m*i+u];

endfor

endfor

圖2 BLAS2子程序GEMV多級(jí)展開(kāi)循環(huán)體

Fig.2 Multi-level loop unrolling for subroutine GEMV of BLAS2

BLAS 3(Level 3)包括矩陣與矩陣操作子程序，包含GEMM,SYMM,SYRK,SYR2K,TRMM和TRSM等子程序.矩陣乘GEMM子程序是BLAS3的核心，其它BLAS3子程序可以基于GEMM來(lái)實(shí)現(xiàn).下面以GEMM子程序?yàn)槔齺?lái)進(jìn)行說(shuō)明.所要求解的矩陣乘形式如下：

C=βC+αA×B

其中A∈Rm×k,B∈Rk×n,C∈Rm×n，α和β是實(shí)數(shù).

矩陣乘在多級(jí)存儲(chǔ)結(jié)構(gòu)上獲得高性能的基本方法是分塊算法，將A，B和C劃分成如下子矩陣：

其中Aip∈Rmc×kc,Bpj∈Rkc×nr，Cij∈Rmc×nr.

原矩陣乘算法轉(zhuǎn)化為多個(gè)子矩陣相乘，根據(jù)Cache和TLB的大小來(lái)選擇子矩陣的大小.具體實(shí)現(xiàn)時(shí)選取Aip的大小為L(zhǎng)2Cache大小的一半，同時(shí)保證不發(fā)生TLB失效，并駐留在二級(jí)Cache中，直到不再使用為止，即完成如下操作：

計(jì)算時(shí)，子矩陣Bp,j(j=1,…,N)以流水的方式進(jìn)入L1DCache.Ci,j的數(shù)據(jù)不重用，不必長(zhǎng)時(shí)間保存在L1D和L2Cache中.

把寄存器總數(shù)的一半分配給矩陣C使用，另一半寄存器被矩陣A，矩陣B和臨時(shí)變量平均使用.針對(duì)FT1000微處理器，雙精度寄存器數(shù)目為32，對(duì)于GEMM子程序，數(shù)組C至多使用16個(gè)寄存器，數(shù)組A和B分別使用4個(gè)寄存器，剩下的作為臨時(shí)變量寄存器使用.圖3給出了BLAS3子程序GEMM按照4*4*4展開(kāi)的循環(huán)體.

通過(guò)指令的合理調(diào)度、數(shù)據(jù)的預(yù)取和寄存器的合理使用，當(dāng)Ai,p，Bp,j和Ci,j在Cache中時(shí)，可以發(fā)揮CPU的峰值性能.

GEMM_kernel使用寄存器的循環(huán)展開(kāi)體:

浮點(diǎn)寄存器數(shù)目為q，令t=4

循環(huán)展開(kāi)因子為m*t、n*t和k*t

數(shù)組D(*,*)、a(*)和b(*)表示不同的寄存器

for j ← 0 to n*t-1 step t do

for i ← 0 to m*t-1 step t do

for k ← 0 to q*t-1 step t do

for u ←0 to t-1 do

for v ←0 to t-1 do

for w ←0 to t-1 do

ldd A[t*i+u, t*j+w], a[w];

ldd B[t*i+w, t*j+w], b[w];

fmuld, a[w], b[w] , tmp;

faddd M[u,v], tmp, M[u,v];

endfor

endfor

endfor

endfor

endfor

endfor

圖3 BLAS3子程序GEMM按照4*4*4展開(kāi)循環(huán)體

Fig.3 Multi-level loop unrolling for subroutine GEMM of BLAS3

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

矩陣和向量是存放于內(nèi)存中的，計(jì)算過(guò)程中，通過(guò)多級(jí)存儲(chǔ)結(jié)構(gòu)，將數(shù)據(jù)取到寄存器中開(kāi)始運(yùn)算，數(shù)據(jù)的預(yù)取可以很好地實(shí)現(xiàn)數(shù)據(jù)傳輸和計(jì)算的重疊，是提高數(shù)據(jù)空間局部性的性能優(yōu)化方法.開(kāi)始執(zhí)行計(jì)算時(shí)，參與運(yùn)算的矩陣和向量存放在內(nèi)存中，開(kāi)始的取數(shù)不會(huì)命中Cache，F(xiàn)T1000上從內(nèi)存中取數(shù)到二級(jí)Cache的延遲大于100拍(即100個(gè)CPU時(shí)鐘周期)，Cache塊大小為32字節(jié)，每次內(nèi)存訪(fǎng)問(wèn)，同時(shí)會(huì)把內(nèi)存中連續(xù)的32字節(jié)分別存放到二級(jí)Cache和一級(jí)數(shù)據(jù)Cache中.

數(shù)據(jù)預(yù)取就是將后面要用的數(shù)據(jù)提前取到L2 Cache中，將訪(fǎng)存的數(shù)據(jù)傳輸操作和乘加計(jì)算操作重疊起來(lái)，如果計(jì)算時(shí)間大于數(shù)據(jù)傳輸時(shí)間，那么整個(gè)計(jì)算就可以得到CPU的峰值性能，如矩陣相乘.通過(guò)參數(shù)可以調(diào)整數(shù)據(jù)預(yù)取的間隔，獲得最佳性能.

BLAS1，BLAS2和BLAS3均用到數(shù)據(jù)預(yù)取，本文僅以雙精度AXPY為例來(lái)說(shuō)明數(shù)據(jù)預(yù)取方法，圖4給出了帶數(shù)據(jù)預(yù)取的展開(kāi)循環(huán)體.

AXPY帶數(shù)據(jù)預(yù)取的展開(kāi)循環(huán)體：

令k=8，PREFECHSIZE=8，SIZE=8

循環(huán)展開(kāi)因子為m*8

數(shù)組a(*)、b(*)和c(*)表示不同的寄存器

for i ← 0 to m*8 step 8 do

for j ← 0 to 4 do

ldd x[m*i+j], a[j];

fmuld a[j], alpha, a[j];

ldd y[m*i+j], b[j];

faddd a[j], b[j], b[j];

std b[j], y[m*i+j];

endfor

prefetch [&x[m*i]+PREFETCHSIZE * SIZE], 0

prefetch [&y[m*i]+PREFETCHSIZE * SIZE], 0

for j ← 5 to 8 do

ldd x[m*i+j], a[j];

fmuld a[j], alpha, a[j];

ldd y[m*i+j], b[j];

faddd a[j], b[j], b[j];

std b[j], y[m*i+j];

endfor

prefetch [&x[m*i+4]+PREFETCHSIZE * SIZE], 0

prefetch [&y[m*i+4]+PREFETCHSIZE * SIZE], 0

endfor

圖4 BLAS1子程序AXPY帶數(shù)據(jù)預(yù)取的展開(kāi)循環(huán)體

Fig.4 Multi-level loop unrolling with the prefetching data for subroutine AXPY of BLAS1

FT1000通過(guò)預(yù)取指令對(duì)內(nèi)存數(shù)據(jù)進(jìn)行操作，SIZE表示一個(gè)數(shù)據(jù)所占的內(nèi)存字節(jié)數(shù)，PREFECHSIZE表示提前預(yù)取數(shù)據(jù)間隔，可以根據(jù)需要改變大小.圖4中，預(yù)取的數(shù)據(jù)提供給下一次循環(huán)體使用，將數(shù)據(jù)從內(nèi)存中，取到二級(jí)Cache中，預(yù)取數(shù)的時(shí)間和圖4中的兩個(gè)小循環(huán)進(jìn)行時(shí)間重疊，用計(jì)算重疊數(shù)據(jù)傳輸.

3 多線(xiàn)程并行計(jì)算方法

對(duì)BLAS1采用向量一維平均劃分的方式組織并行計(jì)算.對(duì)BLAS2中涉及的二維矩陣采用按行或按列一維劃分.對(duì)BLAS3中涉及的二維矩陣采用按行和按列二維劃分.

下面重點(diǎn)對(duì)BLAS3中的矩陣乘并行計(jì)算方法展開(kāi)討論.

P個(gè)線(xiàn)程按pr×pc劃分成二維拓?fù)浣Y(jié)構(gòu), 滿(mǎn)足P=pr×pc.假設(shè)P(r,c) (0≤r

多線(xiàn)程矩陣乘并行算法：

全局?jǐn)?shù)組

局部數(shù)組 r，c，

每一個(gè)線(xiàn)程P(r, c)，其中 ，執(zhí)行：

for i→r×mrto r×mr+mrdo

for l←0 to k-1 do

第r行線(xiàn)程中的任一線(xiàn)程將子矩陣 拷貝到

for j←c×ncto c×nc+ncdo

end for

end for

end for

圖5 避免重復(fù)拷貝A子矩陣的多線(xiàn)程矩陣乘并行算法

Fig.5 Multi-thread parallel matrix multiplication algorithm avoiding the redundant packing ofA

BLAS3中的SYMM，SYR2K，SYRK，TRMM和TRSM子程序并行算法均基于GEMM實(shí)現(xiàn)，參考文獻(xiàn)[12].

4 性能測(cè)試與分析

測(cè)試環(huán)境的硬件平臺(tái)為8核FT1000微處理器，主頻為1 GHz，雙精度浮點(diǎn)性能8Gflops.操作系統(tǒng)為銀河麒麟操作系統(tǒng)，編譯器為gcc-4.5.1.

圖6給出了BLAS1的雙精度性能測(cè)試結(jié)果，固定向量的長(zhǎng)度n不變，統(tǒng)一取為256 000 000.從圖6可以看出，從1線(xiàn)程到8線(xiàn)程各BLAS1程序具有明顯的加速效果，各程序在16或32線(xiàn)程時(shí)，達(dá)到最高性能，64線(xiàn)程時(shí)性能略有下降.造成32或64線(xiàn)程性能下降的主要原因是：1)BLAS1程序的計(jì)算訪(fǎng)存比是1或2，受限于訪(fǎng)存；2)BLAS1程序訪(fǎng)存有空間局部性，沒(méi)有時(shí)間局部性，也就是不存在數(shù)據(jù)的復(fù)用；3)FT1000微處理器的訪(fǎng)存帶寬在16或32線(xiàn)程時(shí)達(dá)到最大.

No of threads

圖7給出了BLAS2的雙精度性能測(cè)試結(jié)果，固定矩陣的長(zhǎng)度n不變，統(tǒng)一取為16 000.從圖7可以看出，浮點(diǎn)性能明顯高于BLAS1的性能，性能最好的SYMV在64線(xiàn)程時(shí)達(dá)到最高性能，3.11Gflops/s，是峰值性能的38.8%，SYMV的計(jì)算訪(fǎng)存比是6，存在空間和時(shí)間局部性.BLAS2中的其他子程序的計(jì)算訪(fǎng)存比是3，只有空間局部性，不存在時(shí)間局部性.

No of threads

BLAS3雙精度浮點(diǎn)性能測(cè)試結(jié)果由圖8給出，由于BLAS3的計(jì)算訪(fǎng)存比大，對(duì)于階為N的方陣，計(jì)算量為2N3，訪(fǎng)存量為3N2，在64線(xiàn)程時(shí)獲得最高性能，對(duì)不同的計(jì)算規(guī)模展開(kāi)測(cè)試.從圖8中可以看出，矩陣乘GEMM的最高性能達(dá)到6.91Gflops，是峰值性能的86.4%.SYMM，SYR2K，SYRK，TRMM和TRSM子程序的最高性能分別達(dá)到6.75，6.73，6.74，6.75和6.69Gflops，和矩陣乘GEMM的性能接近.

影響矩陣乘GEMM性能的主要因素包括：1)數(shù)據(jù)拷貝開(kāi)銷(xiāo).為了充分發(fā)揮多級(jí)Cache的性能，采用了分塊算法，分塊子矩陣在內(nèi)存中的存儲(chǔ)是不連續(xù)的，為了減少對(duì)TLB沖突，需要將矩陣A和B的數(shù)據(jù)預(yù)先拷貝到一個(gè)連續(xù)的內(nèi)存空間.2)數(shù)據(jù)延遲時(shí)間.每個(gè)子矩陣塊相乘前，每個(gè)線(xiàn)程需要進(jìn)行數(shù)據(jù)的填充，數(shù)據(jù)需要從內(nèi)存中傳輸?shù)絃2 Cache，從L2 Cache到L1 Cache，從L1 Cache到寄存器，每個(gè)數(shù)據(jù)大概需要150拍，150拍過(guò)后，每拍可以取一個(gè)雙精度數(shù)據(jù).3)數(shù)據(jù)對(duì)存儲(chǔ)帶寬的競(jìng)爭(zhēng).64線(xiàn)程并行計(jì)算時(shí)，需要同時(shí)從內(nèi)存中取數(shù)，數(shù)據(jù)帶寬成為影響性能的因素.4)計(jì)算不能全覆蓋數(shù)據(jù)的存取.FT1000微處理器中每核啟動(dòng)8線(xiàn)程，8線(xiàn)程共享一套硬件計(jì)算資源，靠多線(xiàn)程的切換來(lái)屏蔽訪(fǎng)存延遲.這種情況主要存在每個(gè)循環(huán)展開(kāi)啟動(dòng)階段，8線(xiàn)程需要同時(shí)取數(shù)，此外循環(huán)展開(kāi)計(jì)算完成后，需要把計(jì)算結(jié)果保存到矩陣C中，此時(shí)存取數(shù)據(jù)量和計(jì)算量處于相同量級(jí).

M=N=K

相對(duì)于其他X86微處理器，F(xiàn)T1000是多核多線(xiàn)體系結(jié)構(gòu)，每8線(xiàn)程共享一個(gè)計(jì)算核，通過(guò)線(xiàn)程切換來(lái)屏蔽訪(fǎng)存延遲，相對(duì)而言更難發(fā)揮其峰值性能，我們認(rèn)為矩陣乘GEMM雙精度能發(fā)揮峰值性能的86.4%已經(jīng)是能達(dá)到的最好浮點(diǎn)性能.

5 結(jié)論與展望

提煉共性基礎(chǔ)數(shù)值算法，研制高性能計(jì)算庫(kù)，統(tǒng)一編程接口，是用戶(hù)充分發(fā)揮微處理器峰值性能的重要手段.BLAS庫(kù)是基礎(chǔ)線(xiàn)性代數(shù)庫(kù)，需要根據(jù)CPU的具體特點(diǎn)，設(shè)計(jì)高性能的算法和數(shù)據(jù)結(jié)構(gòu)，經(jīng)過(guò)高度的手工匯編優(yōu)化，其中的BLAS3可得到接近峰值的浮點(diǎn)性能，滿(mǎn)足用戶(hù)的性能要求.

本文基于國(guó)防科學(xué)技術(shù)大學(xué)研制的多核多線(xiàn)FT1000微處理器，研制了基本線(xiàn)性代數(shù)匯編子程序庫(kù)FitenBLAS，根據(jù)寄存器的數(shù)目，設(shè)計(jì)了向量與向量、矩陣與向量、矩陣與矩陣運(yùn)算的多級(jí)循環(huán)展開(kāi)方法，采用計(jì)算指令流水線(xiàn)調(diào)度、數(shù)據(jù)預(yù)取等通用優(yōu)化技術(shù)，優(yōu)化BLAS庫(kù)串行程序.對(duì)于BLAS3子程序，設(shè)計(jì)了矩陣乘無(wú)冗余數(shù)據(jù)拷貝分塊算法，采用指令重排、訪(fǎng)存與計(jì)算的重疊、分塊等技術(shù)優(yōu)化矩陣乘子程序，基于矩陣乘子程序?qū)崿F(xiàn)了其他BLAS3子程序.其核心子程序矩陣GEMM乘的雙精度計(jì)算性能達(dá)到6.91Gflops，是峰值性能的86.4%.

下一步，面向短向量飛騰微處理器，研制支持向量?jī)?yōu)化運(yùn)算的BLAS庫(kù)，吸收并行編程框架底層調(diào)用的矩陣向量運(yùn)算和稀疏線(xiàn)性數(shù)值算法，完善FitenBLAS庫(kù)，支持更加廣泛的數(shù)值模擬運(yùn)算.

[1] DONGARRA J. Basic linear algebra subprograms technical forum standard [J]. International Journal of High Performance Applications and Supercomputing, 2002, 16(1): 1-111．

[2] Intel MKL homepage. http://software.intel.com/zh-cn/intel-mkl/

[3 AMD ACML homepage. http://developer.amd.com/cpu/ libraries/acml/

[4] IBM ESSL homepage. http://www-03.ibm.com/systems/software/essl/

[5] GotoBLAS homepage. http://www.tacc.utexas.edu/tacc-projects/gotoblas2

[6] GOTO K, VAN DE GEIJN R. Anatomy of high-performance matrix multiplication [J]. ACM Transactions on Mathematical Software, 2008, 34(3): 1-25．

[7] ATLAS homepage. http://math-atlas.sourceforge.net/

[8] WHALEY R, PETITET A, DONGARRA J. Automated empirical optimizations of software and the ATLAS project [J]. Parallel Computing, 2001, 27(1): 3-35．

[9] OpenBLAS homepage. http://xianyi.github.com/ OpenBLAS

[10]張先軼，王茜，張?jiān)迫? OpenBLAS：龍芯3A CPU的高性能BLAS庫(kù) [J]. 軟件學(xué)報(bào), 2011, 22(zk2): 208-216．

ZHANG Xian-yi, WANG Qian, ZHANG Yun-quan. OpenBLAS: a high performance BLAS library on Loongson 3A CPU [J]. Journal of Software, 2011, 22(zk2): 208-216. (In Chinese)

[11]About FT1000 processor, http://www.nscc-tj.gov.cn/ resources/resource_1.asp, 2010-11-17．

[12]GOTO K, VAN DE GEIJIN R. High performance implementation of the level-3 BLAS[J]. ACM Transactions on Mathematical Software, 2008, 35(1): 1-14．

FitenBLAS: High Performance BLAS for a Massively Multithreaded FT1000 Processor

CHI Li-hua?, LIU Jie, YAN Yi-hui, XIE Lin-chuan, GAN Xin-biao,HU Qin-feng, JIANG Jie， LI Sheng-guo

(National Laboratory for Parallel and Distributed Processing, National Univ of Defense Technology,Changsha，Hunan 410073,China)

BLAS library is the fundamental linear algebra library and plays an important role in many large scientific applications. This paper developed a linear algebra library named FitenBLAS on a massively multithreaded FT1000 processor. Based on the hierarchical storage system and the number of registers, multilevel loop unrolling methods were developed for vector-vector, matrix-vector and matrix-matrix linear operations. The codes of FitenBLAS were optimized with instruction layout and data prefetching technology. An avoiding redundant packing method was proposed for parallel matrix-matrix multiplication, and the parallel code was developed. The kernel matrix-matrix multiplication code was optimized with instruction layout, time overlapping of data access and computation, and data blocking. The other BLAS3 subroutines were designed on the matrix multiplication code. The kernel codes of FitenBLAS were developed in assembly language. The performance for the key subroutine of the matrix multiplication reaches 6.91Glops/s, nearly 86.4% of the peak performance of the FT1000.

FT1000 Processor； BLAS； code performance optimization

1674-2974(2015)04-0100-07

2014-10-16

國(guó)家自然科學(xué)基金資助項(xiàng)目(60970033)，National Natural Science Foundation of China(60970033) ；國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)資助項(xiàng)目(2012AA01A301)

遲利華(1970-)，女，山東威海人，國(guó)防科技大學(xué)副研究員，博士

?通訊聯(lián)系人，E-mail:chilihua@nudt.edu.cn

TP332.2

A