吳玫華

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

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

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

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

1 GPU通用計算

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

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

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

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

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

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

式（1）為給定的初始解向量，式（2）為第k+1次迭代中第i個元素的公式。當前后兩次解向量的相對誤差小于某個精度或迭代到一定次數(shù)，則認為算法結束，退出迭代，輸出結果?？梢园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ù)與線程是一一對應的，即第i個線程只訪問對應數(shù)據(jù)序列中的第i個元素。此外，在計算向量點積中還應包含矩陣a第i行的所有元素以及第k次迭代的整個x列向量。

1）基本算法

x_per=x_now= （0， …， 0）； //設置初始解向量，x_per為前一次的解向量；x_now為當前解向量//

do

{

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

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

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

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

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

}while（ture）.

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

2）改進算法

x_per=x_now= （0， …， 0）； //設置初始解向量，x_per為前一次的解向量；x_now為當前解向量//

do{

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

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

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

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

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

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

}

}while（ture）.

3 實驗結果

在 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程序的性能，結果如圖2所示。

圖2 實驗結果性能對比Fig.2 Speedup of algorithms

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

4 結束語

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

[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.