稀疏對稱矩陣的LDLT分解在GPU上的高效實現(xiàn)

陳鑫峰，王武

1.中國科學(xué)院計算機(jī)網(wǎng)絡(luò)信息中心，北京 100190

2.中國科學(xué)院大學(xué)，北京 100049

引 言

稀疏對稱矩陣的LDLT分解對許多科學(xué)工程問題的有效求解非常重要。例如基于有限差分或有限元離散得到的稀疏線性系統(tǒng)，若采用迭代法難以奏效（迭代收斂很慢甚至不收斂）時，直接法LU分解便可有效解決問題。若矩陣具有對稱性，則LU分解轉(zhuǎn)化為LDLT分解。直接法分解過程計算開銷較大，如何提高其性能尤為重要。直接法的步驟可分為符號分解、數(shù)值分解和三角求解，其中數(shù)值分解最耗時。目前在多核或眾核上并行實現(xiàn)稀疏矩陣分解比較困難，因為分解過程中有較多的數(shù)據(jù)依賴和不規(guī)則數(shù)據(jù)訪問，導(dǎo)致難以充分利用這些多核CPU或加速卡，加速效果不佳，甚至無法加速。

基于CPU+GPU的并行異構(gòu)架構(gòu)已成為當(dāng)今高性能計算機(jī)和服務(wù)器的主流體系結(jié)構(gòu)之一。在GPU上稠密矩陣的計算已經(jīng)具備很高的性能，如CUDA的cuBLAS庫[1]，因為稠密問題有利于高效內(nèi)存訪問和存取，以及向量化計算指令集的應(yīng)用。由于非零元素分布不規(guī)則的稀疏方程組難以高效訪存和向量化計算，無法充分利用GPU的計算資源實現(xiàn)高效求解。目前有很多適用于CPU的稀疏矩陣直接法求解器，比如SuperLU、MKL的PARDISO、MUMPS、UMFPACK、HSL[2]的MA57等，但是在GPU上有效實現(xiàn)稀疏矩陣分解的軟件還很少。CUDA的cuSOLVER[3]還沒有提供GPU上稀疏LU分解的接口。加州大學(xué)Peng等[4]開發(fā)了基于GPU的稀疏LU分解算法庫GLU 3.0，并用于集成電路仿真。它的數(shù)值分解較快，但符號分解存在瓶頸。如何有效利用GPU計算能力提升稀疏矩陣分解的性能仍是需要研究的課題[5]。

本文設(shè)計和實現(xiàn)了一種基于GPU的稀疏對稱矩陣的LDLT分解，該方法基于右視（right- looking）的LDLT分解，使得該分解可以實現(xiàn)右視算法的三層循環(huán)并行，這極大的提高了LDLT分解的并行度，并通過GPU動態(tài)并行技術(shù)提升了GPU的計算和訪存效率。本文最后針對一般稀疏對稱矩陣的一個典型測試集測試了該分解的性能。

實驗結(jié)果表明，本文提出的LDLT分解算法在NVIDIA TITAN V GPU測試平臺上達(dá)到較高的計算性能，針對稀疏對稱矩陣測試集，LDLT數(shù)值分解的計算速度相對于UMFPACK最高加速46.2倍，加速倍數(shù)取決于稀疏矩陣在分解過程中的非零元填充比。

本文結(jié)構(gòu)如下：第1部分介紹LDLT的右視分解算法、符號分解和數(shù)值分解；第2部分重點介紹數(shù)值分解在GPU上的并行算法和實現(xiàn)；第3部分給出分解和求解算法在典型測試集上實現(xiàn)的性能，并與UMFPACK、MA57、cuSOLVER和GLU 3.0進(jìn)行了性能對比。最后給出本文結(jié)論。

1 稀疏矩陣的LDLT 分解算法

1.1 LDLT分解的右視算法

階對稱矩陣A的LDLT分解可表示為A=LDLT，其中L是一個對角線元素為1的下三角矩陣，D是一個對角矩陣，而LT表示L的轉(zhuǎn)置。LDLT分解是高斯消去法的一個變型，即：

其中l(wèi)11=1，d11和a11是標(biāo)量，l21和a21是（n?1）×1的列向量，l22和D22以及A22是（n?1）階的子矩陣。由此可得到，d11=a11，l21=a21/d11，L22D22LT22=A22?l21d11lT22。一直重復(fù)該過程，直到求出L和D的所有元素，D的元素存儲在L的主對角線上。

右視算法每次迭代時先求解D的一個元素，然后求解一列L的元素，最后更新剩下的子矩陣。從中可以看到很大的數(shù)據(jù)依賴性，即必須先分解前k?1列，才能分解第k列，這極大的限制了并行性。

為了節(jié)省存儲空間，需要將稀疏矩陣壓縮后存儲，這導(dǎo)致稀疏矩陣非零元素的內(nèi)存訪問極不規(guī)則。由于處理器對內(nèi)存的訪問有空間局部性和時間局部性的特點，而LDLT分解過程中經(jīng)常不規(guī)則跳躍式訪問數(shù)據(jù)，這與內(nèi)存訪問特點相悖，因此極大地降低了性能，所以選擇一個適當(dāng)?shù)南∈杈仃噧Υ娓袷胶苤匾?，本文采用壓縮稀疏列（Compressed Sparse Colum，CSC）格式，它一定程度上減少了不規(guī)則內(nèi)存訪問。

稀疏矩陣A的CSC格式由三部分組成，分別為Ai，Ap和Ax，其中Ai表示稀疏矩陣A的行索引，大小為nnz（A的非零元素個數(shù)），Ap表示A的列指針，大小為n+1，n為A的維數(shù)，Ax表示A的非零元素值，大小為nnz。第k列元素的索引和值分別為Ai[Ap[k]],Ai[Ap[k]+1],…,Ai[Ap[k+1]-1]和Ax[Ap[k]],Ax[Ap[k]+1], …,Ax[Ap[k+1]- 1]。圖1是某個稀疏矩陣的非零元分布結(jié)構(gòu)，●表示非零元素的位置。

圖1 稀疏矩陣分布Fig.1 The pattern of sparse matrix

算法1展示了右視LDLT分解算法的偽代碼，在這個偽代碼中，L表示初始矩陣A被符號分解后的矩陣，包括重排序和非零元填充（fill-in）。k表示當(dāng)前列的索引，j表示當(dāng)前列右邊的列。第k列被計算完成后，第j（j>k）列將會被立即更新。即第k列被計算完以后，該算法會向右看，使用算好的第k列更新第k列所有右邊的列j（j>k）。所以主要的操作變成子矩陣更新，該過程具有較好的并行性。

算法1.基于右視的LDLT分解

分解過程有3層循環(huán)，最外層的循環(huán)表示正在被分解的列，中間一層循環(huán)表示將要被更新的列，而最內(nèi)層的循環(huán)表示待更新列的某個元素。顯而易見，最內(nèi)層的兩個循環(huán)完全可以并行計算，后面介紹如何通過圖劃分實現(xiàn)最外層循環(huán)并行。

1.2 符號分解方法

由于本文待求解的稀疏矩陣是對稱的，并且符號分解不涉及任何具體的數(shù)值，所以本文使用稀疏Cholesky分解的符號分解[6-11]，這一部分時間很短，所以沒有使用GPU加速。這一步結(jié)束以后會將矩陣重新排列，即C=PTAP，其中P是排列矩陣，而且還會生成一個parent數(shù)組，表示列與列的依賴關(guān)系，即parent[i]=k表示第k列依賴第i列，也就表示必須先分解第i列，然后使用分解好的第i列更新第k列，最后才能分解第k列。這個parent數(shù)組也為后續(xù)的最外層循環(huán)并行提供基礎(chǔ)。

符號分解也需要為下一步的數(shù)值分解確定L矩陣的結(jié)構(gòu)，因為對稱稀疏矩陣分解后，原先的零元素可能變成非零元素，即fill-in，所以需要盡可能的減少fill-in的數(shù)量。這就需要重排序算法，即減少fill-in的數(shù)量而提出的算法。符號分解主要包含以下4個步驟：

（1）重排序

該步驟通過對稀疏矩陣的非零元位置置換，即重排序來盡可能減少fill-in的數(shù)量。結(jié)束后會返回置換矩陣P，并使用該矩陣來重排序初始矩陣A，即C=PTAP。常用的重排序有近似最小度[12]（Approximate Minimum Degree, AMD）排序和METIS軟件包的圖劃分排序[13]，我們默認(rèn)調(diào)用METIS重排序，對于規(guī)模較小的矩陣（比如n小于5000），則使用AMD重排序。METIS的排序算法基于多層次遞歸二分切分法、多層K路劃分法以及多約束劃分機(jī)制，執(zhí)行效率更高，填充數(shù)量更少。

（2）尋找依賴圖

該步驟根據(jù)重排序后的矩陣C，尋找列與列之間的依賴關(guān)系，即依賴圖parent，數(shù)組parent[i]=k表示第k列依賴第i列，必須先分解完第i列，然后使用第i列更新第k列，最后才能分解第k列?？赏ㄟ^構(gòu)造etree函數(shù)實現(xiàn)[14]。

圖2為圖1稀疏矩陣符號分解過程的依賴圖：其中節(jié)點表示列，邊表示列與列之間的依賴關(guān)系，且位于上面的節(jié)點所代表的列依賴下面節(jié)點所代表的列。從該依賴圖看出，分解完第1列后才能分解第6列；分解完第6列后才能分解第7列；分解完第7列后才能分解第9列；分解完第9列后才能分解第10列；分解完第10列后才能分解第11列。

圖2 依賴圖Fig.2 The dependency graph

依賴圖為樹結(jié)構(gòu)，對于一個樹節(jié)點，必須先分解完它的孩子節(jié)點所代表的列，然后使用它的孩子節(jié)點所代表的已分解好的列更新它自己所代表的列，最后才能分解它自己所代表的列。可見不同子樹代表的子列分解完全獨立，這進(jìn)一步提升了并行度，即多層圖劃分并行。

（3）根據(jù)依賴圖確定levels

根據(jù)以上分析，我們可以實現(xiàn)層劃分并行。該步驟根據(jù)列與列的依賴關(guān)系，即依賴圖數(shù)組parent確定可并行計算的列，即level[]，同一個level中的列都可以并行分解，不同的level中的列有依賴關(guān)系，只能串行分解。

算法2為確定節(jié)點層劃分的算法偽代碼。其中葉子節(jié)點的層記作level 0，每個節(jié)點的層編號為它的左孩子節(jié)點的level值與右孩子節(jié)點的level值最大值取其較大者加1。level_i[]數(shù)組大小為n，它表示每一層包含的列編號，總層數(shù)為nlevel，數(shù)組level_p[]用于確定每層level_i的起始和終止列，大小為nlevel+1。例如第k層包含的列為level_i[level_p[k]],level_i[level_p[k]+1],…,level_i[level_p[k+1]-1]。這跟儲存稀疏矩陣的CSC格式類似，因此可通過level_i[]和level_p[]索引快速定位每層節(jié)點包含的列編號。

對圖2進(jìn)行多層劃分后的依賴圖見圖3。其中屬于level 0的為第1、2、4、5列，這些列可以被并行分解。屬于level 1的為第3，6列，這些列必須等待level 0中的列被分解完以后才能被分解，同樣這些列也可以被并行分解。屬于level 2的為第7、8列，這些列必須等待level 1中的列被分解完以后才能被分解，同樣這些列可以被并行分解。屬于level 3的為第9列，這些列必須等待level 2中的列被分解完以后才能被分解，同樣這些列可以被并行分解。屬于level 4的為第10列，這些列必須等待level 3中的列被分解完以后才能被分解，同樣這些列可以被并行分解。屬于level 5的為第11列，這些列必須等level 4中的列被分解完以后才能被分解，同樣這些列可以被并行分解。

圖3 多層劃分后的依賴圖Fig.3 The dependency graph after level partition

（4）根據(jù)依賴圖數(shù)組確定L的結(jié)構(gòu)

該步驟根據(jù)列與列之間的依賴關(guān)系，即依賴圖parent數(shù)組，來為下一步的數(shù)值分解確定L的結(jié)構(gòu)。確定L結(jié)構(gòu)的算法基于Cholesky分解的上視分解[6]，算法3為確定填充元位置的偽代碼，其中矩陣L第k行的非零元位置通過ereach函數(shù)[6]得到。圖4為圖1的稀疏矩陣符號分解填充后的L矩陣結(jié)構(gòu)，其中●表示非零元素的位置，○表示填充元的位置。

圖4 填充后的L矩陣Fig.4 The L matrix after fill-in

算法3.確定填充元的位置

17: Lx[p]= x[i];18: x[i]= 0; //為第k+1行清除x 19: end for 20: //存儲Lkk 21: p = c[k]++;22: Li[p]= k;23: Lx[p]= x[k];24: x[k]= 0;25: end for 26: Lp[n]= c[n-1];

1.3 并行數(shù)值分解方法

基于右視的LDLT分解算法見算法1，對應(yīng)的并行數(shù)值分解過程的偽代碼見算法4[14]。

算法4.并行化的LDLT數(shù)值分解

其中第一層循環(huán)是選擇一個level的所有列，這些列可以被并行分解。

在第一層循環(huán)內(nèi)部一共有兩個階段，第一階段為通過使用向量標(biāo)量除法來并行分解當(dāng)前選中的列，可以看出這些被選中的列可以被完全并行分解。第二階段為更新操作，即使用當(dāng)前分解好的列來更新那些依賴該列的子列，這些子列很容易由以下方法定位：

假設(shè)當(dāng)前被分解好的列為k，且lik≠ 0，lik≠ 0(i≥ j)，根據(jù)算法1中的第9行l(wèi)ij= lij- lik* ljk* lkk就可以知道lij≠ 0。所以如果lik≠ 0，則第i 列依賴第k列。為了快速判斷l(xiāng)ik是否等于零，采用臨時非壓縮數(shù)組temp，即temp[i]=lik，如果temp[i]≠ 0，則表示lik≠ 0，這樣就可以快速定位和更新這些子列。從算法1的第9行也可以看到，該更新涉及到兩次乘法，為了提高數(shù)值分解速度，用了兩個臨時數(shù)組：tmpMem1[i]=lik*lkk，tmpMem[i]=lik。而且這些子列只是讀取這些臨時數(shù)組，并沒有更新它們，所以更新某個子列并不會影響其它子列的更新。

從算法1中可以看到，這些子列只更新一次且更新相互獨立，所以這些子列的更新完全可以被并行。第三層循環(huán)本質(zhì)上是一個向量標(biāo)量乘法和向量加法操作，即MAD操作，這個操作用來更新一個子列，算法1中子列每個元素的更新相互獨立，完全可以并行進(jìn)行?？梢钥吹?，算法1的三層循環(huán)通過算法2和算法3實現(xiàn)完全并行，這有利于在GPU上提升性能。

1.4 三角求解方法

完成重排序、符號分解和LDLT數(shù)值分解之后，就可以求解線性方程組Ax = b，其中x可以是向量或矩陣，即同時求解多個有相同系數(shù)矩陣和不同右端向量b的線性方程組。因為C=PTAP，且C=LDLT，基于分解矩陣L、D和置換矩陣P，Ax = b的求解過程為：

(1)x=PTb;

(2)Ly=x;

(3)Dx=y;

(4)LTy=x;

(5)x=Py。

為了提高數(shù)值分解的并行度，沒有采用動態(tài)選主元，而是假設(shè)主對角元不為零，在分解過程中，如果主對角元足夠小，就將其賦值為常量，這會導(dǎo)致主對角元非常小的情況下數(shù)值分解后三角求解的計算結(jié)果有一定的誤差，為減少誤差，三角求解時使用迭代收斂[15-16]方法提高精度。相對于數(shù)值分解，三角求解的計算時間很短，沒有必要使用GPU加速。

2 LDLT數(shù)值分解在GPU上的實現(xiàn)

2.1 GPU架構(gòu)和CUDA編程

GPU在科學(xué)與工程計算的很多應(yīng)用領(lǐng)域發(fā)揮著越來越重要的作用[17-21]，尤其是對于計算密集型的應(yīng)用。GPU由一組線程化的流多處理器（SM）構(gòu)成，它們共享全局內(nèi)存。每個SM含多個標(biāo)量處理器、L1緩存、常量緩存（C-Cache）、多線程指令獲取(MT issue)以及共享內(nèi)存等部件。線程是GPU的基本運算單元，多個線程構(gòu)成線程塊，一個SM上可同時運行多個塊，塊內(nèi)以單指令多線程（SIMT）方式執(zhí)行程序。塊內(nèi)的線程可實現(xiàn)同步。

CUDA作為GPU的編程模型，它支持核調(diào)用和核內(nèi)存分配，從而充分利用GPU的并行計算能力。CUDA假設(shè)CPU和GPU有各自的內(nèi)存空間，分別稱為主機(jī)內(nèi)存和設(shè)備內(nèi)存。同一個塊中所有線程可訪問同一共享內(nèi)存，全局內(nèi)存所有線程均可訪問。CUDA提供了blockIdx.x和threadIdx.x等內(nèi)置變量，調(diào)用核函數(shù)時這些變量被自動確定，指引每個線程和塊的索引，將數(shù)據(jù)映射到對應(yīng)的線程上，從而每個線程可獨立處理對應(yīng)的數(shù)據(jù)。

2.2 LDLT數(shù)值分解在GPU上的實現(xiàn)

LDLT數(shù)值分解的算法偽代碼見偽代碼1，被重排序后的矩陣C=PTAP的列被劃分成了level層，即算法的最外層循環(huán)，每層的列可并行分解。前文已經(jīng)表明，需要兩個臨時非壓縮的內(nèi)存空間來更新子列，但是GPU的全局內(nèi)存有限，所以需要提前判斷當(dāng)前GPU還有多少全局內(nèi)存可用。

在偽代碼1的第3行，free表示剩下可用的全局內(nèi)存，單位是字節(jié)。因為level中的每一列分解都需要nB=2*n*sizeof(數(shù)據(jù)類型)字節(jié)的內(nèi)存（算例中的矩陣數(shù)據(jù)是單精度的，可以是實數(shù)也可以是復(fù)數(shù)），而如果level中的列比free/nB還要大,尤其是n很大且level很小時，無法一次性全部并行分解完level中的所有列，必須分幾次才能完這些列的并行分解。

偽代碼1.基于GPU的并行LDLT數(shù)值分解

一旦確定了要并行分解的列，就要調(diào)用一個核函數(shù)（dynamic）來并行分解這些列，即偽代碼1的第14行，其中restCol表示當(dāng)前需要并行分解的列數(shù)，該核函數(shù)調(diào)用表示每個線程塊負(fù)責(zé)分解一列，其偽代碼如下：

偽代碼2.動態(tài)并行分解核函數(shù)

10: update<<>>(Lp,Li,Lx,tmpMem,tmpMem1,n,k,blockIdx.x);11: cleartmpMem<<<(subColSize+1023)/1024,1024>>>(Lp,Li,tmpMem,tmpMem1,n,k,blockIdx.x);12: end function

在dynamic核函數(shù)中，每個線程首先確定當(dāng)前被分解的列，然后再分別并行調(diào)用三個核函數(shù)factorize、update和clearTmpMem，它們分別表示分解當(dāng)前列，用分解好的當(dāng)前列更新依賴這一列的子列，以及清空分解當(dāng)前列所使用的非壓縮的內(nèi)存空間，它們的偽代碼分別為：

偽代碼3.當(dāng)前列的多線程分解

偽代碼4.更新當(dāng)前列的子列

14: atomicAdd(＆(Lx[Lp[j]+offset+threadIdx.x]), -x);15: end if 16: offset += blockDim.x;17: end while 18: end function

偽代碼5.清空當(dāng)前列多線程分解所用內(nèi)存

本文使用了CUDA C的動態(tài)并行（dynamic parallel）技術(shù)[20,22]。它支持嵌套執(zhí)行內(nèi)核函數(shù)，即從內(nèi)核中啟動內(nèi)核。動態(tài)并行的優(yōu)勢是等內(nèi)核執(zhí)行的時候再配置創(chuàng)建網(wǎng)格和線程塊，這樣就可以動態(tài)的利用GPU硬件調(diào)度器和加載平衡器，通過動態(tài)調(diào)整適應(yīng)負(fù)載。并且在內(nèi)核中啟動內(nèi)核可以減少一部分?jǐn)?shù)據(jù)傳輸?shù)拈_銷。動態(tài)并行技術(shù)尤其適合嵌套調(diào)用和遞歸調(diào)用的內(nèi)核函數(shù)[22]，但只支持Nvidia推出的計算能力大于3.5的GPU設(shè)備。

3 數(shù)值結(jié)果

算法實現(xiàn)基于C/C++和CUDA C，編譯選項為-O3，其中GCC版本為4.8.5，CUDA的版本為11.1。測試平臺包含2個20核的Intel Xeon Gold 6148 CPU和4個NVIDIA TITAN V GPU，CPU主頻為2.4GHz，CPU端的內(nèi)存為384GB，每張GPU卡的顯存為12GB。操作系統(tǒng)為Red Hat Linux 4.8.5-16。

實驗使用的測試集來自佛羅里達(dá)大學(xué)稀疏矩陣集[23]，它被廣泛用于稀疏矩陣基準(zhǔn)測試，評價稀疏矩陣相關(guān)算法的性能。用于對比本文求解器性能的UMFPACK-5.6.2是稀疏數(shù)學(xué)庫SuiteSparse[24]的一個子集，它同時支持復(fù)數(shù)和實數(shù)矩陣的數(shù)值分解和直接法求解。

用于性能對比的稀疏矩陣直接法求解器還有HSL的MA57、CUDA自帶的cuSOLVER和加州大學(xué)Peng等開發(fā)的GLU。cuSPARSE只提供不完全LU分解，cuSOLVER的cusolverSp模塊雖然實現(xiàn)了稀疏矩陣的LU分解、Cholesky分解和QR分解，但目前的版本還沒有提供設(shè)備端的接口，因此無法通過GPU加速。GLU3.0在GPU上實現(xiàn)了稀疏矩陣的LU分解，并用于電路仿真，但只支持實數(shù)類型的稀疏矩陣。

利用CUDA實現(xiàn)的稀疏對稱矩陣的LDLT分解求解器在測試平臺上分解和求解各階段的運行時間見表1（其中數(shù)值分解在GPU上進(jìn)行）；在CPU上采用UMFPACK、MA57和GPU上的LDLT數(shù)值分解的時間見圖5；采用UMFPACK、MA57、cuSOLVER和GPU上的LDLT求解器求解各階段的整體時間見圖6；本文求解器相比UMFPACK的LDLT數(shù)值分解的加速比和整個方程組求解時間的加速比見表2和圖7。圖表中的時間單位為毫秒（下同）。

圖5 不同軟件數(shù)值分解的時間（毫秒）Fig.5 Runtime of LDLT (in ms)

圖6 不同求解器的整體求解時間（毫秒）Fig.6 Total runtime of different solvers (in ms)

圖7 數(shù)值分解和整體求解的加速比Fig.7 Speedup of LDLT and total solver

表中的稀疏對稱矩陣以維數(shù)n排序的，考慮到并行度和內(nèi)存限制，排除了1萬階以下和25萬階以上的矩陣集。其中matrix表示矩陣名稱，n表示矩陣的維數(shù)，nz表示初始矩陣的非零元個數(shù)（只包含了下三角部分），nnz表示本文求解器符號分解后的L矩陣的非零元個數(shù)。

nnz/nz表示fill-in的倍數(shù)。symbolic表示符號分解運行時間，numeric表示數(shù)值分解運行時間，solve表示三角求解運行時間，total表示本文求解器總的運行時間，為符號分解和數(shù)值分解以及三角求解運行時間之和，T(umf)表示UMFPACK總的運行時間，即符號分解和數(shù)值分解以及三角求解運行時間之和，Sp(ldlt)表示數(shù)值分解的加速比，Sp(total)表示總運行時間的加速比。

針對表1中列舉的矩陣，從表2可以看出，三角求解的運行時間遠(yuǎn)少于數(shù)值分解的運行時間，因此表1中只對LDLT分解進(jìn)行GPU加速。從表1和表2還可以看出，非零元填充比越高（即nnz/nz越大），LDLT數(shù)值分解的性能提升越明顯，因為更高的非零元填充比可以更好地利用右視分解算法的三層循環(huán)并行性，而且計算訪存比更高。

表1 基于本文的直接法求解、LDLT分解和求解階段的時間（ms）Table 1 Runtime (in ms)of LDLT decomposition and solving phases using our direct solver

從圖5、圖6可以看出，GPU上的LDLT求解器性能優(yōu)于UMFPACK、MA57和cuSolver求解器，雖然對少數(shù)矩陣（如bcsstk37、bcsstk35），LDLT數(shù)值分解部分與MA57相比沒有加速，但求解總時間是最快的，因為其符號分解策略更佳。從圖7和表2可以看出，和UMFPACK相比，LDLT求解器數(shù)值分解階段最高可加速46.2倍，分解和求解總時間最高加速25.2倍。我們將LDLT與GLU 3.0的LU分解進(jìn)行對比， GPU上的數(shù)值分解部分運行時間見圖8。從圖中可以看出，對于大多數(shù)算例，LDLT比GLU快2-8倍，僅對bcsstk37、bcsstk35兩個矩陣慢了10%。LDLT比GLU 3.0快的原因是采用了層次劃分和動態(tài)并行技術(shù)，且LDLT針對的是對稱矩陣，GLU3.0針對普通稀疏矩陣（對稱矩陣還要補(bǔ)齊另一半非零元數(shù)據(jù)）。但GLU在CPU端的符號分解部分的性能比2.2節(jié)的策略慢10倍以上，因此LDLT整體矩陣分解時間更快。

圖8 LDLT和GLU3.0在GPU上運行時間（毫秒）Fig.8 Runtime of LDLT and GLU 3.0 on GPU (in ms)

表2 UMFPACK的分解和求解各階段的時間（ms），以及本文求解器相對UMFPACK的加速比Table 2 Runtime (in ms)of LDLT decomposition and solving phases using UMFPACK solver, and the speedup of our solver compared with UMFPACK

4 結(jié)論與展望

本文設(shè)計并實現(xiàn)了一種基于GPU的稀疏對稱矩陣的LDLT求解器，該求解器采用Cholesky符號分解和右視算法的LDLT數(shù)值分解，并且該數(shù)值分解可以對右視分解算法的三層循環(huán)進(jìn)行完全線程級并行計算，并采用依賴圖的層次劃分策略和CUDA的動態(tài)并行技術(shù)，這非常有利于LDLT數(shù)值分解在GPU上的加速。實驗結(jié)果表明，針對稀疏對稱矩陣測試集，GPU上實現(xiàn)的LDLT數(shù)值分解計算速度相對于CPU上的UMFPACK最高加速46.2倍，直接求解整體最高加速25.2倍。而且矩陣分解速度也快于加州大學(xué)Peng等開發(fā)的基于GPU的稀疏矩陣LU分解數(shù)值軟件包GLU3.0。

LDLT分解在GPU上的實現(xiàn)策略可為其它高性能GPU異構(gòu)平臺上開展稀疏方程組的高性能數(shù)值求解器的設(shè)計實現(xiàn)提供借鑒，并有望應(yīng)用于復(fù)雜電路仿真與設(shè)計、結(jié)構(gòu)力學(xué)分析等領(lǐng)域的有限元離散線性方程組在GPU上快速高效求解。

利益沖突聲明

所有作者聲明不存在利益沖突關(guān)系。