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

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

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

1 引言

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

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

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

2 基于CPU的Julia集合

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

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

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

3 基于GPU的Julia集合

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

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

基于GPU的Julia集合算法和CPU的最關鍵區別在于kernel()函數的實現方式。基于GPU的kernel()函數:

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

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

int offset=x+y*gridDim.x;

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

4 兩種方法的實驗結果

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

表1 兩種方法的實驗結果

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

5 結束語

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

［1］陶懋頎.關于復迭代的Julia集的注記［J］.北京工業大學學報，1995(3):5－7.

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

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

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

［5］張艷，代巧蓮，等.高性能運算之CUDA應用分析［J］.電信技術研究，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］.現代計算機，2011(7):3－6.