基于GPU—CUDA的共軛斜量法實現(xiàn)及性能對比

彭土有

摘 要： 偏微分方程數(shù)值解法（包括有限差分法、有限元法）以及大量的數(shù)學(xué)物理方程數(shù)值解法最終都會演變成求解大型線性方程組。因此，探討快速、穩(wěn)定、精確的大型線性方程組解法一直是數(shù)值計算領(lǐng)域不斷深入研究的課題且具有特別重要的意義。在迭代法中，共軛斜量法（又稱共軛梯度法）被公認為最好的方法之一。但是，該方法最大缺點是僅適用于線性方程組系數(shù)矩陣為對稱正定矩陣的情況，而且常規(guī)的CPU算法實現(xiàn)非常耗時。為此，通過將線性方程組系數(shù)矩陣作轉(zhuǎn)換成對稱矩陣后實施基于GPU-CUDA的快速共軛斜量法來解決一般性大型線性方程組的求解問題。試驗結(jié)果表明：在求解效率方面，基于GPU-CUDA的共軛斜量法運行效率高，當(dāng)線性方程組階數(shù)超過3000時，其加速比將超過14；在解的精確性與求解過程的穩(wěn)定性方面，與高斯列主元消去法相當(dāng)。基于GPU-CUDA的快速共軛斜量法是求解一般性大型線性方程組快速而非常有效的方法。

關(guān)鍵詞： GPU； CUDA； 大型線性方程組； 共軛斜量法； 算法； 并行計算

中圖分類號：TP311.1 文獻標(biāo)志碼：A 文章編號：1006-8228（2014）04-04-03

Abstract： The numerical solution for partial differential equations （including finite difference method， and finite element method） and a large number of equations of mathematical physics problems will eventually evolve into solving a large-scale linear equation system. Therefore， studying fast， stable and accurate solutions for large-scale linear equation systems has been a hot topic in the field of numerical calculation for years， which has special significance. Among iterative methods， conjugate gradient method is recognized as one of the best methods. However， this method is only applicable to linear equation systems in which coefficient matrix is symmetric and positive definite. Besides， in conventional CPU implementation， the method for solving a large-scale linear equation system is time-consuming. After the linear equations coefficient matrix A is converted into a symmetric matrix by， the fast conjugate gradient method based on GPU-CUDA is implemented to solve a general large-scale linear equation system. The results show that it is highly efficient. When the rank of the coefficient matrix is over 3000， the speedup will be over 14 times. Besides， it has the same accuracy and stability as Gaussian elimination method. The conjugate gradient method based on GPU-CUDA becomes a fast and effective method for solving large-scale general linear equation systems.

Key words： GPU； CUDA； large-scale linear equation system； conjugate gradient method； algorithm； parallel computation

0 引言

偏微分方程數(shù)值解法（包括有限差分法、有限元法）及大量的數(shù)學(xué)物理方程數(shù)值解法最終都將演變成求解大型線性方程組[4，7]。隨著離散網(wǎng)格數(shù)量的增加，線性方程組的階數(shù)也同步增長。目前，求解實際問題的線性方程組的階數(shù)一般都超過1000階，因此，探討穩(wěn)定、快速、精確的大型線性方程組解法具有特別重要的意義。在迭代法中，共軛斜量法（亦稱共軛梯度法）被公認為最好的方法之一。文獻[5]等所述基于CPU的常規(guī)共軛斜量法已非常成熟，但是，該方法最大缺點是僅適用于線性方程組系數(shù)矩陣為對稱正定矩陣的方程組，而且非常耗時，難以滿足求解1000階以上線性方程組的要求。在應(yīng)用GPU新技術(shù)求解大型線性方程組方面，文獻[6]成功應(yīng)用GPU技術(shù)實現(xiàn)了LU分解和Laplace算法。文獻[1]和文獻[3]分別利用片元著色器編程和預(yù)處理平方根共軛梯度法（Preconditioned Conjugate Gradients Squared，PCGS），實現(xiàn)了基于GPU的共軛梯度法，有效加速了線性方程組的求解，但是，實現(xiàn)方法非常復(fù)雜。本文使用計算統(tǒng)一設(shè)備架構(gòu)CUDA（Compute Unified Device Architecture）技術(shù)，將方程組系數(shù)矩陣A作ATA轉(zhuǎn)換成對稱矩陣后，通過CPU與GPU的協(xié)同來實現(xiàn)基于GPU-CUDA的快速共軛斜量法，實現(xiàn)方法簡單并獲得了較大的加速比。

1 共軛斜量法原理及算法[5，8]

共軛斜量法與最速下降法同屬于一種方法—變分法。首先將線性方程組Ah=b的解h=A-1b（向量）看作是自變量x（向量）的某一個二次函數(shù)的極小點（使函數(shù)f（x）達到極小的點），然后用迭代法求f（x）的極小點，借以求出線性方程組Ah=b的解。

共軛斜量法是在最速下降法基礎(chǔ)上發(fā)展起來的一種快速、有效的線性方程組迭代解法，并且是有限步收斂的迭代法。如下是共軛斜量法的推算過程：

根據(jù)上述共軛斜量法的算法理論，容易設(shè)計出基于CPU的串行算法程序。由于每次迭代都需要進行大量的矩陣與向量乘法、向量與向量的點積運算，所以，當(dāng)線性方程組的階數(shù)很大時，基于CPU的共軛斜量法非常耗時。

2 共軛斜量法的CUDA實現(xiàn)[1-2]

通過分析算法我們發(fā)現(xiàn)，迭代過程中大量的矩陣與向量乘法、向量與向量的點積運算都不需要使用計算完成的結(jié)果，適合于設(shè)計基于GPU-CUDA的并行算法，從而減少計算時間、提高算法的效率。根據(jù)共軛斜量法的推算公式以及CUDA程序設(shè)計思想，即CPU進行數(shù)據(jù)準(zhǔn)備、初始化GPU設(shè)備和執(zhí)行串行代碼，GPU進行并行計算并向CPU返回計算結(jié)果。設(shè)計如下算法完成共軛斜量法基于GPU-CUDA的實現(xiàn)。

⑴ 主機端程序完成數(shù)據(jù)的初始化。

⑵ 申請三塊與系數(shù)矩陣A相等容量的設(shè)備內(nèi)存空間，分別用于存儲系數(shù)矩陣A、系數(shù)矩陣的轉(zhuǎn)置AT及變換結(jié)果ATA；申請六塊與方程組右端向量b相等容量的設(shè)備內(nèi)存空間，用于存儲右端向量b、方程組的解向量x和每次迭代需要的中間向量。

⑶ 將主機端系數(shù)矩陣A及右端向量b拷貝回至設(shè)備端。

⑷ 設(shè)計內(nèi)核并行程序，計算ATA和ATb。

⑸ 調(diào)用內(nèi)核并行程序解線性方程組ATAx=ATb。

⑹ 將方程組的解向量x拷貝回主機。

基于GPU-CUDA的共軛斜量法的關(guān)鍵是并行處理網(wǎng)格的構(gòu)造以及內(nèi)核函數(shù)的設(shè)計，并行處理網(wǎng)格的構(gòu)造及主機內(nèi)核函數(shù)調(diào)用的主要源代碼如下[6，10]。

其中，DIM為線性方程組的階數(shù)。因篇幅所限，部分求解中間向量的內(nèi)核函數(shù)源代碼不再列出。算法實現(xiàn)過程中，為了使每次迭代過程都不進行主機與設(shè)備端的數(shù)據(jù)交換，提高總體并行處理性能，采取了如下三個并行優(yōu)化措施：

⑴ 設(shè)備端一次性分配足夠內(nèi)存空間用于存儲所有矩陣、右端向量以及中間向量；

⑵ 主機端按需要調(diào)用內(nèi)核函數(shù)完成并行處理；

⑶ 并行處理過程中，主機與設(shè)備端不進行任何數(shù)據(jù)交換，一直等待所有并行處理結(jié)束后，才將線性方程組的解x從設(shè)備端拷貝回主機。

3 共軛斜量法的性能測試[9-10]

為了充分驗證算法的正確性、穩(wěn)定性，以及進行算法的性能測試，通過如下源代碼隨機生成任意階線性方程組。

在基于Linux系統(tǒng)的Acer ASPIRE 4736G手提電腦（含512M內(nèi)存、8個流處理器的Geforce G105M圖形卡）上進行了大量共軛斜量法、高斯列主元消去法試驗。試驗結(jié)果表明：并行算法與串行算法求解線性方程組的結(jié)果一致，當(dāng)線性方程組的階數(shù)小于600時，基于GPU-CUDA的并行算法較CPU算法耗時；當(dāng)階數(shù)大于700后，并行算法處理效率明顯高于串行算法，并具有較高的加速比。試驗結(jié)果如表1所示。由于Geforce G105M僅是NVIDIA公司的低端GPU顯卡，所以加速比不高，如果換成具有更多流處理器的GPU顯卡，將會獲得非常大的加速比。

通過與高斯列主元消去法的試算結(jié)果對比，得到如下重要結(jié)果。

⑴ 由于高斯列主元消去法中回代求解過程無法實現(xiàn)并行處理，基于GPU-CUDA的高斯列主元消去法比基于CPU的算法計算速度慢且穩(wěn)定性差。因此，基于CPU的高斯列主元消去法仍然是目前解線性方程組最好的方法之一。

⑵ 對于共軛斜量法，當(dāng)線性方程組階數(shù)大于2000時，基于CPU的算法已相當(dāng)耗時，基于GPU-CUDA的算法仍然具有6以上的加速比。線性方程組階數(shù)越大，加速比越大，階數(shù)超過3000階時，加速比達14以上。

4 結(jié)束語

在分析共軛斜量法求解大型線性方程組原理及算法的基礎(chǔ)上，設(shè)計并實現(xiàn)了基于GPU-CUDA的共軛斜量法的快速算法。通過大量的隨機大型線性方程組的求解試驗，得到以下兩個具有重要意義的結(jié)論。

⑴ 基于GPU-CUDA的共軛斜量法不僅運行效率高，當(dāng)線性方程組階數(shù)超過700階時，算法的加速比超過6，線性方程組階數(shù)超過3000階時，加速比達14以上。與此同時，基于GPU-CUDA的共軛斜量法在求解的精確性方面與列主元消去法相當(dāng)。當(dāng)線性方程組階數(shù)大于3500后，盡管加速比變得更大，但是，其穩(wěn)定性較CPU算法差，主要原因是基于GPU算法的誤差積累所致。

⑵ 通過將一般性線性方程組系數(shù)矩陣A作ATA轉(zhuǎn)換，右端項b作ATb轉(zhuǎn)換后，即可應(yīng)用基于GPU-CUDA的共軛斜量法求解線性方程組，從而拓展了共軛斜量法的應(yīng)用范圍，為求解一般性大型線性方程組提供了快速、有效的方法。

參考文獻：

[1] NVIDIA. CUDA Programming Guide[OL].2012-12-22. http：//www.nvidia.com/object/cuda_home.html.

[2] 夏健明，魏德敏著.共軛梯度法的GPU實現(xiàn)[J].微計算機工程，2009.17：274-276

[3] 李熙銘，歐陽丹彤，白洪濤著.基于GPU的混合精度平方根共軛梯度算法[J].儀器儀表學(xué)報，2012.1：98-104

[4] 李榮華，劉播著.微分方程數(shù)值解法（第四版）[M].高等教育出版社，2010.

[5] 蔣長錦編著.科學(xué)計算與C程序集[M].中國水利水電出版社，2010.

[6] 陳穎，林錦賢，呂暾著.LU分解和Laplace算法在GPU上的實現(xiàn)[J].計算機應(yīng)用，2011.3：851-855

[7] 呂英華編著.計算電磁學(xué)的數(shù)值方法[M].人民教育出版社，2006.

[8] 武漢大學(xué)，山東大學(xué)著.計算方法[M].科學(xué)出版社，1983.

[9] 肖江，胡柯良，鄧元勇著.基于CUDA的矩陣乘法和FFT性能測試[J].計算機工程，2009.10：7-10

[10] 劉麗，沈杰，李洪林著.基于GPU的矩陣求逆性能測試和分析[J].華東理工大學(xué)學(xué)報（自然科學(xué)版），2010.6：812-816