基于多精度SIMD的线性方程组迭代细化求解

张茂全

（上海交通大学电子信息与电气工程学院，上海200240）

0 引言

随着神经网络技术的发展，多精度和低精度加速技术已经成为了十分重要的加速方法[1]。在求解线性方程组的过程中也可以使用多精度计算技术进行加速[2-3]。在求解线性方程组的过程中，主要使用计算机的浮点计算能力，例如单精度浮点（FP32）运算和双精度浮点（FP64）运算。在某些研究领域中，如图像处理、地震资料处理等领域，单精度计算足以满足用户对精度的需要。但对于其他复杂的科学计算，例如通过数值模拟预测天气，双精度浮点运算由于其高精度的优点已经成为了必不可少的组成部分。

在现代处理器架构中，支持多种精度的计算能力，并且较低精度的计算吞吐率通常要比高精度的吞吐率高。例如，单精度浮点乘法的能力通常是双精度浮点乘法的两倍，事实也确实如此，在英特尔最新的通用处理器中单精度浮点运算能力是双精度浮点运算能力的两倍。单精度相比于双精度不仅在计算速度上有优势，由于其位宽更小，其在数据存储、传输以及缓存命中率等方面都有优势。

本文设计了基于多精度SIMD（Single Instruction Multiple Data）[4]的线性方程组迭代细化求解方法，可以在获得相同精度的方程解的同时，减少处理器运算的时间。

1 线性方程组求解方法

1.1 直接求解法

在求解规模较小的线性方程组时，一般选取高斯消元法和后向回代法进行直接求解。高斯消元法其实就是一个初等行变换的过程，以一个具有3个未知数的线性方程组为例：

L2减去2*L1便可以消去L2中的x，L3加上L1便可以消去L3的x。此时的方程组变为：

然后再将L3加上3*L2便可将L3中的y消除。这样可以使整个方程变成一个上三角形的形式，L1中含有x、y、z三个未知数，L2中含有y和z两个未知数，L3中只含有z一个未知数。此时可以通过后向回代法，先根据L3求得z的值，然后再根据L2求得y的值，最后根据L1求得x的值。其中高斯消元法的时间复杂度为O(n3)，后向回代法的时间复杂度为O(n2)，对于规模较小且系数均为整数的方程组，通常只需进行一次高斯消元法和后向回代法便可直接求出最终的结果。

1.2 误差分析

在高斯消元法中有两个潜在的误差来源，第一个误差来源是浮点数的舍入误差[5]。在现代计算机中，使用IEEE754的浮点数格式来表示实数，浮点数是一种有限精度的数字表达方式，即相邻的浮点数之间存在着间隔，随着数字的增大，浮点数之间的间隔越大，处于间隔中的实数只能舍入到离它最近的浮点数，此时便产生了误差。第二个误差来源是由于浮点数的淹没，一个非常大的浮点数和一个正常大小的浮点数相加，会导致后者被淹没。例如：(1 +2100)-2100=0，而交换一下运算顺序，1+( 2100-2100)=1，因此在浮点加法运算中不满足结合律。为了避免淹没问题，在消元的过程中可以进行主元的交换，即改变未知数的消去顺序。

误差的结果主要分为两种：前向误差和后向误差。前向误差表示已经得到的解x'与真实解x之间的差异，其大小采用无穷范数表示‖x'-x‖∞；后向误差为残差r的无穷范数，其中残差r的计算方法为所示：

r=b-Ax'，一般认为，当后向误差小于一定数值时，便可认为当前的解为最终的结果。

1.3 迭代求解法

在求解规模较大的线性方程组Ax=b时，基于迭代细化（Iterative Refinement）求解的方法是比较常用的求解方法[6-7]。迭代细化求解方法首先是通过基于高斯消元法的LU分解（DGEFA）将线性方程组的系数矩阵A分解为一个上三角矩阵U和一个下三角矩阵L，原来的方程即变为LUx=b，此时由于三角矩阵的特性，便可以通过后向代入求解（DGESL）求得x的解[8]。但是，在多数情况下，由于浮点数精度的限制，此时求得的方程解可能与精确解之间具有一定的误差，为了尽可能地使误差最小，需要通过多次迭代的方法使得最终的结果非常逼近精确解。

如算法1所示，首先对系数矩阵A进行LU分解，将系数矩阵A分解为上三角矩阵U和下三角矩阵L，分解的时间复杂度为O(n3)。将A分解为L和U的好处是，无论右侧的b向量是否改变，或者有多个b向量时，LU分解只需要进行一次，针对不同的b向量只需要各自进行后向回代操作，这大大减少了运算的步骤，因为后向回代求解过程的时间复杂度为O(n2)，与LU分解的时间复杂度相比是微不足道的。第一次后向回代求解到的x0只是一个初步解，距离精确解还有一定的误差。因此，接下来就需要进行迭代求解来获得更加准确的方程解。在反复迭代求解的过程中，首先计算残差r，接下来便等同于求解纠正方程Ad=r，求解时只需进行后向回代，因为L和U已经得到。将求解得到的d作为修正项，加到已经得到解x'上，即新的解x=x'+d。当新的方程解的后向误差小于一定数量级时，便可以认为新的方程解已经达到了精度的要求，可以作为最终的结果，此时便可以退出循环。

算法1线性方程组的迭代细化求解算法

输入：系数矩阵A，常数项向量b

输出：方程组的解x

在整个迭代循环中，只在循环外进行了一次最耗时的LU分解，在循环内部主要进行的是时间复杂度更小的后向回代和向量的加减法过程。因此，在迭代次数远小于矩阵尺寸n的情况下，整个迭代求解过程中的时间开销主要集中在刚开始的LU分解过程。因此，如果我们能够采取某种方法来加速LU分解的过程，那么也会给整个迭代求解过程带来明显的速度提升。

2 多精度SIMD迭代细化求解

2.1 多精度迭代细化求解

文献[9]的研究表明，可以使用两种浮点精度或三种浮点精度来加速求解线性方程组。在对系数矩阵A进行LU分解时，使用速度更快的单精度浮点乘法，在进行代入求解时，使用精度更高的双精度浮点乘法。由于LU分解在整个求解过程中占据的时间比重更大，对其使用单精度进行计算时，可以使整个求解过程获得显著的加速。因此，可以通过使用高精度和低精度SIMD乘法指令，并通过迭代细化的方法来加速线性方程组Ax=b的求解。

如算法2所示，在使用多精度计算线性方程组的算法中，涉及到三种计算精度：①u是原始数据存储精度的unit-roundoff（单位舍入误差），②uf是LU分解时使用的数据精度的舍入误差，③ur是计算残差时的数据精度的舍入误差。为了保证求解过程最终收敛，这三种精度需要满足以下关系（精度越高、舍入误差越小）：ur≤u≤uf。最简单的情形便是u=uf=ur，即从始至终使用一种精度，例如只使用双精度。也可以使用三种不同的精度ur=FP64，u=FP32，uf=FP16。

算法2基于多精度的线性方程组的迭代细化求解

输入：系数矩阵A(FP64)，常数项向量b(FP64)

输出：方程组的解x(FP64)

2.2 精度选择

如表1所示，双精度（FP64）、单精度（FP32）和半精度（FP16）可以表示的最大值呈指数性下降[11]。因此将双精度（FP64）或单精度（FP32）数据转换为半精度（FP16）时，非常容易发生上溢或下溢。

表1 IEEE754浮点数特性

文献[9]的研究表明在使用多精度求解线性方程组时，系数矩阵A的条件数需要满足一定的条件，才能使整个迭代过程收敛。在ur=FP64、u=FP32、uf=FP16时，只有系数矩阵A的条件数κ∞(A)<104时才能保证收敛，条件数描述了矩阵的病态情况[10]；在ur=FP64、uf=u=FP32时，系数矩阵A的条件数κ∞(A)<108时便能保证收敛。

普通矩阵的条件数有很大的概率超出104，此时如果选取ur=FP64、u=FP32、uf=FP16，可能会因为系数矩阵A的条件数大于104导致迭代求解过程无法收敛。所以，采用ur=FP64、uf=u=FP32进行求解线性方程组时既能保证迭代求解过程收敛，也能获得一定的性能收益。

2.3 多精度SIMD加速LU分解

LU分解是高斯消元法的矩阵形式，通过将高斯消元法转化成矩阵运算的形式，可以更好地进行并行加速。原始的系数矩阵A可以看成下三角矩阵L和上三角矩阵U的乘积，由于L和U均有一半的元素为零元素，因此可以将下三角矩阵L和上三角矩阵U储存在一个n×n的矩阵中，即和原来的系数矩阵的大小相同，而不占用更多的存储空间。

如图1所示，以一个4×4的LU分解为例，四个未知数分别为x,y,z,w。系数矩阵中的每一列表示某个未知数在整个方程组中的系数，因此第一列就表示x在4个方程中的系数。根据上一节介绍的高斯消元法的步骤，首先需要消去第2-4行的x，因此我们需要将a2,1、a3,1、a4,1变为0。保持第一行的数不变，将第一行的数乘以一个特定的标量a并加到第2-4行，为了使a2,1、a3,1、a4,1为0，a应该分别等于。随后保持第二行不变，将第二行的数分别乘以、，依次类推，直至最后一行只有一个未知数，此时便完成LU分解。此时原系数矩阵便成为了最终的上三角矩阵U，其对角线以下的元素均为0。下三角矩阵L其实是记录了每次消元操作的标量a，以消去x为例，a2,1、a3,1、a3,1均会变为0，在后向代入求解时，由于我们知道a2,1、a3,1、a3,1已经为0，不会使用a2,1、a3,1、a3,1的值，所以我们可以将每次的标量a存储在即将要被消去的系数的位置。

图1 LU分解示意图

将一行数乘以一个系数加到另一行数上的操作其实是axpy操作，由于axpy中的各个子操作没有数据依赖，因此可以很容易地使用并行加速方法对其加速。典型的数据级并行加速方法为使用SIMD指令，即单指令多数据（Single Instruction Multiple Data），使用一条指令可以同时对多个元素进行同样的操作。

如图2所示，SIMD指令可以一次对多个操作数进行相同的运算，而普通的指令只能对一组数据运行运算。因此，可以使用单精度SIMD乘法指令和单精度SIMD加法指令来执行LU分解中的axpy操作。单精度是双精度浮点数宽度的一半，因此在不增加向量寄存器的情况下，对一组64位的寄存器进行单精度SIMD乘法或加法运算时，可以同时计算两组单精度乘法或加法，此时axpy需要的乘法和加法指令的数目将减少一半。

图2 SIMD加法指令示意图

3 实验仿真与结果分析

本文的实验主要分为两部分，第一部分为使用MATLAB软件来计算多精度对迭代求解次数的影响。第二部分为使用具有多精度SIMD浮点乘法能力的RISC-V[12]处理器硬件仿真平台来运行具有多精度SIMD指令的汇编程序来观察多精度SIMD指令带来的性能提升。

3.1 多精度迭代次数

在MATLAB软件中使用lu()函数获得一个矩阵的下三角矩阵和上三角矩阵，使用左除符号“”来实现后向代入求解，L其实等价于L-1b，左除符号首先会检查矩阵的形状，如果矩阵是三角矩阵，那么可以直接进行代入求解[13]。首先，只使用双精度浮点数来进行完整的LU分解操作和后向代入求解操作，记录当残差足够小时的迭代次数作为基准数据。随后使用single()函数将双精度的A矩阵转换为单精度浮点格式的矩阵，然后使用lu()函数对其进行LU分解操作，使用双精度浮点数来计算残差和纠正方程，记录残差足够小时的迭代次数。

如表2所示，首先选取了5个在实际应用中产生的数据作为系数矩阵，假定残差的无穷范数小于10-12时，便可认为当前循环中的方程解满足精度要求。可以看到在使用多精度算法进行求解时的迭代次数都要比只使用双精度求解时的迭代次数多，但还是相同的数量级，没有大幅的增加。

表2 多精度求解对迭代次数的影响

3.2 多精度SIMD乘法加速结果

为了验证SIMD乘法指令对LU分解和后向回代操作的加速效果，首先使用Verilog HDL实现了可以计算双精度或同时计算两个单精度（SIMD）的浮点乘法器，该乘法器完全符合IEEE754标准，支持向偶舍入、非规格化数、NaN等。

为了能够运行完整的汇编程序，将具有多精度计算能力的浮点乘法器应用到了一个具有五级流水线单发射的RISC-V处理器的RTL模型中，该处理器支持数据转发、静态分支预测等技术来解决数据冲突和控制冲突，从而减少流水线阻塞情况的发生。为了能够更好地利用数据的局部性[14]，该处理器中配置有16KB大小、两路组相连的L1 Cache，每个cache block的大小为64B。为了支持SIMD乘法指令，需要在译码阶段增加对SIMD乘法指令的译码支持，为了能够对数据进行精度转换，增加了精度转换指令的实现，分别为高精度转为低精度和低精度转为高精度，并且可以将转换后的数据送入目标寄存器。

由于处理器中的cache block大小为64B，可以存储16个单精度数据或8个双精度数据，为了更好地利用时间局部性，减少cache block的替换次数，在进行具体的LU分解操作时，采用了分块的思想：由于LU分解操作实际是将某一行所有数的倍数分别加到剩下的所有行（axpy），同一行的各个数据之间也没有数据相关性，以图3为例，在消去第一列除a1,1外其他的所有数时，可以先在阴影区域内执行axpy操作，因为这样可以保持a1,1-a1,16的数据始终在cache内并命中，如果拥有类似于ARM neon指令集中的向量寄存器，甚至可以确保a1,1-a1,16的数据可以始终驻留在寄存器内，不发生寄存器溢出。

图3 LU分解的分块算法示意图

分别编写了用于LU分解的DGEFA（双精度）和SGEFA（单精度）汇编程序，以及用于后向回代求解的DGESL汇编程序，在SGEFA中使用了SIMD单精度乘法和SIMD单精度加法指令，一次可以进行两个乘法或加法操作。将以上的汇编程序汇编为机器码后加载到RISC-V处理器的instruction-rom中，使用Synopsys VCS软件对RISC-V处理器和多精度浮点乘法器的Verilog文件进行硬件仿真实验。在RISC-V的顶层模块设置有一个计数器用来记录从程序开始到其结束的时钟周期数，在汇编程序结尾添加一条汇编语句AD-DI x31,x0,0作为结束标志，当检测到指令为结束标志时，退出仿真并打印运行周期数。

如表3所示，通过对加载了汇编指令的Verilog文件进行仿真，得到了DGEFA、SGEFA、DGESL的周期数。使用双精度进行求解的总周期数为DGEFA+i×DGESL，使用多精度求解的总周期数为SGEFA+j×DGESL，其中i,j分别为表2中的双精度和多精度的迭代次数。

表3 汇编程序仿真周期结果

如图4所示，使用表3中双精度的周期数除以多精度的周期数便可以得到使用多精度SIMD方法对求解线性方程组带来的加速比，系数矩阵是bcsstk04时的加速比最小，为1.463倍，系数矩阵是steam1时的加速比最大，为1.612倍。

图4 多精度SIMD求解线性方程组加速比

4 结语

针对求解较大规模的线性方程组，本文提出了基于多精度SIMD的加速求解方法，使用低精度的SIMD乘法指令来加速时间复杂度最高的LU分解步骤，使用高精度乘法指令计算残差和纠正方程。通过相应的仿真实验，探讨了多精度计算对迭代次数的影响，以及对求解过程的加速效果。针对不同规模的系数矩阵，该方法的性能可以达到传统方法的1.46-1.61倍，具有较好的加速效果。