基于多精度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倍，具有較好的加速效果。