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

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

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

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

中圖分類號：TP332.2 文獻標識碼：A

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

湖南大學學報（自然科學版）2015年

第4期遲利華等：FitenBLAS：面向FT1000微處理器的高性能線性代數(shù)庫

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

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

1 FT1000微處理器

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

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

每個核內有獨立的一級數(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ù)庫不能直接運行，需要重新設計.

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ù)目進行多級展開，展開因子可以隨意調節(jié)，在展開循環(huán)體內，循環(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)體重復m“*”n/2次.

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

根據(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ù)實際測試情況，進行調整，以找到最佳的循環(huán)展開因子.

BLAS2包括矩陣與向量操作子程序，如矩陣向量乘、rank1和rank2矩陣校正等，涉及二維數(shù)組A和一維數(shù)組x和y間的計算操作，計算結果為一個一維數(shù)組.以矩陣向量乘GEMV子程序為例進行說明，為了復用一維數(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ù)實際測試情況，進行調整，以找到最佳的循環(huán)展開因子.圖2給出了BLAS2子程序GEMV多級展開循環(huán)體.

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ù).

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

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

Ci，0…Ci，N-1=Ai，p×

Bp，0，…，Bp，N-1.

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

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

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

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

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

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

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

AXPY帶數(shù)據(jù)預取的展開循環(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ù)預取的展開循環(huán)體

Fig.4 Multilevel loop unrolling with the prefetching

data for subroutine AXPY of BLAS1

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

3 多線程并行計算方法

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

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

P 個線程按 pr×pc 劃分成二維拓撲結構， 滿足 P=pr×pc.假設 P（r， c） （0≤r

多線程矩陣乘并行算法：

全局數(shù)組

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

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

for i→r×mr to r×mr+mr do

for l←0 to k-1 do

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

for j←c×nc to c×nc+nc do

將ij拷貝到

Call GEMM_kernel（bm，bn，bk，，，ij， LDC）

end for

end for

end for

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

Fig.5 Multithread parallel matrix multiplication

algorithm avoiding the redundant packing of A

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

No of threads

圖6 BLAS1不同線程數(shù)性能測試結果（n=256 000 000）

Fig.6 Computation performance for BLAS1

on difference number of threads（n=256 000 000）

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

No of threads

圖7 BLAS2不同線程數(shù)性能測試結果（n=16 000）

Fig.7 Computation performance for BLAS2

on difference number of threads（n=16 000）

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

M=N=K

圖8 BLAS3不同計算規(guī)模性能測試結果（64線程）

Fig.8 Computation performance

of difference scale of matrixes on 64 threads

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

5 結論與展望

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

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

下一步，面向短向量飛騰微處理器，研制支持向量優(yōu)化運算的BLAS庫，吸收并行編程框架底層調用的矩陣向量運算和稀疏線性數(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/zhcn/intelmkl/

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

[4] IBM ESSL homepage. http：//www03.ibm.com/systems/software/essl/

[5] GotoBLAS homepage. http：//www.tacc.utexas.edu/taccprojects/gotoblas2

[6] GOTO K， VAN DE GEIJN R. Anatomy of highperformance matrix multiplication [J]. ACM Transactions on Mathematical Software， 2008， 34（3）： 1-25.

[7] ATLAS homepage. http：//mathatlas.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]. 軟件學報， 2011， 22（zk2）： 208-216.

ZHANG Xianyi， WANG Qian， ZHANG Yunquan. 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.nscctj.gov.cn/ resources/resource_1.asp， 2010-11-17.

[12]GOTO K， VAN DE GEIJIN R. High performance implementation of the level3 BLAS[J]. ACM Transactions on Mathematical Software， 2008， 35（1）： 1-14.