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

高家全，王志超

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

?

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

高家全，王志超

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

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

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

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

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

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

(1)

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

(2)

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

2GPU實現(xiàn)

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

圖1　計算框架Fig.1　The computing architecture

2.1初始化

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

2.2計算階段

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

(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)存來提高計算速度

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)存訪問不對齊，保證全局內(nèi)存合并，減少內(nèi)存交易 (Memory transaction) 次數(shù)，本研究采用了128對齊填充模式.假定A的前四行非零元素個數(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實驗和分析

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

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

表1　實驗使用的矩陣類型

3.1實驗分析

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

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

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

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

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

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

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

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

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

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

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

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

表5　CPU和GPU時間比較

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

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

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

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

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

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

4結(jié)論

提出了一種GPU加速的SSOR稀疏近似逆預(yù)條件子，調(diào)研了一階和二階近似逆的有效性.為保證其稀疏性，提出了一種解決方案.為提高算法的性能，利用共享內(nèi)存、內(nèi)存訪問合并及填充等策略進(jìn)行優(yōu)化.實驗證明設(shè)計的SSOR稀疏近似逆預(yù)條件子是有效的，能顯著提高PCG算法的性能，而且對比在單個CPU上使用Python實現(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é)報，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].浙江工業(yè)大學(xué)學(xué)報，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上稀疏矩陣與矢量乘積運算的一種改進(jìn)[J].計算機(jī)系統(tǒng)應(yīng)用，2010，19(5)：1906-1912.

[11]周洪，樊曉椏，趙麗麗.基于CUDA的稀疏矩陣與矢量乘法的優(yōu)化[J].計算機(jī)測量與控制，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

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

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

中圖分類號：TP338.6

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

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