蘇 輝，邱夏青，馬文鵬

(信陽師范學院 a.網(wǎng)絡信息與計算中心;b. 計算機與信息技術(shù)學院,河南 信陽 464000)

0 引言

對于很多大規(guī)模數(shù)值模擬來說，稀疏線性系統(tǒng)的求解占用了大部分的計算時間和資源.隨著科學與工程問題的規(guī)模和復雜程度的提高，對稀疏性系統(tǒng)的求解規(guī)模和速度提出了更高的要求，利用高性能計算技術(shù)是當前解決大規(guī)模稀疏線性系統(tǒng)求解的重要途徑[1].在圖形處理器(GPU)硬件架構(gòu)的迅猛發(fā)展和NVIDIA CUDA[2]并行計算平臺的不斷完善之下，GPU的可編程能力也在不斷增強.隨著CUDA運算平臺的推出，GPU的并行計算已經(jīng)滲透到各個學科，如在生物科學研究領域中蛋白質(zhì)與基因序列串的匹配等；在計算金融方面利用GPU加速債券定價、風險分析以及算法交易[3]；在化學計算方面進行量子化學的計算加速；在數(shù)據(jù)挖掘、分析學方面的信息高效排序等，其本質(zhì)就是利用GPU多線程的特點來輔助CPU計算，以此來提高計算機整體的運算速度.

目前，針對GPU的有限元分析在求解大規(guī)模稀疏線性方程組的研究較多，主要是由于在求解線性方程組時數(shù)據(jù)量大并且計算比較集中，因此主要研究稀疏矩陣的數(shù)據(jù)結(jié)構(gòu)和稀疏矩陣向量乘(Sparse Matrix Vector Multiplication, SpMV)的算法[4].同時，很多研究學者也研究了針對GPU的基于CUDA平臺的共軛梯度(Conjugate Gradient Method, CG)算法及各種改進算法，如BICG、BiCGSTAB算法.

本文所述算法是采用行壓縮存儲(Compressed Row Storage, CSR)模式存儲大規(guī)模的稀疏矩陣，在稀疏矩陣向量乘(SpMV)的基礎上采用預處理共軛梯度(Preconditioned Conjugate Gradient Method, PCG)算法，對大規(guī)模稀疏矩陣組成的線性方程組迭代求解.采用CUDA和Matlab混合編程，實現(xiàn)有限元方程的數(shù)值求解，通過多線程求解來實現(xiàn)GPU加速，算法性能有明顯的優(yōu)越性.

1 求解向量方程組的方法與技術(shù)

1.1 有限元方法

有限元方法(Finite Element Analysis, FEA)是一種高效常用的數(shù)值計算方法，由于有限元方法計算精度較高，具有適應大多數(shù)復雜情況的能力，使其成為力學工程分析問題的有效解決手段.化整為零、積整為零是有限元方法的基本思想，這種解決問題的思想與并行計算中的分治思想相協(xié)調(diào)[5].因此，對于規(guī)模較大的物理工程問題，可以借助于有限元方法將大規(guī)模問題分解，使用GPU多線程結(jié)構(gòu)將運算效率低下的串行方法改進成為并行分析方法，使得計算機的執(zhí)行效率變得更高.

1.2 稀疏矩陣行壓縮存儲

矩陣的類型有很多，按照存儲數(shù)據(jù)的稠密程度分類，可以分為稠密矩陣和稀疏矩陣.在稀疏矩陣中，其存儲的零元素較多，而其存儲的有效值較少；而稠密矩陣所存儲的有效值較多，非有效值較少.因此本文采用一種新型的存儲方式，即將矩陣信息壓縮存儲，也是將矩陣中的零元素都剔除，只存儲矩陣中的有效值，對于內(nèi)存空間的利用將會更加充分.

1.3 稀疏矩陣向量乘算法

SpMV算法是用來求解矩陣和向量相乘的函數(shù)，通過SpMV算法，可以求得矩陣與向量相乘的結(jié)果，此算法主要是為后面的預處理共軛梯度(PCG)算法做鋪墊、以供給其調(diào)用來求解線性方程組.

1.4 預處理共軛梯度算法

據(jù)科學統(tǒng)計，預條件共軛梯度(PCG)算法的迭代速度隨著方程的階數(shù)的增加而下降，從而可以推斷出：對于高階稀疏矩陣的線性方程組，無論其階數(shù)有多高，PCG算法迭代的次數(shù)都是可以被接受的[6]，這種方法使求解大型或者超大型方程組成為可能.同時，對于PCG算法，還允許在計算的過程中保留矩陣的稀疏性，允許矩陣采取更加高效的方式來存儲它.此特點使得計算機在處理拓撲形狀復雜的稀疏矩陣的時候節(jié)省大量的內(nèi)存空間.同時，PCG算法適合于求解大型對稱的矩陣，而本次實驗所測試的矩陣正好符合這一特點.

1.5 Matlab軟件的優(yōu)勢

Matlab是一款通常用于數(shù)據(jù)可視化、數(shù)據(jù)計算、算法開發(fā)以及數(shù)據(jù)分析數(shù)學軟件[4].同時，Matlab可以將數(shù)值計算、科學數(shù)據(jù)可視化、矩陣運算以及圖形圖像仿真都集成在一個視圖中，使得工程設計和科學研究中的有效數(shù)值計算問題能夠得到更加有效的解決.Matlab最大的優(yōu)點就是能夠使數(shù)值計算結(jié)果和程序?qū)崿F(xiàn)過程可視化，其強大的圖形圖像處理能力使得靜態(tài)的數(shù)據(jù)轉(zhuǎn)化為生動的圖像，使得數(shù)學規(guī)律更加直觀明了，使人從枯燥無味的數(shù)學計算中解脫出來.其次，在Matlab平臺可以調(diào)用CUDA函數(shù)，實現(xiàn)復雜計算調(diào)用GPU中的核函數(shù)完成，從而實現(xiàn)整體運算加速[7].

2 二維拉普拉斯方程的有限元求解

2.1 求解應用方程

理想流體流動的連續(xù)方程為：

(1)

(2)

對其可描述為理想流體流動的速度勢方程.該方程的形式與拉普拉斯方程的形式一致.

2.2 插值函數(shù)及相關計算

三角形線性單元插值函數(shù)為：

Φi=ai+bix+ciy,i=1,2,3,

(3)

其中

(4)

在計算公式(3)和(4)時，當i=1時，j=2，k=3；當i=2時，j=3，k=1；當i=3時，j=1，k=2.其中Δ為單元面積，可計算如下：

(5)

2.3 加權(quán)余量方程

在伽遼金有限元方法中，權(quán)函數(shù)等于插值函數(shù).式(2)與插值函數(shù)相乘后再整個計算區(qū)域內(nèi)積分，可以得到：

(6)

式(6)的左邊用格林公式進一步推算和展開，得到系統(tǒng)(2)的加權(quán)余量方程：

(7)

2.4 單元方程的建立和總體方程的組合

單元方程的形式如下：

(8)

但是，由于存在第二類邊界條件，則在建立單元方程時，應該分為內(nèi)部單元和邊界單元，這里不再敘述.

總體方程的組合分兩步完成，首先裝配CSR格式的剛度矩陣；最終裝配目標剛度矩陣和裝配好邊界向量.

2.5 求解線性方程組

在裝配完剛度矩陣和求解出邊界向量之后，就可以使用預處理共軛梯度(PCG)求解線性方程組了.在求解線性方程組時，需要先求解出CSR存儲格式的系數(shù)矩陣和由邊界條件計算的向量，然后再調(diào)用PCG算法實現(xiàn)線性方程組的求解.由此計算出來的時間可以進行與后面的并行效率對比分析.

3 基于GPU的二維拉普拉斯方程并行加速

網(wǎng)格是GPU線程的組織方式，若干線程塊組成了一個網(wǎng)格，最多512個線程能組成一個線程塊，屬于同一個線程塊中線程的指令地址相同，這些線程具有很強的可并行性，同時還能夠通過Share memory和barrier實現(xiàn)塊內(nèi)通信，可以形成細粒度并行[8].在線程中通過粗細粒度的雙層并行可以實現(xiàn)GPU運算的多方面加速.

由于GPU多線程的特點，使其在大規(guī)模并行計算方面具有很大的優(yōu)勢，所以對于此次設計中的裝配剛度矩陣和預處理共軛梯度算法來說，采用在GPU端細粒度并行方式可以實現(xiàn)計算提速[9].

3.1 以格點為主序裝配非零元素

在本次剛度矩陣的填充中，填充的剛度矩陣格式為CSR模式，直接將每行的數(shù)值填入value數(shù)組，填充位置的確定需要根據(jù)index數(shù)組確定每行非零元素在value數(shù)組中的起始地址和每行的非零元素的序號來確定.最終的value裝配形式是按照每行的數(shù)據(jù)段從前到后順序執(zhí)行.

3.2 基于Matlab平臺的GPU加速

3.2.1 并行裝配剛度矩陣

在Matlab中調(diào)用MEX文件極為方便，因為其調(diào)用方式和Matlab的內(nèi)建函數(shù)幾乎完全相同，只需要在命令窗口輸入函數(shù)名和對應的參數(shù)即可.因此，本文在Matlab平臺上編寫一個接口函數(shù)mexFunction，其格式如下：

void mexFunction(int nlhs,mxArray*plhs[],int nrhs,const mxArray*prhs[])

nlhs:輸出參數(shù)的數(shù)目

plhs:指向輸出參數(shù)的指針

nrhs:輸入?yún)?shù)的數(shù)目

prhs:指向輸入?yún)?shù)的指針

mexFunction函數(shù)只是提供了一個Matlab和CUDA的接口，其實現(xiàn)的功能就是獲取參數(shù)指針，并且調(diào)用FillValueArray函數(shù)來啟動核函數(shù)，實現(xiàn)多線程的計算[9].

在FillValueArray.cu文件中，面臨的主要任務就是用C語言實現(xiàn)對value數(shù)組的按格點填充.對于每個格點分配一個線程，使其并行完成對value數(shù)組的填充.在GPU的每個thread中，執(zhí)行的操作都相同，區(qū)別是每個線程執(zhí)行時的數(shù)據(jù)不同.

實現(xiàn)CUDA和Matlab的混合編程，把.cu文件用CUDA的編譯器nvcc編譯成.obj文件，作為一個庫函數(shù)供.cpp調(diào)用，然后.cpp與.obj文件共同用MEX接口編譯成供Matlab調(diào)用的函數(shù).在此部分最核心的部分就是用CUDA語言編寫.cu文件，實現(xiàn)對CSR存儲格式的系數(shù)矩陣的在GPU上并行填充.完成在device端分配顯存、數(shù)值計算、寫回host端和顯存釋放等操作.

3.2.2 并行求解線性方程組

在核函數(shù)中，首先，為要使用的變量開辟臨時顯存空間，聲明調(diào)cusparseDcsrmv函數(shù)中的各個參數(shù)類型;其次，使用CUBLAS庫中的一些固有的函數(shù)進行向量的操作;最后，分配線程執(zhí)行函數(shù)，將device端計算的結(jié)果寫回host端，釋放在顯存中開辟的臨時空間.

對于CUDA中的可執(zhí)行文件，主要是獲取PCG函數(shù)的參數(shù)，獲取變量的指針，在顯存中開辟顯存空間，然后再調(diào)用CUDA編寫的.cu函數(shù)，完成并行操作.

對于PCG算法的并行，也是在Matlab平臺用CUDA語言編寫.cu文件與.h文件用nvcc命令編寫成.obj文件，執(zhí)行在GPU上[10].然后.cpp與.obj文件共同用mex接口編譯成供Matlab調(diào)用的函數(shù).在Matlab端對生成的函數(shù)調(diào)用，實現(xiàn)在Matlab平臺實現(xiàn)對大規(guī)模線性方程組的并行求解操作.

4 計算結(jié)果及分析

4.1 測試環(huán)境與測試算例

本次實驗是在Windows7，64位操作系統(tǒng)上運行，處理器的型號為Intel(R)Core(TM)i7-4790，8核，主頻為3.60 GHz，安裝內(nèi)存(RAM)為8.00 GB，GPU型號為NVIDIA GeForce GT 705.本次實驗運行的軟件運行平臺是Matlab 2015a、Visual Stdio 2010和CUDA8.0版本.

本次實驗的測試算例是以流水通過兩平板之間的速度勢變化為物理背景.測試的初始條件為寬H為5 m，高L為10 m，水流速度v為2 m/s，其中水平方向分割段數(shù)Nx為500，垂直方向分割段數(shù)Ny為400，則通過計算可以得出劃分的網(wǎng)格數(shù)為40萬，網(wǎng)格格點數(shù)目為20萬.

4.2 實驗結(jié)果與效率分析

通過計算大規(guī)模線性方程組為載體，來對比通過GPU加速前后程序的執(zhí)行時間，計算GPU在大規(guī)模工程運算中的提速效率.下面通過測試各個程序段的時間來對比使用GPU加速的效果.

1)實驗的并行提速效果

裝配剛度矩陣(fillvalue)提速：3.222÷0.401=8.035倍；

求解線性方程組(PCG)提速：2.810÷0.411=6.84倍；

算法總體提速：6.957÷1.737=4.01倍.

2)并行加速效率分析

通過上述數(shù)據(jù)，可以看出在裝配剛度矩陣中使用GPU中的線程計算可以提速8倍，在求解線性方程組時使用GPU中線程可以提速7倍，對于求解整體的方程提速為4倍，對于剛度矩陣的填充和線性方程組的求解來說,效果比較理想.

5 總結(jié)

本文是在Matlab平臺上實現(xiàn)Matlab和CUDA的混合編程，借助于GPU的多線程結(jié)構(gòu)，使得GPU發(fā)揮細粒度方面的并行優(yōu)勢，實現(xiàn)整體算法的并行加速.綜上所述，本次設計中通過在GPU中的兩次提速，在大規(guī)模數(shù)值計算問題中，能節(jié)省較多的計算時間，實驗效果比較理想，達到在大規(guī)模數(shù)值計算的GPU提速.