Julia集合CPU和GPU方法的分析比較

李改紅，劉金義，謝 陽(yáng)，馬 梁

(遼寧石油化工大學(xué)計(jì)算機(jī)與通訊工程學(xué)院，撫順 113001)

1 引言

Julia集合是一個(gè)著名的分形集，它是復(fù)數(shù)經(jīng)過(guò)迭代得到的。它的定義是，在復(fù)平面上，對(duì)于復(fù)數(shù)Z和C，如果變換Z=Z2+C不會(huì)使Z向無(wú)窮逃逸，那么所有這些初始的復(fù)數(shù)Z所構(gòu)成的集合稱為Julia集，它隨著C的變化而變化。它的特點(diǎn)是經(jīng)過(guò)迭代后，最后的Z值有三種可能:①Z值沒(méi)有界限增加(趨向無(wú)窮);②Z值衰減(趨向Z0，使Z0=Z20+C);③Z值是變化的，即非1或非2。Julia集的形狀基本上分三種:像塵埃一樣的結(jié)構(gòu)、穩(wěn)定的固態(tài)型或樹枝狀。

Julia集的基本算法是，通過(guò)一個(gè)簡(jiǎn)單的迭代等式對(duì)復(fù)平面中的點(diǎn)求值。如果計(jì)算某個(gè)點(diǎn)時(shí)，迭代等式的計(jì)算結(jié)果是發(fā)散的，那么這個(gè)點(diǎn)就不屬于Julia集合。更明確的說(shuō)，就是在迭代等式中計(jì)算得到的一系列值都朝著無(wú)窮大的方向增加，那么這個(gè)點(diǎn)就不屬于Julia集合。相反，如果在迭代等式中計(jì)算得到的一系列值都在某個(gè)邊界范圍之內(nèi)，那么這個(gè)點(diǎn)就屬于 Julia 集合［1－2］。

近年來(lái)，圖形處理器GPU在可編程能力、并行計(jì)算能力和應(yīng)用范圍方面不斷提升和擴(kuò)展，其發(fā)展速度已經(jīng)遠(yuǎn)遠(yuǎn)超過(guò)了CPU。而且GPU在通用計(jì)算領(lǐng)域得到了廣泛的運(yùn)用，并且很多算法都得到了性能提升。本文就分析比較了Julia集合的CPU方法和GPU方法的性能。

2 基于CPU的Julia集合

該方法是使用標(biāo)準(zhǔn)C語(yǔ)言進(jìn)行編程的。首先，在main()函數(shù)中通過(guò)工具庫(kù)創(chuàng)建了一個(gè)大小合適的位圖圖像bitmap(DIM，DIM)，又聲明了一個(gè)指向位圖數(shù)據(jù)的指針ptr，然后調(diào)用核函數(shù)，并把指針傳遞給核函數(shù)。

核函數(shù)就是對(duì)將要繪制的所有點(diǎn)進(jìn)行迭代，并且在每次迭代時(shí)都調(diào)用julia()來(lái)判斷該點(diǎn)是不是屬于Julia集合。如果該點(diǎn)屬于集合，那么函數(shù)julia()將返回1，相反就返回0。如果返回1，就把該點(diǎn)的顏色設(shè)置為紅色，如果是0，顏色設(shè)置為黑色。變量offset是計(jì)算點(diǎn)在位圖中的索引。核函數(shù)如下:

在核函數(shù)中調(diào)用的julia()函數(shù)是整個(gè)算法的核心。Julia()函數(shù)首先將像素坐標(biāo)轉(zhuǎn)換為復(fù)數(shù)空間的坐標(biāo)。為了把復(fù)平面的原點(diǎn)定位到圖像的中心，在這里把像素位置移動(dòng)了DIM/2。而且，為了確保圖像范圍在(－1.0，1.0)之間，又把圖像的坐標(biāo)縮小了DIM/2倍。這樣給定一個(gè)圖像點(diǎn)(i，j)，就可以計(jì)算出復(fù)平面空間中的一個(gè)((DIM/2－i)/(DIM/2)，(DIM/2－j)/(DIM/2))。在計(jì)算出復(fù)空間的點(diǎn)之后，就需要判斷該點(diǎn)是不是屬于Julia集合。這是通過(guò)計(jì)算迭代等式Zn+1=Z2n+C來(lái)判斷的。C是一個(gè)任意的復(fù)數(shù)常量，可以任意賦值。Julia()函數(shù):

3 基于GPU的Julia集合

基于GPU的Julia集合的程序和CPU的很相似，但是執(zhí)行原理不同，GPU是并行計(jì)算，能同時(shí)啟動(dòng)很多線程并行運(yùn)行，執(zhí)行效率很高。

該方法的main()函數(shù)和CPU的執(zhí)行流程是一樣的。首先也是通過(guò)工具庫(kù)創(chuàng)建一個(gè)DIM*DIM大小的位圖圖像，也聲明了一個(gè)指針dev_bitmap，用來(lái)存儲(chǔ)設(shè)備上的數(shù)據(jù)副本，并用cudaMalloc()為它分配內(nèi)存。然后就是調(diào)用內(nèi)核函數(shù)，計(jì)算結(jié)果，并且把結(jié)果傳回到主機(jī)內(nèi)存。值得注意的是，在調(diào)用內(nèi)核函數(shù)kernel()時(shí)，程序指定了多個(gè)并行線程塊同時(shí)執(zhí)行它。因?yàn)槊總€(gè)像素點(diǎn)的計(jì)算和其他像素點(diǎn)的計(jì)算是相互獨(dú)立的，所以為每個(gè)計(jì)算的像素點(diǎn)都運(yùn)行核函數(shù)的一個(gè)副本［3－4］。

基于GPU的Julia集合算法和CPU的最關(guān)鍵區(qū)別在于kernel()函數(shù)的實(shí)現(xiàn)方式。基于GPU的kernel()函數(shù):

在上面的kernel()函數(shù)中不是使用for()循環(huán)生成像素索引，在CUDA運(yùn)行時(shí)通過(guò)變量blockIdx生成像素的索引，并傳遞給julia()函數(shù)。因?yàn)槁暶髁艘粋€(gè)DIM*DIM的線程格，線程格的每一維與圖像的每一維大小是相同的，所以在(0，0)和(DIM－1，DIM－1)之間的每個(gè)像素點(diǎn)都會(huì)得到一個(gè)線程塊［5－6］。

在計(jì)算出像素點(diǎn)的索引之后，就需要得到輸出緩沖區(qū)ptr中的線性偏移量。這個(gè)值是通過(guò)內(nèi)置gridDim計(jì)算的，gridDim是一個(gè)常量，保存的是線程格每一維的大小［7］。在這里，gridDim的值是(DIM，DIM)，所以，線性偏移量是行索引乘以線程格的寬度，然后加上列索引。即:

int offset=x+y*gridDim.x;

最后，定義了一個(gè)判斷某個(gè)點(diǎn)是不是屬于julia集合的方法julia()，這個(gè)方法的代碼和CPU的相同，只是該方法前有個(gè)device修飾符，這表示該方法是在GPU上執(zhí)行的。

4 兩種方法的實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)使用的CPU是AMD Athlon(速龍)II X2 B24雙核，它的一級(jí)數(shù)據(jù)緩存是64KB，一級(jí)代碼緩存也是64KB，二級(jí)緩存是 2MB，速度為 3.00GHz。GPU使用的是GeForce gts 450，它有192個(gè)流處理器，128bit顯存控制器，同時(shí)具有16個(gè)光柵單元和32個(gè)紋理單元，內(nèi)存帶寬為64GB/s，理論計(jì)算能力是1.008TFLOPs。CUDA 版本為2.3。實(shí)驗(yàn)結(jié)果如表1所示。

表1 兩種方法的實(shí)驗(yàn)結(jié)果

由上表可以看出，GPU方法的執(zhí)行時(shí)間比CPU的快了將近10倍。在CPU方法中每次只對(duì)位圖中的一個(gè)點(diǎn)計(jì)算它是不是屬于julia集合，遍歷完整個(gè)位圖中的所有點(diǎn)需要DIM*DIM次，而在GPU方法中使用了并行線程塊，在調(diào)用內(nèi)核函數(shù)時(shí)同時(shí)啟動(dòng)了DIM*DIM個(gè)線程塊，每個(gè)線程塊中有一個(gè)并行線程，那么同時(shí)就有DIM*DIM個(gè)線程在運(yùn)行，每個(gè)線程計(jì)算一個(gè)位圖點(diǎn)，這樣同時(shí)就有一個(gè)活動(dòng)塊的線程在執(zhí)行，所以執(zhí)行時(shí)間就縮短了，這就是并行計(jì)算的優(yōu)勢(shì)。由此可以看出GPU的并行計(jì)算能力很強(qiáng)，很適合大規(guī)模的數(shù)據(jù)計(jì)算。

5 結(jié)束語(yǔ)

從以上的實(shí)例比較可以看出GPU的計(jì)算能力很強(qiáng)，隨著應(yīng)用程序的不斷開拓與發(fā)展，基于GPU的通用計(jì)算將越來(lái)越成熟。同時(shí)，應(yīng)用CUDA的高性能計(jì)算，使數(shù)據(jù)并行處理的能力加速了。因?yàn)镚PU自身的通用計(jì)算能力和CUDA的開發(fā)平臺(tái)，相信在不久的未來(lái)，GPU在數(shù)據(jù)庫(kù)、天文計(jì)算、導(dǎo)航識(shí)別等等領(lǐng)域都將得到更好更廣泛的應(yīng)用。

［1］陶懋頎.關(guān)于復(fù)迭代的Julia集的注記［J］.北京工業(yè)大學(xué)學(xué)報(bào)，1995(3):5－7.

［2］詹棠森，李海林，史華平.分形理論中的Julia集的迭代函數(shù)系統(tǒng) IFS算法的改進(jìn)［J］.福建電腦，2005(8):43－44.

［3］吳恩華，柳有權(quán).基于圖形處理器的通用計(jì)算［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2004，16(5):603－604.

［4］劉進(jìn)鋒，郭雷.CPU與GPU上幾種矩陣乘法的比較與分析［J］.計(jì)算機(jī)工程與應(yīng)用，2011，47(19):9 －11.

［5］張艷，代巧蓮，等.高性能運(yùn)算之CUDA應(yīng)用分析［J］.電信技術(shù)研究，2011(3):43－45.

［6］Barcellos B，Coutinho S，et al.GPU based fluid animation over elastic surface models［J］.Brazilian Symposium on Games and Digital Entertainment，2009(18):119 －129.

［7］呂亞飛，賈 陽(yáng).基于CUDA的快速中值濾波算法［J］.現(xiàn)代計(jì)算機(jī)，2011(7):3－6.