黃敏 丁萍, 羅海飚

(1.華南理工大學(xué) 軟件學(xué)院， 廣東 廣州 510006； 2.廣州中國科學(xué)院軟件應(yīng)用技術(shù)研究所 智能視頻實驗室， 廣東 廣州 511458)

共軛梯度法在GPU及Xeon Phi下的并行優(yōu)化及比較*

黃敏1丁萍1,2羅海飚2

(1.華南理工大學(xué) 軟件學(xué)院， 廣東 廣州 510006； 2.廣州中國科學(xué)院軟件應(yīng)用技術(shù)研究所 智能視頻實驗室， 廣東 廣州 511458)

為了充分利用多核處理器的強大計算能力并滿足具有高并行度應(yīng)用的需求，提出一種基于大規(guī)模稀疏矩陣特征問題求解的并行共軛梯度算法.對圖形處理器(GPU)上的計算，有效利用GPU多層次的存儲器體系，采用線程與矩陣映射、數(shù)據(jù)合并訪問、數(shù)據(jù)復(fù)用等優(yōu)化手段，并通過高效的線程調(diào)度來隱藏全局存儲器的高延遲訪問；對Xeon Phi處理器上的計算，有效利用Xeon Phi的高并行度計算對數(shù)據(jù)通信/傳遞、減少數(shù)據(jù)依賴、向量化、異步計算等進行優(yōu)化，并通過高效的線程調(diào)度來隱藏全局存儲器的高延遲訪問.文中還通過實驗驗證了算法的可行性和正確性，并對比了不同方式下的運行效率，發(fā)現(xiàn)共軛梯度法在GPU下比在Xeon Phi下的加速效果更好.

共軛梯度法；圖形處理器；Xeon Phi；并行優(yōu)化；稀疏矩陣向量乘

城市大氣污染的湍流模擬[1- 2]、視頻處理、流體和力學(xué)計算、生物醫(yī)藥分析等科學(xué)計算中牽涉到很多大規(guī)模數(shù)值模擬，傳統(tǒng)的串行計算方式無法滿足大規(guī)模數(shù)值模擬對計算速度的要求，如何提高求解速度是科研人員關(guān)注的焦點[3].在城市空氣污染的湍流模擬中，如何高效求解預(yù)條件方程是提高迭代求解共軛梯度算法性能的關(guān)鍵[4].本研究的應(yīng)用背景是模擬城市區(qū)域多尺度的百萬網(wǎng)格的計算[5- 6]，在建立有限元的數(shù)值模型后，采取雅可比迭代法[7- 8]求解大型稀疏矩陣線性方程組[9- 10].已有實驗結(jié)果表明，稀疏矩陣向量乘[11](SpMV)占了算法計算總時間的80%以上[12]，成為整個計算過程的瓶頸[13].

目前國內(nèi)外高性能的并行計算主要基于多個中央處理器(CPU)[14]，作為底層計算的稀疏矩陣向量乘和線性方程組也不例外.國內(nèi)外的很多研究都是基于多CPU處理器下的SpMV和線性方程組的并行優(yōu)化[15].針對不同的稀疏矩陣的特征，有學(xué)者提出了適宜的存儲格式，這些存儲格式可以節(jié)省存儲空間或加速讀取[16].此外，數(shù)據(jù)預(yù)取、分塊技術(shù)、數(shù)據(jù)復(fù)用和矩陣重組等優(yōu)化技術(shù)也可以應(yīng)用于其中[17].袁娥等[18]通過將稀疏矩陣分成一個個小的稠密塊，使之適宜用塊壓縮存儲(BCSR)格式進行存儲，同時采用了寄存器分塊和啟發(fā)式分塊大小選擇算法，使得該SpMV計算內(nèi)核的性能大為提高.

由于圖形處理器(GPU)的性能得到較大提高，一些研究者開始探討稀疏矩陣向量乘和線性方程組在GPU上的優(yōu)化，目前這方面的研究已經(jīng)比較成熟，例如，Bell等[19]檢驗了并列格式(COO格式)、行壓縮存儲格式(CSR格式)和ELLPACK格式(ELL格式)，提出了新的混合格式(HYB格式)，且在公共測試集上驗證了所提格式相對于其他格式的有效性和優(yōu)越性.Yuan和孫相征等[20- 21]在研究非零元素值被限制在少量的矩陣對角線上的情況時，針對同時具有“帶狀”性質(zhì)的稀疏矩陣，提出和實現(xiàn)了一種高效的矩陣向量乘存儲格式——bDIA格式.這種方法適宜于每行至多有K個非零值元素的M×N矩陣的情況.Vazquez等[22]發(fā)現(xiàn)ELLPACK格式更適合存儲.劉杰等[23]在求解稀疏線性方程組時，研究了共軛梯度法在GPU上的并行實現(xiàn)，發(fā)現(xiàn)該方法本身不適宜直接放到GPU上計算，需要進行優(yōu)化和改進，最終提出了高效求解的方法并加以實現(xiàn).白洪濤等[24- 25]結(jié)合矩陣格式探討了如何在GPU上優(yōu)化稀疏矩陣向量乘.吳洋和王迎瑞等[26- 27]研究了基于GPU的針對巨型稀疏矩陣的SpMV并行算法.在稀疏矩陣方面，李佳佳等[28]就如何針對不同稀疏矩陣選擇存儲格式展開了研究.在目前所研究的格式中，CSR格式[29]是一種流行的、具有一般用途的稀疏矩陣范例，適合于矩陣本身沒有明顯規(guī)律的情形.Li等[30]研究了新的矩陣格式——JAD格式，這種格式將矩陣行進行變換，使之類似于對角線格式，結(jié)果顯示，采用不完全上下三角矩陣(LU)分解的廣義最小殘量法求解線性方程組，在NVIDIA TESLA C1060上可達到近4倍的加速.

近幾年來，Xeon Phi[31]日益普及.它將多個計算核心整合在一起，形成超多核處理器[32]，使計算性能大幅度提升.相比SpMV和線性方程組在GPU上的優(yōu)化，對已有CPU程序不需要做大的改動，就可利用因特爾眾核架構(gòu)(MIC)的計算資源進行優(yōu)化，這樣就不會造成對原有軟件資源的浪費[33].目前關(guān)于Xeon Phi的研究主要集中在Xeon Phi的架構(gòu)和底層研究上，有關(guān)共軛梯度法和SpMV在Xeon Phi上應(yīng)用的研究還很少.Liu等[34]在Xeon Phi上做了有關(guān)SpMV并行算法的研究，提出了一種新的稀疏矩陣格式ESB(ELLPACK Sparse Block，這種格式是在ELLPACK格式的基礎(chǔ)上進行改進，適用于結(jié)構(gòu)規(guī)律的矩陣)，并使該矩陣適用于Xeon Phi的架構(gòu)，達到負載均衡，算例顯示SpMV在Xeon Phi Coprocessor Codenamed Knights上相較于在Intel Xeon Processor E5-2680和CUDA(NVIDIA Tesla K20X)上分別達到了3.52和1.32的加速比[35].

文中根據(jù)應(yīng)用問題的非結(jié)構(gòu)網(wǎng)格所產(chǎn)生的稀疏矩陣的特點，針對現(xiàn)代異構(gòu)模式多層次的高性能計算機體系架構(gòu)，探討新的稀疏矩陣存儲格式和高效的稀疏矩陣向量乘并行算法，綜合采用分塊技術(shù)、矩陣重組、預(yù)先提取數(shù)據(jù)等針對CPU處理器的優(yōu)化手段，同時，對GPU處理器上的計算采用線程與矩陣映射、數(shù)據(jù)合并訪問、數(shù)據(jù)復(fù)用等優(yōu)化手段，有效利用GPU多層次的存儲器體系，并通過高效的線程調(diào)度來隱藏全局存儲器的高延遲訪問；對Xeon Phi處理器上的計算，有效利用Xeon Phi適于解決邏輯比較復(fù)雜且并行度高的問題的優(yōu)勢.文中主要從以下幾個方面對Xeon Phi上的并行程序進行優(yōu)化：數(shù)據(jù)通信/傳遞、減少數(shù)據(jù)依賴、向量化、異步計算等；在此基礎(chǔ)上，使用其高效的線程調(diào)度來隱藏全局存儲器的高延遲訪問；最終通過對比GPU、OpenMP和Xeon Phi這些不同方式下的計算效率，結(jié)合文中的應(yīng)用背景，提出適合求解SpMV和線性方程組的并行算法.

1 問題描述和算法改進

文中所研究的應(yīng)用問題的實驗數(shù)據(jù)集源自城市大氣污染項目[5]，采集到的數(shù)據(jù)構(gòu)成的稀疏矩陣規(guī)模約為310 000×310 000，數(shù)據(jù)類型為雙精度浮點型.由于數(shù)據(jù)量大，若一次性傳輸這些數(shù)據(jù)，所耗費的時間和空間對于實驗分析而言過于龐大，且在傳輸過程還會出現(xiàn)數(shù)據(jù)遺漏的情況.為了實驗的準(zhǔn)確性和高效性，文中采用基于消息傳遞接口(MPI)的網(wǎng)格預(yù)處理方法對網(wǎng)格方程組進行區(qū)域分裂并行計算[1].在對網(wǎng)格方程組進行區(qū)域分裂后，將巨型稀疏矩陣劃分為若干進程下較小的稀疏矩陣，在此基礎(chǔ)上，文中需要解決快速求解大型線性方程組的問題.共扼梯度法是求解大型稀疏線性方程組的有效方法之一[7]，文中采用共軛梯度法[8]并通過求解泊松方程組來驗證所提方法的有效性.

算法描述如下：

Compute r(0)=b-Ax(0)for some initial guess x(0)

for i=1,2,…

solve Mz(i-1)=r(i-1)

ρi-1=r(i-1)Tz(i-1)

if i=1

p(1)=z(0)

else

p(i)=z(i-1)+βi-1p(i-1)

end if

q(i)=Ap(i)

x(i)=x(i-1)+αip(i)

r(i)=r(i-1)-αiq(i)

check convergence;continue if necessary

end

結(jié)合文中的應(yīng)用背景，優(yōu)化后的共軛梯度法在GPU上的求解流程如下.

步驟1 CPU端讀取矩陣文件，GPU從CPU端獲取稀疏矩陣數(shù)據(jù)和向量文件.

步驟2 在CPU端將不完整CSR格式存儲的數(shù)據(jù)提取出來，轉(zhuǎn)換為完整CSR格式，從而解除計算中的原子操作.

步驟3 計算殘差r(i)=b-x(i)(i=0,1,2,…).

步驟4 交換不同節(jié)點邊界的殘差.

步驟5 得到雅克比系數(shù)diagonal_diff=M-1.

步驟6 在CPU運算平臺求出迭代矩陣.

步驟7 計算SpMV：q=A×p，其中q表示稀疏矩陣和向量乘積的結(jié)果.

步驟8 計算全局誤差：Error=sum(rr×rr)，其中rr為標(biāo)準(zhǔn)差.

步驟9 求解x的值，檢驗是否符合精度要求或者是否已經(jīng)達到最大迭代次數(shù)，從而確定迭代是否終止.若條件成立，x即為方程的解向量，否則轉(zhuǎn)步驟3.

步驟1讀取矩陣數(shù)據(jù)和格式轉(zhuǎn)化，運算一次即可，故在CPU端執(zhí)行該步驟運算.步驟5和6涉及多次迭代，也是整個迭代運算的核心，需要在GPU上進行優(yōu)化.步驟7中的矩陣向量乘是整個方程求解中運算量最大的部分，該步驟需要在GPU上重點優(yōu)化.

最后將獲得的解向量由GPU傳送至主機內(nèi)存中.

上述步驟共涉及稀疏矩陣與向量相乘2次、對角矩陣與向量相乘1次、向量與向量相乘2次、向量更新3次、標(biāo)量相除2次，其中矩陣與向量相乘、向量與向量相乘、向量更新都可以利用數(shù)據(jù)級的并行操作，從而可以實現(xiàn)在CUDA和XeonPhi上的并行計算[7,35].上述流程中，SpMV的計算時間占了總計算時間的80%以上[12].SpMV計算性能與所計算的矩陣有直接的關(guān)系，矩陣的規(guī)模、分布情況、稀疏程度以及非零元素個數(shù)等都會影響加速效率.文中所研究項目來源于針對城市空氣污染模型的大渦模擬，其通過網(wǎng)格并行劃分后產(chǎn)生的稀疏矩陣的非零元素值大致走勢的Matlab仿真圖如圖1所示.可見，稀疏矩陣呈對稱正定線性分布，而非分塊或條形對角線分布，除了對稱外，不具有某種特殊格式.文中采用更加靈活且易于操作的CSR格式，并在此格式基礎(chǔ)上進行改進和優(yōu)化.

圖1 稀疏矩陣非零元素分布圖Fig.1 Distribution diagram of non-zero elements of sparse matrix

CSR[29]用3個相關(guān)向量來表示稀疏矩陣，其中向量A用于存儲稀疏矩陣的非零元素，向量Col_A用于標(biāo)識非零元素對應(yīng)的列值，向量Row_A用于存儲每一行第1個非零元素的偏移量，即Row_A[i]表示第0行到第i-1行的非零元素個數(shù)和，初始Row_A[0]賦值為0.對于一個M×N的稀疏矩陣，Row_A向量的長度為M+1，其作用是把每一行的偏移量存到Row_A[i]中，其最后一個元素存儲矩陣中非零元素的個數(shù).圖2示出了如何用CSR格式中的3個一維數(shù)組表示稀疏矩陣的過程.

圖2 CSR范例

1.1 基于GPU的SpMV的Kernel函數(shù)

計算機GPU體系架構(gòu)在不斷發(fā)展，文中以NVIDIA Tesla K20[16]體系架構(gòu)為代表對GPU的并行層次進行分析.NVIDIA Tesla K20由可伸縮流處理器陣列(SP)和存儲器系統(tǒng)組成，并由一個片上的互聯(lián)網(wǎng)絡(luò)連接.GPU利用大量功能簡單的流多處理器協(xié)同運算來提高數(shù)據(jù)的吞吐率和計算速度.與之不同的是，CPU則采用相對復(fù)雜的分支預(yù)測和邏輯控制來迅速獲得數(shù)據(jù)和指令，從而減少了延遲帶來的耗時.GPU在計算時，是以線程(Thread)為單位，每個線程有本地內(nèi)存和多個寄存器(Register).其一個線程塊(Block)由若干個線程組成，各線程塊中的線程可以通過共享內(nèi)存進行通信，多個線程塊組成了一個線程網(wǎng)格(Grid).針對GPU的體系架構(gòu)可以分析出：基于GPU平臺的稀疏矩陣向量乘方法存在的挑戰(zhàn)主要包括線程映射和存儲層次的管理、數(shù)據(jù)訪問和數(shù)據(jù)復(fù)用、非零元素分布不均.

常見的基于CSR格式、利用CUDA加速的SpMV有3種實現(xiàn)方法：①使用一個線程來處理一行，一個warp包含32個線程，那么一次可以處理32行,但因每行非零元素個數(shù)的不同，線程空轉(zhuǎn)的問題會非常嚴重，同時，在線程空轉(zhuǎn)情況下，無法聯(lián)合訪問GPU的顯存，降低了訪存效率；②每個線程塊計算結(jié)果向量的一個元素，同樣也會存在每個線程塊處理非零元素不均而造成線程空等的情況；③若使用一個warp來計算矩陣的一行與向量乘，雖然在一定程度上減輕了線程空轉(zhuǎn)的問題，并且warp內(nèi)可聯(lián)合訪問，使得運行效率有所提高，但是空轉(zhuǎn)問題依然比較嚴重，訪問效率仍然存在提高的空間.用CSR進行存儲時，基于GPU平臺下的SpMV還存在著以下問題：第一，線程間等待的問題可能是由于同一個warp內(nèi)的線程間的計算量負載不均衡而引起的；第二，線程不滿足對全局存儲器的合并訪問導(dǎo)致訪存延遲.這兩個問題都可能影響到稀疏矩陣方程組的求解性能.結(jié)合文中稀疏矩陣的特點，每一行非零元素值不超過32個，故選擇第3種方式并在此基礎(chǔ)上進行優(yōu)化.

由于文獻[5- 6]中初始程序里的計算方式不是完整CSR格式，而是按照對角線和上三角矩陣來進行計算，若是分別進行存儲，再放到GPU上計算，會產(chǎn)生原子操作，從而影響計算速度，所以開始需要掃描一遍稀疏矩陣，生成完整CSR格式，解除CUDA上的原子操作(1).偽代碼如下：

for i<-THREAD_ID to nde ∥nde表示矩陣維度

for k<-nnz[i] to nnz[i+1] ∥數(shù)組nnz[]存儲每行非零元素個數(shù)

tmp_kk<-colume[k] ∥數(shù)組colume[]存儲非零元素的列下表值

dbl_tmp<-delta_t*diff[k]∥dbl_tmp代表矩陣A的非零元素值

qq_i<-qq_i+dbl_tmp*pp[tmp_kk]∥A*p

atomicAdd((double*)&(temp_qq[tmp_kk])，double(dbl_tmp * pp[i]))∥原子操作(1)，避免同一線程塊內(nèi)不同線程對同一元素進行加操作以及同一線程網(wǎng)格內(nèi)不同線程塊內(nèi)不同線程對同一元素進行加操作

temp_qq[tmp_kk]<-temp_qq[tmp_kk] + dbl_tmp*pp[i]∥數(shù)組temp_qq存儲A*p的值

end

針對GPU架構(gòu)上SpMV的計算資源沒有充分利用以及負載不平衡等問題，文中對CSR格式進行改良，即采用非零元素個數(shù)升序排序，把非零元素個數(shù)相似的矩陣行盡可能多地放在一個warp或者多個warp中進行處理，以避免計算資源浪費和解決負載不平衡的問題，進而提高稀疏矩陣向量乘的計算性能.

基于CUDA的稀疏矩陣向量乘的具體優(yōu)化算法流程如下，其中“2”、“3”-“7”采用異步操作“cudaEventCreate(&start);cudaEventCreate(&stop);”，使申請GPU變量空間和計算誤差等步驟同步進行，這樣可以隱藏申請GPU變量空間的時間，具體的算法步驟如下：

1 cudaEventCreate(&start)

2 malloc variable space，initial the whole CSR array col_A、row_A、val_A，transfer array col_A、row_A、val_A、q to GPU

3 poisson_solver_residual()∥calculate the residual error

4 exchange the Boundary elements()

5 From initial value x，calculate r=b-Ax

6 p=M-1r and transfer the value p to GPU

7 ρ=rTp

8 cudaEventCreate(&stop)

9 for i=1，2，…{

10 GPU kernel

11 {q=Ap}

12 transfer p and q to CPU

13 CPU {

14 exchange the elements of each boundary();

15 α_tmp=pTq

16 α=ρ/α

17 ρold=ρ

18 r=r-αq

19 e_tmp=rTr

20 q=M-1r

21 ρ_tmp=rTq

22 β=ρ/ρold

23 }

24 exchange the error of each boundary();

25 transfer p and q to GPU

26 GPU kernel {

27 p=q+βp

28 q=val_diag×p

29 x=x+αp

30 }

31 transfer p，q to CPU

32 Check the convergence，end if satisfied

33 end of for

34 end

根據(jù)以上描述，文中的算法流程可以用圖3來描述.

圖3 文中優(yōu)化算法的流程

1.2 基于GPU的SpMV的Half-Warp優(yōu)化

采用線程負責(zé)計算結(jié)果向量的一個元素時，因CUDA的體系架構(gòu)，處理器是將線程塊劃分為多個warp，并且以warp為執(zhí)行單元，同一warp內(nèi)的線程是天然同步的，所以當(dāng)一個warp內(nèi)的線程對應(yīng)處理行的非零元素個數(shù)相差較大時，每個線程的進度會有很大的差異，這將造成許多線程計算資源的空轉(zhuǎn)等待.與DIA和ELL格式不同，CSR存儲格式每行擁有可變數(shù)量的非零元素個數(shù)，而這種差異導(dǎo)致了線程間的差異.例如，當(dāng)稀疏矩陣每行包含高度可變數(shù)目的非零元素時，可能會出現(xiàn)的情況是：在包含非零元素個數(shù)最多的行的線程繼續(xù)重復(fù)執(zhí)行操作的同時，同一warp中的其他線程卻保持空閑狀態(tài).另外，線程的同步會產(chǎn)生額外的開銷.若每個線程塊計算結(jié)果向量的一個元素，當(dāng)稀疏矩陣某一行的非零元素個數(shù)不能被block內(nèi)的線程總數(shù)整除時，也會出現(xiàn)部分線程等待的情況.而矢量內(nèi)核的有效執(zhí)行需要矩陣的行包含的非零元素個數(shù)大于warp的大小(一般為32).所以，當(dāng)所計算行的非零元素個數(shù)小于GPU的執(zhí)行單元warp的線程數(shù)32時，利用warp計算的策略仍然會存在空轉(zhuǎn)等待的情況.因此，矢量內(nèi)核的性能與矩陣每一行的非零元素值的個數(shù)是緊密相關(guān)的.文中研究對象所產(chǎn)生的稀疏矩陣各行中最大的非零元素個數(shù)不超過32，且超過70%的矩陣行的非零元素個數(shù)不超過16，進一步地結(jié)合GPU的全局存儲器合并訪問的需要，可知文中適合采用Half-Warp的策略.

以Half-Warp為單位循環(huán)計算的優(yōu)化策略即每個Half-Warp分別計算稀疏矩陣中某一行元素與向量的乘積.首先，根據(jù)block和grid線程個數(shù)的分配以及矩陣行非零元素的個數(shù)，動態(tài)地將數(shù)據(jù)以Half-Warp為單位進行分割，每個單位計算稀疏矩陣的一行元素值與向量的乘積，然后對Half-Warp的中間結(jié)果進行reduction歸約求和，由于中間結(jié)果存在數(shù)據(jù)復(fù)用，為了有效地提高寄存器的訪問速度，文中使用了共享存儲器SharedMemory來存儲.文中SpMV的Kernel函數(shù)的程序偽代碼如下(優(yōu)化編譯選項為-arch sm_35)：

t<-threadIdx.x∥Each block’s Thread ID

id<-t &(warpSize-1)∥Thread ID within warp

warpsPerBlock<-blockDim.x/warpSize

myRow<-(blockIdx.x*warpsPerBlock)+(t/warpSize)∥One row per warp

Create shared memory array val[NUM_THREADS]

∥以上為賦初始值

if myRow

warpStart<-nnz[myRow]

warpEnd<-nnz[myRow+1]

qq_i<-0.0

for jj<-warpStart+id to warpEnd

tmp_kk<-colume[jj]

qq_i<-value[jj]*pp[tmp_kk]+qq_i

jj<-warpSize+jj

val[t]<-qq_i

∥Reduce the partial sums

if id

val[t]<-val[t+mun]+val[t]

end for

∥Write the result

if id==0

temp_qq[myRow]<-val[t]+temp_qq[myRow]

end

1.3 數(shù)據(jù)的傳輸

在進行GPU通用計算時，有時會出現(xiàn)數(shù)據(jù)從CPU端傳到GPU端的時間比kernel的運行時間要長的情況，而數(shù)據(jù)的傳輸又是基于GPU通用計算之外的開銷，因此會影響整體的運行速度，進而降低加速比.文中算法采用了分頁鎖定(Page-Locked)內(nèi)存空間的方式，用cudaMallocHost()申請一塊CPU內(nèi)存，再將數(shù)據(jù)拷貝到顯存,這樣將會減少一定的數(shù)據(jù)傳輸時間，從而提高計算速度.但是當(dāng)數(shù)據(jù)規(guī)模非常大時，在CPU上鎖定太多的內(nèi)存同樣會影響kernel的運行時間，所以文中僅用cudaMallocHost()申請需要在CPU和GPU之間頻繁傳輸?shù)闹虚g計算結(jié)果pp和qq.

1.4 其他優(yōu)化

其他優(yōu)化如下：①兩向量點乘的操作和向量更新通過自編kernel函數(shù)，采用并行規(guī)約求和方法，比直接調(diào)用CUBLAS庫中的CublasDdot實現(xiàn)點乘運算和調(diào)用CublasDaxpy更新向量要快1.6倍；②為了加速數(shù)據(jù)讀取、減少讀取延遲，文中用紋理寄存器綁定CSR格式中的val值.

1.5 Xeon Phi優(yōu)化

1.5.1 并行優(yōu)化

首先，為了減少開銷，文中只進行一次創(chuàng)建和結(jié)束OpenMP并行區(qū)域的操作.其次，將r、e_tmp、q和ρ_tmρ放在同一個OpenMP循環(huán)內(nèi)，并且把p、q和x放在另外一個循環(huán)內(nèi)，以便使接下來的計算可以從緩存中直接讀取數(shù)據(jù)，文中把r和向量在for循環(huán)計算開始前讀入高速緩存中.例如，q=M-1r可以從緩存中預(yù)先讀取下一次計算所需的r值，直到最后一次循環(huán)結(jié)束.最后，用變量q替換變量z并進行分組計算，這樣可以省去變量z在算法中所帶來的開銷.優(yōu)化流程如下(所使用的優(yōu)化選項為-O3)：

1 #PARAGMA OMP PARALLEL {

2 for i=1，2，…

3 # PARAGMA OMP PARALLEL FOR

4 q=Ap

5 #pragma omp single {exchange boundary elements}

6 # PARAGMA OMP FOR REDUCTION (+:α_tmp)

7 α_tmp=pTq

8 # PARAGMA OMP SINGLE{

9 α=ρ/α

10 ρold=ρ

11 }

12 # PARAGMA OMP FOR REDUCTION(+:ρ_tmp，error_tmp){

13 r=r-αq

14 e_tmp=rTr

15 q=M-1r

16 ρ_tmp=rTq

17 }

18 # PARAGMA OMP MASTER {

19

20 β=ρ/ρold

21 }

22 # PARAGMA OMP SINGLE {

23 exchange the error of boundary

24 }

25 # PARAGMA OMP FOR NOWAIT{

26 p=q+βp

27 q=val_diag×p

28 x=x+αp

29 }

30 Check the convergence，end if satisfied

31 end of for

32 end of #PRAGMA

1.5.2 減少數(shù)據(jù)依賴

將上三角格式分別計算、轉(zhuǎn)換為完整CSR格式來計算SpMV，可去掉內(nèi)層循環(huán)的數(shù)據(jù)依賴關(guān)系，即去掉了多線程下的原子操作.對于Xeon Phi平臺，由于要開啟上百個線程，減少線程開啟對平臺性能的影響更為重要.在內(nèi)層計算之前，與上述GPU優(yōu)化一樣，先轉(zhuǎn)換為完整CSR格式并進行排序，在排序的同時將內(nèi)層循環(huán)需要隨機調(diào)度的變量按順序存放，以便進行向量化優(yōu)化，這樣可以隱藏線程間的交互時間.

1.5.3 向量化優(yōu)化

計算誤差和更新向量結(jié)果值的操作在調(diào)整程序循環(huán)結(jié)構(gòu)(比如循環(huán)嵌套和交換次序等)并使之沒有數(shù)據(jù)依賴關(guān)系后，使用自動化向量化優(yōu)化循環(huán).

1.5.4 計算異步

計算進程從主進程獲取初始化后的數(shù)據(jù)以及稀疏矩陣所要乘的向量的值.申請Xeon Phi變量空間和計算誤差、轉(zhuǎn)換數(shù)據(jù)格式等步驟同步進行，這樣可以隱藏在Xeon Phi上申請變量的時間.最后釋放Xeon Phi上已申請的空間.為了隱藏數(shù)據(jù)傳輸時間，文中以流水線方式執(zhí)行傳輸數(shù)據(jù)和計算，具體如圖4所示.

圖4 異步執(zhí)行示意圖

2 實驗分析與總結(jié)

文中測試平臺如下：CPU為Intel(R) Core(TM) i5-3470 CPU 3.20 GHz；GPU為Nvidia Tesla K20M；Xeon Phi為Intel(R) Xeon(R) CPU E5-2650 0 @2.00 GHz，7 110P，8 GB GDDR5，61 cores；操作系統(tǒng)為64位Linux Red Hat，1個CPU，1個GPU，1個MIC.所測試的數(shù)據(jù)來源于大渦流模擬有限元計算產(chǎn)生的稀疏矩陣[5- 6]，稀疏矩陣用WRF_PCG表示，選取主進程所產(chǎn)生的稀疏矩陣數(shù)據(jù)作為代表進行測試，它的主要特征參數(shù)見表1.

表1 參數(shù)設(shè)置

首先，在CPU+GPU平臺單進程下測試SpMV的運行時間，測試結(jié)果見圖5.

從圖5可看出：格式優(yōu)化后WRF_PCG在GPU上的運算速度大幅度提升，耗時從串行程序的0.879 s加速至0.214 s，說明將不完整CSR格式轉(zhuǎn)換為完整CSR格式后，所解除的原子操作以及優(yōu)化數(shù)組排序的操作使得WRF_PCG更適合在GPU上進行計算；Half-Warp及數(shù)據(jù)異步優(yōu)化、數(shù)據(jù)傳輸?shù)绕渌麅?yōu)化所帶來的加速效果分別對應(yīng)于0.148、0.140和0.111s的耗時.以Half-Warp為單位的計算適合文中稀疏矩陣各行非零元素個數(shù)不超過16的情況，在提高寄存器訪問速度的同時降低了線程空轉(zhuǎn)的可能；數(shù)據(jù)異步優(yōu)化隱藏了在GPU上開辟空間傳輸數(shù)據(jù)的時間，數(shù)據(jù)傳輸優(yōu)化和其他優(yōu)化增強了數(shù)據(jù)復(fù)用；從耗時加速情況來看，文中各種策略是有效的.

圖5 WRF_PCG在CPU+GPU下優(yōu)化后的SpMV運行時間

Fig.5 SpMV run-time of WRF_PCG after CPU+GPU optimization

然后，在CPU下用OpenMP優(yōu)化后測試SpMV的運行時間，結(jié)果如圖6所示.

圖6 WRF_PCG在OpenMP優(yōu)化后的SpMV運行時間

從圖6可看出，在1、2、4個線程下，OpenMP優(yōu)化后的加速比基本是隨著線程數(shù)的增多而接近線性增長，由于測試環(huán)境的CPU是4核，每個核2個線程，所以在8個線程開啟時，運行加速效果最佳.

最后，在CPU+Xeon Phi下用Xeon Phi優(yōu)化后測試SpMV的運行時間，結(jié)果如圖7所示.

從圖7可看出，在其他條件相同的情況下，1個線程時Xeon Phi的優(yōu)化效果并不理想，比CPU單線程下運行慢很多.進一步測試發(fā)現(xiàn)，在Xeon Phi上開辟變量空間和數(shù)據(jù)傳輸時間對運行速度影響更大.當(dāng)線程數(shù)達到Xeon Phi cores的兩倍時，運算時間最短，達0.259 s，與初始程序相比，加速比達3.421，加速效果較為明顯.由于啟動Xeon Phi卡的耗時過長和線程切換耗時，在Xeon Phi下的優(yōu)化效果不如NVIDIA Tesla K20.

圖7 WRF_PCG在CPU+Xeon Phi優(yōu)化后的SpMV運行時間

Fig.7 SpMV run-time of WRF_PCG after CPU+Xeon Phi optimization

在其他條件相同，1個GPU卡、1個MIC卡和2個CPU核四線程開啟的情況下，WRF_PCG在各優(yōu)化方式下的SpMV和PCG的加速比情況如圖8、9所示.可以看出，對WRF_PCG而言，SpMV和PCG加速效果最好的是CPU+CUDA方式，其次是CPU+Xeon Phi，最后是OpenMP.在百萬網(wǎng)格產(chǎn)生的數(shù)據(jù)量下，CPU+CUDA異構(gòu)模式在高性能計算機上用共軛梯度求解稀疏線程方程比僅在CPU下求解加速3.2倍，而在整個計算最耗時的計算SpMV部分，加速比約達7.9.隨著數(shù)據(jù)量的增大，CUDA的計算速度會得到更大程度的發(fā)揮，所以該算法優(yōu)化程序的擴展性良好.

圖8 WRF_PCG在各優(yōu)化方式下的SpMV加速比

Fig.8 SpMV speedup ratio of WRF_PCG after each optimi-zation

為了證明文中所選存儲格式是最適合的，分別在DIA、COO、ELL、HYB和文中格式下測試文中數(shù)據(jù)集(相應(yīng)參數(shù)列于表1).

8個進程下稀疏矩陣的非零元素值的分布見圖2，測試結(jié)果如圖10所示.

圖9 WRF_PCG在各優(yōu)化方式下的PCG加速比

Fig.9 PCG speedup ratio of WRF_PCG after each optimization

圖10 不同存儲格式下運行時間的比較

Fig.10 Comparison of run-time under different storage formats

從實驗結(jié)果來看，不同存儲格式的計算性能有很大的差異，文中所討論的應(yīng)用問題所產(chǎn)生的稀疏矩陣的非零元素分布是對稱的，但不存在與對角線平行的直線，且有很多離散的點分布，這種分布與對角線格式的要求差異最大，在用對角線格式存儲時會填充很多零值，既占據(jù)大量存儲空間又不能給計算帶來便利，故該方法的優(yōu)化效果最差，較串行算法僅加速10%.COO格式具有普遍應(yīng)用性，但行和列指針都被顯式地存儲，開辟的存儲空間較其他格式要多，所以其優(yōu)化效果并不是最佳，較串行算法其加速比約為2；當(dāng)矩陣每行的非零元素值的最大數(shù)目與平均數(shù)目相差不大時，ELL格式是比較好的選擇，這種格式更適合于規(guī)則的矩陣.文中的矩陣形狀并不是由半結(jié)構(gòu)化和結(jié)構(gòu)化的網(wǎng)格得到的矩陣抽象圖，也不是最適合的格式，較串行算法的加速比約為2.6.HYB格式是將ELL和COO兩種格式進行混合，這種格式適合于當(dāng)稀疏矩陣每行的非零元素值的數(shù)目變化特別大時，相比兩種單一的格式，其加速效果明顯，但文中矩陣的每行非零元素的個數(shù)變化不大，所以加速效果沒有預(yù)期顯著，加速比約為3.文中的優(yōu)化算法是針對應(yīng)用問題產(chǎn)生的矩陣進行的優(yōu)化，從GPU的多層次存儲器體系、線程與矩陣映射、數(shù)據(jù)合并訪問、數(shù)據(jù)復(fù)用等方面逐一進行分析和優(yōu)化，并最后達到比較理想的加速效果.雖然不同進程下產(chǎn)生的稀疏矩陣的數(shù)據(jù)不同，但是其非零元素分布圖大同小異，實驗結(jié)果也進一步證明，在分布情況如圖2所示，且矩陣本身沒有特別的格式特征，用其他格式也達不到比較好的加速效果的情況下，采用文中的優(yōu)化算法最佳，其相較于串行算法的加速比約為8.

為了驗證文中所提算法的拓展性和有效性，選取與文獻[18]相同的10個測試矩陣來測試稀疏矩陣向量乘的性能，測試集來源于佛羅里達稀疏矩陣集[36]，該測試集中的數(shù)據(jù)可以真實地反應(yīng)SpMV在真實應(yīng)用中的性能.所選取的10個稀疏矩陣的結(jié)構(gòu)信息如表2所示(按照行列數(shù)排序).

表2 稀疏矩陣參數(shù)

文獻[18]和[28]與文中研究內(nèi)容相似，因此，在相同測試環(huán)境和測試集下，分別根據(jù)文獻[18]和[28]所提供的算法并將參數(shù)調(diào)試到其最優(yōu)情況，測試表2中稀疏矩陣的加速比，結(jié)果如圖11所示.

圖11 3種SpMV算法對10種稀疏矩陣的加速比

Fig.11 Speedup ratios of three SpMV algorithms on ten sparse matrices

從圖11可看出:對于3種基于CSR格式的算法而言，若矩陣的結(jié)構(gòu)由單一分塊組成(如矩陣1、2、4、5)，SpMV的加速比會更大;若矩陣的結(jié)構(gòu)不適宜分塊(如矩陣3和7)，則SpMV的加速比相對較小.文獻[28]和文中算法都采用異步計算的方式，隱藏了轉(zhuǎn)換和優(yōu)化格式的時間，而文獻[18]中算法沒有隱藏，其轉(zhuǎn)換BCSR格式的時間使得加速比不如其他兩種算法.文中算法在選取參數(shù)和應(yīng)用數(shù)據(jù)上都是以實際應(yīng)用為背景，故在應(yīng)用性上較文獻[28]中算法更具有普遍性和有效性，在測試10個不同矩陣時，文中算法從加速比結(jié)果來看是最優(yōu)的.

對文中提出的優(yōu)化方法在公共稀疏矩陣方程集上進行測試，選用的是文獻[36]中的11種不同規(guī)模的矩陣.這些矩陣選取自科學(xué)和工程計算中常用的矩陣數(shù)據(jù)，規(guī)模由千量級至百萬量級，跨度很大.測試矩陣數(shù)據(jù)來源于流體力學(xué)、熱能計算、空氣動力學(xué)、計算機圖形學(xué)、生物和環(huán)境科學(xué)等領(lǐng)域，具有不同的非零元素個數(shù)和矩陣特征.表3列出了測試矩陣的具體參數(shù).

表3 測試矩陣參數(shù)

在相同測試環(huán)境和測試數(shù)據(jù)下，根據(jù)文獻[15]所提供的算法并將參數(shù)調(diào)試到其最優(yōu)情況，測試表3中的稀疏矩陣，文獻[15]和文中優(yōu)化算法在CPU情況下的加速比對比如圖12所示.

從圖12可看出，兩種算法的加速比均大于10，最高達70以上，但從整體來看，文中所提算法的加速比優(yōu)于文獻[15]算法，且矩陣規(guī)模越大加速效果越明顯.當(dāng)矩陣規(guī)模到達1 585 478×1 585 478時，文中算法的加速比約為文獻[15]算法的1.4倍，加速效果顯著.另外，文中算法在數(shù)據(jù)處理時采用分塊的方法，且對存儲器的優(yōu)化適用于實際應(yīng)用數(shù)據(jù)，故在應(yīng)用性上較文獻[15]算法更具普遍性和有效性.

圖12 文中算法和文獻[15]算法對測試矩陣的加速比

Fig.12 Speedup ratios of the proposed algorithm and the algorithm proposed in Literature [15] on test matrices

為驗證文中提出的基于Xeon Phi的優(yōu)化算法在其他應(yīng)用上的可拓展性，測試了該算法在表3中所給稀疏矩陣上的加速比，并與文獻[34]中提出的基于Xeon Phi上的稀疏線性系統(tǒng)的優(yōu)化算法進行比較，結(jié)果見圖13.

圖13 文中算法和文獻[34]算法在公共集方程上的加速比

Fig.13 Speedup ratios of the proposed algorithm and the algorithm proposed in Literature [34] on the public equation matrix

從測試結(jié)果來看，在Xeon Phi上的整體加速比較在GPU上的略低，這是因為SpMV本身計算的邏輯性并不復(fù)雜，而Xeon Phi中開啟線程和線程間的開銷比計算所需時間要多.文中算法整體加速比均大于10，最優(yōu)情況達到68.14，總體而言，文中算法優(yōu)于文獻[34]算法，同時隨著矩陣的規(guī)模增大，文中算法的加速效果愈發(fā)明顯.當(dāng)矩陣規(guī)模為4 307×4 307時，文中算法較文獻[34]算法的加速比增加約6.4%；當(dāng)矩陣規(guī)模達到1 585 478×1 585 478時，文中算法加速比約為文獻[34]算法的1.5倍，加速效果顯著.在所測試的11個不同矩陣上，文中算法均獲得了更好的加速效果.

為了說明文中算法的可擴展性，在廣州超級計算中心的天河二號超級計算機上，根據(jù)應(yīng)用問題和超級計算機硬件的特性設(shè)計計算參數(shù)，用現(xiàn)代并行理論評測工具測評并行程序，找出問題瓶頸，實現(xiàn)程序的大規(guī)模測試優(yōu)化.不同節(jié)點數(shù)時文中算法在GPU上的運行時間和加速比如圖14所示.可見，在2 064核上仍有加速，加速比近20，表明文中算法具有良好的可擴展性.

圖14 文中算法在天河二號超級計算機上的超多核測試結(jié)果

Fig.14 Test results of the proposed algorithm on the Tianhe-2 supercomputer under many-core condition

3 結(jié)語

文中針對并行共軛梯度應(yīng)用問題中的存儲器瓶頸類算法SpMV進行研究，將問題的特征與GPU的特殊體系結(jié)合在一起，提出了CSR格式下的優(yōu)化策略，實驗結(jié)果驗證了文中方法的有效性；同時，進一步分析OpenMP版本的性能結(jié)果，證實了在Xeon Phi上優(yōu)化的可能性.由于OpenMP版本的性能結(jié)果良好，移植到Xeon Phi上并進行各種優(yōu)化后，已達到預(yù)期的優(yōu)化效果.對比CUDA和Xeon Phi的實驗結(jié)果不難發(fā)現(xiàn)，對于隨機存儲性大且并行度不高、邏輯簡單的應(yīng)用而言，CUDA比Xeon Phi更具優(yōu)勢.目前，關(guān)于稀疏矩陣向量乘的優(yōu)化在CUDA上已有一些研究成果，但在Xeon Phi上的算法優(yōu)化還處于初期探索階段，因此，文中所得實驗數(shù)據(jù)和研究成果對科研和工程實踐都將有所裨益.

目前的工作只是在單GPU和單Xeon Phi上進行，未來擬將優(yōu)化算法移植到多GPU和MIC上，充分挖掘CPU與GPU、MIC的協(xié)同計算能力，通過實驗歸納出適于兩者異步計算的策略.另外，文中算法在MIC上的優(yōu)化效果并不如在GPU上顯著，其并行計算潛力還有待繼續(xù)挖掘，今后將進一步尋找、優(yōu)化適宜于在MIC上并行的算法.后續(xù)研究中，還可以針對邏輯復(fù)雜、并行度高的算法進行優(yōu)化，并與GPU上的優(yōu)化進行對比、總結(jié);同時，還可針對解決稀疏線性系統(tǒng)的其他算法進行優(yōu)化，形成解決這一類問題的計算策略、規(guī)則和模式.在Xeon Phi上求解文中問題的做法是來回通信，在后續(xù)研究中，將考慮將算法全部在Xeon Phi上求解，以達到更好的優(yōu)化效果.

[1] Mielikainen Jarno,Huang Bormin,Huang Allen Hung-Lung,et al.Improved GPU/CUDA based parallel weather and research forecast(WRF)simple moment 5-class(WSM5)cloud Xeon phirophysic [J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2012,5(4):1256- 1265.

[2] 方維,孫廣中,吳超,等.一種三維快速傅里葉變換并行算法 [J].計算機研究與發(fā)展,2011,48(3):440- 446. Fang Wei,Sun Guangzhong,Wu Chao,et al.A parallel algorithm of three-dimensional fast Fourier transform [J].Journal of Computer Research and Development,2011,48(3):440- 446.

[3] 蔡勇,李光耀,王琥.GPU通用計算平臺上中心差分格式顯式有限元并行計算 [J].計算機研究與發(fā)展,2013,50(2):412- 419. Cai Yong,Li Guangyao,Wang Hu.Parallel computing of central difference explicit finite element based on GPU general computing platform [J].Journal of Computer Research and Development,2013,50(2):412- 419.

[4] 唐亮,駱祖塋,趙國興,等.利于GPU計算具有線性并行度的P/G網(wǎng)SOR求解算法 [J].計算機研究與發(fā)展,2013,50(7):1491- 1500. Tang Liang,Luo Zuying,Zhao Guoxing,et al.SOR-based P/G solving algorithm of linear parallelism for GPU computing [J].Journal of Computer Research and Development,2013,50(7):1491- 1500.

[5] Cheng W C,Liu C H.Large-eddy simulation of flow and pollutant transports in and above two-dimensional idea-lized street canyons [J].Boundary-Layer Meteorol,2011,139(4):411- 437.

[6] Liu C H,Barth M C.Large-eddy simulation of flow and scalar transport in a modeled street canyon [J].Journal of Appl Meteor,2002,41(5):660- 673.

[7] Maringanti A,Athavale V,Patkar S.Acceleration of conjugate gradient method for circuit simulation using CUDA [C]∥Proceedings of HiPC 2009.Kochi:IEEE,2009:438- 444.

[8] Matam K K,Kothapalli K.Accelerating sparse matrix vector multiplication in iterative methods using GPU [C]∥Proceedings of 2011 International Conference on Parallel Processing.Taipei:ACM,2011:612- 621.

[9] Li Xiangge.Preconditioned conjugate gradient solver for structural problems [D].Columbia:Department of Computer Science,University of Missouri-Columbia,2013.

[10] Naumov M.Parallel solution of sparse triangular linear systems in the preconditioned iterative methods on the GPU,NVR-2011- 001 [R].Santa Clara:NVIDIA,2011.

[11] Reguly István,Giles Mike.Efficient sparse matrix-vector multiplication on cache-based GPUs [C]∥Proceedings of Innovative Parallel Computing.San Jose:IEEE,2012.

[12] Dziekonski A,Lamecki A,Mrozowski M.A memory efficient and fast sparse matrix vector product on a GPU [J].Progress in Electromagnetics Research,2011,116(1):49- 63.

[13] El Zein A,HRendell A P.Generating optimal CUDA sparse matrix-vector product implementations for evolving GPU hardware [J].Concurrency and Computation:Practice and Experience,2012,24(1):3- 13.

[14] Monakov Alexander,Lokhmotov Anton,Avetisyan Arutyun.Automatically tuning sparse matrix-vector multiplication for GPU architectures [J].Lecture Notes in Computer Science,2010,59(52):111- 125.

[15] 張建飛,沈德飛.基于GPU的稀疏線性系統(tǒng)的預(yù)條件共軛梯度法 [J].計算機應(yīng)用,2013,33(3):825- 829. Zhang Jianfie,Shen Defei.GPU-based preconditioned conjugate gradient method for solving sparse linear systems [J].Journal of Computer Applications,2013,33(3):825- 829.

[16] NVIDIA.Fermi compute architecture whitepaper [R].Santa Clara:NVIDIA,2009.

[17] NVIDIA.CUSPARSE library [Z].Santa Clara:NVIDIA,2011.

[18] 袁娥,張云泉,劉芳芳,等.SpMV的自動性能優(yōu)化實現(xiàn)技術(shù)及其應(yīng)用研究 [J].計算機研究與發(fā)展,2009,46(7):1117- 1126. Yuan E,Zhang Yunquan,Liu Fangfang,et al.Automatic performance tuning of sparse matrix-vector multiplication:implementation techniques and its application research [J].Journal of Computer Research and Development,2009,46(7):1117- 1126.

[19] Bell N,Garland M.Efficient sparse matrix-vector multiplication on CUDA [R].Santa Clara:NVIDIA,2008.

[20] Yuan Liang,Zhang Yunquan,Sun Xiangzheng,et al.Optimizing sparse matrix vector multiplication using diagonal storage matrix format [C]∥Proceedings of the 12th IEEE International Conference on High Performance Computing and Communications.Cidade de Goa:IEEE,2010:585- 590.

[21] 孫相征,張云泉,王婷,等.對角稀疏矩陣的SpMV自適應(yīng)性能優(yōu)化 [J].計算機研究與發(fā)展,2013,50(3):648- 656. Sun Xiangzheng,Zhang Yunquan,Wang Ting,et al.Optimizing auto-tuning of SpMV for diagonal sparse matrices [J].Journal of Computer Research and Development,2013,50(3):648- 656.

[22] Vazquez F,Ortega G,Fernandez J J,et al.Improving the performance of the sparse matrix vector product with GPUs [C]∥Proceedings of the 10th IEEE International Conference on Computer and Information Technology(CIT 2010).Bradford:IEEE,2010:1146- 1151.

[23] 劉杰,劉興平,遲利華,等.一種改進的適合并行計算的共軛剩余算法 [J].計算機學(xué)報,2006,29(3):495- 499. Liu Jie,Liu Xing-ping,Chi Li-hua,et al.An improved conjugate residual algorithm for large symmetric linear systems [J].Chinese Journal of Computer,2006,29(3):495- 499.

[24] 白洪濤,歐陽丹彤,李熙銘,等.基于GPU的稀疏矩陣向量乘優(yōu)化 [J].計算機科學(xué),2010,37(8):168- 181. Bai Hong-tao,Ouyang Dan-tong,Li Xi-ming,et al.Optimizing sparse matrix-vector multiplication based on GPU [J].Computer Science,2010,37(8):168- 181.

[25] 白洪濤.基于GPU的高性能并行算法研究 [D].長春:吉林大學(xué)計算機科學(xué)與技術(shù)學(xué)院,2010.

[26] 吳洋,趙永華,紀國良.一類大規(guī)模稀疏矩陣特征問題求解的并行算法 [J].數(shù)值計算與計算機應(yīng)用2013,34(2):137- 146. Wu Yang,Zhao Yonghua,Ji Guoliang.Parallel solving large-scale sparse matrix eigenvalue problem [J].Journal on Numerical Methods and Computer Applications,2013,34(2):137- 146.

[27] 王迎瑞,任江勇,田榮.基于GPU的高性能稀疏矩陣向量乘及CG求解器優(yōu)化 [J].計算機科學(xué),2013,40(3):46- 49. Wang Ying-rui,Ren Jiang-yong,Tian Rong.Efficient sparse matrix-vector multiplication and CG solver optimization on GPU [J].Computer Science,2013,40(3):46- 49.

[28] 李佳佳,張秀霞,譚光明,等.選擇稀疏矩陣乘法最優(yōu)存儲格式的研究 [J].計算機研究與發(fā)展,2014,51(4):882- 894. Li Jiajia,Zhang Xiuxia,Tan Guangming,et al.Study of choosing the optimal storage format of sparse matrix vector multiplication [J].Journal of Computer Research and Development,2014,51(4):882- 894.

[29] 趙加強.基于OpenCL的稀疏矩陣向量乘優(yōu)化 [D].長春:吉林大學(xué)軟件學(xué)院,2012.

[30] Li Ruipeng,Saad Yousef.GPU-accelerated preconditioned iterative linear solvers [J].The Journal of Supercom-puting,2013,63(2):443- 466.

[31] Larry Meadows.Experiments with WRF on Intel?many integrated core(Intel MIC)architecture [C]∥Proceedings of the IWOMP 2012.Rome:Springer-Verlag Berlin Heidelberg,2012:130- 139.

[32] Intel.Intel Xeon Phi software architecture KNC technical bootcamp [EB/OL].(2014- 03- 15)[2015- 02- 03].http:∥www.intel.cn/content/www/cn/zh/processors/xeon/xeon-phi-coprocessor-system-software-developers-guide.html.

[33] 英特爾亞太研發(fā)有限公司,北京并行科技公司.釋放多核潛能:英特爾Parallel Studio并行開發(fā)指南 [M].北京:清華大學(xué)出版社,2010.

[34] Liu X,Smelyanskiy M,Chow E,et al.Efficient sparse matrix-vector multiplication on x86-based many-core processors [C]∥Proceedings of the 27th International Conference on Supercomputing(ICS).Eugene:Association for Computing Machinery,2013:4996- 5013.

[35] Liu X,Patel A,Chow E.A new scalable parallel algorithm for Fock matrix construction [C]∥Proceedings of the 28th IEEE International Parallel and Distributed Processing Symposium(IPDPS).Phoenix:IEEE,2014:902- 914.

[36] Davis T A,Hu Y.The university of Florida sparse matrix collection [J].ACM Trans on Mathematical Software,2011,38(1):1- 25.

Parallel Optimization and Comparison of Conjugate Gradient Method on GPU and Xeon Phi Coprocessors

HuangMin1DingPing1，2LuoHai-biao2

(1.School of Software Engineering，South China University of Technology，Guangzhou 510006，Guangdong，China;2.Research Center of Parallel Software Research Center，Institute of Software Application Technology，Guangzhou & CAS，Guangzhou 511458，Guangdong，China)

In order to harness the strong horsepower of multi-core processors and meet the demand of high paralle-lism,a new parallel conjugate gradient algorithm is proposed,which focuses on solving the linear equations of large-scale sparse matrices.For the GPU coprocessors,the memory hierarchy of GPU is effectively utilized,optimization methods,such as thread and matrix mappings,data merging and data multiplexing,are adopted,and an effective thread scheduling is conducted to hide the high latency of accessing the global memory of GPU.For Xeon Phi processors,the computation of high parallelism is effectively utilized to optimize data communication and transmission,data dependence reduction,vectorization and asynchronous computation,and effective thread scheduling is also conducted to hide the high latency of accessing global memory of GPU.Finally,the proposed algorithm is proved to be feasible and correct by tests on GPU and Xeon Phi,and its parallel efficiencies in two different ways are compared.It is found that the proposed algorithm on GPU has a better acceleration effect than itself on Xeon Phi.

conjugate gradient method;graphics processing unit;Xeon Phi;parallel optimization;sparse matrix-vector multiplication

2015- 03- 10

廣東省公益研究與能力建設(shè)專項(2014A040401018);廣東省促進科技服務(wù)業(yè)發(fā)展計劃項目(2013B040404009)；廣東省新媒體與品牌傳播創(chuàng)新應(yīng)用重點實驗室資助項目(2013WSYS0002)

黃敏(1976-)，女，博士，副教授，主要從事并行計算和移動云計算研究.E-mail: minh@scut.edu.cn

1000- 565X(2015)11- 0035- 12

TP391.9

10.3969/j.issn.1000-565X.2015.11.006