并行對稱矩陣三對角化算法在GPU集群上的有效實(shí)現(xiàn)

劉世芳 趙永華 于天禹 黃榮鋒

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

2(中國科學(xué)院大學(xué) 北京 100190)

稠密對稱特征值問題存在于科學(xué)計(jì)算和工程模擬中的許多領(lǐng)域，例如計(jì)算材料和計(jì)算化學(xué)等科學(xué)計(jì)算應(yīng)用[1-3].規(guī)約稠密對稱矩陣為三對角形式是求解對稱特征問題的關(guān)鍵一步，且占據(jù)了問題求解的大部分時間.因此對稠密對稱矩陣三對角化算法進(jìn)行加速計(jì)算是非常必要的.雖然在LAPACK(linear algebra package)[4]，ScaLAPACK(scalable linear algebra package)[5]等一些傳統(tǒng)的標(biāo)準(zhǔn)庫中都可以獲得對稱矩陣三對角化算法，但硬件發(fā)展促使我們重新設(shè)計(jì)算法以更有效地使用底層硬件.本文基于GPU和CPU的混合架構(gòu)重新設(shè)計(jì)了對稱矩陣三對角化的并行算法的實(shí)現(xiàn)方法，使其可以有效地利用GPU的并行計(jì)算能力，降低特征值求解的計(jì)算時間.

稠密對稱矩陣三對角化一般分為2種方法：1階段法和2階段法.標(biāo)準(zhǔn)的1階段法是直接將稠密對稱矩陣通過雙邊的Householder正交變換轉(zhuǎn)化為三對角矩陣[4].而2階段法[6-7]包括兩步,第1步通過一系列的Householder變換將矩陣規(guī)約成一個帶狀(band)矩陣，這一步運(yùn)用了許多BLAS(basic linear algebra subroutine)-3(矩陣-矩陣)操作，因此取得了較好的效果[7-9]，第2步將帶狀矩陣轉(zhuǎn)化為三對角矩陣[10].在2階段法中，由于第1步使用了更多的BLAS-3操作，通常2階段法要比1階段法更快.然而，由于2階段法中第2步的計(jì)算不能充分發(fā)揮GPU的性能，并且若要求解原始問題特征向量，2階段法比1階段法多需要1個額外的轉(zhuǎn)化計(jì)算，因此在整體上所需的時間要比1階段法更多[11].因此，在本文中，我們采用1階段法.

在1階段法的三對角化塊算法中，主要分為2部分：條塊內(nèi)規(guī)約與條塊外尾部矩陣的修正.主要的計(jì)算量主要集中在條塊內(nèi)規(guī)約時的對稱矩陣向量積和條塊外尾部矩陣的修正.已經(jīng)有一些工作利用GPU對三對角化中這2部分加速計(jì)算[11-12]，文獻(xiàn)[11-12]采用了僅對尾部對稱矩陣的上或下三角部分進(jìn)行修正的方法，該方法降低了修正過程中的計(jì)算量，但使得修正后的矩陣失去了完整性，導(dǎo)致了條塊內(nèi)規(guī)約過程中對稱矩陣向量積計(jì)算變?yōu)?個不規(guī)則的三角矩陣向量積計(jì)算，這種不規(guī)則的計(jì)算不能充分利用GPU的并行計(jì)算能力.針對GPU集群的架構(gòu)特點(diǎn)，本文給出的算法改進(jìn)包括4個方面：

1) 給出的三對角化算法是對整個矩陣進(jìn)行操作，這避免了條塊內(nèi)規(guī)約過程中不規(guī)則矩陣向量乘運(yùn)算,但這也帶來了一些冗余計(jì)算：在條塊內(nèi)修正為規(guī)約部分時，不僅要修正列條塊內(nèi)為規(guī)約部分，也要修正行條塊內(nèi)未規(guī)約部分.但從整體上減少了算法的用時，增加了性能.

2) 增加計(jì)算密集度.在條塊內(nèi)規(guī)約中將2個小的長條矩陣-向量積合并成1個更大的矩陣-向量積.并且，在條塊外尾部矩陣的修正中，將調(diào)用2次矩陣-矩陣乘的秩-b修正合并成1個調(diào)用1次矩陣-矩陣乘秩-2b修正.加大計(jì)算密集度能更加充分發(fā)揮GPU的計(jì)算能力，增加GPU的利用率，從而提升算法的性能.

3) 利用多個CUDA流同時啟動多個CUDA內(nèi)核，使得在三對角化算法中獨(dú)立的CUDA操作可以在不同的CUDA流中并發(fā)執(zhí)行.另外，在基于MPI+CUDA的并行算法中，利用了CPU與GPU之間的異步數(shù)據(jù)傳輸操作，使得在不同CUDA流中的數(shù)據(jù)傳輸操作與GPU內(nèi)的核函數(shù)同時執(zhí)行，隱藏了數(shù)據(jù)傳輸?shù)臅r間，使得并行算法的性能進(jìn)一步提高.

4) 采用完全并行的點(diǎn)-點(diǎn)數(shù)據(jù)通信方式.在并行算法中，由在處理器列上的向量獲得對應(yīng)的在處理器行上的向量操作，若使用全局通信，則需要將屬于在處理器列上的各個列處理的部分向量進(jìn)行全局求和操作，這是一個全局的非并行的通信操作.本文給出的處理器行列之間的向量轉(zhuǎn)置算法是一個完全并行的點(diǎn)-點(diǎn)的并行通信形式：處理器列上的各個處理器并行地將屬于本處理器的部分向量發(fā)送到對應(yīng)的處理器行上的處理器，處理器行上的各個處理器并行地接收發(fā)送來的數(shù)據(jù).這有效地改善了三對角化算法實(shí)現(xiàn)的通信性能.

本文三對角化算法的實(shí)現(xiàn)是在我們開發(fā)的PSEPS(parallel symmetric eigenproblem solver)特征求解器軟件包的基礎(chǔ)上完成的.我們從多個方面優(yōu)化了三對角化算法實(shí)現(xiàn)，使其主要計(jì)算部分并行地在多個GPU上進(jìn)行加速計(jì)算.計(jì)算部分通過調(diào)用NVIDIA的cuBLAS(CUDA basic linear algebra subroutine)庫完成，從而有效地保證了多GPU實(shí)現(xiàn)的性能.在中國科學(xué)院超級計(jì)算機(jī)系統(tǒng)“元”上，使用Nvidia Tesla K20 GPGPU卡(峰值雙精度浮點(diǎn)性能1.17 TFLOPS)測試了算法的實(shí)現(xiàn)性能，取得了較好的加速效果與可擴(kuò)展性.在對30 000×30 000階矩陣?yán)肎PU進(jìn)行加速計(jì)算時，文獻(xiàn)[8]中利用4個NDIVIA Tesla C2050 GPUs(峰值雙精度浮點(diǎn)性能515.2GFLOPS)加速時，性能約為115GFLOPS，約達(dá)到峰值性能的22.32%.我們利用4個Nvidia Tesla K20 GPUs時性能約為376.60 GFLOPS，約達(dá)到峰值性能的32.19%.

1 三對角化串行塊算法

T=Hn-2…H2H1AH1H2…Hn-2

為三對角對稱矩陣.

(I-βuuT)A(I-βuuT)=A-uwT-wuT，

Fig. 1 Panel reduction and trailing submatrix correction圖1 條塊規(guī)約與尾部矩陣修正

塊算法將尾部矩陣的秩-2修正聚合為秩-b修正.A22的延遲修正使向量v的生成變得復(fù)雜，因?yàn)榇鎯υ诰仃嘇22中的元素在條塊規(guī)約過程中沒有及時更新.因此，當(dāng)求解w時，v的計(jì)算需要替換為

v=β(A22-UWT-WUT)u，

其中,UWTu與WUTu對應(yīng)著尾部矩陣的修正.進(jìn)一步：

同樣，在尾部矩陣A22的修正中：

這將需要2次矩陣乘法的秩-b修正合并成了只需要1次矩陣乘法的秩-2b修正.在實(shí)現(xiàn)中，我們可以將向量w和u交叉存儲在矩陣M和N中：

M=(u1,w1,u2,w2,…,ub,wb)，
N=(w1,u1,w2,u2,…,wb,ub).

將這些運(yùn)算合并成一個更大的計(jì)算，這使得計(jì)算更加密集，更適合GPU的架構(gòu).密集計(jì)算可以更加充分地利用GPU的計(jì)算能力，從而提升了三對角化算法的性能.

算法1.利用GPU加速計(jì)算的三對角化串行塊算法.

劃分矩陣,其中dATL是0×0階矩陣，置i=1;

do untildABR是0×0階矩陣{

確定列條塊的大小：bc=min{b,dABR的列寬}；

重新劃分矩陣,其中dA01，dA11，dA21的列數(shù)為bc:

/*條塊內(nèi)規(guī)約*/

/*將u從GPU傳到CPU*/

Synchronize(stream1);

Step1. 計(jì)算(u,β,μ)=House(u):

β→φi;μ→zi;i=i+1;

(stream1)u→du;

/*將向量u傳到GPU*/

/*將du保存在dA中*/

Step2. 計(jì)算矩陣向量積：

Step3. 修正矩陣向量積:

Step3.1. (stream2)dvt=dNTdu;

Synchronize(stream1)；

Step3.2. (stream2)dv=dv-βdMdvt;

Step4. 計(jì)算修正向量:

Synchronize(stream2);

Step5. 將du和dw交叉保存到矩陣dM和dN中：

Step5.1. (stream1)dM=[dM|du|*];

Step5.2. (stream2)dM=[dM|*|dw];

Step5.3. (stream3)dN=[dN|dw|*];

Step5.4. (stream4)dN=[dN|*|du];

Step6. 更新條塊內(nèi)未規(guī)約部分：

Step6.1. 更新列條塊：

Step6.2. 更新行條塊：

Synchronize(stream1&stream2&stream3

&stream4&stream5&stream6);

對下列矩陣形式繼續(xù)上述過程：

}end while

Step7. 對尾部矩陣dA22執(zhí)行秩-2b修正：

(stream1)dA22=dA22-dNdMT;

對下列矩陣形式繼續(xù)上述過程:

}end do

在算法1內(nèi)的條塊內(nèi)規(guī)約過程中，可以發(fā)現(xiàn)Step2和Step3.1以及Step5.1～Step5.4，Step6.1～Step6.2這些操作互不依賴，每個操作可以看成是1個獨(dú)立的任務(wù)，這些任務(wù)可以并發(fā)執(zhí)行.在利用GPU進(jìn)行加速計(jì)算的三對角化算法中，可以將每一個CUDA流視為GPU上的1個任務(wù)，這樣通過啟動多個CUDA流，使得這些獨(dú)立的任務(wù)在不同的CUDA流中并發(fā)執(zhí)行，使時間進(jìn)一步壓縮，從而提升三對角化算法的性能.注意在利用多個CUDA流并發(fā)啟動的三對角化算法中，在有數(shù)據(jù)依賴的位置需要進(jìn)行同步操作，以保證算法的正確性，如圖2所示:

Fig. 2 Parallel execution with multiple CUDA streams in the panel reduction圖2 在條塊規(guī)約內(nèi)利用多個CUDA流并發(fā)執(zhí)行

2 基于MPI+CUDA的三對角化并行塊算法

2.1 數(shù)據(jù)在處理器上的分布方式

MPI并行算法的效率較大地依賴數(shù)據(jù)在處理器間的移動代價(jià),這個代價(jià)又與數(shù)據(jù)在處理器上的初始映射(分布)和處理器間通信所使用的算法有關(guān),并且數(shù)據(jù)在處理器的分布方式極大地影響著算法的負(fù)載平衡和通信特性.在稠密線性代數(shù)中，當(dāng)前普遍使用的是2維塊循環(huán)分布.2維塊循環(huán)具有支持可拓展的高性能并行稠密線性代數(shù)編碼的特性[14-16]，并被HPF(high performance fortran handbook)[17]，ScaLAPCK[5]和其他并行稠密線性代數(shù)庫選為首要的數(shù)據(jù)分布.因此在基于MPI+CUDA的并行三對角化塊算法中，我們采用了2維塊循環(huán)的數(shù)據(jù)分布方法.將稠密對稱矩陣A采用2維塊循環(huán)方式分布到p=pr×pc處理器網(wǎng)格中，存儲塊的大小為b.這里，假設(shè)矩陣A的元素A(1,1)映射到處理器p(0,0)中的本地矩陣元素為Aloc(1,1)，則全局元素A(i,j)映射到處理器

中的本地矩陣元素為Aloc(iloc,jloc)，其中

2.2 MPI通信方式

在2維通信域中，處理器間的通信通常發(fā)生在處理器網(wǎng)格中同行處理器之間或同列處理器之間.在行、列映射方式下矩陣的整行或整列被播送.基于MPI+CUDA并行三對角化算法中的通信主要操作包括行播送、列播送、行聚集求和、轉(zhuǎn)置操作.

1) 行播送.2維通信域中處理器行內(nèi)某處理器將數(shù)據(jù)廣播到該處理器所在行的所有處理器，在算法2中此操作用Rbroadcast表示.

2) 列播送.2維通信域中處理器列內(nèi)某處理器將數(shù)據(jù)廣播到該處理器所在列的所有處理器，在算法2中用Cbroadcast表示.

3) 行聚集求和.在處理器行上對數(shù)據(jù)進(jìn)行聚集求和，在算法2中用Rsum表示.

4) 轉(zhuǎn)置操作.將處理器行/列上的向量轉(zhuǎn)置到對應(yīng)的處理器列/行上的向量.在算法2中用Transpose表示.

2.3 基于MPI+CUDA的并行塊算法

基于串行的利用GPU加速計(jì)算的的三對角化塊算法1，算法2給出了基于MPI+CUDA的并行三對角化塊算法在GPU集群上的實(shí)現(xiàn).在算法2中，(GC-數(shù)字)表示數(shù)據(jù)從GPU傳到CPU，(CG-數(shù)字)表示數(shù)據(jù)從CPU傳到GPU.(CC-數(shù)字)表示CPU之間的通信操作.(G-數(shù)字)表示在GPU端的操作.算法2中的變量前邊加“d”表示在GPU端對應(yīng)的變量.

在基于MPI+CUDA的并行算法2中，除了將獨(dú)立的CUDA操作在不同的CUDA流中進(jìn)行運(yùn)算，還利用了CPU與GPU之間的異步數(shù)據(jù)傳輸操作，使得在不同CUDA流中的數(shù)據(jù)傳輸操和GPU內(nèi)的核函數(shù)同時執(zhí)行： (G-3)GPU內(nèi)的拷貝操作與(GC-3)數(shù)據(jù)傳輸操作，(CG-4)數(shù)據(jù)傳輸操作與(G-4) GPU內(nèi)的計(jì)算操作.這可以隱藏CPU與GPU之間數(shù)據(jù)傳輸?shù)臅r間，使得基于MPI+CUDA的并行算法的性能進(jìn)一步提高.

算法2.基于MPI+CUDA的并行三對角化塊算法.

/*算法中的dot表示向量的點(diǎn)乘運(yùn)算，scal表示向量的數(shù)乘運(yùn)算，copy表示向量的拷貝操作，gemv表示矩陣-向量乘運(yùn)算，axpy表示向量的加法運(yùn)算，ger表示矩陣的一階修正運(yùn)算，gemm表示矩陣-矩陣乘運(yùn)算*/

for eachbt∈[1,bn]

as=bt×b;/*條塊內(nèi)第一列*/

ae=min(as+b-1,n-2);/*最后一列*/

for eachk∈[as,ae] /*條塊內(nèi)規(guī)約*/

Step1. 計(jì)算Householder向量及其范數(shù):

ifk∈CSthenCS=CS(k) endif

ifk∈RSthen

Step1.1. (GC-1) (stream1)dA(k,CS)→uc;

Synchronize(stream1);

Step1.2. (CC-1)Cbroadcast(uc);

Step1.3. (CG-2) (stream1)uc→duc;

endif

/*dot*/

Synchronize(stream1);

Step1.5. (CC-2)t=Rsum(g);

ifj+1∈CSthen /*行播送*/

Step1.6. (CC-3)Rbroadcast(u(k+1));

endif

Step1.7. 計(jì)算次對角元素μ和范數(shù)β

Step1.8. (G-2) (stream1)duc=βduc;

/*scal*/

ifk∈RSthen

Synchronize(stream1);

Step1.9. (G-3) (stream2)duc→dA(k,CS);

/*copy*/

endif

ifk∈RSthenRS=RS(k) endif

Step1.10. (GC-3) (stream1)duc→uc;

Synchronize(stream1);

Step1.11. (CC-4)ur=Transpose(uc);

ifk+1∈CSthen

Step1.12. (CC-5)Rbroadcast(ur);

Step1.13. (CG-4) (stream2)ur→dur;

endif

Step2.計(jì)算矩陣向量積(gemv):

(G-4) (stream1)dxr=βdA(RS,CS)duc;

Step3. 修正矩陣向量積:

Step3.1. (G-5) 計(jì)算部分向量dxr(gemv×2)

Synchronize(stream1&stream2);

(G-5.2) (stream0)dxr-=dMrdxt;

Step3.2. (GC-5) (stream0)dxr→xr;

Synchronize(stream0);

Step3.3. (CC-6)wr=Rsum(xr);

Step3.4. (CG-6) (stream2)wr→dwr;

Step3.5. (CC-7)wc=Transpose(wr);

ifk∈RSthen

Step3.6. (CC-8)Cbroadcast(wc);

Step3.7. (CG-7) (stream1)wc→dwc;

endif

Step4. 修正矩陣向量積:

/*dot*/

Synchronize(stream1);

Step4.2. (CC-9)σ=Rsum(s);

Step4.3. (G-7) 計(jì)算dwc,dwr(axpy×2);

Synchronize(stream1&stream2);

Step5.(G-8) 將向量dwr，dwc，dur，duc交叉保存到矩陣dMr,dNc的相應(yīng)位置(copy×4):

(G-8.1) (stream1)dMr=[dMr|dur|*];

(G-8.2) (stream2)dMr=[dMr|*|dwr];

(G-8.3) (stream3)dNc=[dNc|dwc|*];

(G-8.4) (stream4)dNc=[dNc|*|duc];

Step6. (G-9)更新條塊內(nèi)未規(guī)約部分:

(G-9.1) (stream5) /*ger×2*/

(G-9.2) (stream6) /*ger×2*/

Synchronize(stream1&stream2&stream3&

stream4&stream5&stream6);

endfor

Step7. 對尾部矩陣執(zhí)行秩-2b修正(gemm):

Synchronize(stream0);

endfor

Fig. 3 Operations of Algorithm 2圖3 算法2的操作

圖4描述了將處理器列上的向量vc轉(zhuǎn)置到對應(yīng)的處理器行上的向量vr的過程，每一個列處理器并行的將部分向量發(fā)送到對應(yīng)的行處理器上，行處理器并行地接收部分向量：pse,0→pre,0，pse,1→pre,1，pse,2→pre,0，pse,3→pre,3.將處理器行上的向量轉(zhuǎn)置到對應(yīng)的處理器列上的操作類似，這里不再討論.

Fig. 4 Transpose operation圖4 轉(zhuǎn)置操作

CPU+GPU的異構(gòu)編程需要CPU與GPU之間的數(shù)據(jù)傳輸，要想利用GPU進(jìn)行加速計(jì)算，需要將計(jì)算數(shù)據(jù)上傳到GPU端.因?yàn)槭褂玫亩鄠€GPU分布在多個節(jié)點(diǎn)上，則GPU間的通信需要在CPU中通過MPI進(jìn)行，首先將要通信的數(shù)據(jù)從GPU端傳回到CPU端，之后通過MPI通信，通信完成后，再將數(shù)據(jù)傳回到相應(yīng)GPU中.算法2中CPU與GPU的數(shù)據(jù)傳輸操作如圖2所示，共有7次.

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

因?yàn)镕rank矩陣

A=(aij),aij=n-max(i,j)-1,
i,j=1,2,…,n

具有解析特征值：

所以為了驗(yàn)證算法實(shí)現(xiàn)的準(zhǔn)確性和測試算法的性能，我們使用了Frank矩陣作為測試數(shù)據(jù).

我們在中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)中心超級計(jì)算機(jī)系統(tǒng)“元”上實(shí)現(xiàn)了本文描述的算法，并在GPU隊(duì)列上進(jìn)行測試.每臺GPGPU計(jì)算節(jié)點(diǎn)配置2塊Nvidia Tesla K20 GPGPU卡(峰值雙精度浮點(diǎn)性能(板卡)1.17 TFLOPS)，2顆Intel E5-2680 V2(Ivy Bridge|10C|2.8 GHz)處理器.使用的cuBLAS版本是6.0.37，使用的MPI是openmpi-1.8.4.在測試時，為了使得多個GPU并行加速計(jì)算，每個節(jié)點(diǎn)申請2個進(jìn)程，1個進(jìn)程對應(yīng)1個GPU.測試的數(shù)據(jù)類型是雙精度浮點(diǎn)數(shù).

3.1 基于MPI+CUDA的三對角化并行塊算法性能分析

表1展示了20 000×20 000階矩陣三對角化塊算法在CPU調(diào)用上MKL(math kernel library) 2013的BLAS庫和通過多個CUDA流并發(fā)調(diào)用GPU的cuBLAS庫使用不同進(jìn)程數(shù)(1，2，4，8，16，32)時的時間對比.從表1中可以看出，首先利用GPU進(jìn)行加速計(jì)算的三對角化算法取得了明顯的加速效果.其次，隨著使用GPU數(shù)目的增多，加速倍數(shù)在不斷減小，這是因?yàn)殡S著隨著使用GPU數(shù)目的增多，在單個GPU上負(fù)載減小，并且隨著進(jìn)程數(shù)目的增多，通信時間、CPU與GPU之間的數(shù)據(jù)時間等開銷占比增大，因此計(jì)算加速效果不如負(fù)載大的好.總的來說，利用GPU對三對角化算法進(jìn)行加速取得了較好的加速效果(這里b=64).

Fig. 5 Peak performance of tridiagonal parallel algorithm圖5 并行三對角化算法的峰值性能

圖5展示了當(dāng)b=64時不同規(guī)模的矩陣使用不同GPU個數(shù)(1，2，4，8，16，32)時達(dá)到的性能.從圖5中可以發(fā)現(xiàn)，在使用GPU的數(shù)目一樣時，矩陣規(guī)模越大達(dá)到的性能越高，這是因?yàn)閱栴}規(guī)模越小，GPU的負(fù)載越小，而GPU適合計(jì)算密集型任務(wù)，這樣加速效果就不如規(guī)模大的好.并且規(guī)模越大，計(jì)算時間所占的比例越大，加速效果越明顯，則達(dá)到的性能越高.從圖5可以看出，基于MPI+CUDA的三對角化并行塊算法取得了較好的性能.

圖6展示了不同矩陣規(guī)模(10 000×10 000，15 000×10 000，20 000×20 000，25 000×25 000，30 000×30 000)使用不同個數(shù)的GPU進(jìn)行加速計(jì)算時相對于使用1個進(jìn)程/GPU時加速比曲線.

Fig. 6 Speedup curves for different matrix size圖6 不同矩陣規(guī)模的加速比

從圖6中可以看出當(dāng)矩陣規(guī)模一定時，隨著使用GPU數(shù)目增加，加速比增加的幅度在不斷降低，這是由于隨著使用GPU個數(shù)增加，GPU與CPU之間的數(shù)據(jù)傳輸、通信等開銷所占的比例增加導(dǎo)致的.當(dāng)使用的GPU的數(shù)目一樣時，較大規(guī)模矩陣的加速比大于較小規(guī)模矩陣的加速比，加速比隨著矩陣規(guī)模的增大而增大，這是因?yàn)榇笠?guī)模的矩陣在GPU上的負(fù)載較大，比小規(guī)模矩陣更能夠充分地利用GPU的計(jì)算能力.并且，隨著使用GPU數(shù)目的增多，較大規(guī)模的加速比與較小規(guī)模的加速比相差越來越大，這是因?yàn)樵趩蝹€GPU上的負(fù)載越大，取得的性能越好.總的來說，基于MPI+CUDA的三對角化并行塊算法具有良好的可擴(kuò)展性.

3.2 加大計(jì)算密集度的算法性能分析

不同的條塊(矩陣存儲塊)的大小會影響三對角化算法的性能，當(dāng)三對角化算法只啟動1個CUDA流時，分別測試了如表2所示的當(dāng)b分別取值32，64，128時的不同規(guī)模的矩陣(10 000×10 000，20 000×20 000，30 000×30 000)利用不同GPU個數(shù)(2，4，8，16)加速時的用時.可以看出，當(dāng)b=64時用時最短，因此在只啟動1個CUDA流的三對角化算法的分析中，取b=64.

表3展示了在三對角化算法中啟動1個CUDA流時20 000×20 000階矩陣使用1個GPU加速時各個計(jì)算部分在GPU上的用時.其中，算法2中(G-4)gemv計(jì)算條塊內(nèi)規(guī)約的矩陣向量積v=Au的計(jì)算量約為4n3/3，占據(jù)了在GPU上計(jì)算總時間的約86.62%，因此，采用完整的矩陣進(jìn)行運(yùn)算提升這部分的性能，將會提升整個程序的性能.

表4展示20 000×20 000階矩陣的使用不同GPU個數(shù)(4，8，16，32)時分別采用對上三角和下三角直接調(diào)用2次cublasDgemv與使用完整矩陣只調(diào)用1次cublasDgemv的時間對比，可以發(fā)現(xiàn)采用完整矩陣的方式這部分時間減少了約一半，又因?yàn)檫@部分占了三對角計(jì)算時間的大部分，這部分性能的提升將會大幅度提升整個算法的性能.但采用完整矩陣的方式將增加(G-9)ger的計(jì)算量，表5展示了增加前后的時間對比，雖然這部分的用時增加了，但是相對于(G-4)gemv時間的減少這部分的增加的時間非常少，因此整個程序的性能相對提升了.

對條塊內(nèi)第(G-5)gemv步進(jìn)行優(yōu)化，將原來需要4次矩陣向量積的運(yùn)算合并成只需要2次矩陣向量積的運(yùn)算.表6展示了20 000×20 000階矩陣的使用不同GPU個數(shù)(1，2，4，8，16，32)時合并前后的時間對比，從表6中可以發(fā)現(xiàn)，將矩陣-向量運(yùn)算合并后，這部分的時間減少了大約42%.這是因?yàn)閷⑿∫?guī)模的矩陣向量乘運(yùn)算合并成較大規(guī)模的矩陣向量乘運(yùn)算后，更加充分地利用了GPU的資源，因此程序的性能得到了進(jìn)一步地提升.

(G-10)gemm條塊外尾部矩陣的秩-2b修正計(jì)算量約為4n3/3，雖然這部分的計(jì)算量與(G-4)gemv的計(jì)算量一樣，但占據(jù)的時間僅為7.27%(20 000×20 000階矩陣使用1個GPU加速時)，這2部分時間占比相差非常大是因?yàn)锽LAS-2與BLAS-3對GPU的利用率不同導(dǎo)致的.在(G-10)gemm中我們將需要2次gemm的秩-b修正合并成調(diào)用1次gemm的秩-2b修正，表7展示了20 000×20 000階矩陣使用不同GPU個數(shù)(1，2，4，8，16，32)時優(yōu)化前后的時間對比，從表7可以發(fā)現(xiàn)，將秩-b修正合并為秩-2b修正后，此部分的時間約減少了23%.合并計(jì)算后使得計(jì)算更加密集，更加充分地利用了GPU的計(jì)算能力，提高了GPU的利用率，因此獲得了更高的性能.

Fig. 7 The time of each part of tridiagonal parallel algorithm圖7 并行三對角化算法各部分的時間

圖7展示了在三對角化算法中啟動1個CUDA流時20 000×20 000階矩陣使用不同GPU數(shù)目(2，4，8，16，32)加速時三對角化算法各個部分所用占用的時間(compu_gpu表示在GPU上的計(jì)算時間；memcpy表示CPU與GPU之間的數(shù)據(jù)傳輸時間；ctime表示CPU之間的通信時間；others表示其他時間).從圖7中可以看出,在整個三對角化的過程中，計(jì)算時間占據(jù)了大部分，因此利用GPU進(jìn)行加速計(jì)算能夠明顯提升三對角化算法的性能.并且，隨著GPU的數(shù)目增多，計(jì)算時間compu_gpu逐漸遞減.同時，通信時間ctime隨著處理器數(shù)目(從2個開始)的增加也在逐漸遞減.總的來說，整個基于MPI+CUDA的三對角化算法的用時隨著使用GPU數(shù)目的增多逐漸遞減，整個算法表現(xiàn)出了良好的強(qiáng)擴(kuò)展性.

3.2 利用多個CUDA流的算法性能分析

在通過多個CUDA流并發(fā)調(diào)用多個GPU內(nèi)核的三對角化算法中，我們同樣測試了條塊(矩陣存儲塊)b取不同值(32，64，128)時對不同規(guī)模的矩陣(10 000×10 000，20 000×20 000，30 000×30 000)利用不同GPU個數(shù)(2，4，8，16)加速時的影響，如表8所示，當(dāng)b=64時用時最短，因此在啟動多個CUDA流的三對角化算法的分析中，取b=64.

表9展示了在基于MPI+CUDA的并行三對角化塊算法中啟動1個CUDA流與啟動多個CUDA流時，20 000×20 000階矩陣使用不同進(jìn)程數(shù)(1，2，4，8，16，32)時的時間對比，可以發(fā)現(xiàn)啟動多個CUDA流相比于只啟動1個CUDA流用時更短，取得了較好的加速效果.另外隨著進(jìn)程數(shù)目的增加，加速效果越來越好，這是因?yàn)殡S著進(jìn)程數(shù)目的增加，在單個GPU上的負(fù)載不斷減小，GPU中可用的SM(streaming multiprocessor)等資源越多，多個流并發(fā)執(zhí)行的可能性越大，因此規(guī)模越小，加速效果越好.

Table 9 Time Comparison Between Starting a CUDA Stream

Table 1 Time Comparison of Using CPU or GPU表1 使用CPU或GPU的耗時比較

Table 2 Time Comparison for Different Matrix Size with b Taking Different Size表2 不同規(guī)模的矩陣當(dāng)b取不同大小時的耗時比較 s

Table 3 Time and Proportion of Each Calculation Part表3 各個部分計(jì)算部分用時和占比

Table 4 Time Comparison of Calling cublasDgemv Twice and Once表4 調(diào)用2次cublasDgemv與1次cublasDgemv的 耗時比較 s

Table 5 Time Comparison of Before and After (G-9)ger Increase表5 (G-9)ger增加前后耗時比較 s

Table 6 Time Comparison Before and After Merge of the(G-5)gemv Within the Panel Reduction表6 條塊內(nèi)規(guī)約(G-5)gemv合并前后耗時比較

Table 8 Time Comparison for Different Matrix Size with b Taking Different Size表8 不同規(guī)模的矩陣當(dāng)b取不同大小時時間對比 s

Table 10 Time Comparison Between Starting a CUDA Stream and Multiple CUDA Streams of Different Scale Matrices When Using 32 GPUs表10 不同規(guī)模的矩陣使用32個GPU加速時啟動1個 CUDA流與多個CUDA流的時間比較

and Multiple Streams of 20 000×20 000 Order Matrix

表9 20 000×20 000階矩陣啟動1個CUDA流與多個CUDA流的時間比較

Number ofGPUsTime∕sOne StreamMultiple StreamsReduction∕%1213.21205.513.612117.49109.916.45468.6863.487.57843.0139.597.961628.7625.5811.073221.6219.1011.66

表10展示了不同規(guī)模的矩陣(10 000×10 000，20 000×20 000，30 000×30 000，40 000×40 000，50 000×50 000)利用32個GPU進(jìn)行加速計(jì)算時啟動1個CUDA流與啟動多個CUDA流時的時間比較，可以發(fā)現(xiàn)隨著矩陣規(guī)模的增大，加速效果不斷減小.這是因?yàn)椋m然從軟件的角度來看，CUDA操作在不同的流中并發(fā)運(yùn)行；而從硬件上來看，不一定總是如此.根據(jù)每個SM資源的可用性，完成不同的CUDA流可能仍然需要互相等待.矩陣規(guī)模的越大，GPU上的負(fù)載越大，而GPU中SM資源是有限的，規(guī)模越大，需要的資源越多，因此多個內(nèi)核同時并發(fā)執(zhí)行的可能性就越小.并且，多個CUDA流并發(fā)執(zhí)行，所用的總時間主要去取決于用時最多的那個任務(wù)，利用多個CUDA流并發(fā)調(diào)用多個內(nèi)核，并沒有減少用時最長的任務(wù)所用的時間，在三對角化算法中，主要是條塊內(nèi)規(guī)約中的矩陣向量積(G-4)gemv操作，矩陣規(guī)模越大，這部分的計(jì)算時間占的比例越多，因此對于整個三對角化算法，利用多個CUDA流并發(fā)調(diào)用多個內(nèi)核，問題規(guī)模越大加速效果反而不如規(guī)模小的明顯.總的來說，在單個GPU上的負(fù)載越小，利用多個CUDA流進(jìn)行并發(fā)操作提升的性能越明顯.

4 結(jié)束語

本文描述了基于MPI+CUDA的稠密對稱矩陣三對角化并行塊算法在GPU集群上實(shí)現(xiàn)，主要的計(jì)算部分由GPU完成，而通信部分在CPU上通過MPI來完成，從而充分發(fā)揮CPU+GPU的異構(gòu)平臺的優(yōu)勢.為了適應(yīng)GPU架構(gòu)，我們重構(gòu)了算法，在條塊內(nèi)規(guī)約過程中的對稱矩陣向量相乘時使用了完整的矩陣，雖然帶來了額外的計(jì)算，但比直接分別對上三角和下三角矩陣進(jìn)行2次矩陣-向量相乘時時間明顯減少.同時，在條塊內(nèi)規(guī)約中合并了一些矩陣-向量積運(yùn)算，這部分的時間減少了約42%.而且，在條塊外尾部矩陣的修正過程中，將矩陣的秩-b修正合并為秩-2b修正，此部分的時間約減少了23%.本文提出的三對角化算法通過加大計(jì)算密集度來發(fā)揮GPU的計(jì)算能力，增加GPU的利用率，從而提升了算法的性能.另外，在三對角化算法中，對在GPU中互不依賴的操作，利用了多個CUDA流同時啟動多個內(nèi)核，使這些獨(dú)立的CUDA操作可以在不同的流中并發(fā)執(zhí)行，使得時間進(jìn)一步壓縮.并且，在基于MPI+CUDA的并行三對角化算法中，利用了CPU與GPU之間的異步數(shù)據(jù)傳輸，使得在不同流中的數(shù)據(jù)傳輸與核函數(shù)同時執(zhí)行，隱藏了數(shù)據(jù)傳輸?shù)臅r間，進(jìn)一步提升了并行算法的性能.在三對角化并行算法中，塊的大小將會影響通信次數(shù)、負(fù)載平衡以及利用GPU計(jì)算的性能，在使用時需要權(quán)衡.總的來說，本文中的異構(gòu)三對角化算法取得了較好的性能與加速效果，并且具有良好的可擴(kuò)展性.