基于多精度SIMD的線性方程組迭代細(xì)化求解

張茂全

（上海交通大學(xué)電子信息與電氣工程學(xué)院，上海200240）

0 引言

隨著神經(jīng)網(wǎng)絡(luò)技術(shù)的發(fā)展，多精度和低精度加速技術(shù)已經(jīng)成為了十分重要的加速方法[1]。在求解線性方程組的過程中也可以使用多精度計(jì)算技術(shù)進(jìn)行加速[2-3]。在求解線性方程組的過程中，主要使用計(jì)算機(jī)的浮點(diǎn)計(jì)算能力，例如單精度浮點(diǎn)（FP32）運(yùn)算和雙精度浮點(diǎn)（FP64）運(yùn)算。在某些研究領(lǐng)域中，如圖像處理、地震資料處理等領(lǐng)域，單精度計(jì)算足以滿足用戶對(duì)精度的需要。但對(duì)于其他復(fù)雜的科學(xué)計(jì)算，例如通過數(shù)值模擬預(yù)測(cè)天氣，雙精度浮點(diǎn)運(yùn)算由于其高精度的優(yōu)點(diǎn)已經(jīng)成為了必不可少的組成部分。

在現(xiàn)代處理器架構(gòu)中，支持多種精度的計(jì)算能力，并且較低精度的計(jì)算吞吐率通常要比高精度的吞吐率高。例如，單精度浮點(diǎn)乘法的能力通常是雙精度浮點(diǎn)乘法的兩倍，事實(shí)也確實(shí)如此，在英特爾最新的通用處理器中單精度浮點(diǎn)運(yùn)算能力是雙精度浮點(diǎn)運(yùn)算能力的兩倍。單精度相比于雙精度不僅在計(jì)算速度上有優(yōu)勢(shì)，由于其位寬更小，其在數(shù)據(jù)存儲(chǔ)、傳輸以及緩存命中率等方面都有優(yōu)勢(shì)。

本文設(shè)計(jì)了基于多精度SIMD（Single Instruction Multiple Data）[4]的線性方程組迭代細(xì)化求解方法，可以在獲得相同精度的方程解的同時(shí)，減少處理器運(yùn)算的時(shí)間。

1 線性方程組求解方法

1.1 直接求解法

在求解規(guī)模較小的線性方程組時(shí)，一般選取高斯消元法和后向回代法進(jìn)行直接求解。高斯消元法其實(shí)就是一個(gè)初等行變換的過程，以一個(gè)具有3個(gè)未知數(shù)的線性方程組為例：

L2減去2*L1便可以消去L2中的x，L3加上L1便可以消去L3的x。此時(shí)的方程組變?yōu)椋?/p>

然后再將L3加上3*L2便可將L3中的y消除。這樣可以使整個(gè)方程變成一個(gè)上三角形的形式，L1中含有x、y、z三個(gè)未知數(shù)，L2中含有y和z兩個(gè)未知數(shù)，L3中只含有z一個(gè)未知數(shù)。此時(shí)可以通過后向回代法，先根據(jù)L3求得z的值，然后再根據(jù)L2求得y的值，最后根據(jù)L1求得x的值。其中高斯消元法的時(shí)間復(fù)雜度為O(n3)，后向回代法的時(shí)間復(fù)雜度為O(n2)，對(duì)于規(guī)模較小且系數(shù)均為整數(shù)的方程組，通常只需進(jìn)行一次高斯消元法和后向回代法便可直接求出最終的結(jié)果。

1.2 誤差分析

在高斯消元法中有兩個(gè)潛在的誤差來源，第一個(gè)誤差來源是浮點(diǎn)數(shù)的舍入誤差[5]。在現(xiàn)代計(jì)算機(jī)中，使用IEEE754的浮點(diǎn)數(shù)格式來表示實(shí)數(shù)，浮點(diǎn)數(shù)是一種有限精度的數(shù)字表達(dá)方式，即相鄰的浮點(diǎn)數(shù)之間存在著間隔，隨著數(shù)字的增大，浮點(diǎn)數(shù)之間的間隔越大，處于間隔中的實(shí)數(shù)只能舍入到離它最近的浮點(diǎn)數(shù)，此時(shí)便產(chǎn)生了誤差。第二個(gè)誤差來源是由于浮點(diǎn)數(shù)的淹沒，一個(gè)非常大的浮點(diǎn)數(shù)和一個(gè)正常大小的浮點(diǎn)數(shù)相加，會(huì)導(dǎo)致后者被淹沒。例如：(1 +2100)-2100=0，而交換一下運(yùn)算順序，1+( 2100-2100)=1，因此在浮點(diǎn)加法運(yùn)算中不滿足結(jié)合律。為了避免淹沒問題，在消元的過程中可以進(jìn)行主元的交換，即改變未知數(shù)的消去順序。

誤差的結(jié)果主要分為兩種：前向誤差和后向誤差。前向誤差表示已經(jīng)得到的解x'與真實(shí)解x之間的差異，其大小采用無窮范數(shù)表示‖x'-x‖∞；后向誤差為殘差r的無窮范數(shù)，其中殘差r的計(jì)算方法為所示：

r=b-Ax'，一般認(rèn)為，當(dāng)后向誤差小于一定數(shù)值時(shí)，便可認(rèn)為當(dāng)前的解為最終的結(jié)果。

1.3 迭代求解法

在求解規(guī)模較大的線性方程組Ax=b時(shí)，基于迭代細(xì)化（Iterative Refinement）求解的方法是比較常用的求解方法[6-7]。迭代細(xì)化求解方法首先是通過基于高斯消元法的LU分解（DGEFA）將線性方程組的系數(shù)矩陣A分解為一個(gè)上三角矩陣U和一個(gè)下三角矩陣L，原來的方程即變?yōu)長(zhǎng)Ux=b，此時(shí)由于三角矩陣的特性，便可以通過后向代入求解（DGESL）求得x的解[8]。但是，在多數(shù)情況下，由于浮點(diǎn)數(shù)精度的限制，此時(shí)求得的方程解可能與精確解之間具有一定的誤差，為了盡可能地使誤差最小，需要通過多次迭代的方法使得最終的結(jié)果非常逼近精確解。

如算法1所示，首先對(duì)系數(shù)矩陣A進(jìn)行LU分解，將系數(shù)矩陣A分解為上三角矩陣U和下三角矩陣L，分解的時(shí)間復(fù)雜度為O(n3)。將A分解為L(zhǎng)和U的好處是，無論右側(cè)的b向量是否改變，或者有多個(gè)b向量時(shí)，LU分解只需要進(jìn)行一次，針對(duì)不同的b向量只需要各自進(jìn)行后向回代操作，這大大減少了運(yùn)算的步驟，因?yàn)楹笙蚧卮蠼膺^程的時(shí)間復(fù)雜度為O(n2)，與LU分解的時(shí)間復(fù)雜度相比是微不足道的。第一次后向回代求解到的x0只是一個(gè)初步解，距離精確解還有一定的誤差。因此，接下來就需要進(jìn)行迭代求解來獲得更加準(zhǔn)確的方程解。在反復(fù)迭代求解的過程中，首先計(jì)算殘差r，接下來便等同于求解糾正方程Ad=r，求解時(shí)只需進(jìn)行后向回代，因?yàn)長(zhǎng)和U已經(jīng)得到。將求解得到的d作為修正項(xiàng)，加到已經(jīng)得到解x'上，即新的解x=x'+d。當(dāng)新的方程解的后向誤差小于一定數(shù)量級(jí)時(shí)，便可以認(rèn)為新的方程解已經(jīng)達(dá)到了精度的要求，可以作為最終的結(jié)果，此時(shí)便可以退出循環(huán)。

算法1線性方程組的迭代細(xì)化求解算法

輸入：系數(shù)矩陣A，常數(shù)項(xiàng)向量b

輸出：方程組的解x

在整個(gè)迭代循環(huán)中，只在循環(huán)外進(jìn)行了一次最耗時(shí)的LU分解，在循環(huán)內(nèi)部主要進(jìn)行的是時(shí)間復(fù)雜度更小的后向回代和向量的加減法過程。因此，在迭代次數(shù)遠(yuǎn)小于矩陣尺寸n的情況下，整個(gè)迭代求解過程中的時(shí)間開銷主要集中在剛開始的LU分解過程。因此，如果我們能夠采取某種方法來加速LU分解的過程，那么也會(huì)給整個(gè)迭代求解過程帶來明顯的速度提升。

2 多精度SIMD迭代細(xì)化求解

2.1 多精度迭代細(xì)化求解

文獻(xiàn)[9]的研究表明，可以使用兩種浮點(diǎn)精度或三種浮點(diǎn)精度來加速求解線性方程組。在對(duì)系數(shù)矩陣A進(jìn)行LU分解時(shí)，使用速度更快的單精度浮點(diǎn)乘法，在進(jìn)行代入求解時(shí)，使用精度更高的雙精度浮點(diǎn)乘法。由于LU分解在整個(gè)求解過程中占據(jù)的時(shí)間比重更大，對(duì)其使用單精度進(jìn)行計(jì)算時(shí)，可以使整個(gè)求解過程獲得顯著的加速。因此，可以通過使用高精度和低精度SIMD乘法指令，并通過迭代細(xì)化的方法來加速線性方程組Ax=b的求解。

如算法2所示，在使用多精度計(jì)算線性方程組的算法中，涉及到三種計(jì)算精度：①u是原始數(shù)據(jù)存儲(chǔ)精度的unit-roundoff（單位舍入誤差），②uf是LU分解時(shí)使用的數(shù)據(jù)精度的舍入誤差，③ur是計(jì)算殘差時(shí)的數(shù)據(jù)精度的舍入誤差。為了保證求解過程最終收斂，這三種精度需要滿足以下關(guān)系（精度越高、舍入誤差越小）：ur≤u≤uf。最簡(jiǎn)單的情形便是u=uf=ur，即從始至終使用一種精度，例如只使用雙精度。也可以使用三種不同的精度ur=FP64，u=FP32，uf=FP16。

算法2基于多精度的線性方程組的迭代細(xì)化求解

輸入：系數(shù)矩陣A(FP64)，常數(shù)項(xiàng)向量b(FP64)

輸出：方程組的解x(FP64)

2.2 精度選擇

如表1所示，雙精度（FP64）、單精度（FP32）和半精度（FP16）可以表示的最大值呈指數(shù)性下降[11]。因此將雙精度（FP64）或單精度（FP32）數(shù)據(jù)轉(zhuǎn)換為半精度（FP16）時(shí)，非常容易發(fā)生上溢或下溢。

表1 IEEE754浮點(diǎn)數(shù)特性

文獻(xiàn)[9]的研究表明在使用多精度求解線性方程組時(shí)，系數(shù)矩陣A的條件數(shù)需要滿足一定的條件，才能使整個(gè)迭代過程收斂。在ur=FP64、u=FP32、uf=FP16時(shí)，只有系數(shù)矩陣A的條件數(shù)κ∞(A)<104時(shí)才能保證收斂，條件數(shù)描述了矩陣的病態(tài)情況[10]；在ur=FP64、uf=u=FP32時(shí)，系數(shù)矩陣A的條件數(shù)κ∞(A)<108時(shí)便能保證收斂。

普通矩陣的條件數(shù)有很大的概率超出104，此時(shí)如果選取ur=FP64、u=FP32、uf=FP16，可能會(huì)因?yàn)橄禂?shù)矩陣A的條件數(shù)大于104導(dǎo)致迭代求解過程無法收斂。所以，采用ur=FP64、uf=u=FP32進(jìn)行求解線性方程組時(shí)既能保證迭代求解過程收斂，也能獲得一定的性能收益。

2.3 多精度SIMD加速LU分解

LU分解是高斯消元法的矩陣形式，通過將高斯消元法轉(zhuǎn)化成矩陣運(yùn)算的形式，可以更好地進(jìn)行并行加速。原始的系數(shù)矩陣A可以看成下三角矩陣L和上三角矩陣U的乘積，由于L和U均有一半的元素為零元素，因此可以將下三角矩陣L和上三角矩陣U儲(chǔ)存在一個(gè)n×n的矩陣中，即和原來的系數(shù)矩陣的大小相同，而不占用更多的存儲(chǔ)空間。

如圖1所示，以一個(gè)4×4的LU分解為例，四個(gè)未知數(shù)分別為x,y,z,w。系數(shù)矩陣中的每一列表示某個(gè)未知數(shù)在整個(gè)方程組中的系數(shù)，因此第一列就表示x在4個(gè)方程中的系數(shù)。根據(jù)上一節(jié)介紹的高斯消元法的步驟，首先需要消去第2-4行的x，因此我們需要將a2,1、a3,1、a4,1變?yōu)?。保持第一行的數(shù)不變，將第一行的數(shù)乘以一個(gè)特定的標(biāo)量a并加到第2-4行，為了使a2,1、a3,1、a4,1為0，a應(yīng)該分別等于。隨后保持第二行不變，將第二行的數(shù)分別乘以、，依次類推，直至最后一行只有一個(gè)未知數(shù)，此時(shí)便完成LU分解。此時(shí)原系數(shù)矩陣便成為了最終的上三角矩陣U，其對(duì)角線以下的元素均為0。下三角矩陣L其實(shí)是記錄了每次消元操作的標(biāo)量a，以消去x為例，a2,1、a3,1、a3,1均會(huì)變?yōu)?，在后向代入求解時(shí)，由于我們知道a2,1、a3,1、a3,1已經(jīng)為0，不會(huì)使用a2,1、a3,1、a3,1的值，所以我們可以將每次的標(biāo)量a存儲(chǔ)在即將要被消去的系數(shù)的位置。

圖1 LU分解示意圖

將一行數(shù)乘以一個(gè)系數(shù)加到另一行數(shù)上的操作其實(shí)是axpy操作，由于axpy中的各個(gè)子操作沒有數(shù)據(jù)依賴，因此可以很容易地使用并行加速方法對(duì)其加速。典型的數(shù)據(jù)級(jí)并行加速方法為使用SIMD指令，即單指令多數(shù)據(jù)（Single Instruction Multiple Data），使用一條指令可以同時(shí)對(duì)多個(gè)元素進(jìn)行同樣的操作。

如圖2所示，SIMD指令可以一次對(duì)多個(gè)操作數(shù)進(jìn)行相同的運(yùn)算，而普通的指令只能對(duì)一組數(shù)據(jù)運(yùn)行運(yùn)算。因此，可以使用單精度SIMD乘法指令和單精度SIMD加法指令來執(zhí)行LU分解中的axpy操作。單精度是雙精度浮點(diǎn)數(shù)寬度的一半，因此在不增加向量寄存器的情況下，對(duì)一組64位的寄存器進(jìn)行單精度SIMD乘法或加法運(yùn)算時(shí)，可以同時(shí)計(jì)算兩組單精度乘法或加法，此時(shí)axpy需要的乘法和加法指令的數(shù)目將減少一半。

圖2 SIMD加法指令示意圖

3 實(shí)驗(yàn)仿真與結(jié)果分析

本文的實(shí)驗(yàn)主要分為兩部分，第一部分為使用MATLAB軟件來計(jì)算多精度對(duì)迭代求解次數(shù)的影響。第二部分為使用具有多精度SIMD浮點(diǎn)乘法能力的RISC-V[12]處理器硬件仿真平臺(tái)來運(yùn)行具有多精度SIMD指令的匯編程序來觀察多精度SIMD指令帶來的性能提升。

3.1 多精度迭代次數(shù)

在MATLAB軟件中使用lu()函數(shù)獲得一個(gè)矩陣的下三角矩陣和上三角矩陣，使用左除符號(hào)“”來實(shí)現(xiàn)后向代入求解，L其實(shí)等價(jià)于L-1b，左除符號(hào)首先會(huì)檢查矩陣的形狀，如果矩陣是三角矩陣，那么可以直接進(jìn)行代入求解[13]。首先，只使用雙精度浮點(diǎn)數(shù)來進(jìn)行完整的LU分解操作和后向代入求解操作，記錄當(dāng)殘差足夠小時(shí)的迭代次數(shù)作為基準(zhǔn)數(shù)據(jù)。隨后使用single()函數(shù)將雙精度的A矩陣轉(zhuǎn)換為單精度浮點(diǎn)格式的矩陣，然后使用lu()函數(shù)對(duì)其進(jìn)行LU分解操作，使用雙精度浮點(diǎn)數(shù)來計(jì)算殘差和糾正方程，記錄殘差足夠小時(shí)的迭代次數(shù)。

如表2所示，首先選取了5個(gè)在實(shí)際應(yīng)用中產(chǎn)生的數(shù)據(jù)作為系數(shù)矩陣，假定殘差的無窮范數(shù)小于10-12時(shí)，便可認(rèn)為當(dāng)前循環(huán)中的方程解滿足精度要求?？梢钥吹皆谑褂枚嗑人惴ㄟM(jìn)行求解時(shí)的迭代次數(shù)都要比只使用雙精度求解時(shí)的迭代次數(shù)多，但還是相同的數(shù)量級(jí)，沒有大幅的增加。

表2 多精度求解對(duì)迭代次數(shù)的影響

3.2 多精度SIMD乘法加速結(jié)果

為了驗(yàn)證SIMD乘法指令對(duì)LU分解和后向回代操作的加速效果，首先使用Verilog HDL實(shí)現(xiàn)了可以計(jì)算雙精度或同時(shí)計(jì)算兩個(gè)單精度（SIMD）的浮點(diǎn)乘法器，該乘法器完全符合IEEE754標(biāo)準(zhǔn)，支持向偶舍入、非規(guī)格化數(shù)、NaN等。

為了能夠運(yùn)行完整的匯編程序，將具有多精度計(jì)算能力的浮點(diǎn)乘法器應(yīng)用到了一個(gè)具有五級(jí)流水線單發(fā)射的RISC-V處理器的RTL模型中，該處理器支持?jǐn)?shù)據(jù)轉(zhuǎn)發(fā)、靜態(tài)分支預(yù)測(cè)等技術(shù)來解決數(shù)據(jù)沖突和控制沖突，從而減少流水線阻塞情況的發(fā)生。為了能夠更好地利用數(shù)據(jù)的局部性[14]，該處理器中配置有16KB大小、兩路組相連的L1 Cache，每個(gè)cache block的大小為64B。為了支持SIMD乘法指令，需要在譯碼階段增加對(duì)SIMD乘法指令的譯碼支持，為了能夠?qū)?shù)據(jù)進(jìn)行精度轉(zhuǎn)換，增加了精度轉(zhuǎn)換指令的實(shí)現(xiàn)，分別為高精度轉(zhuǎn)為低精度和低精度轉(zhuǎn)為高精度，并且可以將轉(zhuǎn)換后的數(shù)據(jù)送入目標(biāo)寄存器。

由于處理器中的cache block大小為64B，可以存儲(chǔ)16個(gè)單精度數(shù)據(jù)或8個(gè)雙精度數(shù)據(jù)，為了更好地利用時(shí)間局部性，減少cache block的替換次數(shù)，在進(jìn)行具體的LU分解操作時(shí)，采用了分塊的思想：由于LU分解操作實(shí)際是將某一行所有數(shù)的倍數(shù)分別加到剩下的所有行（axpy），同一行的各個(gè)數(shù)據(jù)之間也沒有數(shù)據(jù)相關(guān)性，以圖3為例，在消去第一列除a1,1外其他的所有數(shù)時(shí)，可以先在陰影區(qū)域內(nèi)執(zhí)行axpy操作，因?yàn)檫@樣可以保持a1,1-a1,16的數(shù)據(jù)始終在cache內(nèi)并命中，如果擁有類似于ARM neon指令集中的向量寄存器，甚至可以確保a1,1-a1,16的數(shù)據(jù)可以始終駐留在寄存器內(nèi)，不發(fā)生寄存器溢出。

圖3 LU分解的分塊算法示意圖

分別編寫了用于LU分解的DGEFA（雙精度）和SGEFA（單精度）匯編程序，以及用于后向回代求解的DGESL匯編程序，在SGEFA中使用了SIMD單精度乘法和SIMD單精度加法指令，一次可以進(jìn)行兩個(gè)乘法或加法操作。將以上的匯編程序匯編為機(jī)器碼后加載到RISC-V處理器的instruction-rom中，使用Synopsys VCS軟件對(duì)RISC-V處理器和多精度浮點(diǎn)乘法器的Verilog文件進(jìn)行硬件仿真實(shí)驗(yàn)。在RISC-V的頂層模塊設(shè)置有一個(gè)計(jì)數(shù)器用來記錄從程序開始到其結(jié)束的時(shí)鐘周期數(shù)，在匯編程序結(jié)尾添加一條匯編語(yǔ)句AD-DI x31,x0,0作為結(jié)束標(biāo)志，當(dāng)檢測(cè)到指令為結(jié)束標(biāo)志時(shí)，退出仿真并打印運(yùn)行周期數(shù)。

如表3所示，通過對(duì)加載了匯編指令的Verilog文件進(jìn)行仿真，得到了DGEFA、SGEFA、DGESL的周期數(shù)。使用雙精度進(jìn)行求解的總周期數(shù)為DGEFA+i×DGESL，使用多精度求解的總周期數(shù)為SGEFA+j×DGESL，其中i,j分別為表2中的雙精度和多精度的迭代次數(shù)。

表3 匯編程序仿真周期結(jié)果

如圖4所示，使用表3中雙精度的周期數(shù)除以多精度的周期數(shù)便可以得到使用多精度SIMD方法對(duì)求解線性方程組帶來的加速比，系數(shù)矩陣是bcsstk04時(shí)的加速比最小，為1.463倍，系數(shù)矩陣是steam1時(shí)的加速比最大，為1.612倍。

圖4 多精度SIMD求解線性方程組加速比

4 結(jié)語(yǔ)

針對(duì)求解較大規(guī)模的線性方程組，本文提出了基于多精度SIMD的加速求解方法，使用低精度的SIMD乘法指令來加速時(shí)間復(fù)雜度最高的LU分解步驟，使用高精度乘法指令計(jì)算殘差和糾正方程。通過相應(yīng)的仿真實(shí)驗(yàn)，探討了多精度計(jì)算對(duì)迭代次數(shù)的影響，以及對(duì)求解過程的加速效果。針對(duì)不同規(guī)模的系數(shù)矩陣，該方法的性能可以達(dá)到傳統(tǒng)方法的1.46-1.61倍，具有較好的加速效果。