基于GPU—CUDA的共軛斜量法實現及性能對比

彭土有

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

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

中圖分類號：TP311.1 文獻標志碼：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 引言

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

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

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

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

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

2 共軛斜量法的CUDA實現[1-2]

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

⑴ 主機端程序完成數據的初始化。

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

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

⑷ 設計內核并行程序，計算ATA和ATb。

⑸ 調用內核并行程序解線性方程組ATAx=ATb。

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

基于GPU-CUDA的共軛斜量法的關鍵是并行處理網格的構造以及內核函數的設計，并行處理網格的構造及主機內核函數調用的主要源代碼如下[6，10]。

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

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

⑵ 主機端按需要調用內核函數完成并行處理；

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

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

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

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

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

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

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

4 結束語

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

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

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

參考文獻：

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

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

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

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

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

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

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

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

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

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