李改紅，劉金義，謝 陽，馬 梁

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

1 引言

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

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

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

2 基于CPU的Julia集合

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

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

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

3 基于GPU的Julia集合

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

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

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

在上面的kernel()函數(shù)中不是使用for()循環(huán)生成像素索引，在CUDA運行時通過變量blockIdx生成像素的索引，并傳遞給julia()函數(shù)。因為聲明了一個DIM*DIM的線程格，線程格的每一維與圖像的每一維大小是相同的，所以在(0，0)和(DIM－1，DIM－1)之間的每個像素點都會得到一個線程塊［5－6］。

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

int offset=x+y*gridDim.x;

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

4 兩種方法的實驗結(jié)果

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

表1 兩種方法的實驗結(jié)果

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

5 結(jié)束語

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

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

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

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

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

［5］張艷，代巧蓮，等.高性能運算之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］呂亞飛，賈 陽.基于CUDA的快速中值濾波算法［J］.現(xiàn)代計算機，2011(7):3－6.