FitenBLAS：面向FT1000微处理器的高性能线性代数库

迟利华　刘杰　晏益慧　谢林川　甘新标　胡庆丰　蒋杰　李胜国

摘 要：BLAS库是基本线性代数子程序库，是许多大型科学与工程计算的核心计算程序，FitenBLAS库是在多核多线FT1000微处理器上开发的基本线性代数库，其研制对FT1000微处理器在科学与工程计算中的应用具有重要意义.根据多级存储结构和寄存器的数目，设计了向量与向量、矩阵与向量和矩阵与矩阵运算的多级循环展开方法，采用指令调度、数据预取等通用优化技术，优化BLAS库串行程序.对于BLAS3子程序，设计了矩阵乘无冗余数据拷贝分块算法，采用指令重排、访存与计算的重叠、分块等技术优化矩阵乘子程序，基于矩阵乘子程序实现了其他BLAS3子程序.研制了汇编线性代数程库FitenBLAS，其核心子程序矩阵乘的双精度计算性能达到6.91Gflops，是峰值性能的86.4%.

关键词：FT1000微处理器；BLAS库；性能优化

中图分类号：TP332.2 文献标识码：A

基本线性代数子程序BLAS（Basic Linear Algebra Subroutines）库，提供最基本的线性代数函数接口[1]，分为三级：BLAS 1（Level 1）包括向量与向量操作子程序，如点积、向量相加等.BLAS 2（Level 2）包括矩阵与向量操作子程序，如矩阵向量相乘等.BLAS 3（Level 3）包括矩阵与矩阵操作子程序，如矩阵与矩阵相乘等.

BLAS库是每款微处理器要移植和优化的数学库，是许多大型科学与工程计算的核心计算模块，同时BLAS库子程序可以反映许多应用程序的计算特点，如果BLAS库可以在微处理器上获得高性能，同样的应用程序也可以获得好的性能，BLAS库程序可以验证微处理器的功能和计算性能.因此各个厂家在新型号微处理器推出时，都会配套针对微处理器特点研制、优化和推出高性能BLAS库，BLAS库已经成为微处理器的必备数学库之一.

湖南大学学报（自然科学版）2015年

第4期迟利华等：FitenBLAS：面向FT1000微处理器的高性能线性代数库

Intel公司针对通用CPU开发了MKL基本数学运算库[2]，包含采用多线程进行并行计算的函数库，可以在结点内实现高性能；AMD公司针对通用CPU开发了ACML基本数学运算库[3]，具有MKL类似的特点；IBM公司开发了ESSL基本数学运算库[4].上述厂商开发的基本数学运算库，均包含了使用汇编语言手工优化的高性能BLAS库.GotoBLAS[5-6]是开源的针对不同类型的微处理器开发的BLAS库，提供了OpenMP多线程并行版本，是性能最好的BLAS库之一，2008后不再更新，不支持最新推出的微处理器；ATLAS[7-8]针对不同的微处理器提供可以自动优化的BLAS库.OpenBLAS是针对龙芯等微处理器开发的高性能BLAS库[9-10].

FT1000是由国防科学技术大学研制的单芯片多线程（CMT）处理器，是天河1A计算结点采用的微处理器之一[11].BLAS库在FT1000上获得高性能是需要研究的重要问题，目前没有可以在FT1000上运行的高性能BLAS库，本文结合FT1000多核多线微处理器特点，设计了并行计算方法和数据结构，研制了手工汇编子程序，进行了针对性的性能优化，研制了高性能线性代数库FitenBLAS.

1 FT1000微处理器

FT1000微处理器包含8个处理器核，每核包含8套硬件现场，支持8个线程.每个线程有1个完整的寄存器文件，大部分ASI，ASR和特权寄存器都是每个线程1份.

每个核包含2条整数流水线、1条浮点流水线和1条存储流水线.8个线程共享浮点流水线和存储流水线.8个线程分成2组，每组4个线程，共享1条整数流水线.虽然8个线程同时运行，但是在任意时刻，最多两个线程是活跃的，这两个线程发射的指令只可能是下面的组合：1对整数操作、1个整数和1个浮点、1个整数和1个存储、1个浮点和1个存储.每个组内的线程按照LRU算法每个周期进行切换.

每个核内有独立的一级数据cache和一级指令cache.L1I Cache大小为 16 kB，8路组相联，块大小32字节.L1D Cache大小为8 kB，4路组相联，块大小32字节.指令TLB为64项全相联，数据TLB为128项全相联.8个线程共享L1I，L1D和TLB，通过TLB中的自动释放机制使多个线程更新TLB.

FT1000微处理器采用SPARC V9指令集，现有的BLAS线性代数库不能直接运行，需要重新设计.

2 程序优化方法

2.1 循环展开方法

循环展开可以减少循环体的循环次数，减少分支执行的时间，为流水线提供更多的并行机会，是一种通用的程序优化方法，是手工汇编优化BLAS子程序的主要方法.在循环展开中，展开因子的选择是核心问题，目前的编译器一般只对循环体很小的循环进行循环展开，且使用的展开因子是很小的常数，这可能损失一些计算性能，特别是对于多层嵌套循环，编译器的优化效果不明显，只能采用手工的循环展开优化.

BLAS 1（Level 1）包括向量与向量操作子程序，包含SUM，DOT，AXPY，COPY，SCAL等子程序.当跨步为1时，BLAS1程序中的计算操作针对一维数据进行连续访问，具有良好的空间数据局部性，可以采用手工循环展开来优化性能.采用汇编编程时寄存器的数目限制了展开次数，为此，本文根据寄存器的数目进行多级展开，展开因子可以随意调节，在展开循环体内，循环使用寄存器.以AXPY为例来说明多级展开方法，图1给出了多级展开循环体，双精度浮点寄存器的数目为n，涉及2个向量运算，平均使用向量寄存器，每个向量使用的寄存器最大数目为n/2，循环展开因子可以取为m“*”n/2.ldd，fmuld，faddd和std是汇编指令，分别表示取数、相乘、相加和存数.为了表示方便，图1中使用了循环控制变量，在实际展开的循环体中，没有循环控制变量，而将图1中的循环体重复m“*”n/2次.

AXPY使用寄存器的多级展开循环体：

根据寄存器数目，图1中给出的是AXPY多级循环展开的一般方法，BLAS1其它子程序中可以根据图1中给出的展开方法，进行多级循换展开.针对FT1000微处理器，双精度寄存器数目为32，对于AXPY子程序，每个向量至多使用16个寄存器.使用16个寄存器时，循环展开因子就是16的倍数，同样地循环展开因子可以是12，13，14和15等数的倍数，可以根据实际测试情况，进行调整，以找到最佳的循环展开因子.

BLAS2包括矩阵与向量操作子程序，如矩阵向量乘、rank1和rank2矩阵校正等，涉及二维数组A和一维数组x和y间的计算操作，计算结果为一个一维数组.以矩阵向量乘GEMV子程序为例进行说明，为了复用一维数组x，对数组x进行分段，对数组A进行分块.在分配寄存器时，二维数组A，一维数组x和y以及保存临时变量的一维数组z给分配寄存器总数的1/4.针对FT1000微处理器，双精度寄存器数目为32，对于GEMV子程序，数组A至多使用8个寄存器，一维数组x和y分别使用8个寄存器.使用8个寄存器时，循环展开因子就是8的倍数，同样地循环展开因子可以是4，5，6和7等数的倍数，可以根据实际测试情况，进行调整，以找到最佳的循环展开因子.图2给出了BLAS2子程序GEMV多级展开循环体.

BLAS 3（Level 3）包括矩阵与矩阵操作子程序，包含GEMM，SYMM，SYRK，SYR2K，TRMM和TRSM等子程序.矩阵乘GEMM子程序是BLAS3的核心，其它BLAS3子程序可以基于GEMM来实现.下面以GEMM子程序为例来进行说明.所要求解的矩阵乘形式如下：

C=βC+αA×B

其中A∈Rm×k，B∈Rk×n，C∈Rm×n，α和β是实数.

矩阵乘在多级存储结构上获得高性能的基本方法是分块算法，将A，B和C划分成如下子矩阵：

原矩阵乘算法转化为多个子矩阵相乘，根据Cache和TLB的大小来选择子矩阵的大小.具体实现时选取Aip的大小为L2 Cache大小的一半，同时保证不发生TLB失效，并驻留在二级Cache中，直到不再使用为止，即完成如下操作：

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

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

计算时，子矩阵Bp，j（j=1，…，N）以流水的方式进入L1D Cache.Ci，j的数据不重用，不必长时间保存在L1D和L2 Cache中.

把寄存器总数的一半分配给矩阵C使用，另一半寄存器被矩阵A，矩阵B和临时变量平均使用.针对FT1000微处理器，双精度寄存器数目为32，对于GEMM子程序，数组C至多使用16个寄存器，数组A和B分别使用4个寄存器，剩下的作为临时变量寄存器使用.图3给出了BLAS3子程序GEMM按照4*4*4展开的循环体.

通过指令的合理调度、数据的预取和寄存器的合理使用，当Ai，p，Bp，j和Ci，j在Cache中时，可以发挥CPU的峰值性能.

2.2 数据预取

矩阵和向量是存放于内存中的，计算过程中，通过多级存储结构，将数据取到寄存器中开始运算，数据的预取可以很好地实现数据传输和计算的重叠，是提高数据空间局部性的性能优化方法.开始执行计算时，参与运算的矩阵和向量存放在内存中，开始的取数不会命中Cache，FT1000上从内存中取数到二级Cache的延迟大于100拍（即100个CPU时钟周期），Cache块大小为32字节，每次内存访问，同时会把内存中连续的32字节分别存放到二级Cache和一级数据Cache中.

数据预取就是将后面要用的数据提前取到L2 Cache中，将访存的数据传输操作和乘加计算操作重叠起来，如果计算时间大于数据传输时间，那么整个计算就可以得到CPU的峰值性能，如矩阵相乘.通过参数可以调整数据预取的间隔，获得最佳性能.

BLAS1，BLAS2和BLAS3均用到数据预取，本文仅以双精度AXPY为例来说明数据预取方法，图4给出了带数据预取的展开循环体.

AXPY带数据预取的展开循环体：

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

循环展开因子为m*8

数组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带数据预取的展开循环体

Fig.4 Multilevel loop unrolling with the prefetching

data for subroutine AXPY of BLAS1

FT1000通过预取指令对内存数据进行操作，SIZE表示一个数据所占的内存字节数，PREFECHSIZE表示提前预取数据间隔，可以根据需要改变大小.图4中，预取的数据提供给下一次循环体使用，将数据从内存中，取到二级Cache中，预取数的时间和图4中的两个小循环进行时间重叠，用计算重叠数据传输.

3 多线程并行计算方法

对BLAS1采用向量一维平均划分的方式组织并行计算.对BLAS2中涉及的二维矩阵采用按行或按列一维划分.对BLAS3中涉及的二维矩阵采用按行和按列二维划分.

下面重点对BLAS3中的矩阵乘并行计算方法展开讨论.

P 个线程按 pr×pc 划分成二维拓扑结构， 满足 P=pr×pc.假设 P（r， c） （0≤r

多线程矩阵乘并行算法：

全局数组

局部数组 r，c，

每一个线程P（r， c），其中 ，执行：

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实现，参考文献＼[12＼].

4 性能测试与分析

测试环境的硬件平台为8核FT1000微处理器，主频为1 GHz，双精度浮点性能8Gflops.操作系统为银河麒麟操作系统，编译器为gcc-4.5.1.

图6给出了BLAS1的双精度性能测试结果，固定向量的长度n不变，统一取为256 000 000.从图6可以看出，从1线程到8线程各BLAS1程序具有明显的加速效果，各程序在16或32线程时，达到最高性能，64线程时性能略有下降.造成32或64线程性能下降的主要原因是：1）BLAS1程序的计算访存比是1或2，受限于访存；2）BLAS1程序访存有空间局部性，没有时间局部性，也就是不存在数据的复用；3）FT1000微处理器的访存带宽在16或32线程时达到最大.

No of threads

图6 BLAS1不同线程数性能测试结果（n=256 000 000）

Fig.6 Computation performance for BLAS1

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

图7给出了BLAS2的双精度性能测试结果，固定矩阵的长度n不变，统一取为16 000.从图7可以看出，浮点性能明显高于BLAS1的性能，性能最好的SYMV在64线程时达到最高性能，3.11Gflops/s，是峰值性能的38.8%，SYMV的计算访存比是6，存在空间和时间局部性.BLAS2中的其他子程序的计算访存比是3，只有空间局部性，不存在时间局部性.

No of threads

图7 BLAS2不同线程数性能测试结果（n=16 000）

Fig.7 Computation performance for BLAS2

on difference number of threads（n=16 000）

BLAS3双精度浮点性能测试结果由图8给出，由于BLAS3的计算访存比大，对于阶为N的方阵，计算量为2N3，访存量为3N2，在64线程时获得最高性能，对不同的计算规模展开测试.从图8中可以看出，矩阵乘GEMM的最高性能达到6.91Gflops，是峰值性能的86.4%.SYMM，SYR2K，SYRK，TRMM和TRSM子程序的最高性能分别达到6.75，6.73，6.74，6.75和6.69Gflops，和矩阵乘GEMM的性能接近.

影响矩阵乘GEMM性能的主要因素包括：1）数据拷贝开销.为了充分发挥多级Cache的性能，采用了分块算法，分块子矩阵在内存中的存储是不连续的，为了减少对TLB冲突，需要将矩阵A和B的数据预先拷贝到一个连续的内存空间.2）数据延迟时间.每个子矩阵块相乘前，每个线程需要进行数据的填充，数据需要从内存中传输到L2 Cache，从L2 Cache到L1 Cache，从L1 Cache到寄存器，每个数据大概需要150拍，150拍过后，每拍可以取一个双精度数据.3）数据对存储带宽的竞争.64线程并行计算时，需要同时从内存中取数，数据带宽成为影响性能的因素.4）计算不能全覆盖数据的存取.FT1000微处理器中每核启动8线程，8线程共享一套硬件计算资源，靠多线程的切换来屏蔽访存延迟.这种情况主要存在每个循环展开启动阶段，8线程需要同时取数，此外循环展开计算完成后，需要把计算结果保存到矩阵C中，此时存取数据量和计算量处于相同量级.

M=N=K

图8 BLAS3不同计算规模性能测试结果（64线程）

Fig.8 Computation performance

of difference scale of matrixes on 64 threads

相对于其他X86微处理器，FT1000是多核多线体系结构，每8线程共享一个计算核，通过线程切换来屏蔽访存延迟，相对而言更难发挥其峰值性能，我们认为矩阵乘GEMM双精度能发挥峰值性能的86.4%已经是能达到的最好浮点性能.

5 结论与展望

提炼共性基础数值算法，研制高性能计算库，统一编程接口，是用户充分发挥微处理器峰值性能的重要手段.BLAS库是基础线性代数库，需要根据CPU的具体特点，设计高性能的算法和数据结构，经过高度的手工汇编优化，其中的BLAS3可得到接近峰值的浮点性能，满足用户的性能要求.

本文基于国防科学技术大学研制的多核多线FT1000微处理器，研制了基本线性代数汇编子程序库FitenBLAS，根据寄存器的数目，设计了向量与向量、矩阵与向量、矩阵与矩阵运算的多级循环展开方法，采用计算指令流水线调度、数据预取等通用优化技术，优化BLAS库串行程序.对于BLAS3子程序，设计了矩阵乘无冗余数据拷贝分块算法，采用指令重排、访存与计算的重叠、分块等技术优化矩阵乘子程序，基于矩阵乘子程序实现了其他BLAS3子程序.其核心子程序矩阵GEMM乘的双精度计算性能达到6.91Gflops，是峰值性能的86.4%.

下一步，面向短向量飞腾微处理器，研制支持向量优化运算的BLAS库，吸收并行编程框架底层调用的矩阵向量运算和稀疏线性数值算法，完善FitenBLAS库，支持更加广泛的数值模拟运算.

参考文献

[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.