實戰(zhàn)：使用MATLAB進行GPU高級編程

在GPU上執(zhí)行能夠加快我的應(yīng)用程序嗎?

GPU能夠?qū)Ψ弦韵聵藴实膽?yīng)用程序進行加速：

大規(guī)模并行—計算能夠被分割成上百個或上千個獨立的工作單元。

計算密集型—計算消耗的時間顯著超過了花費轉(zhuǎn)移數(shù)據(jù)到GPU內(nèi)存以及從GPU內(nèi)存轉(zhuǎn)移出數(shù)據(jù)的時間。

不滿足上述標準的應(yīng)用程序在GPU上運行時可能會比CPU要慢。

使用MATLAB進行GPU編程

FFT，IFFT以及線性代數(shù)運算超過了100個內(nèi)置的MATLAB函數(shù)，通過提供一個類型為GPUArray(由并行計算工具箱提供的特殊數(shù)組類型)的輸入?yún)?shù)，這些函數(shù)就能夠直接在GPU上運行。這些啟用GPU的函數(shù)都是重載的，換句話說，這些函數(shù)根據(jù)傳遞的參數(shù)類型的不同而執(zhí)行不同的操作。

例如，以下代碼使用FFT算法查找CPU上偽隨機數(shù)向量的離散傅里葉變換：

A = rand(2^16,1);

B = fft (A);

為在GPU上執(zhí)行相同的操作，我們首先使用gpuArray命令將數(shù)據(jù)從MATLAB工作空間轉(zhuǎn)移至GPU設(shè)備內(nèi)存。然后我們能夠運行重載函數(shù)fft：

A = gpuArray(rand(2^16,1));

B = fft (A);

fft操作在GPU上而不是在CPU上執(zhí)行，因為輸入?yún)?shù)(GPUArray)位于GPU的內(nèi)存中。

結(jié)果B存儲在GPU當中。然而，B在MATLAB工作空間中依舊可見。通過運行class(B)，我們看到B是一個GPUArray。

class(B)

ans =

parallel.gpu.GPUArray

我們能夠使用啟用GPU的函數(shù)繼續(xù)對B進行操作。例如，為可視化操作結(jié)果，plot命令自動處理GPUArrays。

plot(B);

為將數(shù)據(jù)返回至本地的MATLAB工作集，你可以使用gather命令。例如

C = gather(B);

C現(xiàn)在是MATLAB中的double，能夠被處理double變量的所有MATLAB函數(shù)操作。

在這個簡單的例子當中，執(zhí)行單個FFT函數(shù)節(jié)省的時間通常少于將向量從MATLAB工作集移動到設(shè)備內(nèi)存的時間。一般來說是這樣的但是也取決于硬件和陣列規(guī)模。數(shù)據(jù)傳輸開銷可能變得異常顯著以至于降低了應(yīng)用的總體性能，尤其是當你重復(fù)地在CPU和GPU之間交換數(shù)據(jù)，執(zhí)行相對來說很少的計算密集型操作時。更有效率的方式是當數(shù)據(jù)處于GPU當中時對數(shù)據(jù)進行一些操作，只在必要的情況下才將數(shù)據(jù)返回至CPU。

需要指出的是，和CPU類似，GPU的內(nèi)存也是有限的。然而，與CPU不同，GPU不能在內(nèi)存和硬盤之間交換數(shù)據(jù)。因此，你必須核實你希望保留在GPU當中的數(shù)據(jù)不會超出內(nèi)存的限制，尤其是當用到大規(guī)模矩陣時。通過運行g(shù)puDevice命令，可以查詢GPU卡，獲取信息比如名稱，總內(nèi)存以及可用內(nèi)存。

采用MATLAB解波動方程

為將上述例子應(yīng)用到具體的環(huán)境中，我們在一個實際的問題中實現(xiàn)GPU的功能。計算目標是解二階波動方程。

當u=0時到達臨界值。我們使用基于波譜法的算法解空間方程，使用基于二階中心有限差分法的算法解時間方程。

波譜法通常用于解決偏微分方程。采用波譜法的解決方案接近連續(xù)基函數(shù)比如正弦和余弦的線性組合。在這個例子中，我們應(yīng)用了切比雪夫波譜法，使用切比雪夫多項式作為基函數(shù)。

我們在每一個時間步長使用切比雪夫波普法計算當前解決方案的在x象限和y象限的二次導(dǎo)數(shù)。我們同時使用這些中間數(shù)值與舊的解決方案和新的解決方案，應(yīng)用二階中心有限差分法(也稱為蛙跳法)計算新的解決方案。我們選擇了保持蛙跳法穩(wěn)定性的時間步長。

MATLAB算法是計算密集型的，當網(wǎng)格中元素的數(shù)目超過了計算解決方案的增長，算法的執(zhí)行時間將顯著增加。當在單個CPU上使用2048x2048的網(wǎng)格執(zhí)行時，完成50個時間步長需要一分多鐘。需要指出的是我們計算的時間已經(jīng)包括了MATLAB內(nèi)在的多線程性能優(yōu)勢。自從R2007a起，MATLAb的一些函數(shù)就支持多線程計算。這些函數(shù)自動在多線程上執(zhí)行，并不需要在代碼中顯示指定命令去創(chuàng)建線程。

當考慮如何使用并行計算工具箱加速計算時，我們將關(guān)注每個時間步長所執(zhí)行的計算指令代碼。圖3距離說明了為獲取在GPU上運行的算法需要做出的改變。需要指出的是涉及MATLAB操作的計算指令、啟用GPU的重載函數(shù)可以從并行計算工具箱獲取。這些操作包括FFT，IFFT，矩陣乘法，以及各種元素明智(element-wise)操作。因此，我們不必改變算法就能夠在GPU執(zhí)行。只需要在進入每個時間步長計算結(jié)果的循環(huán)前使用gpuArray將數(shù)據(jù)轉(zhuǎn)移到GPU當中。  #p#page_title#e#

 
圖 3. 代碼對比工具顯示了CPU版本和GPU版本的差異。

CPU和GPU版本共享的代碼超過了84%(在111行當中有94行)。

計算指令在GPU上執(zhí)行后，我們將計算結(jié)果從GPU轉(zhuǎn)移至CPU。被啟用GPU的函數(shù)所引用的每個變量必須在GPU上創(chuàng)建或者在使用前轉(zhuǎn)移到GPU上。

為將用于光譜分化的一個權(quán)重轉(zhuǎn)變?yōu)镚PUArray變量，我們使用

W1T = gpuArray(W1T);

某些類型的數(shù)組能夠直接在GPU上構(gòu)造，不用從MATLAB工作集轉(zhuǎn)移。例如，為直接在GPU上創(chuàng)建全零矩陣，我們使用

uxx = parallel.gpu.GPUArray.zeros(N+1,N+1);

我們使用gather函數(shù)將數(shù)據(jù)從GPU中轉(zhuǎn)移回MATLAB工作集;例如：

vvg = gather(vv);

需要指出的是這只是將一個數(shù)據(jù)轉(zhuǎn)移至GPU，然后從GPU轉(zhuǎn)移回MATLAB工作集。每個時間步長的所有計算指令都是在GPU上執(zhí)行的。