在GPU上實現(xiàn)Jacobi迭代法的分析與設(shè)計

吳玫華

（福建省凱特科技有限公司 福建 福州 350002）

隨著 GPU（Graphic Processing Unit，圖形處理器）技術(shù)的快速發(fā)展，GPU的浮點運算能力迅速提升，并大幅度超過了通用CPU（Central Processing Unit）。早期GPU主要用于圖形處理，硬件設(shè)計上主要考慮數(shù)據(jù)的并行，一般擁用大量的計算單元，而缺少CPU中用于分支預(yù)測、亂序執(zhí)行等機(jī)制的邏輯控制單元，因此可利用GPU上的大量計算單元提升計算性能。Nvidia公司在2008年推出的GeForce GTX 280顯卡的浮點計算峰值達(dá)到950Gflop/s，相對于當(dāng)前主流多核CPU的性能提升了許多[1]。GPU的強(qiáng)大計算性能使得大批研究者投入到將大量計算需要和復(fù)雜的問題映射成GPU可處理的問題，即利用GPU來做的通用計算GPGPU（General Purpose Computation on GPU）[2]。

Jacobi迭代法[3]是科學(xué)計算領(lǐng)域中常用的計算方法，其核心計算是矩陣乘向量。Jacobi迭代法廣泛用于求解線性方程組，是許多數(shù)值計算問題的核心，用于構(gòu)成更高級的計算方法。研究這個經(jīng)典算法在GPU上的實現(xiàn)對在GPU平臺上開發(fā)相關(guān)科學(xué)計算程序具有重要的借鑒意義。

筆者在一個含有 Intel Core（TM）Duo E7200 CPU和Nvidia GeForce 9500GT GPU的計算機(jī)上實現(xiàn)Jacobi迭代法，并對算法性能進(jìn)行測試。實驗結(jié)果表明把算法映射到GPU上執(zhí)行能夠獲得較好的加速比，最高得到10.2的加速比。

1 GPU通用計算

最早的GPGPU直接使用了圖形學(xué) API（Application Program Interface）編程。這種開發(fā)方式要求編程人員將數(shù)據(jù)打包成紋理，將計算任務(wù)映射為對紋理的渲染過程，用匯編或者高級著色語言（如 GLSL、Cg、HLSL）編寫 Shader程序，然后通過圖形學(xué) API（如 Direct3D、OpenGL）執(zhí)行。這種方式要求編程人員不僅要熟悉自己需要實現(xiàn)的計算和并行算法，還要對圖形學(xué)硬件和編程接口有深入的了解。由于開發(fā)難度大，傳統(tǒng)的GPGPU沒有被廣泛地應(yīng)用[4-8]。

盡管在這期間許多非視覺的科學(xué)計算問題被成功地映射到顯卡上[9-10]，但是Nvidia公司發(fā)布的CUDA和開放運算語言O(shè)penCL[11]框架才真正引領(lǐng)高性能計算進(jìn)入了通用GPU計算時代。Nvidia公司于2007年正式發(fā)布的CUDA，是一種將GPU作為數(shù)據(jù)并行計算設(shè)備的軟硬件體系。

Nvidia的CUDA架構(gòu)通過對硬件的重新組織，使GPU能夠應(yīng)用于更一般的領(lǐng)域。CUDA試圖通過提供一般的高層和低層的API來訪問GPU的并行元素，以減輕將問題映射到GPU上的不便[12]。現(xiàn)在的GPU特別適合計算密集型、高度并行化的高性能計算[13]。CUDA提供了對顯卡的抽象，把顯卡作為能夠同時執(zhí)行大量輕量級線程的設(shè)備。這些線程組織成塊，塊中的每個線程能夠訪問它們自身的寄存器和塊的共享存儲器，每個線程也能夠與其相鄰的線程進(jìn)行同步。一個內(nèi)核函數(shù)的代碼由一個或多個這樣的塊執(zhí)行。由于具有CUDA能力的設(shè)備允許DRAM的類屬存?。ň奂头稚ⅲ?，所以每個線程都能訪問GPU板卡上的顯存和紋理存儲器。

2 Jacobi迭代的實現(xiàn)

Jacobi迭代是用于求n元一次方程組ax=b的解。其中a為系數(shù)矩陣，x為解向量，b為常數(shù)向量。

Jacobi迭代的數(shù)學(xué)表示為：

式（1）為給定的初始解向量，式（2）為第k+1次迭代中第i個元素的公式。當(dāng)前后兩次解向量的相對誤差小于某個精度或迭代到一定次數(shù)，則認(rèn)為算法結(jié)束，退出迭代，輸出結(jié)果?？梢园l(fā)現(xiàn)式（2）中包含了矩陣乘向量以及向量的四則運算，這些都是可以并行的。

圖1 Jacobi迭代的實現(xiàn)過程Fig.1 Basically implement of Jacobi

根據(jù)CUDA的編程模型，把每次解向量的計算放到GPU中運行，即kernel函數(shù)中實現(xiàn)。每個線程計算解向量的一個元素，因此可以得出每個線程上的數(shù)據(jù)和計算，如圖1所示。映射到第i個線程的數(shù)據(jù)包括：矩陣a對角線的第i個元素aii；常數(shù)向量b的第i個元素bi；第k次迭代的解向量x中第i個元素xki。這3個數(shù)據(jù)與線程是一一對應(yīng)的，即第i個線程只訪問對應(yīng)數(shù)據(jù)序列中的第i個元素。此外，在計算向量點積中還應(yīng)包含矩陣a第i行的所有元素以及第k次迭代的整個x列向量。

1）基本算法

x_per=x_now= （0， …， 0）； //設(shè)置初始解向量，x_per為前一次的解向量；x_now為當(dāng)前解向量//

do

{

count++；//統(tǒng)計迭代次數(shù)//

Jacobi_kernel＜＜＜grid， thread＞＞＞（ a， b， x_per， x_now ）；//調(diào)用GPU內(nèi)核函數(shù)處理，每個內(nèi)核函數(shù)處理一個解向量中的一個元素//

CopyDataDeviceToHost（x_per，x_now ）； //將 GPU 中的數(shù)據(jù)復(fù)制回內(nèi)存//

if（distance（x_per，x_now ）＜eps||count＞MaxSteps）

Break；//兩次的近似相對誤差足夠小或達(dá)到迭代次數(shù)上限則退出循環(huán)//

}while（ture）.

以上為Jacobi迭代的基本實現(xiàn)方法。注意到GPU處理浮點數(shù)的速度較快，但GPU與CPU之間的數(shù)據(jù)傳輸會消耗大量時間。因此對基本實現(xiàn)方法進(jìn)行一些修改，不用每次迭代后就將數(shù)據(jù)傳輸?shù)紺PU進(jìn)行相對誤差和迭代次數(shù)的判斷，可以執(zhí)行若干次迭代后再將數(shù)據(jù)傳輸回CPU進(jìn)行判斷，這樣減少傳輸時間，從而提高執(zhí)行效率。在理想情況下，在存儲器傳輸進(jìn)行的同時，GPU中的各線程也在計算。

2）改進(jìn)算法

x_per=x_now= （0， …， 0）； //設(shè)置初始解向量，x_per為前一次的解向量；x_now為當(dāng)前解向量//

do{

count++；//統(tǒng)計迭代次數(shù)//

Jacobi_kernel＜＜＜grid， thread＞＞＞（ a， b， x_per， x_now ）；//調(diào)用GPU內(nèi)核函數(shù)處理，每個內(nèi)核函數(shù)處理一個解向量中的一個元素//

if（count%10==0） //每執(zhí)行10次迭代后進(jìn)行判斷//{

CopyDataDeviceToHost（ x_per， x_now ）； //將 GPU 中的數(shù)據(jù)復(fù)制回內(nèi)存//

if（distance（x_per，x_now ）＜eps||count＞MaxSteps）

Break；//兩次的近似相對誤差足夠小或達(dá)到迭代次數(shù)上限則退出循環(huán)//

}

}while（ture）.

3 實驗結(jié)果

在 Intel Core2DuoE7200CPU（2core@2.53GHz）的 OpenMP 2.0計算平臺和 NVIDIA GeForce 9500GT GPU （650 MHz）的CUDA SDK 2.0計算平臺上分別測試由算法改寫的程序在GPU上的運行情況和原始程序在CPU上的運行情況，其中CPU上的程序采用C++語言編寫，并使用OpenMP在雙核上并行執(zhí)行。

在不同矩陣規(guī)模下測試了在GPU下的兩種實現(xiàn)策略與雙核CPU程序的性能，結(jié)果如圖2所示。

圖2 實驗結(jié)果性能對比Fig.2 Speedup of algorithms

圖2 分別給出了Jacobi迭代在各種規(guī)模下雙核CPU程序、GPU基本實現(xiàn)、GPU改進(jìn)實現(xiàn)的性能對比。從圖中可以看出，GPU程序加速比隨著數(shù)據(jù)規(guī)模的增大而增大。GPU處理單元的頻率（650 MHz）遠(yuǎn)小于 CPU（2.53 GHz）的頻率。 GPU的性能優(yōu)勢來自于眾多的計算單元，以龐大的吞吐率同時進(jìn)行多點計算，因此要保證GPU性能的發(fā)揮，必須提供足夠多的可以并行的計算量。從圖2中還可以看出，在矩陣規(guī)模為256×256時，GPU的性能低于雙核CPU程序。一方面是因為計算量較小，無法發(fā)揮GPU的性能，另一方面是因為調(diào)用kernel函數(shù)也存在一定開銷，如果計算量太小則計算的開銷無法抵消這部分調(diào)用的開銷。

4 結(jié)束語

將數(shù)值計算有效地映射到GPU上執(zhí)行必須考慮多方面的因素，包括訪問的效率、分支語句的處理等。本文將科學(xué)計算領(lǐng)域中常用的Jacobi迭代在GPU上實現(xiàn)，對數(shù)據(jù)訪問進(jìn)行優(yōu)化，獲得了良好的加速比。從實驗結(jié)果也可以看出，隨著矩陣規(guī)模的增大GPU程序的加速比也在增大，說明算法的可擴(kuò)展性良好。

[1]NVIDIA.NVIDIA CUDA programming guid version 2.0[S].2008.

[2]JD Owens，M Houston，D Luebke， et al.GPU Computing[J].Proceedings of the IEEE，2008，96（5）:879-899.

[3]Jeffery JL.Numerical Analysis and Scientific Computation[M].Addison Wesley，2004.

[4]Kessenich J， Baldwin D， Rost R.The OpenCL Shading Language[R].2003.

[5]Mark W R，Glanville R S，Akeley K，et al.Cg:A system for programming graphics hardware in a clike language[J].ACM Transactions on Graphics-Proceedings of ACM SIGGRAPH，2003，22（3）:897-907.

[6]Microsoft.High2Level Shader Language[S].2003.

[7]Pike A.DirectX 8 Tutorial[R].Retrieved March 15，2006.

[8]Shreiner D，Woo M，Neider J，etal.OpenGLProgramming Manual[M].OpenGL ARB， Boston.5th ed.Addison Wes2ley，2005.

[9]Svetlin A M.CUDA compatible GPU as an efficient hard-ware accelerator for AEScryptography[C]//In Proceedings of IEEE International Conference on Signal Processing and Communication，2007:65-68.

[10]Szerwinski R，Guneysu T.Exploiting the power of GPUs for Asymmetric Cryptography[J].Lecture Notes in Computer Science，2008:79-99.

[11]KHRONOS:OpenCL overview web page[EB/OL]. （2009）http://www.khronos.org/opencl/，2009.

[12]Halfhill T.Parallel Processing with CUDA[M].Microprocessor，2008.

[13]Gummaraju J，Rosenblum M. Stream processing in generalpur-pose processors[S].2005.