面向FT-M7002的一種Cholesky分解向量處理算法

李慧祥,張會(huì)福

(1.湖南科技大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院，湖南 湘潭，411201;2.湖南科技大學(xué) 服務(wù)計(jì)算與軟件服務(wù)新技術(shù)湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 湘潭，411201)

Cholesky分解是正定對(duì)稱矩陣常用的矩陣三角分解方法，相比于通用的LU分解，其計(jì)算量約減少一半。因此，Cholesky分解廣泛應(yīng)用于正定對(duì)稱線性方程組求解[1]、加速矩陣求逆[2]和Kalman濾波并行[3]等問題中。為了提高Cholesky分解算法的執(zhí)行效率，研究人員對(duì)其進(jìn)行了廣泛的研究。吳榮騰[4]在多處理器平臺(tái)上研究了Cholesky分解動(dòng)態(tài)調(diào)度算法，解決各處理器之間的負(fù)載均衡問題。文獻(xiàn)[5-7]在GPU(graphic processing unit)平臺(tái)上采用小尺寸分塊并行的方法加速Cholesky分解，減少CPU與GPU之間的數(shù)據(jù)通信時(shí)間。文獻(xiàn)[8-10]基于FPGA(field programmable gate array)平臺(tái)分析了Cholesky分解的數(shù)據(jù)依賴關(guān)系，設(shè)計(jì)了Cholesky分解的細(xì)粒度流水并行結(jié)構(gòu)。由于處理器的硬件架構(gòu)千差萬別，Cholesky分解已有的優(yōu)化方法并不能充分利用SIMD(single instruction multiple data)DSP的硬件資源，需要研究進(jìn)一步的優(yōu)化方法。

在SIMD體系結(jié)構(gòu)的DSP中，由于多個(gè)并行單元共享一套取指、譯碼、派發(fā)邏輯，控制結(jié)構(gòu)簡(jiǎn)單，可以有效降低硬件代價(jià)，能夠在較低的功耗下獲得較高的性能。近年來，國(guó)內(nèi)的高性能DSP也取得了較大的成就，例如飛騰、龍芯、魂芯、申威等多個(gè)系列處理器。飛騰FT-M7002是國(guó)防科技大學(xué)微電子所自主研發(fā)的高性能SIMD DSP，其片上集成1個(gè)RISC CPU核和2個(gè)FT-MT2 DSP核，在1 GHz運(yùn)行時(shí)，雙精度浮點(diǎn)的峰值性能高達(dá)100 GFlops，單精度浮點(diǎn)的峰值性能高達(dá)200 GFlops[11]。硬件不斷進(jìn)步帶來性能提升的同時(shí)，也對(duì)上層軟件庫(kù)提出了更高的要求。因此，根據(jù)硬件特點(diǎn)優(yōu)化相應(yīng)的算法軟件，充分發(fā)揮出國(guó)產(chǎn)高性能DSP的硬件架構(gòu)優(yōu)勢(shì)具有重要意義。

本文將結(jié)合FT-M7002處理器的硬件資源，研究Cholesky分解算法的優(yōu)化。

1 Cholesky分解算法的基本原理

Cholesky分解的定義是：給定一個(gè)正定對(duì)稱矩陣A，則存在唯一的主對(duì)角線元素全為正數(shù)的下三角矩陣L滿足：

A=LLT

(1)

根據(jù)式(1)，對(duì)于任意一個(gè)對(duì)稱正定矩陣都可以得到式(2)：

(2)

根據(jù)式(2)中矩陣乘法的運(yùn)算規(guī)律，可以得到下三角矩陣L第j列更新的公式，見式(3)和式(4)。

(3)

(4)

2 FT-M7002處理器介紹

FT-M7002是一款40 nm工藝的高性能SIMD DSP，整個(gè)處理器使用三級(jí)存儲(chǔ)結(jié)構(gòu)。單個(gè)DSP核內(nèi)擁有32 kB的標(biāo)量存儲(chǔ)空間(scalar memory，SM)和512 kB的向量存儲(chǔ)空間(vector memory，VM)，所有的計(jì)算核共享2 MB的全局共享Cache，核外擁有32 GB的大容量DDR存儲(chǔ)空間[12-13]。本文的相關(guān)研究基于單個(gè)DSP核。

FT-M7002的DSP核基于VLIW結(jié)構(gòu)，包含一個(gè)五流出的標(biāo)量處理器單元(scalar process unit,SPU)和一個(gè)六流出的向量處理單元(vector process unit,VPU)，兩個(gè)處理單元以緊耦合的方式工作[13]，具體結(jié)構(gòu)見圖1。SPU只包含一個(gè)處理單元，主要負(fù)責(zé)串行任務(wù)處理和程序控制。VPU由16個(gè)向量計(jì)算引擎(vector process engine,VPE)構(gòu)成，最多支持16路32位數(shù)據(jù)進(jìn)行向量運(yùn)算，主要面向密集型的計(jì)算提供并行處理。DMA(direct memory access)為內(nèi)核提供了高速數(shù)據(jù)傳輸通路，實(shí)現(xiàn)核外DDR與SM和VM的快速數(shù)據(jù)交換。VM是VPU獨(dú)享的數(shù)據(jù)存儲(chǔ)器，在訪存不沖突時(shí)，可以同時(shí)支持2個(gè)向量讀/寫，2個(gè)DMA讀/寫共4個(gè)并行請(qǐng)求。通過合理的數(shù)據(jù)排布，可以實(shí)現(xiàn)DMA傳輸和向量訪存并行。

圖1 FT-MT2 內(nèi)核結(jié)構(gòu)Fig.1 FT-MT2 kernel structure

FT-M7002具有大容量的VM、豐富的指令集、大寬度的向量運(yùn)算和高效的數(shù)據(jù)傳輸?shù)忍攸c(diǎn)，在計(jì)算密集型的矩陣算法中能發(fā)揮出較大的優(yōu)勢(shì)。

3 基于FT-M7002的Cholesky分解實(shí)現(xiàn)與優(yōu)化

3.1 Cholesky分解串行算法分析

根據(jù)式(3)、式(4)可以得到Cholesky分解生成下三角矩陣L的串行算法見算法1。分解結(jié)束后矩陣L存儲(chǔ)在矩陣A的下三角位置，下面對(duì)矩陣L的更新將描述成矩陣A相應(yīng)位置的更新。

算法1

Input:矩陣A, 矩陣A的階數(shù)order

Output:矩陣A

Begin

1.forj=0:1:order

2./*對(duì)角線元素A[j,j]的更新*/

3.fork=0:1:j

4.A[j,j]=A[j,j]-A[j,k]×A[j,k]

5.endfor

6.A[j,j]=sqrt(A[j,j])

7. /*第j列其他元素A[i,j](i>j)的更新*/

8.fori=j+1∶1∶order

9.fork=0∶1∶j

10.A[i,j]=A[i,j]-A[i,k]×A[j,k]

11.endfor

12.A[i,j]=A[i,j]/A[j,j]

13.endfor

14.endfor

End

算法2

Input:矩陣A, 矩陣A的階數(shù)order

Output:矩陣A

Begin

1.forj=0∶1∶order

2. /*對(duì)角線元素A[j,j]的更新*/

3.fork=0∶1∶j

4.A[j,j]=A[j,j]-A[k,j]×A[k,j]

5.endfor

6.A[j,j]=sqrt(A[j,j])

7. /*第j行其他元素A[j,i](i>j)的更新*/

8.fori=j+1∶1∶order

9.fork=0∶1∶j

10.A[j,i]=A[j,i]-A[k,i]×A[k,j]

11.endfor

12.A[j,i]=A[j,i]/A[j,j]

13.endfor

14.endfor

End

從算法1可以看出，矩陣A第j列的更新分為兩個(gè)任務(wù)：1)對(duì)角線元素A[j,j]的更新；2)第j列其他元素A[i,j](i>j)的更新。第二個(gè)任務(wù)與第一個(gè)任務(wù)有數(shù)據(jù)依賴關(guān)系，二者只能串行執(zhí)行。

算法1中共有3層嵌套循環(huán)，下文相關(guān)描述中最外層循環(huán)指的是forj=0∶1∶order。第一個(gè)中間層循環(huán)指的是fork=0∶1∶j，完成對(duì)角線元素A[i,j]的更新任務(wù)。第二個(gè)中間層循環(huán)指的是fori=j+1∶1∶order，完成第j列其他元素A[i,j]更新的任務(wù)。兩個(gè)任務(wù)之間有數(shù)據(jù)依賴關(guān)系，二者只能串行執(zhí)行。

第一個(gè)中間層循環(huán)是矩陣A第j行部分元素的平方和，循環(huán)結(jié)束后更新對(duì)角線元素A[j,j]。第二個(gè)中間層循環(huán)是一個(gè)二重循環(huán)，它外層循環(huán)的所有迭代更新矩陣A第j列其他元素A[i,j](i>j)，內(nèi)層循環(huán)是矩陣A第i行和第j行部分元素的點(diǎn)積?？梢钥闯觯喝绻麑?duì)它的內(nèi)層循環(huán)進(jìn)行展開，讓多次循環(huán)迭代遍布在SIMD單元并行執(zhí)行，那么在內(nèi)層循環(huán)結(jié)束后需要對(duì)多個(gè)SIMD單元中的結(jié)果規(guī)約求和，一次完整的內(nèi)層循環(huán)只能得到矩陣A第j列的一個(gè)標(biāo)量結(jié)果，這種計(jì)算方式在SIMD DSP上的執(zhí)行效率較低。如果對(duì)它的外層循環(huán)進(jìn)行展開，讓多次循環(huán)迭代遍布在SIMD單元并行執(zhí)行，這樣就可以將內(nèi)層循環(huán)的點(diǎn)積轉(zhuǎn)換成不同SIMD單元內(nèi)部的乘累加操作，但是這種方式需要對(duì)矩陣A進(jìn)行列訪問，在并行存儲(chǔ)器中這種非連續(xù)訪存效率非常低下，僅為連續(xù)訪存效率的1/15[14]。

綜上所述，算法1直接映射到SIMD DSP上并不能充分利用處理器的硬件資源?？紤]到矩陣對(duì)稱的幾何特征，可以按照同樣的計(jì)算原理生成上三角矩陣來代替下三角矩陣L，在使用時(shí)只需把訪問下三角矩陣的L[i,j]轉(zhuǎn)換成訪問上三角矩陣的LT[j,i]。經(jīng)過上述轉(zhuǎn)換，生成下三角矩陣L計(jì)算過程中對(duì)矩陣A的按列訪問被轉(zhuǎn)換成生成上三角矩陣LT計(jì)算過程中對(duì)矩陣A的按行訪問，SIMD單元間的規(guī)約求和轉(zhuǎn)換成SIMD單元內(nèi)部的乘累加操作。根據(jù)算法1可以得到Cholesky分解生成上三角矩陣LT的串行算法見算法2，在算法2中，兩個(gè)中間層循環(huán)分別完成兩個(gè)任務(wù)：1)對(duì)角線元素A[j,j]的更新；2)第j行其他元素A[j,i](i>j)的更新。

3.2 基于FT-M7002的Cholesky分解優(yōu)化

3.2.1 向量化方法設(shè)計(jì)

在對(duì)角線元素A[j,j]的更新任務(wù)中，為了減少對(duì)A[j,j]的重復(fù)訪問，用一個(gè)局部變量vSum存儲(chǔ)A[k,j]平方和的累加值，在循環(huán)結(jié)束之后再進(jìn)行計(jì)算得到更新后的A[j,j]。需要注意的是，在存儲(chǔ)A[j,j]時(shí)只開啟第1個(gè)VPE，防止覆蓋矩陣的有效數(shù)據(jù)。

在第j行其他元素A[j,i](i>j)的更新任務(wù)中，根據(jù)3.1節(jié)的算法分析，將第二個(gè)中間層循環(huán)的外層循環(huán)展開后，被展開的循環(huán)派發(fā)到所有的VPE上并行執(zhí)行。對(duì)于多次循環(huán)中A[k,j]的重復(fù)訪問問題，通過混洗單元將A[k,j]復(fù)制到所有的VPE中參與運(yùn)算。同時(shí)為了減少內(nèi)層循環(huán)對(duì)A[j,i]的重復(fù)訪問，將0～j-1行的乘累加結(jié)果存儲(chǔ)在局部變量vSum中。經(jīng)過上述優(yōu)化之后的偽代碼見算法3，為了便于理解，算法3中沒有描述出計(jì)算長(zhǎng)度與處理器SIMD寬度不一致的代碼。

在算法3中，mov_to_vlr(N)表示開啟N個(gè)VPE，在FT-M7002上N最大為16。vec_mula()是乘累加函數(shù)，vSum=vec_mula(A[k,i],vtkj,vSum)的含義是vSum=A[k,j]×A[k,j]+vSum。shuffle()是混洗函數(shù)，通過配置特定的混洗模式，將向量變量中第一個(gè)值復(fù)制到所有的VPE中。

3.2.2 對(duì)角線元素更新任務(wù)的循環(huán)合并

在算法3中，第j行其他元素A[j,i](i>j)的更新任務(wù)已經(jīng)充分向量化。但是在對(duì)角線元素A[j,j]的更新任務(wù)中，由于計(jì)算過程中只使用了一個(gè)VPE，浪費(fèi)了大量的計(jì)算能力，需要對(duì)其進(jìn)行優(yōu)化。從外層循環(huán)分析，當(dāng)j=0時(shí)，對(duì)角線A[j,j]的更新任務(wù)中只需要對(duì)對(duì)角線元素A[j,j]開根號(hào)，不需要累加第j列部分元素的平方和。當(dāng)j=n+1時(shí)，對(duì)角線A[j,j]的更新任務(wù)中需要累加的第j列部分元素會(huì)在第j-1行其他元素A[j,i](i>j)的更新任務(wù)中訪問到，見圖2框中數(shù)據(jù)。

定義一個(gè)中間變量vDiag存儲(chǔ)第j-1行其他元素A[j,i](i>j)更新任務(wù)中訪問到的第一個(gè)值，這樣可以直接消除對(duì)角線元素A[j,j]更新任務(wù)的循環(huán)，實(shí)現(xiàn)循環(huán)合并。循環(huán)合并后的偽代碼見算法4。

算法3

Input: 矩陣A, 矩陣A的階數(shù)order

Output:矩陣A

Begin

1.forj=0∶1∶order

2. /*對(duì)角線元素A[j,j]的更新*/

3.vSum=0

4.mov_to_vlr(1)

5.fork=0∶1∶j

6.vSum=vec_mula(A[k,j],A[k,j],vSum)

7.endfor

8.A[j,j]=sqrt(A[j,j]-vSum)

9.mov_to_vlr(N)

10. /*第j行其他元素A[j,i](i>j)的更新*/

11.fori=j+1:N:order

12.vSum=0

13.fork=0∶1∶j

14.vtkj=shuffle(A[k,j])

15.vSum=vec_mula(A[k,i],vtkj,vSum)

16.endfor

17.A[j,i]=(A[j,i]-vSum)/A[j,j]

18.endfor

19.endfor

End

算法4

Input:矩陣A, 矩陣A的階數(shù)order

Output:矩陣A

Begin

1.vDiag=0

2.forj=0∶1∶order

3. /*對(duì)角線元素A[j,j]的更新*/

4.mov_to_vlr(1)

5.A[j,j]=sqrt(A[j,j]-vDiag)

6.mov_to_vlr(N)

7.vDiag=0

8. /*第j行其他元素A[j,i](i>j)的更新*/

9.fori=j+1∶N∶order

10.vSum=0

11.fork=0∶1∶j

12.vtkj=shuffle(A[k,j])

13.vSum=vec_mula(A[k,i],vtkj,vSum)

14.if(i=j+1)

15.vDiag=vec_mula(A[k,i],A[k,i],vDiag)

16.endif

17.endfor

18.A[j,i]=(A[j,i]-vSum)/A[j,j]

19.if(i=j+1)

20.vDiag=vec_mula(A[j,i],A[j,i],vDiag)

21.endif

22.endfor

23.endfor

End

算法4中相關(guān)函數(shù)的含義與算法3一致。在算法4中，通過在第j行其他元素A[j,i](i>j)的更新任務(wù)中以增加兩個(gè)if分支為代價(jià)消除了對(duì)角線元素A[j,j]更新任務(wù)的循環(huán)，在后續(xù)的匯編代碼中這兩個(gè)if分支使用條件執(zhí)行代替，消除跳轉(zhuǎn)指令。

4 結(jié)果分析

為了體現(xiàn)優(yōu)化后的算法在發(fā)揮處理器體系結(jié)構(gòu)性能方面的表現(xiàn)，引入德州儀器(Texas Instruments,TI)的TMS320C6678處理器的庫(kù)函數(shù)DSPF_sp_cholesky的性能進(jìn)行對(duì)比，分析兩類處理器的單核計(jì)算能力。

4.1 平臺(tái)運(yùn)行參數(shù)

算法2、算法3和算法4的性能是通過在FT-M7002開發(fā)板上運(yùn)行手工匯編函數(shù)測(cè)得，TI平臺(tái)的性能是通過在CCS5.5.0模擬器中選用TMS320C6678處理器運(yùn)行庫(kù)函數(shù)DSPF_sp_cholesky測(cè)得。平臺(tái)運(yùn)行參數(shù)見表1。

表1 實(shí)驗(yàn)平臺(tái)參數(shù)Table 1 Experimental platform parameters

4.2 性能分析

實(shí)驗(yàn)測(cè)試了64、128、192、256、320和384階矩陣的Cholesky分解，測(cè)得的節(jié)拍數(shù)見表2，其中算法2、算法3、算法4的匯編代碼沒有使用循環(huán)展開和軟件流水等循環(huán)優(yōu)化方法。

表2 Cholesky分解算法節(jié)拍數(shù)Table 2 The cycle of cholesky decomposition Algorithm

根據(jù)表2數(shù)據(jù)得到算法4相對(duì)于算法3、算法4相對(duì)于TI庫(kù)函數(shù)、算法3相對(duì)于算法2的加速比,具體的折線圖見圖3，其中橫坐標(biāo)表示矩陣的階數(shù)，縱坐標(biāo)表示加速比。

從圖3中算法3相對(duì)于算法2的加速比可以看出，在對(duì)算法2進(jìn)行向量?jī)?yōu)化之后，優(yōu)化后的算法(算法3)相對(duì)于串行算法(算法2)獲得了7.21～12.81的加速比。從趨勢(shì)上看，隨著矩陣規(guī)模的增大，加速比先增加，最后趨于穩(wěn)定。因?yàn)镕T-M7002處理器共有16個(gè)VPE，理論上向量?jī)?yōu)化后的算法相對(duì)于串行算法最高能達(dá)到16倍的加速比。但是，由于本文的實(shí)驗(yàn)是將源數(shù)據(jù)放入核外DDR空間，在計(jì)算過程中由于數(shù)據(jù)傳輸耗時(shí)以及Cholesky分解算法本身會(huì)存在數(shù)據(jù)計(jì)算長(zhǎng)度與SIMD寬度不對(duì)齊的情況，因此加速比達(dá)不到理論值，最后加速比穩(wěn)定在13左右。

算法4在算法3的基礎(chǔ)上消除了對(duì)角線元素A[j,j]更新任務(wù)的循環(huán)，對(duì)算法結(jié)構(gòu)進(jìn)行了優(yōu)化。在圖3的6組數(shù)據(jù)中，算法4相對(duì)于算法3獲得了1.93～2.29的加速比。在算法4中，矩陣規(guī)模越大，對(duì)角線元素A[j,j]更新任務(wù)的循環(huán)耗時(shí)占比越小。因此，隨著矩陣規(guī)模的增大，算法4相對(duì)于算法3的加速比緩慢減小。

圖3中算法4相對(duì)于TI庫(kù)函數(shù)只達(dá)到了1.90～2.82的加速比，并不是Cholesky分解算法在FT-M7002處理器上最佳的性能。循環(huán)展開和軟件流水是常用的循環(huán)優(yōu)化方法，二者搭配使用可以增加循環(huán)內(nèi)指令的并行度，提高循環(huán)的執(zhí)行效率。循環(huán)優(yōu)化后的算法4相對(duì)于TI庫(kù)函數(shù)獲得了3.29～7.01的加速比，見圖4，加速效果較為明顯。

圖3 算法加速比Fig.3 Algorithm speedup

圖4 循環(huán)優(yōu)化之后算法4相對(duì)于TI庫(kù)函數(shù)的加速比Fig.4 Acceleration ratio of algorithm 4 relative to TI library function after cyclic optimization

5 結(jié)語(yǔ)

針對(duì)FT-M7002處理器向量SIMD處理的特點(diǎn)，研究了Cholesky分解算法的優(yōu)化。根據(jù)矩陣對(duì)稱的幾何特征，提出生成上三角矩陣LT代替下三角矩陣L，避免了內(nèi)存的非連續(xù)訪問。分析了最外層循環(huán)相鄰兩次循環(huán)間的數(shù)據(jù)訪問規(guī)律，消除了對(duì)角線元素更新的循環(huán)。結(jié)果表明：最終優(yōu)化后的算法相對(duì)于TI庫(kù)函數(shù)DSPF_sp_cholesky獲得了3.29～7.01的加速比，加速效果較為明顯。研究針對(duì)向量DSP單核的Cholesky分解算法的優(yōu)化，有助于進(jìn)一步研究針對(duì)多核情況相應(yīng)的算法優(yōu)化。