高家全，王志超

(浙江工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，浙江 杭州 310023)

?

基于GPU的SSOR稀疏近似逆預(yù)條件研究

高家全，王志超

(浙江工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，浙江 杭州 310023)

摘要：由于SSOR預(yù)條件共軛梯度算法中預(yù)條件方程求解需要前推和回代，導(dǎo)致算法遷移到GPU平臺(tái)上并行效率不高.為此，基于諾依曼多項(xiàng)式分解技術(shù)，提出了一種GPU加速的SSOR稀疏近似逆預(yù)條件子(GSSORSAI).它不僅保持了原線性系統(tǒng)系數(shù)矩陣的稀疏和對(duì)稱正定特性，而且預(yù)條件方程求解僅需一次稀疏矩陣矢量乘運(yùn)算，避免了前推和回代過程.實(shí)驗(yàn)結(jié)果表明：在NVIDIA Tesla C2050 GPU上，對(duì)比使用Python在單個(gè)CPU上SSOR稀疏近似逆預(yù)條件子實(shí)現(xiàn)方法，GSSORSAI平均快將近100倍；應(yīng)用到并行的PCG算法中，相比無預(yù)條件的CG算法，平均提高了算法的3倍的收斂速度.

關(guān)鍵詞：SSOR預(yù)條件子；預(yù)條件共軛梯度算法；稀疏近似逆；GPU

求解稀疏對(duì)稱正定的線性系統(tǒng)Ax=b是科學(xué)計(jì)算和工程應(yīng)用的核心問題.一般采用迭代法求解，其中共軛梯度 (Conjugate gradient, CG) 算法是目前流行的解此類線性系統(tǒng)的迭代法之一[1-2].與其他迭代法一樣，CG算法的性能和效率依賴于線性系統(tǒng)的系數(shù)矩陣條件數(shù)的好壞，壞的系數(shù)矩陣條件數(shù)將會(huì)導(dǎo)致CG算法收斂慢或不收斂[1-4].為解決這方面的問題，求解線性系統(tǒng)之前需進(jìn)行預(yù)處理.預(yù)處理技術(shù)就是構(gòu)造一個(gè)預(yù)條件子M，使得Max=Mb或者AMy=b,x=My.選擇合適的預(yù)條件子能夠減少算法的迭代次數(shù)進(jìn)而減少其執(zhí)行時(shí)間.SSOR (Symmetric successive over-relaxation) 預(yù)條件子用于CG算法的預(yù)處理已證明是有效的，然而由于SSOR預(yù)條件方程求解需要前推和回代過程，很難進(jìn)行并行化處理，這很大程度上降低了SSOR預(yù)條件CG (PCG) 算法解決大規(guī)模線性系統(tǒng)的性能和效率[3-4].

基于諾依曼(Neumann)多項(xiàng)式分解技術(shù)，提出了一種SSOR稀疏近似逆預(yù)條件子，并在GPU (Graphic processing unit)上并行實(shí)現(xiàn).提出的SSOR稀疏近似逆預(yù)條件子的主要運(yùn)算由矩陣求和與兩矩陣相乘組成，易于并行實(shí)現(xiàn)，而且應(yīng)用于PCG算法，預(yù)條件方程求解僅需一個(gè)稀疏矩陣矢量乘運(yùn)算，避免了SSOR預(yù)條件子的前推和回代過程.筆者主要做了如下幾個(gè)貢獻(xiàn)：首先提出了一種SSOR稀疏近似逆預(yù)條件子，不僅保持了稀疏性，而且還保證了其對(duì)稱正定特性；然后探討了SSOR一階和二階稀疏近似逆的性能，并給出了高效的GPU并行實(shí)現(xiàn)算法；最后應(yīng)用于PCG算法，驗(yàn)證了其有效性.

1SSOR近似逆預(yù)條件子

(1)

為保證其稀疏性，給出了如下方案：如果

(2)

為簡便，僅討論SSOR一階近似逆和二階近似逆兩種情況.

2GPU實(shí)現(xiàn)

基于CUDA，設(shè)計(jì)了一個(gè)有效的SSOR稀疏近似逆預(yù)條件子計(jì)算框架(圖1).框架分初始化、計(jì)算和應(yīng)用三個(gè)階段.

圖1　計(jì)算框架Fig.1　The computing architecture

2.1初始化

首先復(fù)制矩陣A到GPU.對(duì)于稀疏矩陣A，為提高運(yùn)算效率和節(jié)省存儲(chǔ)空間，采用CSR的壓縮存儲(chǔ)方式[9-10].

2.2計(jì)算階段

首先計(jì)算對(duì)角矩陣D.對(duì)角矩陣D由矩陣A的對(duì)角元素組成，核函數(shù)設(shè)計(jì)如下：讓每個(gè)線程 (Thead)負(fù)責(zé)尋找A的一行中的對(duì)角元素值，并存儲(chǔ)到全局內(nèi)存中.共需要Row(A的行數(shù))個(gè)線程，如果每個(gè)線程塊 (Block) 包含256個(gè)線程，則需要Row/256個(gè)線程塊.

(3)

Input: A, K, M;

Output: M;

__shared__ P[];

//利用GPU中的線程、線程塊和網(wǎng)格屬性

tid←blockDim.x* blockIdx.x+ threadIdx.x;

lane←tid & 31;wid←threadIdx.x>> 5;row←wid;

forifrom K.mPtr[row] to K.mPtr[row+1] withi++

P[i] = K.mIndex[i];

end

forjfrom A.mPtr[row]+ lane to A.mPtr[row+1]withj+= 32

// Each element of M← one row of K * one column of K

column←A.mIndex[j];

forkfrom K.mPtr[column] to K.mPtr[column] withk++

forlfrom K.mPtr[row] to K.mPtr[row+1] withl++

//使用共享內(nèi)存來提高計(jì)算速度

if(K.mIndex[k] == P[l])

dot←dot + K.mData[k] * K.mData[l];

break;

end

end

M.mData[j]←M.mData[j] + dot;

M.mIndex[j]←A.mIndex[j];

end

為解決內(nèi)存訪問不對(duì)齊，保證全局內(nèi)存合并，減少內(nèi)存交易 (Memory transaction) 次數(shù)，本研究采用了128對(duì)齊填充模式.假定A的前四行非零元素個(gè)數(shù)分別為12，26，8，16，并且精度為雙精度 (Double)類型.從圖2可以看出：填充前，共需7 次內(nèi)存交易，然而填充后，前4 行的交易次數(shù)可以從7 次減小為5 次，進(jìn)而提高了將近30%左右的內(nèi)存訪問效率.

圖2　填充模式Fig.2　Padding scheme

2.3應(yīng)用階段

3實(shí)驗(yàn)和分析

這節(jié)將從如下幾個(gè)方面來測(cè)試筆者提出的SSOR預(yù)條件子的一階SSOR稀疏近似逆 (表示為M1) 和二階SSOR稀疏近似逆 (表示為M2) 的性能： 1) 分析參數(shù)w和α對(duì)M1和M2性能的影響；2) 驗(yàn)證設(shè)計(jì)的SSOR稀疏近似逆算法的并行性；3) 應(yīng)用到PCG算法中，來驗(yàn)證M1和M2的有效性.

在下面實(shí)驗(yàn)里， PCG算法的終止條件為最大迭代次數(shù)不多于5 000或殘差小于10-12.所有的實(shí)驗(yàn)都執(zhí)行在NVIDIA Tesla C2050 GPU上，測(cè)試矩陣來自Florida大學(xué)的稀疏矩陣集[12]，表1給出了它們簡要的描述.

表1　實(shí)驗(yàn)使用的矩陣類型

3.1實(shí)驗(yàn)分析

首先，測(cè)試稀疏性參數(shù)α對(duì)一階SSOR稀疏近似逆M1和二階稀疏SSOR近似逆M2性能的影響.為簡便，在這個(gè)實(shí)驗(yàn)里，取w=0.8.表2給出了實(shí)驗(yàn)結(jié)果，每一欄對(duì)應(yīng)α值的M1/ M2應(yīng)用到PCG中算法收斂的迭代次數(shù).迭代次數(shù)越少，表明獲得的M1/ M2的性能越好.

從表2可以看出：對(duì)M1和M2來說，PCG收斂迭代次數(shù)一般隨著α值的增加而減少.這意味著，一般情況下，只要取α=1.0時(shí)，保持M1和M2與A有著相同的稀疏結(jié)構(gòu)時(shí)，M1和M2的性能最好.

表2　參數(shù)對(duì)M1和M2性能的影響1)

注：1) N/A表示迭代次數(shù)>5 000.

其次，調(diào)查w對(duì)M1和M2性能的影響.不失一般性，在這個(gè)實(shí)驗(yàn)里，取α=1.0，表3，4分別列出了w變化對(duì)應(yīng)M1和M2的PCG算法收斂的迭代次數(shù).同樣，迭代次數(shù)越少，表明獲得的M1/M2的性能越好.

從表3，4可以看出：類似于，w值的選取對(duì)M1和M2性能有影響，而且不同的w值，獲得的M1/M2性能差異明顯.因此，為得到最優(yōu)性能的M1和M2，對(duì)給定的稀疏矩陣，最優(yōu)w值的選取實(shí)驗(yàn)是必要的.

表3　參數(shù)w對(duì)M1性能的影響1)

注：1) N/A表示迭代次數(shù)>5 000.

表4　參數(shù)w對(duì)M2性能的影響1)

注：1) N/A表示迭代次數(shù)>5 000.

3.2SSOR稀疏近似逆算法的并行性

此節(jié)來測(cè)試提出的SSOR稀疏近似逆算法的并行性.從表3，4很容易計(jì)算出：為使用的6個(gè)測(cè)試矩陣，當(dāng)α=1.0，w=0.9時(shí)，PCG平均迭代次數(shù)最小.不失一般性，這節(jié)實(shí)驗(yàn)取w=1.0 和α=0.9，分別計(jì)算一階/二階SSOR稀疏近似逆算法的GPU執(zhí)行時(shí)間和對(duì)應(yīng)串行CPU的執(zhí)行時(shí)間 (串行的實(shí)現(xiàn)是基于編程模型Python).表5給出CPU和GPU執(zhí)行時(shí)間的比較.

表5　CPU和GPU時(shí)間比較

從表5可以看出：一階稀疏近似逆情況下，CPU與GPU的執(zhí)行時(shí)間比率介于54.1～227之間，平均比率大約101；二階稀疏近似逆情況下，CPU與GPU的執(zhí)行時(shí)間比率介于43.4～220之間，平均比率大約90；總平均比率95.這表明本研究設(shè)計(jì)的GPU加速的一階/二階SSOR稀疏近似逆并行算法有著高度并行性.

3.3SSOR稀疏近似逆算法的性能

將M1和M2應(yīng)用到PCG算法中，測(cè)試其性能.為每一個(gè)測(cè)試矩陣，α=1.0，w取表3，4中對(duì)應(yīng)此測(cè)試矩陣平均迭代次數(shù)最小的w值.例如，對(duì)矩陣ecology2，取w=1.1.圖3，4分別給出了GPU加速的CG算法、一階SSOR稀疏近似逆PCG算法及二階SSOR稀疏近似逆PCG算法的迭代次數(shù)和執(zhí)行時(shí)間.時(shí)間單位為毫秒(ms).

圖3　執(zhí)行時(shí)間Fig.3　The execution time

圖4　迭代次數(shù)Fig.4　The number of iterations

從圖3，4可以看出：相對(duì)于直接CG算法，使用SSOR一階/二階稀疏近似逆預(yù)條件子后，PCG算法的迭代次數(shù)顯著減少，而且執(zhí)行時(shí)間也明顯少于直接CG算法.這表明提出的SSOR一階/二階稀疏近似逆預(yù)條件子是有效的.相對(duì)于一階稀疏近似逆預(yù)條件子M1，雖然取得二階稀疏近似逆預(yù)條件子M2花費(fèi)的時(shí)間略多(3.2節(jié))，但大大減少了PCG算法的迭代次數(shù)，進(jìn)而減少了其執(zhí)行時(shí)間.因此，可以斷言，一般情況下，M2優(yōu)于M1.

4結(jié)論

提出了一種GPU加速的SSOR稀疏近似逆預(yù)條件子，調(diào)研了一階和二階近似逆的有效性.為保證其稀疏性，提出了一種解決方案.為提高算法的性能，利用共享內(nèi)存、內(nèi)存訪問合并及填充等策略進(jìn)行優(yōu)化.實(shí)驗(yàn)證明設(shè)計(jì)的SSOR稀疏近似逆預(yù)條件子是有效的，能顯著提高PCG算法的性能，而且對(duì)比在單個(gè)CPU上使用Python實(shí)現(xiàn)的SSOR稀疏近似逆預(yù)條件子，提出的GPU加速的SSOR稀疏近似逆預(yù)條件子算法平均快將近100 倍.

參考文獻(xiàn)：

[1]GOLUB G H, VAN LOAN C F. Matrix computations[M]. Baltimore: The John Hopkins University Press，1989.

[2]王廠文，張有正.正定矩陣和行列式不等式[J].浙江工業(yè)大學(xué)學(xué)報(bào)，2006，34(3)：351-355.

[3]KAASSCHIETER E F. Preconditioned conjugate gradients for solving singular systems[J]. CAM journal，1988(24)：265-275.

[4]龍愛芳.避免二階導(dǎo)數(shù)計(jì)算的迭代法[J].浙江工業(yè)大學(xué)學(xué)報(bào)，2005，33(5)：602-604.

[5]AMENT M, KNITTEL G, WEISKOPF D. A parallel preconditioned conjugate gradient solver for the Poisson problem on a multi-GPU[J]. Parallel distributed and network-based computing，2010(6)：583-592.

[6]SAAD Y. Iterative Methods for Sparse Linear Systems[M]. Philadelphia: SIAM，2003.

[7]CHOW E, SAAD Y. Approximate inverse preconditioners via sparse-sparse iterations[ J]. SIAM journal，1998(19)：995-1023.

[8]GORTE M J, HUCKLE T. Parallel preconditioning with sparse approximate inverses[J]. SIAM journal，1997(18)：838-853.

[9]BELL N, GARLAND M. Efficient sparse matrix-vector multiplication on CUDA[R]. Santa Clara: NVIDIA,2008.

[10]馬超，韋剛，裴頌文，等.GPU上稀疏矩陣與矢量乘積運(yùn)算的一種改進(jìn)[J].計(jì)算機(jī)系統(tǒng)應(yīng)用，2010，19(5)：1906-1912.

[11]周洪，樊曉椏，趙麗麗.基于CUDA的稀疏矩陣與矢量乘法的優(yōu)化[J].計(jì)算機(jī)測(cè)量與控制，2010，18(8)：116-120.

[12]DAVIS T A, HU Y. The University of Florida sparse matrix collection [J]. ACM transactions on mathematical software，2011，38(1)：1-25.

(責(zé)任編輯：陳石平)

Research of the SSOR sparse approximate inverse preconditioner on GPUs

GAO Jiaquan, WANG Zhichao

(College of Computer Science and Technology Engineering, Zhejiang University of Technology, Hangzhou 310023, China)

Abstract:For the SSOR preconditioned conjugate gradient algorithm, the preconditioner equation solving needs the forward/backward substitutions, which greatly prevents parallelizing SSOR PCG algorithms on the GPU platform due to their strong serial processing. Thus, based on the Neumann series approximation, a GPU accelerated SSOR sparse approximate inverse preconditioner is proposed. For GSSORSAI, it preserves the sparse and symmetric positive characteristics of the original coefficient matrix in the linear system, and the preconditioner equation solving only needs a sparse matrix-vector multiplication operation, which avoids the forward/backward substitutions. Experiments results show on the NVIDIA Tesla C2050 GPU, GSSORSAI is generated on average 100 times faster than the implementation by Python on single CPU. Compared to the convergence of the CG algorithm, the PCG algorithm with GSSORSAI has on average 3 times faster convergent rate.

Keywords:SSOR preconditioner; preconditioned conjugate gradient algorithm; sparse approximate inverse; GPU

收稿日期：2015-09-21

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目(61379017)

作者簡介：高家全(1972—)，男，河南固始人，副教授，博士，研究方向?yàn)楦咝阅苡?jì)算，大數(shù)據(jù)分析與可視化，智能算法及應(yīng)用，E-mail：gaojq@zjut.edu.cn.

中圖分類號(hào)：TP338.6

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1006-4303(2016)02-0140-06