非靜壓海洋數(shù)值模式加速計(jì)算的CUDA實(shí)現(xiàn)*

王春暉，苗春葆，沈 飆**

（1.中國(guó)海洋大學(xué)物理海洋實(shí)驗(yàn)室，山東 青島266100；2.國(guó)家海洋局北海環(huán)境監(jiān)測(cè)中心海洋溢油監(jiān)別與損害評(píng)估技術(shù)重點(diǎn)實(shí)驗(yàn)室，山東 青島266033）

自從NVIDIA提出CUDA框架以后，使用GPU進(jìn)行通用計(jì)算變得相當(dāng)簡(jiǎn)單［1］。CUDA以擴(kuò)展的C語(yǔ)言為基礎(chǔ)，取代了以前基于GPU的程序設(shè)計(jì)模式，采用CUDA技術(shù)，程序員不需要將計(jì)算轉(zhuǎn)化成圖形函數(shù)，提高了程序開(kāi)發(fā)的靈活性［2］。在CUDA編程模型中，GPU被視為1個(gè)協(xié)處理器，通過(guò)大量線(xiàn)程的并行執(zhí)行來(lái)加快計(jì)算速度。在CUDA架構(gòu)下，1個(gè)程序分為兩部分：主機(jī)（Host）端和設(shè)備（Device）端。主機(jī)端是指在CPU上執(zhí)行的部分，而設(shè)備端是在GPU上執(zhí)行的部分，在GPU端執(zhí)行的程序稱(chēng)為“核”（Kernel）。CUDA程序執(zhí)行時(shí)，驅(qū)動(dòng)程序會(huì)將GPU端代碼編譯成可執(zhí)行程序，并傳送到GPU運(yùn)行。在GPU上會(huì)產(chǎn)生大量的線(xiàn)程，每一個(gè)線(xiàn)程都會(huì)去執(zhí)行核程序，雖然程序有同一份，但由于每個(gè)線(xiàn)程的索引值不同，因而各個(gè)線(xiàn)程就能對(duì)不同的數(shù)據(jù)進(jìn)行計(jì)算［3］。

在自然界發(fā)生的很多流體運(yùn)動(dòng)中，表面水波是最容易被觀察到的，由于地球表面大約三分之二的面積被水覆蓋，因此對(duì)水波的研究就具有重要的實(shí)際意義。對(duì)淺水波進(jìn)行數(shù)值模擬時(shí)，若垂向壓力一直滿(mǎn)足靜壓假設(shè)，則Navier－Stokes方程可以得到簡(jiǎn)化，不用求解泊松方程。但在自由表面有較大變化、地形突變、或垂向摻混比較嚴(yán)重的區(qū)域，采用這種假設(shè)的模擬結(jié)果會(huì)有很大誤差；若使用非靜壓模式進(jìn)行模擬，計(jì)算量會(huì)有明顯的增加，難以在較短的時(shí)間內(nèi)完成計(jì)算任務(wù)。因此，本項(xiàng)目使用CUDA實(shí)現(xiàn)對(duì)表面水波數(shù)值模擬的GPU加速，實(shí)現(xiàn)非靜壓表面水波的較快的數(shù)值模擬。

與CPU串行程序相比，可以在不犧牲模擬精度的前提下（相對(duì)誤差不超過(guò)2×10－3）大大提高計(jì)算效率；計(jì)算速度的顯著提高展現(xiàn)了海洋數(shù)值模式應(yīng)用GPU進(jìn)行加速的廣闊前景。

1 數(shù)值模型

1.1 基本方程

流體作為物質(zhì)的一種形態(tài)，遵循自然界中關(guān)于物質(zhì)運(yùn)動(dòng)的普遍規(guī)律，例如質(zhì)量守恒、動(dòng)量守恒和能量守恒等。為簡(jiǎn)單起見(jiàn)，本項(xiàng)目采用描述海水運(yùn)動(dòng)的二維垂向切片模型［4］，其控制方程為：

其中：u代表水平方向的流速；w代表垂向流速；ρo代表海水密度；P代表壓強(qiáng)；η表示海面水位起伏。

1.2 差分格式

參 考 Jochen Kampf［4］， 對(duì) 上 述 控 制 方 程 在Arakawa C網(wǎng)格上進(jìn)行差分計(jì)算（見(jiàn)圖1），垂向采用z坐標(biāo)。

動(dòng)壓強(qiáng)P可以分為兩部分：

圖1 Arakawa C網(wǎng)格中各變量位置分布Fig.1 Arakawa C－grid for a vertical ocean slice

其中：p代表平靜海平面的靜壓強(qiáng)；q由兩部分組成，一部分是海平面傾斜造成的壓強(qiáng)變化，第二部分是非靜壓效應(yīng)造成的壓強(qiáng)變化。若海水密度是一個(gè)常量，方程（5）可以簡(jiǎn)寫(xiě)為：P＝q。把q分解為當(dāng)前時(shí)刻n的壓強(qiáng)加上下一時(shí)刻n＋1的變化值，即：

相應(yīng)的，方程（1）～（4）的求解過(guò)程可以分為2步：第一步，在時(shí)刻n顯式求出初始猜測(cè)流速u(mài)＊跟w＊；第二步，隱式求出Δq。

流速的初始猜測(cè)值可以由時(shí)間前差方案獲得，即：

其中：i跟k分別代表網(wǎng)格的水平和垂向索引。動(dòng)量方程的有限差分格式為：

將方程（8）、（9）代入方程（3），兩邊同時(shí)乘以ΔzΔx，可以得到：

該方程在數(shù)學(xué)上被稱(chēng)為泊松方程，其中：

圖2給出了這些系數(shù)在Arakawa－C網(wǎng)格中的位置關(guān)系。方程右邊包含了初始猜測(cè)流速場(chǎng)（u＊，w＊）的散度，即：

方程（10）的解意味著流場(chǎng)沒(méi)有散度，即滿(mǎn)足流場(chǎng)的不可壓縮條件，新的壓強(qiáng)場(chǎng)則可以表示為：

圖2 ae，aw，at，ab 的位置關(guān)系圖Fig.2 Location used to define the coefficients ae，aw ，atand ab.

1.3 泊松方程的求解算法

方程（10）可以采用高斯 －賽德?tīng)柕℅auss－Seidel，簡(jiǎn)稱(chēng) G－S方法）或超松弛迭代法（Successive Over－Relaxation，簡(jiǎn)稱(chēng)SOR方法）進(jìn)行求解［4］。用方程表示如下：

其中：r＝0，1，2，…代表迭代次數(shù)；上標(biāo)l由Δq是否更新決定，即l＝r＋1或l＝r；ω的取值與采用哪種方法有關(guān)：若采用G－S方法，則ω取值為1，若采用SOR方法，則ω的典型取值范圍為1.2～1.4。SOR方法與G－S方法相比其優(yōu)勢(shì)在于：通過(guò)選擇松弛因子ω，可以使得迭代過(guò)程收斂較快。

G－S方法及SOR方法中Δq的初始值可以設(shè)置為0，但是如果初始值采取前一步的取值，會(huì)使計(jì)算迭代次數(shù)明顯減少，提高計(jì)算效率。即：

海表面壓力在每次進(jìn)行迭代時(shí)都需要進(jìn)行計(jì)算。具體計(jì)算過(guò)程如下所

（1）對(duì)流速進(jìn)行更新

（2）在垂向上積分水平流速

（3）忽略海平面變化引起的水體厚度的輕微變化，海表壓強(qiáng)的變化可用下面公式表示：

重復(fù)G－S方法或SOR方法，直到相鄰2次迭代的差滿(mǎn)足要求，即：

其中：ε是2次迭代之間壓強(qiáng)的允許誤差。根據(jù)方程（11）可知，此時(shí)意味著：

其中：上標(biāo)為r＋1的值即為G－S方法或SOR方法的解。圖3給出了G－S／SOR方法計(jì)算流程圖。雖然SOR方法收斂速度更快，但其各個(gè)網(wǎng)格的計(jì)算具有嚴(yán)格的順序性，無(wú)法實(shí)現(xiàn)并行計(jì)算，因此這里采用G－S算法對(duì)方程（10）進(jìn)行求解。

圖3 G－S／SOR方法計(jì)算流程圖Fig.3 Flow chart of the G－S／SOR method

1.4 串行程序的實(shí)現(xiàn)

用C語(yǔ)言實(shí)現(xiàn)以上求解算法，計(jì)算流程圖如下：

計(jì)算過(guò)程中所用到的主要函數(shù)及其功能為：

（1）init函數(shù)：初始化變量的值，為變量賦初值；

（2）dyn函數(shù)：動(dòng)力過(guò)程的計(jì)算，計(jì)算下一個(gè)時(shí)間步各變量的值；

（3）store＿surf＿dp函數(shù)：將海表面壓強(qiáng)變化保存到另一個(gè)數(shù)組；

圖4 C語(yǔ)言程序計(jì)算流程圖Fig.4 Flow chart of the C language procedure code

2 基于CUDA的并行程序的實(shí)現(xiàn)與優(yōu)化

通過(guò)對(duì)各子程序進(jìn)行并發(fā)性（concurrency）分析可知，每個(gè)子程序均可以采用并行算法進(jìn)行計(jì)算加速，但由于init函數(shù)計(jì)算量較小，且僅調(diào)用一次，因此不需要采用GPU進(jìn)行加速；但對(duì)于動(dòng)力方程計(jì)算dyn的7個(gè)函數(shù)，計(jì)算量均較大，因而可以將這些程序改為kernel函數(shù)，用GPU進(jìn)行計(jì)算。程序的優(yōu)化可以分為以下5個(gè)部分：

2.1 全局存儲(chǔ)訪(fǎng)問(wèn)的程序優(yōu)化

全局存儲(chǔ)器（global memory）即普通的顯存，GPU與CPU都可以進(jìn)行讀寫(xiě)訪(fǎng)問(wèn)，網(wǎng)格中的任意線(xiàn)程都可以對(duì)全局存儲(chǔ)器的任意位置進(jìn)行讀寫(xiě)。全局存儲(chǔ)器具有較高的訪(fǎng)問(wèn)延遲，比較容易成為性能瓶頸［3］。能否滿(mǎn)足合并訪(fǎng)問(wèn)在很多情況下會(huì)使CUDA程序的速度產(chǎn)生高達(dá)一個(gè)數(shù)量級(jí)的差異。CUDA中實(shí)際執(zhí)行單元是以warp為單位的，在Fermi架構(gòu)中，1個(gè)warp由32個(gè)線(xiàn)程組成，每次GPU調(diào)度1個(gè)warp里的32個(gè)線(xiàn)程執(zhí)行同一命令，同在1個(gè)warp的線(xiàn)程以不同數(shù)據(jù)資源執(zhí)行相同的指令。如果同一個(gè)warp中的線(xiàn)程對(duì)全局存儲(chǔ)的訪(fǎng)問(wèn)滿(mǎn)足合并訪(fǎng)問(wèn)條件，則可以達(dá)到對(duì)存儲(chǔ)帶寬的完全充分利用［5］。

C語(yǔ)言的二維數(shù)組是行序存放的，使用CUDA訪(fǎng)問(wèn)global memory時(shí)，要使訪(fǎng)問(wèn)效率最高，應(yīng)該從256字節(jié)對(duì)齊的地址（即addr＝0，256，512…）開(kāi)始進(jìn)行連續(xù)訪(fǎng)問(wèn)，因此，為提高內(nèi)存訪(fǎng)問(wèn)的效率，有2種方法：第1種是使用cudaMallocPitch函數(shù)進(jìn)行自動(dòng)對(duì)齊，該函數(shù)分配的內(nèi)存中，數(shù)組每一行中第1個(gè)元素開(kāi)始的地址都保證是對(duì)齊的。由于每行有多少個(gè)數(shù)據(jù)并不確定，即每行的元素不一定是256的倍數(shù)，因此，為了保證數(shù)組每一行中第1個(gè)元素開(kāi)始的地址是對(duì)齊的，cudaMallocPitch函數(shù)在分配內(nèi)存時(shí)，每行都會(huì)多分配出一些字節(jié)，從而保證每行中的元素加上多分配的字節(jié)恰好是256的倍數(shù)（對(duì)齊）；第2種方法是手動(dòng)對(duì)齊，即定義數(shù)組時(shí)保證是數(shù)組第二維是256的倍數(shù)，在本項(xiàng)目中作者采用第二種方法，保證數(shù)組的第二維，即nx＋2為256的倍數(shù)，此時(shí)有效帶寬達(dá)到最大。

2.2 使用共享存儲(chǔ)器進(jìn)行程序優(yōu)化

共享存儲(chǔ)器（shared memory）位于GPU片內(nèi)，是一塊可讀寫(xiě)高速存儲(chǔ)器，能被同一block中的所有線(xiàn)程訪(fǎng)問(wèn)。對(duì)共享存儲(chǔ)器變量的聲明通過(guò)＿＿shared＿＿來(lái)實(shí)現(xiàn)。在并行訪(fǎng)問(wèn)時(shí)為了能夠獲得高帶寬，其被劃分為能被同時(shí)訪(fǎng)問(wèn)且大小相等的存儲(chǔ)器模塊（bank）。不同的bank可以互不干擾的同時(shí)進(jìn)行工作，因此可以同時(shí)訪(fǎng)問(wèn)位于n個(gè)bank上的n個(gè)地址，有效帶寬為僅有一個(gè)bank時(shí)的n倍。但假如warp請(qǐng)求訪(fǎng)問(wèn)的很多個(gè)地址位于同一個(gè)bank中，此時(shí)bank在1個(gè)時(shí)刻不能響應(yīng)多個(gè)請(qǐng)求，因此這些請(qǐng)求就必須被串行地完成，即出現(xiàn)存儲(chǔ)器沖突，這時(shí)硬件會(huì)將造成存儲(chǔ)器沖突的訪(fǎng)存請(qǐng)求進(jìn)行劃分，分成幾次不存在沖突的獨(dú)立請(qǐng)求，同時(shí)有效帶寬會(huì)降低數(shù)倍，降低的倍數(shù)即拆分后得到的不存在沖突情況的請(qǐng)求個(gè)數(shù)請(qǐng)求［6］。

計(jì)算核函數(shù)void calc＿dp和calc＿u＿w時(shí)多次用到Δq，因此可以在程序中將數(shù)組Δq加載到共享存儲(chǔ)器，提高有效帶寬，但由于程序中用到的頻率并不是太高，有效帶寬提高不大，計(jì)算時(shí)間只是略有減少。

2.3 使用常數(shù)存儲(chǔ)器進(jìn)行程序優(yōu)化

常數(shù)存儲(chǔ)器（constant memory）是只讀的地址空間，其數(shù)據(jù)位于顯存，且擁有緩存機(jī)制，用以節(jié)約帶寬，加快訪(fǎng)問(wèn)速度［7］。常數(shù)存儲(chǔ)器與全局存儲(chǔ)器相比，在讀取或者存入數(shù)據(jù)時(shí)都具有很小的延遲、很高的速度。聲明常數(shù)存儲(chǔ)器通過(guò)＿＿constant＿＿來(lái)實(shí)現(xiàn)。

在本項(xiàng)目中，除邊界處外，方程（11）中的ae，aw，at，ab在大多數(shù)網(wǎng)格中都是相同的，因此將其由數(shù)組改為常數(shù)，然后使用常數(shù)存儲(chǔ)器進(jìn)行優(yōu)化。該方法有效減少了對(duì)全局存儲(chǔ)的訪(fǎng)問(wèn)，大大提高了計(jì)算效率。

2.4 使用多個(gè)流進(jìn)行程序優(yōu)化

程序管理并發(fā)需要通過(guò)流（stream）來(lái)實(shí)現(xiàn)，流是按順序執(zhí)行的一系列命令，但不同流之間則沒(méi)有這一順序性，在支持多個(gè)流并行的設(shè)備上是并行執(zhí)行的。這樣，可以使一個(gè)流的計(jì)算與另一個(gè)流的數(shù)據(jù)傳輸同時(shí)進(jìn)行，從而提高GPU中資源的利用率。

流的定義方法是創(chuàng)建一個(gè)cudaStream＿t對(duì)象，并在啟動(dòng)內(nèi)核和進(jìn)行memcpy時(shí)將該對(duì)象作為參數(shù)傳入，參數(shù)相同的屬于同一個(gè)流，參數(shù)不同的屬于不同的流。執(zhí)行參數(shù)沒(méi)有流參數(shù)，或使用0作為流參數(shù)時(shí)，則使用默認(rèn)的流。當(dāng)使用默認(rèn)的流進(jìn)行任何內(nèi)核啟動(dòng)、內(nèi)存設(shè)置或拷貝函數(shù)時(shí)，只有在之前設(shè)備上所有的操作均已完成后才會(huì)開(kāi)始。

在本項(xiàng)目中，設(shè)置參數(shù)isconverge作為判斷循環(huán)是否終止的條件，若滿(mǎn)足方程（16），則isconverge＝1，迭代終止；相反，若不滿(mǎn)足方程（16），則isconverge＝0，迭代過(guò)程繼續(xù)。由于每次迭代是否收斂的判斷需要在設(shè)備端進(jìn)行，而是否退出整個(gè)迭代循環(huán)的操作需要在主機(jī)端進(jìn)行，所以在每次迭代時(shí)都會(huì)把isconverge的值由設(shè)備端拷貝到主機(jī)端，這會(huì)花費(fèi)一定的時(shí)間。而isconverge的拷貝是可以與核函數(shù)執(zhí)行同時(shí)進(jìn)行的，因此，這里創(chuàng)建2個(gè)流，將核函數(shù)調(diào)用放在stream0中，而設(shè)備端到主機(jī)端的數(shù)據(jù)拷貝放在stream1中，從而使用計(jì)算時(shí)間掩蓋了主機(jī)設(shè)備通信時(shí)間，顯著提高了程序的效率。

2.5 網(wǎng)格及線(xiàn)程塊的維度設(shè)計(jì)

CUDA在執(zhí)行內(nèi)核函數(shù)時(shí)，block會(huì)被分配到SM（流多處理器）中執(zhí)行。對(duì)于block的尺寸應(yīng)該優(yōu)先考慮，而對(duì)于grid的尺寸，應(yīng)該根據(jù)問(wèn)題的規(guī)模予以分配。

block的尺寸與數(shù)據(jù)劃分緊密相關(guān)。較小的block使用資源較少，一般在1個(gè)SM中能夠有更多的active block；而block較大時(shí)，意味著可以進(jìn)行通信的線(xiàn)程更多，此時(shí)指令流效率也更高。為了有效利用執(zhí)行單元，應(yīng)該讓每個(gè)block中的線(xiàn)程數(shù)量是32的整數(shù)倍，且最好為128～256之間，因此在本項(xiàng)目核函數(shù)中block的維度，均取為256。

為了提高訪(fǎng)問(wèn)全局存儲(chǔ)器和共享存儲(chǔ)器的效率，blockDim.x在Telsa架構(gòu)中應(yīng)設(shè)計(jì)為16的整數(shù)倍，在Fermi架構(gòu)中，blockDim.x 應(yīng)為32的整數(shù)倍［8］。在block的尺寸及維度確定以后，就可以按照實(shí)際問(wèn)題的規(guī)模對(duì)grid中各個(gè)維度的block數(shù)量進(jìn)行確定。一般的，grid在某個(gè)維度上的block數(shù)量為：

grid在某方向的block數(shù)量＝（問(wèn)題在該方向上的尺寸－1）／每個(gè)block在該方向的尺寸＋1。

因?yàn)檎麛?shù)除法只會(huì)取結(jié)果的整數(shù)部分，這樣可能使問(wèn)題的邊界得不到處理，引起錯(cuò)誤，因此需要在后面加上1，以保證block數(shù)量滿(mǎn)足實(shí)際要求；若問(wèn)題在某方向的尺寸恰好為每個(gè)block在該方向尺寸的整數(shù)倍，此時(shí)后面＋1反而會(huì)導(dǎo)致block數(shù)量多于實(shí)際要求的數(shù)量，因此需要在前面的括號(hào)內(nèi)減去1。

3 模擬結(jié)果及性能分析

3.1 測(cè)試平臺(tái)

進(jìn)行測(cè)試時(shí)，所有不必要程序全部關(guān)閉，測(cè)試的軟硬件環(huán)境配置如下：

出自《圣經(jīng)》。上帝對(duì)人類(lèi)所犯下的罪孽非常憂(yōu)傷，決定用洪水消滅人類(lèi)。諾亞是個(gè)正直的人，上帝吩咐他造船避災(zāi)。經(jīng)過(guò)40個(gè)晝夜的洪水，除諾亞一家和部分動(dòng)物外，其他生物都被洪水吞沒(méi)。

CPU AMD Phenom（tm）II x4 965 內(nèi)存：4G

GPU NVIDIA Geforce GTX 560Ti 顯存：1G

主板 技嘉GA－MA790GP－UD3H

主板芯片組 AMD 790GX

總線(xiàn)接口標(biāo)準(zhǔn) PCI－Expressx16 2.0

操作系統(tǒng) Windows HPC Server 2008R2 64位

3.2 模擬結(jié)果

模式中水平網(wǎng)格大小為5m，垂向網(wǎng)格大小為2m，時(shí)間步長(zhǎng)為0.05s。圖5為模擬的各時(shí)刻的水位起伏，可以看到，隨著時(shí)間的推移，重力波從左逐漸向右傳播，重力波造成的水位起伏在海面處最大，隨著深度的增加逐漸減小，到海底處為0。

下面對(duì)摸擬結(jié)果進(jìn)行誤差分析，保留計(jì)算結(jié)果中每個(gè)時(shí)刻、每個(gè)網(wǎng)格dp的值，僅對(duì)dp不為0的點(diǎn)進(jìn)行統(tǒng)計(jì)誤差分析，CUDA程序與CPU串行程序的相對(duì)誤差計(jì)算公式為：

其中：dpcuda代表CUDA程序計(jì)算結(jié)果；dpc代表CPU串行程序的計(jì)算結(jié)果；M代表結(jié)果不為0的網(wǎng)格點(diǎn)總數(shù)。統(tǒng)計(jì)結(jié)果見(jiàn)表1。由表1可以看出，網(wǎng)格數(shù)較小時(shí)，CUDA程序與CPU串行程序計(jì)算結(jié)果的相對(duì)誤差較小，隨著網(wǎng)格數(shù)的增大，相對(duì)誤差也隨之增大，并逐漸趨向穩(wěn)定；不管是采用單精度計(jì)算，還是采用雙精度計(jì)算，相對(duì)誤差均未超過(guò)2×10－3，這說(shuō)明模擬結(jié)果穩(wěn)定，可靠。

表1 不同網(wǎng)格數(shù)CUDA程序模擬誤差Table 1 Simulation errors of CUDA code with different grids

3.3 性能分析

C語(yǔ)言的time庫(kù)函數(shù)中包含很多時(shí)間函數(shù)，其中的clock函數(shù)可以對(duì)程序運(yùn)行時(shí)間進(jìn)行檢測(cè)。為盡可能減小誤差，對(duì)每個(gè)程序運(yùn)行10次，將平均值作為最終結(jié)果。單精度的計(jì)算時(shí)間對(duì)比及加速比見(jiàn)圖6和7。

圖6 單精度計(jì)算時(shí)間對(duì)比Fig.6 Comparison of computation times of GPU code and CPU code with single precision

圖7 單精度計(jì)算加速比Fig.7 Speedup of the GPU code relative to the serial CPU code with single precision

圖8 雙精度計(jì)算時(shí)間對(duì)比Fig.8 Comparison of computation times of GPU code and CPU code with double precision

由于網(wǎng)格點(diǎn)數(shù)較少時(shí)計(jì)算時(shí)間較少，所以將y軸采用對(duì)數(shù)坐標(biāo)，從圖7中可以看出，當(dāng)網(wǎng)格點(diǎn)數(shù)較少（256×16）時(shí)，GPU運(yùn)行速度與CPU運(yùn)行速度相當(dāng)，加速比僅為1.82，這是由于數(shù)據(jù)傳入或傳出設(shè)備（GPU）、啟動(dòng)kernel等都需要浪費(fèi)一定的時(shí)間，當(dāng)網(wǎng)格數(shù)較少時(shí)，這部分所占時(shí)間比重較大，因此總的計(jì)算時(shí)間相差不大；隨著網(wǎng)格數(shù)的增加，CPU計(jì)算時(shí)間均急劇上升，而GPU計(jì)算時(shí)間上升較慢，即此時(shí)加速比（CPU計(jì)算時(shí)間／GPU計(jì)算時(shí)間）也急劇增大，但當(dāng)網(wǎng)格數(shù)增大到2 048×128以后，加速比增速明顯放緩，當(dāng)網(wǎng)格數(shù)為8 192×512時(shí)，加速比可以達(dá)到231.8倍。這說(shuō)明問(wèn)題規(guī)模越大，加速比也越大，越適合應(yīng)用CUDA程序進(jìn)行GPU并行加速計(jì)算，而問(wèn)題規(guī)模較小時(shí)，加速比較小，使用CUDA程序進(jìn)行計(jì)算意義不大。

當(dāng)采用雙精度計(jì)算時(shí)，可以得到與采用單精度計(jì)算時(shí)大致相同的結(jié)論（見(jiàn)圖8和9），只是加速比比用單精度計(jì)算要小，網(wǎng)格數(shù)為8 192×512時(shí)的加速比為141.7?？梢钥闯霰M管Fermi架構(gòu)相對(duì)Tesla架構(gòu)雙精度浮點(diǎn)處理能力已有很大的提升，但仍有很大的改進(jìn)余地。

圖9 雙精度計(jì)算加速比Fig.9 Speedup of the GPU code relative to the serial CPU code with double precision

4 結(jié)論

自從CUDA在2007年首次出現(xiàn)以來(lái)，已在醫(yī)學(xué)圖像、視頻播放、信號(hào)處理、生物計(jì)算等領(lǐng)域獲得了廣泛的應(yīng)用。這是因?yàn)镚PU相對(duì)CPU擁有更高帶寬的獨(dú)立顯存、適合處理并行計(jì)算任務(wù)且能夠大幅降低系統(tǒng)成本，從而為這些問(wèn)題提供了新的解決方案。

在物理海洋研究中，數(shù)值模擬是一個(gè)非常重要的研究手段。而高分辨率海洋數(shù)值模擬的計(jì)算量是很大的，為了盡快取得模擬結(jié)果，需要借助高性能計(jì)算機(jī)來(lái)減少程序的運(yùn)行時(shí)間。而CPU的計(jì)算能力無(wú)法滿(mǎn)足海洋數(shù)值模擬對(duì)計(jì)算能力的要求。海洋數(shù)值模式具有非常高的計(jì)算密度，非常適合使用GPU來(lái)進(jìn)行加速計(jì)算，而目前大部分海洋數(shù)值模式的計(jì)算都基于CPU為核心，還沒(méi)有移植到GPU平臺(tái)，因而無(wú)法利用GPU強(qiáng)大的計(jì)算能力。本項(xiàng)目使用CUDA實(shí)現(xiàn)了一個(gè)非靜壓海洋數(shù)值模式的GPU加速計(jì)算，在海洋數(shù)值模式的GPU加速方面進(jìn)行了初步的嘗試。與CPU串行程序相比，基于GPU的程序可以在不犧牲模擬精度的前提下（相對(duì)誤差不超過(guò)2×10－3）大大提高計(jì)算效率；單精度計(jì)算的加速比最高可以達(dá)到232，雙精度計(jì)算的加速比最高可以達(dá)到142。本項(xiàng)目為海洋數(shù)值模式向GPU平臺(tái)的移植積累了豐富的經(jīng)驗(yàn)，計(jì)算速度的顯著提高展現(xiàn)了海洋數(shù)值模式應(yīng)用GPU進(jìn)行加速的廣闊前景。

同時(shí)值得注意的是，CUDA中GPU的編程模型和CPU有很大的差別。數(shù)據(jù)的輸入輸出以及很多初始化工作需要由CPU完成，高密度的計(jì)算部分需要在GPU上完成，而CPU和GPU能訪(fǎng)問(wèn)的存儲(chǔ)空間是隔離的，這樣需要在主機(jī)端和設(shè)備端進(jìn)行數(shù)據(jù)的顯式傳輸，而且主機(jī)端代碼和設(shè)備端代碼有明顯的區(qū)分，因而將傳統(tǒng)的大量的海洋數(shù)值模式完全移植到GPU平臺(tái)是非常耗時(shí)的。除了直接使用CUDA外，還可以使用OpenACC或者一些GPU函數(shù)庫(kù)來(lái)實(shí)現(xiàn)GPU上的加速計(jì)算，但這些方式是否能夠既方便地將海洋數(shù)值模式移植到GPU平臺(tái)，又充分發(fā)揮GPU的計(jì)算能力，還需要進(jìn)行實(shí)際的使用和檢驗(yàn)。

［1］ Folkert B，Rob H B，Henk A D.Accelerating a barotropic ocean model using a GPU［J］.Ocean Modelling，2012，41：16－21.

［2］ 周 洪，樊曉椏，趙麗麗.基于CUDA的稀疏矩陣與矢量乘法的優(yōu)化［J］.計(jì)算機(jī)測(cè)量與控制，2010，18（8）：1906－1908.

［3］ 張 舒，褚艷利，趙開(kāi)勇，等.GPU高性能運(yùn)算之CUDA［M］.北京：中國(guó)水利水電出版社，2009：14－81.

［4］ Jochen K.Advanced Ocean Modelling：Using Open－Source Software［M］.Berlin：Germany Springer，2010：21－35.

［5］ CUDA C Best Practices Guide Version 4.1［EB／OL］.［2012－01－11］，［2012－12－03］.http：／／nvidia.com.

［6］ NVIDIA CUDA C Programming Guide Version 4.2［EB／OL］.［2012－04－16］，［2012－12－03］.http：／／nvidia.com.

［7］ Jason S，Edward K.CUDA by Example：An Introduction to General－Purpose GPU Programming［M］.USA：Addison－Wesley，2010：37－57.

［8］ Tuning CUDA Applications for Fermi Version 1.5［EB／OL］.［2011－05－03］，［2012－12－03］.http：／／nvidia.com.