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.