陳鑫峰，王武

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

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

引 言

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

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

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

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

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

1 稀疏矩陣的LDLT 分解算法

1.1 LDLT分解的右視算法

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

其中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的主對角線上。

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

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

稀疏矩陣A的CSC格式由三部分組成，分別為Ai，Ap和Ax，其中Ai表示稀疏矩陣A的行索引，大小為nnz（A的非零元素個(gè)數(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是某個(gè)稀疏矩陣的非零元分布結(jié)構(gòu)，●表示非零元素的位置。

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

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

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

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

1.2 符號分解方法

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

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

（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ù)實(shí)現(xiàn)[14]。

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

圖2 依賴圖Fig.2 The dependency graph

依賴圖為樹結(jié)構(gòu)，對于一個(gè)樹節(jié)點(diǎn)，必須先分解完它的孩子節(jié)點(diǎn)所代表的列，然后使用它的孩子節(jié)點(diǎn)所代表的已分解好的列更新它自己所代表的列，最后才能分解它自己所代表的列?？梢姴煌訕浯淼淖恿蟹纸馔耆?dú)立，這進(jìn)一步提升了并行度，即多層圖劃分并行。

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

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

算法2為確定節(jié)點(diǎn)層劃分的算法偽代碼。其中葉子節(jié)點(diǎn)的層記作level 0，每個(gè)節(jié)點(diǎn)的層編號為它的左孩子節(jié)點(diǎn)的level值與右孩子節(jié)點(diǎn)的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é)點(diǎn)包含的列編號。

對圖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)是選擇一個(gè)level的所有列，這些列可以被并行分解。

在第一層循環(huán)內(nèi)部一共有兩個(gè)階段，第一階段為通過使用向量標(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í)非壓縮數(shù)組temp，即temp[i]=lik，如果temp[i]≠ 0，則表示lik≠ 0，這樣就可以快速定位和更新這些子列。從算法1的第9行也可以看到，該更新涉及到兩次乘法，為了提高數(shù)值分解速度，用了兩個(gè)臨時(shí)數(shù)組：tmpMem1[i]=lik*lkk，tmpMem[i]=lik。而且這些子列只是讀取這些臨時(shí)數(shù)組，并沒有更新它們，所以更新某個(gè)子列并不會影響其它子列的更新。

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

1.4 三角求解方法

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

(1)x=PTb;

(2)Ly=x;

(3)Dx=y;

(4)LTy=x;

(5)x=Py。

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

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

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

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

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

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

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

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

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

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

偽代碼2.動(dòng)態(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ù)中，每個(gè)線程首先確定當(dāng)前被分解的列，然后再分別并行調(diào)用三個(gè)核函數(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的動(dòng)態(tài)并行（dynamic parallel）技術(shù)[20,22]。它支持嵌套執(zhí)行內(nèi)核函數(shù)，即從內(nèi)核中啟動(dòng)內(nèi)核。動(dòng)態(tài)并行的優(yōu)勢是等內(nèi)核執(zhí)行的時(shí)候再配置創(chuàng)建網(wǎng)格和線程塊，這樣就可以動(dòng)態(tài)的利用GPU硬件調(diào)度器和加載平衡器，通過動(dòng)態(tài)調(diào)整適應(yīng)負(fù)載。并且在內(nèi)核中啟動(dòng)內(nèi)核可以減少一部分?jǐn)?shù)據(jù)傳輸?shù)拈_銷。動(dòng)態(tài)并行技術(shù)尤其適合嵌套調(diào)用和遞歸調(diào)用的內(nèi)核函數(shù)[22]，但只支持Nvidia推出的計(jì)算能力大于3.5的GPU設(shè)備。

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

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

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

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

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

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

圖6 不同求解器的整體求解時(shí)間（毫秒）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表示初始矩陣的非零元個(gè)數(shù)（只包含了下三角部分），nnz表示本文求解器符號分解后的L矩陣的非零元個(gè)數(shù)。

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

針對表1中列舉的矩陣，從表2可以看出，三角求解的運(yùn)行時(shí)間遠(yuǎn)少于數(shù)值分解的運(yùn)行時(shí)間，因此表1中只對LDLT分解進(jìn)行GPU加速。從表1和表2還可以看出，非零元填充比越高（即nnz/nz越大），LDLT數(shù)值分解的性能提升越明顯，因?yàn)楦叩姆橇阍畛浔瓤梢愿玫乩糜乙暦纸馑惴ǖ娜龑友h(huán)并行性，而且計(jì)算訪存比更高。

表1 基于本文的直接法求解、LDLT分解和求解階段的時(shí)間（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相比沒有加速，但求解總時(shí)間是最快的，因?yàn)槠浞柗纸獠呗愿?。從圖7和表2可以看出，和UMFPACK相比，LDLT求解器數(shù)值分解階段最高可加速46.2倍，分解和求解總時(shí)間最高加速25.2倍。我們將LDLT與GLU 3.0的LU分解進(jìn)行對比， GPU上的數(shù)值分解部分運(yùn)行時(shí)間見圖8。從圖中可以看出，對于大多數(shù)算例，LDLT比GLU快2-8倍，僅對bcsstk37、bcsstk35兩個(gè)矩陣慢了10%。LDLT比GLU 3.0快的原因是采用了層次劃分和動(dòng)態(tài)并行技術(shù)，且LDLT針對的是對稱矩陣，GLU3.0針對普通稀疏矩陣（對稱矩陣還要補(bǔ)齊另一半非零元數(shù)據(jù)）。但GLU在CPU端的符號分解部分的性能比2.2節(jié)的策略慢10倍以上，因此LDLT整體矩陣分解時(shí)間更快。

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

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

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

利益沖突聲明

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