使用GPU加速計算矩陣的Cholesky分解

沈　聰　高火濤

(武漢大學(xué)電子信息學(xué)院　湖北 武漢 430072)

?

使用GPU加速計算矩陣的Cholesky分解

沈聰高火濤

(武漢大學(xué)電子信息學(xué)院湖北 武漢 430072)

針對大型實對稱正定矩陣的Cholesky分解問題，給出其在圖形處理器(GPU)上的具體實現(xiàn)。詳細(xì)分析了Volkov計算Cholesky分解的混合并行算法，并在此基礎(chǔ)上依據(jù)自身計算機(jī)的CPU以及GPU的計算性能，給出一種更為合理的三階段混合調(diào)度方案，進(jìn)一步減少CPU的空閑時間以及避免GPU空閑情況的出現(xiàn)。數(shù)值實驗表明，當(dāng)矩陣階數(shù)超過7000時，新的混合調(diào)度算法相比標(biāo)準(zhǔn)的MKL算法獲得了超過5倍的加速比，同時對比原Volkov混合算法獲得了顯著的性能提升。

圖形處理器喬里斯基分解加速比混合算法

0　引　言

近年來，隨著計算機(jī)技術(shù)的發(fā)展，圖形處理器(GPU)越來越強(qiáng)大，并且在計算能力上已經(jīng)超過了通用CPU。使用GPU計算可以以低廉的價格獲得巨大的計算性能，因此成為了科學(xué)計算領(lǐng)域的一個應(yīng)用熱點。自2007年NVIDIA公司推出了CUDA運算平臺，并使用C語言為CUDA構(gòu)架編寫程序以來，GPU計算技術(shù)已經(jīng)廣泛應(yīng)用于諸如信號處理、圖像處理、信息安全等熱門領(lǐng)域中。在數(shù)值線性代數(shù)中，GPU可以用于加速大規(guī)模的矩陣計算問題，包括矩陣的分解和求特征值以及特征向量等。許多組織和機(jī)構(gòu)在GPU上實現(xiàn)了LAPACK庫中的函數(shù)，并發(fā)布了相關(guān)軟件包以供科研人員使用。比較著名的軟件包有CULA和MAGMA等。

Cholesky分解是矩陣計算中的一個基本分解問題，常常用于求解系數(shù)矩陣為對稱正定矩陣的線性方程組，由于其計算量比使用LU分解求解一般線性方程組的算法約減少一半。因此在科學(xué)計算中也有廣泛的應(yīng)用。此外，Cholesky分解也應(yīng)用于Kalman 濾波[1]以及Monte Carlo仿真[2]等算法中。隨著基本矩陣算法在科學(xué)計算中廣泛使用以及其處理的數(shù)據(jù)矩陣越來越大，研究使用GPU加速矩陣計算問題是十分必要的。

CUBLAS是在GPU上實現(xiàn)的BLAS。利用CUBLAS庫可以很方便地移植CPU代碼到GPU上進(jìn)行加速計算。自2008年Volkov等人提出CPU-GPU混合計算矩陣的LU、Cholesky和QR分解算法[5]以來，混合算法思想被廣泛應(yīng)用于各個計算領(lǐng)域，如流體仿真計算[6]、光線跟蹤算法[7]等?；旌纤惴ㄋ枷胍惨驗閂olkov等人的工作正逐步應(yīng)用于各基本矩陣算法中，如矩陣的Hessenberg約化、二對角化以及矩陣的特征值求解中。使用CPU-GPU混合算法可以通過減少CPU的空閑時間，進(jìn)一步加速矩陣計算。然而混合算法的性能很大程度上取決于調(diào)度算法的優(yōu)劣。好的調(diào)度方案可以最小化CPU與GPU的空閑時間，達(dá)到混合最優(yōu)的目的；不適當(dāng)?shù)恼{(diào)度策略可能導(dǎo)致計算錯誤或者計算時間的延長。基于此，本文將在最新的Kerpler GPU構(gòu)架上重新考察Volkov的混合算法，并給出一種新的調(diào)度方案，以達(dá)到混合更優(yōu)的目的。數(shù)值實驗表明基于新的三階段混合調(diào)度算法的函數(shù)New_Sche相比MKL中的dpotrf函數(shù)最高獲得了5.2倍的加速比，而且其計算性能顯著優(yōu)于原Volkovde 混合算法。文中給出的針對Cholesky分解的調(diào)度策略同樣適用于矩陣的其他分解及約化算法。

1　使用CPU計算Cholesky分解

給定一個對稱正定矩陣A，則存在主對角元素全為正數(shù)的下三角矩陣L滿足：

A=LLT

(1)

通常稱下三角矩陣L為矩陣A的Cholesky因子。下面給出Cholesky分解的直接算法[3]及基于塊的下遞推更新算法[4]。

1.1直接Cholesky分解

對于對稱正定矩陣A直接由式(1)得到如下等式：

(2)

則對于第j個對角元素ajj，有：

(3)

對于位置為(i,j)(i>j)的元素，我們可得如下等式：

aij=ri1rj1+ri2rj2+…+rijrjj

(4)

由此可得到求解Cholesky分解的公式如下：

(5)

(6)

若令j=1,2,…,n，循環(huán)使用式(5)、式(6)計算L的下三角元素，這樣便得到計算Cholesky分解的直接算法。由于該直接算法主要是Level 2 BLAS操作，不適合移植到GPU上計算。

1.2基于塊的Cholesky分解

BLAS庫內(nèi)部實現(xiàn)了各種線性代數(shù)的基本運算。由于其代碼的高效性和完整性，因此一直為科研人員廣泛使用。其中Level3BLAS是依據(jù)現(xiàn)代計算機(jī)的內(nèi)存等級進(jìn)行了大量的優(yōu)化處理，使得其計算性能遠(yuǎn)高于Level2BLAS。類似的結(jié)論也在GPU上成立。

為了能利用矩陣的子塊進(jìn)行計算，下面使用Cholesky分解的塊遞推更新算法。對矩陣A進(jìn)行分塊，如下：

A=A11BTBA＾[]　A11∈

則由式(1)有：

A=A11BTBA＾[]=L110SL＾[]LT11ST0L＾T[]

于是得到如下遞推公式：

L11=cholesky(A11)

(7)

(8)

(9)

對矩陣按式(7)-式(9)依次遞推更新下去，便可得到Cholesky分解。

具體的遞推過程如圖1所示。

圖1　基于塊的Cholesky分解

1.3算法及分析

基于塊的Cholesky分解算法如下：

Step1使用IntelMKL庫中的dpotrf函數(shù)計算第一個子塊A1的Cholesky分解；

Step2使用BLAS庫中的Level3函數(shù)trsm求解矩陣方程X×A1T=A2；

Step3使用BLAS庫中的Level3函數(shù)syrk更新尾部矩陣A3=A3-A2×A2T；

Step4對尾部矩陣A3重新分塊并針對A3重復(fù)執(zhí)行以上三步，直到計算完最后一個子塊的Cholesky分解則結(jié)束。

在余部矩陣A3的階數(shù)較大時，Step3占據(jù)一次循環(huán)約80%的計算時間。而A1由于分塊大小的固定，計算量不會改變，Step1所占的計算時間非常小。當(dāng)A2的行遠(yuǎn)大于分塊值時，Step2幾乎占據(jù)了剩余20%的計算時間。

2　使用GPU計算Cholesky分解

2.1CUBLAS庫

CUBLAS是由NVIDIA公司提供的并在CUDA構(gòu)架的GPU上實現(xiàn)BLAS的函數(shù)庫。利用CUBLAS庫可以很方便地移植矩陣計算的CPU代碼到GPU上進(jìn)行計算。由于矩陣算法中主要涉及矩陣與矩陣相乘等的Level 3操作以及矩陣和向量相乘等的Level 2操作，且BLAS和CUBLAS庫中的函數(shù)是經(jīng)過高度優(yōu)化的，并被全世界廣泛認(rèn)可，因此在計算Cholesky分解中，可直接使用該庫中的函數(shù)在CPU或GPU上實現(xiàn)矩陣的Level 2或Level 3操作。這樣，可以將主要精力放在Cholesky分解算法的任務(wù)劃分及任務(wù)調(diào)度上。

2.2使用CUBLAS計算Cholesky分解

圖2　使用CUBLAS加速Cholesky分解

由于CUBLAS庫中提供了trsm和syrk函數(shù)，但沒有提供potrf函數(shù)。對此，可以在每次循環(huán)時先將需要計算Cholesky分解的子矩陣塊傳回到CPU上并利用MKL中的potrf函數(shù)計算其Cholesky分解，然后在將結(jié)果傳輸至GPU對應(yīng)的位置。再依次使用CUBLAS庫函數(shù)在GPU上更新矩陣的剩余部分。然而這樣做的結(jié)果是當(dāng)GPU(或CPU)執(zhí)行計算任務(wù)時，CPU(或GPU)處于閑置狀態(tài)，這導(dǎo)致了硬件資源的浪費。使用CUBLAS計算Cholesky分解的具體過程如圖2所示。

這樣直接使用CUBLAS庫函數(shù)加速實對稱正定矩陣的Cholesky分解的效果也是十分明顯的。首先由于在GPU上的操作全為Level 3的BLAS操作；其次，相對較小的分塊值nb，如nb<512，在本實驗環(huán)境下，使用MKL庫函數(shù)計算512階對稱正定矩陣的Cholesky分解耗時約為4 ms，而來回傳輸512階矩陣總共耗時約2 ms。相比較大的對稱正定矩陣而言，直接使用CUBLAS庫函數(shù)已經(jīng)可以稱之為充分利用GPU進(jìn)行計算了。

2.3Volkov的混合算法及分析

Volkov等人于2008年給出了CPU-GPU混合計算Cholesky分解的算法[5]，并在NVIDIA的GTX 280上獲得了良好的加速效果。該算法隨后被MAGMA軟件包采納，成為計算Cholesky分解的標(biāo)準(zhǔn)算法之一。圖3-圖6以4×4的矩陣塊為例，簡單描述Volkov的混合算法。

圖3　第1步更新

圖4　第2步更新

圖5　第3步更新

圖6　第4步更新

圖7　某次剩余矩陣

Volkov的混合算法巧妙地移動了矩陣循環(huán)的更新順序，這樣余部矩陣的上次更新和本次子塊的Cholesky分解及其所在列的更新可以同時進(jìn)行。

考察Volkov算法中的混合并行部分，對GPU上的某次迭代更新后，有如圖7所示的剩余矩陣：

其中矩陣A1、A2、A3為剩余待更新矩陣，子塊B1、B2為已求出的Cholesky分解部分。之后，在GPU上按如下公式對A1與A2進(jìn)行更新：

A1=A1-B1×B1T

(10)

A2=A2-B2×B1T

(11)

然后將A1與A2異步傳輸?shù)紺PU。同時，在GPU上對余部矩陣按如下公式進(jìn)行更新：

A3=A3-B2×B2T

(12)

在CPU上再次對A1和A2更新：

A1=potrf(A1)

(13)

A2=trsm(A1,A2)

(14)

其混合并行的示意如圖8所示。

圖8　原Volkov的混合調(diào)度策略

可以發(fā)現(xiàn)在GPU上更新矩陣A1、A2與更新A3是相互獨立的。而且當(dāng)剩余待更新矩陣階數(shù)較大時，可以考慮將A1的首次更新移到CPU上進(jìn)行。注意，對A2的兩次更新不能同時移到CPU上進(jìn)行，因為對A2的兩次更新計算量都比較大，在CPU上計算時，耗時很多，極容易超過GPU上對A3的更新耗時，從而導(dǎo)致GPU的空閑。

此外，當(dāng)剩余待更新的矩陣較小時，在GPU上更新A3的時間較短。此時不適合再將子矩陣A2轉(zhuǎn)移到CPU上進(jìn)行更新，否則會使得GPU出現(xiàn)長時間空閑。

整個混合計算過程中對數(shù)據(jù)傳輸采用異步方式，使得每次數(shù)據(jù)傳輸時間被GPU計算或CPU計算隱藏。

對此，我們綜合考慮了計算的CPU以及GPU的計算性能，以Volkov的混合算法為基礎(chǔ)，給出了一種新的調(diào)度策略。

2.4新的混合調(diào)度算法

記剩余待更新的矩陣A3的階數(shù)為ns。將新的混合調(diào)度算法分為以下三個階段：

階段1當(dāng)ns>x時，更新子矩陣A3耗時較長，因此CPU可以做更多的工作以減少空閑時間。由于對子矩陣塊A2的兩次更新只能有一次移到CPU上進(jìn)行，故有以下2種調(diào)度策略：

1) 將A2的第一次更新放在GPU上執(zhí)行，如圖9所示。

圖9　階段1的混合調(diào)度策略1

2) 將A2的第二次更新放在GPU上執(zhí)行,如圖10所示。

圖10　階段1的混合調(diào)度策略2

通過分析A2兩次更新的計算量，可以發(fā)現(xiàn)A2的前次更新計算量為A2再次更新計算量的兩倍。為防止GPU出現(xiàn)空閑，只有當(dāng)ns>x1(x1>x)時，GPU上對A3的更新任務(wù)能完全隱藏CPU上的計算任務(wù)，采用混合調(diào)度策略2；否則，x

階段2當(dāng)y

階段3當(dāng)ns

圖11　階段3的混合調(diào)度策略

其中x與y的值需要依據(jù)各自計算的CPU和GPU的性能來大致確定，以保證GPU不會出現(xiàn)空閑時間并且CPU的空閑時間最小。這樣的混合調(diào)度策略則是非常合適的。

同時，要對Volkov混合算法的第一步和最后一步進(jìn)行調(diào)整。在第一步中，對第一主列塊的更新與余部矩陣傳輸?shù)紾PU上可以異步進(jìn)行，這樣其中的一個耗時會被隱藏。最后一步更新時，可以先將最后一個子塊先傳回至CPU，然后將更新最后一個子塊與傳回剩余矩陣數(shù)據(jù)異步執(zhí)行。

對于分塊大小nb，可以參考MAGMA軟件包中關(guān)于分塊值的設(shè)定。針對Kerpler構(gòu)架的GPU，可將分塊值設(shè)定如下：

(15)

3　數(shù)值實驗及結(jié)果分析

3.1實驗環(huán)境

本文中描述的GPU計算過程均是在NVIDIA的Kerpler構(gòu)架的GPU上實現(xiàn)的。有關(guān)數(shù)值試驗的硬件和軟件平臺見表1所示。

表1　實驗環(huán)境

3.2實驗結(jié)果與分析

實驗測試1000～10 000階對稱正定矩陣，而且使用雙精度數(shù)據(jù)類型，以保證計算結(jié)果的精確性。

首先分別給出本實驗環(huán)境下MKL中執(zhí)行dgemm以及dtrsm的耗時以及CUBLAS中一執(zhí)行dgemm,dtrsm以及syrk的耗時的一個參考值，如表2所示。該數(shù)據(jù)用于估計劃分階段的x(包括x1)與y值。

表2　各函數(shù)測試結(jié)果(單位：毫秒)

其中針對dgemm函數(shù)，m為表格的第一列，n=k=512。而dtrsm與dsyrk函數(shù)的n為表格第一列，k=512。

表3給出了數(shù)據(jù)矩陣在CPU與GPU之間的傳輸耗時的一個參考值。

表3　數(shù)據(jù)傳輸測試(單位：毫秒)

下面考慮x=6000時，使用階段1的混合調(diào)度策略1是合適的。首先依據(jù)表3可得到，異步傳輸A1,B1(矩陣規(guī)模均為512×512) 以及在CPU上更新A1的時間被GPU上更新A2隱藏。而異步傳輸A2(矩陣規(guī)模6000×512)耗時t1<20 ms，在CPU上再次更新A2耗時t2<80 ms，異步傳回A1,A2耗時t3<20 ms，而GPU上對矩陣A2(規(guī)模6000階)的更新耗時為t≈150 ms>(t1+t2+t3)，故可取x=6000。類似，我們可以由此得到比較合理的x1，x以及y的值，分別如下：

x1=8000x=6000y=4000

最后，對1000～10000階對稱正定矩陣分別使用4種方法進(jìn)行計算。表4給出了4種方法計算1000～10 000階對稱正定矩陣的Cholesky分解的結(jié)果。其中可以看到基于新調(diào)度策略的函數(shù)New_Sche明顯優(yōu)于使用CUBLAS直接計算Cholesky分解以及基于Volkov的混合算法的函數(shù)Volkov_S。

表4　各算法測試結(jié)果(單位：秒)

4　結(jié)　語

科學(xué)計算中使用CPU與GPU混合并行計算的例子層出不窮，然而最困難的在于設(shè)計一種負(fù)載均衡的混合調(diào)度方案，從而使得CPU的空閑時間盡可能小，同時又不能讓GPU空閑。糟糕的混合調(diào)度方案會使得整體計算性能下降，或者使得計算不正確，從而失去了混合計算的目的，這都是需要避免的。合適的混合調(diào)度策略往往需要依據(jù)自身計算機(jī)的CPU以及GPU的計算性能來共同確定。首先細(xì)致分析每步操作在CPU以及GPU上的計算速度，然后合理劃分計算任務(wù)到CPU和GPU上執(zhí)行。再者，算法中會大量使用異步傳輸操作，使得這些數(shù)據(jù)傳輸耗時被CPU和GPU上的計算時間隱藏，這也是優(yōu)化加速的一個重要部分。最后，本文中沒有考慮分塊值對混合調(diào)度策略的影響。由實驗的分析可知，分塊值與調(diào)度策略有直接關(guān)系，影響階段值x1，x以及y的確定。如何確定最優(yōu)的分塊值以及相應(yīng)的調(diào)度策略將是后面要進(jìn)行的工作。

[1] Chandrasekar J,Kim I S,Bernstein D S,et al.Reduced-Rank Unscented Kalman filtering using Cholesky-based decomposition[C]//American Control Conference,June,2008:1274-1279.

[2] Yu H,Chung C Y,Wong K P,et al.Probabilistic Load Flow Evaluation With Hybrid Latin Hypercube Sampling and Cholesky Decomposition[J].IEEE Transactions on Power Systems,2009,24 (2):661-667.

[3] David S.Watkins.Fundamentals of Matrix Computations[M].New York:John Wiley and Sons,2013.

[4] Gene H Golub,Charles F，Van Loan.Matrix Computations [M].Baltimore:Johns Hopkins University Press,2013.

[5] Volkov V,Demmel J W.Benchmarking gpus to tune dense linear algebra[C]//Proceedings of the 2008 ACM/IEEE Conference on Supercomputing,Nov,2008:1-11.

[6] 胡鵬飛，袁志勇，廖祥云，等．基于CPU-GPU混合加速的SPH流體仿真方法[J]．計算機(jī)工程與科學(xué)，2014,36(7):1231-1237.

[7] 張健，焦良葆，陳瑞．CPU-GPU混合平臺上動態(tài)場景光線跟蹤的研究[J]．計算機(jī)工程與應(yīng)用，2012,48(21):151-154.

[8] Yaohung M Tsai,Weichung Wang,RayBing Chen.Tunning Block Size for QR Factorization on CPU-GPU Hybrid Systems[C]//Proceedings of the IEEE 6th International Symposium on Embedded Multicore Socs,Sept,2012:205-211.

[9] John Cheng,Max Grossman,Ty McKercher.CUDA C Programming[M].Indianalpois:John Wiley & Sons,2014.

[10] 劉金碩，鄧娟，周崢，等．基于CUDA的并行程序設(shè)計[M]．北京:科學(xué)出版社，2014．

ACCELERATING CALCULATION OF CHOLESKY FACTORISATION OF MATRIX WITH GPU

Shen CongGao Huotao

(School of Electronic Information,Wuhan University,Wuhan 430072,Hubei,China)

A concrete implementation of Cholesky factorisation on graphic processing unit (GPU) for large real symmetric positive definite matrix is described in this article.We analyse the hybrid parallel algorithm presented by Volkov for computing the Cholesky factorisation in detail.On that basis,and according to the computational performances of CPU and GPU on our own computers,we present a more reasonable hybrid three-phase scheduling strategy,which further reduces the idle time of CPU and avoids the occurrence of GPU in idle status.Numerical experiment shows that the new hybrid scheduling algorithm achieves a speedup of more than 5 times compared with the standard MKL algorithm when the order of a matrix is larger than 7000,and it also observably outperforms the performance of original Volkov’s hybrid algorithm.

GPUCholesky factorisationSpeedupHybrid algorithm

2015-04-10。湖北省自然科學(xué)基金重點項目(ZRZ20 14000286)。沈聰，碩士生，主研領(lǐng)域：并行計算。高火濤，教授。

TP361

A

10.3969/j.issn.1000-386x.2016.09.066