張 濤 （長(zhǎng)江大學(xué)信息與數(shù)學(xué)學(xué)院，湖北 荊州434023水資源與水電科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室 （武漢大學(xué)），湖北 武漢430072）廖文軍 （中國(guó)石油集團(tuán)東方地球物理勘探有限責(zé)任公司物探技術(shù)研究中心，河北 涿州072751）

在工程計(jì)算中，經(jīng)常需要對(duì)Hermite方程組進(jìn)行求解，而且由于需要保證求解的精確性，一般都采用直接法進(jìn)行求解，在問(wèn)題規(guī)模較大的情況下，該方程組的求解時(shí)間開(kāi)銷很大。如果單純?cè)贑PU上進(jìn)行計(jì)算，計(jì)算效率往往滿足不了實(shí)際生產(chǎn)的需要。隨著GPU技術(shù)［1］的出現(xiàn)，利用GPU眾核處理能力，采用CUDA計(jì)算［2-3］架構(gòu)，利用GPU多核并行實(shí)現(xiàn)Hermite方程組的求解，將對(duì)解決其計(jì)算瓶頸提供契機(jī)。

1 Hermite方程組解法的GPU并行化分析

筆者所討論的Hermite方程組解法是一種類似LU分解的直接法，本質(zhì)上仍然是Gauss消去法，其主要的計(jì)算過(guò)程是消元操作。無(wú)論Gauss消去法還是LU分解法及其他消去類的直接法，以按列選主元為例，在每一次消元過(guò)程中，各列的消去操作是獨(dú)立的，是可以并行進(jìn)行的，而且在每一列內(nèi)，各行元素的消去操作也是可以并行進(jìn)行的。因此，可以采用GPU來(lái)加速Hermite方程組的求解。

2 Hermite方程組解法的GPU并行算法

2．1 算法設(shè)計(jì)

為了達(dá)到好的加速比，預(yù)測(cè)算子計(jì)算的GPU并行算法應(yīng)遵循GPU的SIMD體系結(jié)構(gòu)特點(diǎn)，并且充分發(fā)揮其眾核計(jì)算能力，整個(gè)算法采用三級(jí)并行模型，其具體設(shè)計(jì)過(guò)程如下：

1）方程組集之間的并行 定義多個(gè)Hermite方程組為一個(gè)方程組集 （每個(gè)集包含相同個(gè)數(shù)的多個(gè)Hermite方程組），把一個(gè)方程組集的求解設(shè)計(jì)為一個(gè)GPU Kernel，由于Fermi架構(gòu)的 GPU支持16個(gè)GPU Kernel同時(shí)并行處理，為充分利用GPU的流處理器 （SM），可以設(shè)計(jì)多個(gè)GPU Kernel使其并行處理多個(gè)方程組集的求解，假設(shè)一個(gè)GPU Kernel只利用了7個(gè)SM資源，一個(gè)Fermi GPU有14個(gè)SM，那么一塊GPU可以并行求解2個(gè)Hermite方程組集，其并行對(duì)應(yīng)關(guān)系如圖1所示。

圖1 Hermite方程組集的并行對(duì)應(yīng)關(guān)系

2）每個(gè)方程組集內(nèi)多個(gè)Hermite方程組并行 每個(gè)方程組集包括N個(gè)Hermite方程組，由于每個(gè)Hermite方程組的求解是獨(dú)立的，彼此之間不存在數(shù)據(jù)依賴性，方程組之間可以并行求解。把每一個(gè)Hermite方程組的求解設(shè)計(jì)為一個(gè)GPU Block，即每一個(gè)GPU Block負(fù)責(zé)求解一個(gè)Hermite方程組，那么對(duì)于一個(gè)GPU Kernel而言，它有N個(gè)GPU Block，由于GPU Block運(yùn)行于SM內(nèi)，每個(gè)Hermite方程組的求解與GPU SM之間的對(duì)應(yīng)關(guān)系如圖2所示。

3）每個(gè)Hermite方程組內(nèi)的并行 筆者所討論的Hermite方程組的解法最耗時(shí)的地方為消元操作，可以利用GPU Thread并行進(jìn)行消元操作。一個(gè)GPU Block包括多個(gè)Warp（32個(gè)GPU Thread），一個(gè)Warp負(fù)責(zé)一列的消元，多個(gè)Warp并行消去多列；一個(gè)Warp內(nèi)一個(gè)Thread負(fù)責(zé)一列中一行的消元，32個(gè)Thread實(shí)現(xiàn)并行多行消元。每一個(gè) Warp的32個(gè)Thread運(yùn)行于SM內(nèi)的GPU Core（一個(gè)SM包括32個(gè)GPU Core）內(nèi)，系數(shù)矩陣的行列與GPU Core并行對(duì)應(yīng)關(guān)系如圖3所示。

圖2 Hermite方程組的求解與GPU SM之間的對(duì)應(yīng)關(guān)系

圖3 系數(shù)矩陣的行列與GPU Core并行對(duì)應(yīng)關(guān)系

2．2 CUDA線程模型設(shè)計(jì)

線程模型是用來(lái)明確CUDA程序內(nèi)核的執(zhí)行配置，根據(jù)GPU硬件資源，如寄存器個(gè)數(shù)、共享內(nèi)存大小等，來(lái)定義網(wǎng)格和線程塊，好的線程模型會(huì)使程序的并發(fā)度最優(yōu)，實(shí)現(xiàn)計(jì)算與訪存之間的相互隱藏，使程序性能達(dá)到最優(yōu)。

1）網(wǎng)格 （Grid） 即如何對(duì)所有線程進(jìn)行分塊，定義線程塊數(shù)和線程塊間的組織方式。筆者把一個(gè)空間窗的一個(gè)頻率當(dāng)成一個(gè)線程塊，定義為dimGrid（N），其中N為一個(gè)方程組集的Hermite方程組個(gè)數(shù)。

2）線程塊 （Block） 即定義一個(gè)線程塊有多少個(gè)線程和線程的組織方式。由于實(shí)際生產(chǎn)中，多個(gè)應(yīng)用的Hermite系數(shù)矩陣的階是不同的，因而GPU Kernel所需要的寄存器數(shù)、共享內(nèi)存數(shù)量也不同。針對(duì)不同的階，為使性能最優(yōu)，每一個(gè)線程塊的線程數(shù)是變化的，其具體定義為dimBlock（M，W），一個(gè)線程塊的總線程數(shù)為M＊W，其中W隨系數(shù)矩陣的階的不同而不同。

3）線程模型 每個(gè)Grid劃分為N個(gè)Block，即計(jì)算一個(gè)方程組集內(nèi)的N個(gè)Hermite方程組；每個(gè)Block求解一個(gè)Hermite方程組，它劃分為W個(gè)Warp。假設(shè)Hermite系數(shù)矩陣為二維矩陣H［n，m］，每一個(gè)Warp負(fù)責(zé)一列的消元，W個(gè)Warp一次并行消去W列，每一個(gè)Warp需要循環(huán)m／W次。一個(gè)Warp內(nèi)每一個(gè)GPU Thread負(fù)責(zé)一列中一行的消元，M個(gè)GPU Thread一次并行消去M行，則每一個(gè)Thread需要循環(huán)n／M次。

2．3 CUDA內(nèi)存使用設(shè)計(jì)

根據(jù)GPU并行算法、數(shù)據(jù)訪問(wèn)特點(diǎn)及GPU內(nèi)存資源特性，選擇不同的內(nèi)存存放不同的數(shù)據(jù)，以達(dá)到性能最優(yōu)。

1）Global memory使用 Fermi GPU是按照Warp方式來(lái)訪存的，為了實(shí)現(xiàn)對(duì)Global Memory合并訪問(wèn)，使其訪存性能達(dá)到最優(yōu)，一個(gè)Warp內(nèi)的32個(gè)Thread應(yīng)同時(shí)訪問(wèn)Global memory內(nèi)的連續(xù)內(nèi)存。由于系數(shù)矩陣H是按照列優(yōu)先方式存放的，每一列內(nèi)元素?cái)?shù)據(jù)是連續(xù)存放的，所以在消元操作時(shí)32個(gè)線程應(yīng)同時(shí)訪問(wèn)同一列內(nèi)的32行數(shù)據(jù)元素，提高訪存性能。

2）Shared memory使用 由于Shared memory為GPU的片上內(nèi)存，訪問(wèn)速度快，對(duì)于一個(gè)Block塊中公共的數(shù)據(jù)，如消元操作時(shí)的公共列元素，可以放入共享內(nèi)存中，將提高訪存性能。

3）L1Cache使用 Fermi GPU提供L1Cache，其為GPU的片上內(nèi)存，由于資源受限，L1Cache與Shared memory大小之和僅為64K，可以把L1Cache動(dòng)態(tài)配置為16K或48K2種方式，這樣可以進(jìn)一步提高訪存性能。

3 性能測(cè)試

測(cè)試環(huán)境包括硬件環(huán)境 （CPU：Intel? CoreI i7CPU 920＠2．67GHZ；內(nèi)存：2．67GHZ，8GB；GPU：Tesla C2050）和軟件環(huán)境 （OS：64位Linux RedHat 4．5）。選取某實(shí)際工業(yè)生產(chǎn)中產(chǎn)生的 Hermite方程組，每個(gè)方程組集包含1000個(gè)Hermite方程組，每個(gè)Hermite系數(shù)矩陣的階分別為350、700、1700，對(duì)這3種Hermite方程組的求解分別運(yùn)行單線程串行程序和GPU并行程序進(jìn)行性能比較測(cè)試。為了保證測(cè)試性能結(jié)果的穩(wěn)定性，每次測(cè)試10個(gè)方程組集，然后取其平均時(shí)間，并且每個(gè)測(cè)試重復(fù)進(jìn)行3次，測(cè)得CPU單線程串行程序平均每個(gè)方程組集計(jì)算時(shí)間 （CPU＿Aver＿Time）和GPU并行程序平均每個(gè)方程組集計(jì)算時(shí)間 （GPU＿Ave＿Time），性能對(duì)比如表1所示。從表1不難看出：①GPU并行程序的性能較CPU單線程串行程序提升明顯，加速比達(dá)到18．55～26．35倍；②由于Hermite方程組的求解計(jì)算量隨著系數(shù)矩陣的階數(shù)增大而增加，筆者設(shè)計(jì)的GPU并行算法性能加速比越來(lái)越高，表明系數(shù)矩陣階越高，Hermite方程組的求解計(jì)算越密集，采用GPU并行算法效果將更好。

表1 2種算法性能結(jié)果比對(duì)表

4 結(jié) 語(yǔ)

針對(duì)Hermite系數(shù)矩陣方程組求解計(jì)算效率低的問(wèn)題，筆者首先對(duì)一種Hermite方程組直接解法進(jìn)行GPU并行可行性分析，再基于GPU利用CUDA技術(shù)對(duì)該解法進(jìn)行并行化設(shè)計(jì)。經(jīng)測(cè)試其性能較原串行算法提升18．55～26．35倍，很好地解決了該直接解法的計(jì)算瓶頸。

［1］張浩，李利軍，林嵐．GPU的通用計(jì)算應(yīng)用研究 ［J］．計(jì)算機(jī)與數(shù)字工程，2005，33（12）：60-62．

［2］郭境峰，蔡偉濤 ．新一代高性能運(yùn)算技術(shù)——CUDA簡(jiǎn)介 ［J］．現(xiàn)代科技，2009，8（6）：29-30．

［3］吳連貴，易瑜，李肯立 ．基于CUDA的地震數(shù)據(jù)相干體并行算法 ［J］．計(jì)算應(yīng)用，2009（3）：294-296．