遲利華,劉 杰,晏益慧,謝林川,甘新標,胡慶豐,蔣 杰,李勝國

(國防科技大學(xué) 并行與分布處理重點實驗室，湖南 長沙 410073)

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

遲利華?,劉 杰,晏益慧,謝林川,甘新標,胡慶豐,蔣 杰,李勝國

(國防科技大學(xué) 并行與分布處理重點實驗室，湖南 長沙 410073)

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

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

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

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

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

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

1 FT1000微處理器

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

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

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

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

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

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

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

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

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

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

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

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

alpha為寄存器，保存標量值

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多級展開循環(huán)體

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

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

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

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

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

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

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

alpha為寄存器，保存標量值

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多級展開循環(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來實現(xiàn).下面以GEMM子程序為例來進行說明.所要求解的矩陣乘形式如下：

C=βC+αA×B

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

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

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

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

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

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

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

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

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

循環(huán)展開因子為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展開循環(huán)體

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

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

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

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

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

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

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

循環(huán)展開因子為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ù)取的展開循環(huán)體

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

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

3 多線程并行計算方法

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

下面重點對BLAS3中的矩陣乘并行計算方法展開討論.

P個線程按pr×pc劃分成二維拓撲結(jié)構(gòu), 滿足P=pr×pc.假設(shè)P(r,c) (0≤r

多線程矩陣乘并行算法：

全局數(shù)組

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

每一個線程P(r, c)，其中 ，執(zhí)行：

for i→r×mrto r×mr+mrdo

for l←0 to k-1 do

第r行線程中的任一線程將子矩陣 拷貝到

for j←c×ncto c×nc+ncdo

end for

end for

end for

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

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

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

4 性能測試與分析

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

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

No of threads

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

No of threads

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

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

M=N=K

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

5 結(jié)論與展望

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

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

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

[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]張先軼，王茜，張云泉. OpenBLAS：龍芯3A CPU的高性能BLAS庫 [J]. 軟件學(xué)報, 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

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

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

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

TP332.2

A