基于龙芯3B的循环规约算法向量化研究

基于龙芯3B的循环规约算法向量化研究
吴淅;黄章进;顾乃杰
【摘要】By analyzing the unique characteristics of reduction cyclic algorithm for solving tri-diagonal linear system, this paper presents an optimized algorithm specially for Godson 3B processor. The reduction cyclic algorithm is implemented on Godson 3B platform based on its vector instructions and vectorization techniques are used to boost performance of the algorithm. Experimental results show that, comparied with non-vectorization algorithm, vectorization algorithm larger enhances the performance.%研究国产CPU龙芯3B的体系结构,分析三对角线性方程组求解中的循环规约算法,并结合算法特性,提出循环规约算法的优化方案.采用向量化级数,利用龙芯3B处理器提供的向量扩展指令对该算法在龙芯3B上进行实现和优化.实验结果表明,与非向量化算法相比,向量化算法的性能提升较大.
【期刊名称】《计算机工程》
【年(卷),期】2013(039)002
【总页数】6页(P293-298)
【关键词】龙芯3B;三对角线性方程;循环规约算法;向量化;泊松方程;离散化
【作者】吴淅;黄章进;顾乃杰
【作者单位】中国科学技术大学计算机科学与技术学院,合肥230027;安徽省计算与通讯软件重点实验室,合肥230027;中国科学技术大学与中国科学院沈阳计算技术研究所网络与通信联合实验室,合肥230027;中国科学技术大学计算机科学与技
术学院,合肥230027;安徽省计算与通讯软件重点实验室,合肥230027;中国科学技
术大学与中国科学院沈阳计算技术研究所网络与通信联合实验室,合肥230027;中
国科学技术大学计算机科学与技术学院,合肥230027;安徽省计算与通讯软件重点
实验室,合肥230027;中国科学技术大学与中国科学院沈阳计算技术研究所网络与
通信联合实验室,合肥230027
【正文语种】中文
【中图分类】TP301.6
1 概述
三对角线性方程组的求解在数学和物理学中有广泛的应用。

它被大量地使用在结构力学、热力学和流体力学的分析中[1]。

同时，对于数值计算中的许多问题，特别
是线性扩散问题，其中的代数方程矩阵通常为三对角，或者最终会演化成三对角矩阵。

因此，怎样高效地实现三对角线性方程组的求解，尤其是大规模数据量的情况，是现在的热点问题之一。

本文针对龙芯3B处理器提供的扩展向量指令[2]，对求解三对角矩阵的循环规约算法[3]进行了向量化。

同时，还将该向量化算法实际运用在求解泊松方程中，提高
方程的求解速度。

2 相关工作
龙芯3B处理器在兼容了MIPS64指令集的同时，实现了针对大规模数据的向量扩展指令。

其设计的向量运算部件，支持SIMD机制，提供了256位向量寄存器，
并实现了包括256位向量访存的向量扩展指令。

这样的结构和指令集设计，使得
龙芯 3B适用于大规模数据处理，比如视频编/解码运算、FFT运算，以及求解大
规模稀疏线性方程组等[4]。

龙芯 3号处理器主要面向服务器和高性能应用[5]，随
着国产龙芯系列不断推出，势必走向产业化。

硬件产品必须配套系统软件和应用软件的支持才能最终被用户所用。

为用户提供面向龙芯体系结构的高性能数学函数库是必不可少的。

求解三对角方程组的基础算法为高斯消去法，该算法首先从前向后逐步更新系数，然后从后向前回代解出未知数。

它是最简单易懂的一种解法，但前后依赖性大，属于基本的串行算法。

从 1960年开始，各种并行度较大的算法被逐渐提出。

其中比较有代表性的有循环规约算法[3]、递归倍增算法[6]和并行循环规约算法等。

文献[3]提出循环规约算法。

文献[7]将该算法应用在GPU上，利用GPU着色语言实现了三对角线性方程组的求解。

文献[8]使用CUDA实现了循环规约算法，并将它应用在GPU上对图形学中的实时浅水流体模拟中。

文献[9]研究了在GPU上实现循环规约算法、并行循环规约算法、递归倍增法以及这3种算法的混合变种算法，并提出一种针对GPU上运行算法的访存问题、计算量以及控制开销等方面的优化方法，同时给出测试和分析方法。

目前为止，有关循环规约算法的研究主要集中在并行实现上，对其的向量化实现较少。

由于龙芯3B平台提供了完善的向量扩展指令，可以将上述算法进行向量化，且在本文所做的工作之前，还没有该算法在该平台上的实现以及优化。

因此，本文以循环规约算法为例，充分发掘平台特性，使用向量指令集，高效地实现该算法，为数学库中的其他算法在龙芯 3B上的向量化实现提供参考。

3 算法向量化的求解
循环规约算法由前向约化和后向回代 2个过程组成[3]。

其基本思想是在约化过程中，每次将待解方程组规模减半，直到规模为 2。

回代过程使用前向约化过程中得到的未知数值，回代解出当前规模下的另一半未知数[8]。

未知数个数为18的循环规约算法如图1所示。

图1 未知数个数为18的循环规约算法
3.1 循环规约算法
三对角线性方程组，形如Ax=d，其中，系数矩阵A是一个三对角矩阵，如下所示:
约定本文中所有方程和所有变量的下标全部从1开始，更新之后的数据以原数据右上方标注“'”表示。

在前向约化过程的每次迭代中，将当前规模方程组中偶数编号方程的所有奇数编号变量消去，只留下偶数编号变量。

其中，更新的偶数编号变量由其前一个和后一个奇数编号变量得到。

更新系数矩阵a[]、b[]、c[]、d[]的方法见式(2)。

其中，a'[i]为a[i]更新之后的值，依次类推：
若编号为i的变量为最后一个，即不存在i+1编号变量，则更新方式为式(3):
在后向回代的每次迭代中，将前面已解的未知数，即前一次迭代中得到的偶数编号未知数的值回代到本次计算过程中，得到本次待解方程组的奇数编号未知数的值。

式(4)为计算编号为i的未知数xi的方法。

根据待求解未知数xi在当前规模下所处的位置:最后一个、第一个和中间的某一位，分别使用式(4)的第一行、第二行和最后一行。

3.2 在龙芯3B上的向量化
本实验中的系数及未知数全部为双精度浮点数，一条指令同时处理4个双精度浮点数据。

对于处于边界的数据，由于其更新计算方式不同，不可与其他数据一起在一条指令中操作，且每次循环只有一个边界数据，因此，不使用向量指令处理。

从图2可较清晰地看出，前向约化每次迭代中，每个待更新的变量只与变量本身、其前一个和后一个变量相关，与当次迭代中其他待更新变量没有任何依赖关系。

在后向回代的每次迭代中，待求解的未知数只依赖自己以及上一轮已求解得到的未知数，而与其他待求解未知数不相关。

图2 向量指令的重整计算操作
假设某一次迭代中待更新数据之间的步长为stride，例如图1上半部分中第一次
迭代时，待更新数据e[1]2和e[1]4之间的步长为stride=2，分别位于原数组的
第2个和第 4个位置，更新 e[1]2需要使用上一轮的数据 e1、e2和e3，更新
e[1]4需要使用上一轮的数据e3、e4和e5。

更新后得到的新数据替换掉数组中更新之前的数据。

使用向量指令可同时操作4个数据，假设首个数据为原数组的第 i 个数据，则一条向量指令同时操作的数据是原数组的第i、i+stride、i+2stride、
i+3stride位置上的数据，假设为e[i]、e[i+ stride]、e[i+2stride]、e[i+3stride]。

这4个数互相之间没有依赖。

e’[i]为e[i]更新后的值，更新完毕后写回数组，覆
盖掉数组中原始数据e[i]。

依次类推。

如图2中上半部分所示，符号泛指计算操作。

因此，更新 e[i]、e[i+stride]、e[i+2stride]、e[i+ 3stride]这4个数时，首先判
断e[i]和e[i+3stride]是否为临界数据。

若是临界数据，则单独处理；若为非临界
数据，则放在同一条向量指令中处理。

一次读取
e[i−1/2stride],e[i],…,e[i+7/2stride]共9个数据，将这9个数据分别读进3个向
量寄存器$z0、$z1和$z2中，然后使用一条指令对向量寄存器$z0、$z1、$z2进行操作:$z3＜—$z0$z1$z2，即可完成上述4个数共涉及9个数据的更新，如图2中下半部分所示。

计算得出的数据存放在向量寄存器$z3中，并从$z3写回到数组
中的相应位置。

前向约化过程中每次迭代更新系数矩阵a、b、c、d这4个数组共需要至少4×3=12个向量寄存器。

后向回代的向量化过程与上述类似，当系统规模大于4时开始进行向量化。

使用
向量指令做计算这一过程一直进行到方程组规模小于4时停止。

此时系统规模过小，且存在边界数据，需要特殊处理，因此，不适合继续使用向量指令操作。

理论上，若向量化部分的数据恰好是4的倍数，则可将每次迭代中的循环次数降
到原来的1/4。

现实中，不可能每次迭代的操作数个数都是4的倍数，比如边界数据的更新方式不同可能会导致连同边界数据在内最多4个数据不能同时在一条指
令中处理，而只能以流水的顺序处理，使得整体的处理速度下降。

此外，随着循环展开后的循环体内指令数目的增多，将不存在操作数相关性和不存在运算部件相关性的指令调度在一起，通过合理的指令调度，从而充分应用龙芯3B的流水线特性，达到提高代码在龙芯3B上的执行速度。

4 算法实现及实验结果分析
本课题在龙芯 3B平台上实现了上节提出的向量化CR算法，并与未向量化的结果进行了比较，给出了不同规模情况下各自的计算时间以及向量化算法的加速比。

最后，根据实验数据给出实验结果分析。

4.1 算法实现
本实验测试用系统规模(即未知数个数)从100～107，实验结果如表1和图3所示。

加速比处于2～3.5之间，加速比的最大值为3.5，出现在系统规模为2×105时。

表1 不同未知数下算法运行时间及加速比未知数个数向量化算法运行时间/ms非
向量化算法运行时间/ms 加速比102 0.035 0.055 1.57103 0.254 0.555 2.19104
2.767 7.407 2.68105 1.104×102
3.765×102 3.41106 1.375×103
4.316×103
3.14107 1.256×104 1.871×104 1.49
图3 不同规模下算法的时间曲线及加速比
4.2 实验结果分析
在本文实验中，可向量化部分在整个算法中所占比例为80%。

在可向量化部分中，每次循环可能都会遇到边界数据，不可使用向量指令处理，而整个可向量化部分的循环中有75%次必须独立处理边界数据，以上的原因都会造成加速比的下降。

(1)访存分析
假设从内存中取数需要L拍，从一级cache中取数需要C1拍，从二级cache取
数需要C2拍，一次乘法运算需要M拍，一次加法运算需要A拍，一次除法运算
需要D拍，一次赋值运算需要E拍，将数据写回内存需要S拍，数据移位整合操
作需要V拍。

龙芯3B各项参数值为L∶100, C1∶2, C2∶5, M∶6, A∶1, D∶20, E∶1, S∶100, V∶18。

在前向约化过程中，可向量化部分的核心为64位浮点运算的d[4]=a[4]b[4]c[4]。

其中，表示计算操作，共包括2个双精度浮点除法、6个浮点乘以及6个浮点加法。

对于向量化方法，4个数组a、b、c、d分别需要从内存中取4个数据进行一次向量计算。

以数组a为例，读取数据时由于一次读取一个Cache行(32 Byte)，正好为4个double型数据，因此只需读取一次，则从L2 Cache中读数到L1 Cache，并从L1 Cache读数到寄存器中需要C2+C1个时钟周期。

其他数组情况类似，则a、b、c、d这 4个数组从 L2 Cache读取数据到寄存器共需要4×(C2+C1)个时
钟周期。

在寄存器中对数据进行乘加除计算和赋值需要6×M+6×A+2×D+E个时钟周期。

将得到的结果数据写回需要C2个时钟周期。

对于非向量化算法，以数组a为例，读取第一个数据到寄存器的时候，由于是按
照Cache行读取，而Cache行大小为32 Byte(相当于4个double型数据)，从
而剩下的 3个数据都存在于 L1 Cache中，则下次直接从 L1 Cache中读取。

因此，对于数组a，其读取4个数据需要C2+3×C1个时钟周期。

其他数组情况类似，则
a、b、c和d这4个数组从内存中取数共需要4×(C2+3×C1)个时钟周期。

4个数据进行乘加和赋值操作需要4×(6×M+6×A+2×D+E)个时钟周期。

4个数据写回
内存需要4×C2。

因此，当数据已在二级Cache中时，非向量化和向量化时间分别如式(5)和式(6)所示。

将龙芯3B指令周期代入公式，得到time_nonv1=524, time_vec1=189，加速比为2.77。

在后向回代过程中，向量化部分的核心为d[4]=a[4]b[4]c[4]，其中，为计算操作。

共包括1个双精度浮点除、2个浮点乘和2个浮点加法。

当数据已在二级 cache中时，非向量化和向量化时间与前向约化过程相似，分别
如式(7)和式(8)所示。

将龙芯 3B指令周期代入公式，得到 time_nonv2=332, time_vec2=141，加速比为2.36。

实验包括前向约化和后向代回2个部分，考虑到计算与访存存在计算重叠部分，
实际加速比最大值为3.41，与理论最大值接近。

(2)峰值分析
在方程组规模n达到2×105(lgn=5.3)时，出现最大的加速比为3.5。

龙芯3B的
二级Cache大小为4 MB=4×1024 KB，当数据全部处于二级Cache中时，此时数据访存时间最少，整体处理速度最大。

本文实验中需要访存4个数组，数组元
素为双精度浮点数，占用4个字节。

假设系统规模为u，即数组元素个数为u，则4个数组占用4×4u=16u Byte。

当数组占用内存16u Byte小于二级Cache大小
4 MB时，所有数据都在二级Cache中，没有从内存读取数据的访存开销。

由
16u＜4×1024×1024得到，u＜262144。

理论上，系统规模超过262144时，加速比开始下降。

当系统规模过小时，向量化操作进行的数量较少，非向量化的比例占的比较高，因此，向量化带来的加速相对有限，不能完美地体现向量化的优越性。

在实际实验过程中，操作系统中其他进程的数据可能占用二级 Cache，造成实际
系统规模不可能达到理论的262144，加速比就已经开始下降。

实验中加速比峰值出现时的系统规模为200000，接近理论值。

5 算法的实际应用
三对角线性方程组的求解是使用有限差分法求解泊松方程的核心步骤。

泊松方程如式(9)所示[10]:
对x域进行离散化均匀采样，每个采样点之间的距离为h，假设总共为n个采样点，则h=1/n。

离散化后的泊松方程可表示成为矩阵向量乘的形式，形如 Av=f，v是一个一维向量v=(v1,v2,…,vn-1)T。

其中，vj表示 u(xj)的逼近值。

f=f(f(x1),f(x2),…,(xn−1))T=(f1,f2,…,fn-1)T 表示各个采样点处的方程右值，A是
一个(n−1)× (n−1)维三对角、正定矩阵。

如式(10)所示:
本文运用上节实现的求解三对角矩阵的向量化方法，求解泊松方程。

方程右值 f(x)及解析解 u(x)分别在式(11)和式(12)中给出。

表2给出不同离散步长h下实验得出的数值解与理论解析解之间的最大误差。

从
表中可看出，步长越小，即采样点越多，最大误差越小。

当步长小于0.2×10-3，即采样点数大于5000时，最大误差为0。

表2 不同离散步长下的最大误差离散步长h 最大误差10−2 3.1×10−3 0.2×10−2 1.2×10−4 10−3 3.0×10−5 0.2×10−3 010−4 0
图4给出了步长h为10−3，即采样点数为1000时，u(x)逼近值的数值解曲线和误差曲线。

其中误差由数值解与解析解相减得到，从图 4可看出，误差曲线接近0，说明两者误差很小。

图4 采样点数为1000的泊松方程数值解及误差曲线
在本文实验中，采样点数为1000时，使用本文提出的求解三对角线性方程组向量化算法求解泊松方程与不使用该向量化算法求解方程的时间分别为231 μs和545 μs，加速比为 2.3。

6 结束语
国产龙芯系列处理器正在走向市场化，为用户提供面向龙芯体系结构特点的高性能数学函数库是必不可少的。

求解三对角矩阵是数学、物理等领域的一个基本问题，优化其性能有着重要的意义。

本文基于龙芯3B处理器，针对其实现了对向量扩展指令支持的特点，对求解三对角矩阵的循环规约算法进行了向量化，并合理安排指令调度。

最后，将此向量化算法运用在求解泊松方程中，并取得了一定的性能提升。

若要进一步提高性能，今后还需要在算法层面做优化。

同时，由于龙芯3B的多核特性，可以考虑使用多核进行方程组求解。

参考文献
[1]Wang H H.A Parallel Method for Tridiagonal Equations[J].ACM Transactions on Mathematical Software, 1981, 7(2)∶170-183.
[2]姜彩霞.基于龙芯的 GCC自动向量化移植与优化[D].北京∶ 中国科学院计算技术研究所, 2009.
[3]Gander W, Golub G H.Cyclic Reduction——History and
Applications[C]//Proc.of the Workshop on Scientific Computing.Hong Kong, China∶ [s.n.], 1997∶ 73-85.
[4]裴晓航, 何颂颂.基于龙芯 3B的 H.264解码器的向量化[J].电子技术, 2010,
37(10)∶ 88-90.
[5]龙芯3B处理器用户手册V1.0版[Z].北京∶ 中国科学院计算技术研究所, 2011.
[6]Stone H S.An Efficient Parallel Algorithm for the Solution of a
Tridiagonal Linear System of Equations[J].Journal of the ACM, 1973, 20(1)∶ 27-38.
[7]Kass M, Lefohn A, Owens J D.Interactive Depth of Field Using Simulated Diffusion[R].Pixar Animation Studios,Technical Report∶ 06-01, 2006. [8]Sengupta S, Harris M, Zhang Yao, et al.Scan Primitives for GPU Computing[C]//Proc.of the 22nd ACM SIGGRAPH/EUROGRAPHICS Symposium on Graphics Hardware.[S.l.]∶ ACM Press, 2007∶ 97-106. [9]Zhang Yao, Cohen J, Owens J D.Fast Tridiagonal Solvers on the
GPU[C]//Proc.of the 15th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming.New York, USA∶ ACM Press, 2010∶127-136.
[10]Briggs W L, Henson V E, McCormick S F.A Multi-grid Tutorial[M].2nd ed.Beijing, China∶ Tsinghua University Press, 2011.。