實(shí)戰(zhàn):使用MATLAB進(jìn)行GPU高級(jí)編程
在GPU上執(zhí)行能夠加快我的應(yīng)用程序嗎?
GPU能夠?qū)Ψ弦韵聵?biāo)準(zhǔn)的應(yīng)用程序進(jìn)行加速:
大規(guī)模并行—計(jì)算能夠被分割成上百個(gè)或上千個(gè)獨(dú)立的工作單元。
計(jì)算密集型—計(jì)算消耗的時(shí)間顯著超過(guò)了花費(fèi)轉(zhuǎn)移數(shù)據(jù)到GPU內(nèi)存以及從GPU內(nèi)存轉(zhuǎn)移出數(shù)據(jù)的時(shí)間。
不滿足上述標(biāo)準(zhǔn)的應(yīng)用程序在GPU上運(yùn)行時(shí)可能會(huì)比CPU要慢。
使用MATLAB進(jìn)行GPU編程
FFT,IFFT以及線性代數(shù)運(yùn)算超過(guò)了100個(gè)內(nèi)置的MATLAB函數(shù),通過(guò)提供一個(gè)類型為GPUArray(由并行計(jì)算工具箱提供的特殊數(shù)組類型)的輸入?yún)?shù),這些函數(shù)就能夠直接在GPU上運(yùn)行。這些啟用GPU的函數(shù)都是重載的,換句話說(shuō),這些函數(shù)根據(jù)傳遞的參數(shù)類型的不同而執(zhí)行不同的操作。
例如,以下代碼使用FFT算法查找CPU上偽隨機(jī)數(shù)向量的離散傅里葉變換:
A = rand(2^16,1);
B = fft (A);
為在GPU上執(zhí)行相同的操作,我們首先使用gpuArray命令將數(shù)據(jù)從MATLAB工作空間轉(zhuǎn)移至GPU設(shè)備內(nèi)存。然后我們能夠運(yùn)行重載函數(shù)fft:
A = gpuArray(rand(2^16,1));
B = fft (A);
fft操作在GPU上而不是在CPU上執(zhí)行,因?yàn)檩斎雲(yún)?shù)(GPUArray)位于GPU的內(nèi)存中。
結(jié)果B存儲(chǔ)在GPU當(dāng)中。然而,B在MATLAB工作空間中依舊可見(jiàn)。通過(guò)運(yùn)行class(B),我們看到B是一個(gè)GPUArray。
class(B)
ans =
parallel.gpu.GPUArray
我們能夠使用啟用GPU的函數(shù)繼續(xù)對(duì)B進(jìn)行操作。例如,為可視化操作結(jié)果,plot命令自動(dòng)處理GPUArrays。
plot(B);
為將數(shù)據(jù)返回至本地的MATLAB工作集,你可以使用gather命令。例如
C = gather(B);
C現(xiàn)在是MATLAB中的double,能夠被處理double變量的所有MATLAB函數(shù)操作。
在這個(gè)簡(jiǎn)單的例子當(dāng)中,執(zhí)行單個(gè)FFT函數(shù)節(jié)省的時(shí)間通常少于將向量從MATLAB工作集移動(dòng)到設(shè)備內(nèi)存的時(shí)間。一般來(lái)說(shuō)是這樣的但是也取決于硬件和陣列規(guī)模。數(shù)據(jù)傳輸開(kāi)銷可能變得異常顯著以至于降低了應(yīng)用的總體性能,尤其是當(dāng)你重復(fù)地在CPU和GPU之間交換數(shù)據(jù),執(zhí)行相對(duì)來(lái)說(shuō)很少的計(jì)算密集型操作時(shí)。更有效率的方式是當(dāng)數(shù)據(jù)處于GPU當(dāng)中時(shí)對(duì)數(shù)據(jù)進(jìn)行一些操作,只在必要的情況下才將數(shù)據(jù)返回至CPU。
需要指出的是,和CPU類似,GPU的內(nèi)存也是有限的。然而,與CPU不同,GPU不能在內(nèi)存和硬盤(pán)之間交換數(shù)據(jù)。因此,你必須核實(shí)你希望保留在GPU當(dāng)中的數(shù)據(jù)不會(huì)超出內(nèi)存的限制,尤其是當(dāng)用到大規(guī)模矩陣時(shí)。通過(guò)運(yùn)行g(shù)puDevice命令,可以查詢GPU卡,獲取信息比如名稱,總內(nèi)存以及可用內(nèi)存。
采用MATLAB解波動(dòng)方程
為將上述例子應(yīng)用到具體的環(huán)境中,我們?cè)谝粋€(gè)實(shí)際的問(wèn)題中實(shí)現(xiàn)GPU的功能。計(jì)算目標(biāo)是解二階波動(dòng)方程。
當(dāng)u=0時(shí)到達(dá)臨界值。我們使用基于波譜法的算法解空間方程,使用基于二階中心有限差分法的算法解時(shí)間方程。
波譜法通常用于解決偏微分方程。采用波譜法的解決方案接近連續(xù)基函數(shù)比如正弦和余弦的線性組合。在這個(gè)例子中,我們應(yīng)用了切比雪夫波譜法,使用切比雪夫多項(xiàng)式作為基函數(shù)。
我們?cè)诿恳粋€(gè)時(shí)間步長(zhǎng)使用切比雪夫波普法計(jì)算當(dāng)前解決方案的在x象限和y象限的二次導(dǎo)數(shù)。我們同時(shí)使用這些中間數(shù)值與舊的解決方案和新的解決方案,應(yīng)用二階中心有限差分法(也稱為蛙跳法)計(jì)算新的解決方案。我們選擇了保持蛙跳法穩(wěn)定性的時(shí)間步長(zhǎng)。
MATLAB算法是計(jì)算密集型的,當(dāng)網(wǎng)格中元素的數(shù)目超過(guò)了計(jì)算解決方案的增長(zhǎng),算法的執(zhí)行時(shí)間將顯著增加。當(dāng)在單個(gè)CPU上使用2048x2048的網(wǎng)格執(zhí)行時(shí),完成50個(gè)時(shí)間步長(zhǎng)需要一分多鐘。需要指出的是我們計(jì)算的時(shí)間已經(jīng)包括了MATLAB內(nèi)在的多線程性能優(yōu)勢(shì)。自從R2007a起,MATLAb的一些函數(shù)就支持多線程計(jì)算。這些函數(shù)自動(dòng)在多線程上執(zhí)行,并不需要在代碼中顯示指定命令去創(chuàng)建線程。
當(dāng)考慮如何使用并行計(jì)算工具箱加速計(jì)算時(shí),我們將關(guān)注每個(gè)時(shí)間步長(zhǎng)所執(zhí)行的計(jì)算指令代碼。圖3距離說(shuō)明了為獲取在GPU上運(yùn)行的算法需要做出的改變。需要指出的是涉及MATLAB操作的計(jì)算指令、啟用GPU的重載函數(shù)可以從并行計(jì)算工具箱獲取。這些操作包括FFT,IFFT,矩陣乘法,以及各種元素明智(element-wise)操作。因此,我們不必改變算法就能夠在GPU執(zhí)行。只需要在進(jìn)入每個(gè)時(shí)間步長(zhǎng)計(jì)算結(jié)果的循環(huán)前使用gpuArray將數(shù)據(jù)轉(zhuǎn)移到GPU當(dāng)中。 #p#page_title#e#
圖 3. 代碼對(duì)比工具顯示了CPU版本和GPU版本的差異。
CPU和GPU版本共享的代碼超過(guò)了84%(在111行當(dāng)中有94行)。
計(jì)算指令在GPU上執(zhí)行后,我們將計(jì)算結(jié)果從GPU轉(zhuǎn)移至CPU。被啟用GPU的函數(shù)所引用的每個(gè)變量必須在GPU上創(chuàng)建或者在使用前轉(zhuǎn)移到GPU上。
為將用于光譜分化的一個(gè)權(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);
需要指出的是這只是將一個(gè)數(shù)據(jù)轉(zhuǎn)移至GPU,然后從GPU轉(zhuǎn)移回MATLAB工作集。每個(gè)時(shí)間步長(zhǎng)的所有計(jì)算指令都是在GPU上執(zhí)行的。