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

張 濤 （長江大學信息與數學學院，湖北 荊州434023水資源與水電科學國家重點實驗室 （武漢大學），湖北 武漢430072）廖文軍 （中國石油集團東方地球物理勘探有限責任公司物探技術研究中心，河北 涿州072751）

在工程計算中，經常需要對Hermite方程組進行求解，而且由于需要保證求解的精確性，一般都采用直接法進行求解，在問題規模較大的情況下，該方程組的求解時間開銷很大。如果單純在CPU上進行計算，計算效率往往滿足不了實際生產的需要。隨著GPU技術［1］的出現，利用GPU眾核處理能力，采用CUDA計算［2-3］架構，利用GPU多核并行實現Hermite方程組的求解，將對解決其計算瓶頸提供契機。

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

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

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

2．1 算法設計

為了達到好的加速比，預測算子計算的GPU并行算法應遵循GPU的SIMD體系結構特點，并且充分發揮其眾核計算能力，整個算法采用三級并行模型，其具體設計過程如下：

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

圖1 Hermite方程組集的并行對應關系

2）每個方程組集內多個Hermite方程組并行 每個方程組集包括N個Hermite方程組，由于每個Hermite方程組的求解是獨立的，彼此之間不存在數據依賴性，方程組之間可以并行求解。把每一個Hermite方程組的求解設計為一個GPU Block，即每一個GPU Block負責求解一個Hermite方程組，那么對于一個GPU Kernel而言，它有N個GPU Block，由于GPU Block運行于SM內，每個Hermite方程組的求解與GPU SM之間的對應關系如圖2所示。

3）每個Hermite方程組內的并行 筆者所討論的Hermite方程組的解法最耗時的地方為消元操作，可以利用GPU Thread并行進行消元操作。一個GPU Block包括多個Warp（32個GPU Thread），一個Warp負責一列的消元，多個Warp并行消去多列；一個Warp內一個Thread負責一列中一行的消元，32個Thread實現并行多行消元。每一個 Warp的32個Thread運行于SM內的GPU Core（一個SM包括32個GPU Core）內，系數矩陣的行列與GPU Core并行對應關系如圖3所示。

圖2 Hermite方程組的求解與GPU SM之間的對應關系

圖3 系數矩陣的行列與GPU Core并行對應關系

2．2 CUDA線程模型設計

線程模型是用來明確CUDA程序內核的執行配置，根據GPU硬件資源，如寄存器個數、共享內存大小等，來定義網格和線程塊，好的線程模型會使程序的并發度最優，實現計算與訪存之間的相互隱藏，使程序性能達到最優。

1）網格 （Grid） 即如何對所有線程進行分塊，定義線程塊數和線程塊間的組織方式。筆者把一個空間窗的一個頻率當成一個線程塊，定義為dimGrid（N），其中N為一個方程組集的Hermite方程組個數。

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

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

2．3 CUDA內存使用設計

根據GPU并行算法、數據訪問特點及GPU內存資源特性，選擇不同的內存存放不同的數據，以達到性能最優。

1）Global memory使用 Fermi GPU是按照Warp方式來訪存的，為了實現對Global Memory合并訪問，使其訪存性能達到最優，一個Warp內的32個Thread應同時訪問Global memory內的連續內存。由于系數矩陣H是按照列優先方式存放的，每一列內元素數據是連續存放的，所以在消元操作時32個線程應同時訪問同一列內的32行數據元素，提高訪存性能。

2）Shared memory使用 由于Shared memory為GPU的片上內存，訪問速度快，對于一個Block塊中公共的數據，如消元操作時的公共列元素，可以放入共享內存中，將提高訪存性能。

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

3 性能測試

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

表1 2種算法性能結果比對表

4 結 語

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

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

［2］郭境峰，蔡偉濤 ．新一代高性能運算技術——CUDA簡介 ［J］．現代科技，2009，8（6）：29-30．

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