GPU上的Mandelbrot集

发布者洛伦舒尔,2011年7月18日

38次观看（过去30天）| 0喜欢| 8评论

今天我欢迎回来的嘉宾博客本托福夫谁上次颠簸了雕刻恐龙.这些天本在办公室工作并行计算工具箱™特别关注GPU。

关于这个演示

Benoit Mandelbrot于去年秋天去世，尽管他发明了“分形”一词作为分形几何学的先驱之一，他留给大多数人的永恒遗产无疑是以他的名字命名的分形。很少有数学家或数学算法像曼德布罗特在80年代末和90年代初的作品那样装饰了青少年卧室的墙壁。

我还怀疑，在不同的平台上，没有其他的分形已经用很多不同的语言实现了。我的第一个版本是一个块状32x24版本，在老化的ZX光谱上有七种亮丽的颜色。现在我们有了更多的图形功能和更多的颜色。事实上，正如我们将要看到的，像Mandelbrot集这样的算法非常适合在GPU（图形处理单元）上运行，而GPU在您的PC中大部分处于空闲状态。

作为一个起点，我将使用来自Clever Moler的松散的Mandelbrot集的版本用MATLAB进行实验电子书。然后，我们将研究并行计算工具箱使用的三种不同方式™ 用于使用GPU加速此过程：

使用现有算法，但使用GPU数据作为输入
使用Arrayfun.独立地在每个元素上执行算法
使用MATLAB / CUDA接口运行一些现有的CUDA / C ++代码

正如您将看到的，这三种方法难度越来越大，但回报的速度却大大提高。选择一种方法来权衡难度和速度是典型的并行问题。

安装程序

下面的值指定了Mandelbrot集合的高度缩放部分，该集合位于主心形线与其左侧的p/q灯泡之间的山谷中。

以下代码通过在这些限制之间创建1000x1000网格和虚部（y）之间的每一个计算，形成每个计算的起点集。对于这个特殊的位置，我碰巧知道500次迭代将足以画出一个漂亮的照片。

最大= 500;gridsize = 1000;XLIM = [-0.748766713922161，-0.7487667077771757];ylim = [0.12364084894862,0.123640851045266];

MATLAB中的Mandelbrot集

下面是使用在CPU上运行的标准MATLAB命令实现的Mandelbrot集。该计算是矢量化的，因此每个位置都会立即更新。

％ 设置t = tic（）;x = Linspace（XLIM（1），XLIM（2），网格化）;y = linspace（ylim（1），ylim（2），gridsize）;[xGrid，ygrid] = meshgrid（x，y）;z0 = Xgrid + 1i * ygrid;count = zeros（尺寸（z0））;%算计z = z0;对于n = 0：maxIterations z = z。* z + z0;内部= abs（z）<= 2;count =内部计数+;终止计数=日志（计数+1）；%展示cputime = toc（t）;设置（GCF，“位置”, [200 200 600 600] ); 图像sc（x，y，计数）；轴形象彩色贴图（[jet（）；flipud（jet（））；0]）；头衔（斯普林特）(“%1.2秒（不带GPU）”，cpuTime））；

使用parallel.gpu.GPUArray-速度提高16倍

parallel.gpu.GPUArray提供可用于创建数据阵列的许多函数的GPU版本，包括邻域,logspace.，和Meshgrid.这里需要的功能。同样，计数使用该函数直接在GPU上初始化数组parallel.gpu.GPUArray.zeros.

除了这些简单的数据初始化的简单变化之外，算法不变（和> 16次更快）：

％ 设置t=tic（）；x=parallel.gpu.GPUArray.linspace（xlim（1），xlim（2），gridSize）；y=parallel.gpu.GPUArray.linspace（ylim（1），ylim（2），gridSize）；[xGrid，yGrid]=meshgrid（x，y）；z0=复合物（xGrid，yGrid）；计数=并行.gpu.GPUArray.zero（大小（z0））；%算计z = z0;对于n = 0：maxIterations z = z。* z + z0;内部= abs（z）<= 2;count =内部计数+;终止计数=日志（计数+1）；%展示天真才= toc（t）;ImageC（x，y，count）轴形象标题（Sprintf（'％1.2fsecs（天真gpu）=％1.1fx更快',......naiveGPUTime，cpuTime/naiveGPUTime））

使用元素操作-速度提高164倍

注意到算法在输入的每个元素上都是平等操作的，我们可以将代码放在一个helper函数中，并使用Arrayfun.. 对于GPU阵列输入，与arrayfun一起使用的函数将编译为本机GPU代码。在本例中，我们将循环放置在processmandelbrotelement.m.如下所示：

功能count = processmandelbrotelement（x0，y0，maxIrlations）z0 =复合物（x0，y0）;z = z0;count = 1;虽然（count <= maxIrtations）... &&（（z）* real（z）+ imag（z）* imag（z））<= 4）count = count + 1;z = z * z + z0;结束count = log（count）;

请注意，由于此函数仅处理单个元素，因此引入了早期中止。对于Mandelbrot集合的大多数视图，大量元素很早就停止，这可以节省大量处理。这个对于循环也被一个尽管循环，因为它们通常效率更高。此函数未提及GPU，也未使用GPU的特定功能-这是标准的MATLAB代码。

使用Arrayfun.这意味着，我们不需要调用数千个单独的GPU优化操作（每次迭代至少6个），而只需要调用一个执行整个计算的并行GPU操作。这大大减少了开销，比以前快了164倍。

％ 设置t = tic（）;x = parallel.gpu.gpuarray.linspace（xlim（1），xlim（2），gridsize）;y = parallel.gpu.gpuarray.linspace（ylim（1），ylim（2），gridsize）;[xGrid，ygrid] = meshgrid（x，y）;%算计count=arrayfun（@processMandelbrotElement、xGrid、yGrid、maxIterations）；%展示gpuarrayfuntime = toc（t）;ImageC（x，y，count）轴形象标题（Sprintf（'％1.2fsecs（gpu arrayfun）=％1.1fx更快',......gpuArrayfunTime，cpuTime/gpuArrayfunTime））；

在CUDA工作-速度提高340倍

在MATLAB中的实验通过将基本算法转换为C-Mex函数，可以提高性能。如果您愿意在C/C++中做一些工作，那么可以使用并行计算工具箱使用MATLAB数据调用预先编写的CUDA内核。这是使用工具箱的parallel.gpu.cudakernel.特征。

元素处理算法的CUDA/C++实现已经用手工编写processmandelbrotelement.cu.. 然后必须使用nVidia的NVCC编译器手动编译，以生成程序集级别processMandelbrotElement.ptx(.ptx代表“并行线程执行语言”）。

CuDA/C++代码比我们目前看到的MATLAB版本要复杂得多，因为C++中缺少复数。但是，该算法的本质没有改变：

__设备uuu无符号整数迭代（双常量realPart0，双常量imagPart0，无符号整数常量maxIters）{//Initialize:z=z0双realPart=realPart0；双imagPart=imagPart0；无符号整数计数=0；//在（（计数<=maxIters）和（（realPart*realPart+imagPart*imagPart）<=4.0）时循环直到转义）{++count；//更新：z=z*z+z0；双常量oldRealPart=realPart；realPart=realPart*realPart-imagPart*realPart+realPart0；imagPart=2.0*oldRealPart*imagPart+imagPart0；}返回计数；}

（完整的源代码在末尾链接的文件交换提交中提供。）

阵列中的每个元素需要一个GPU线程，分为块。内核表示线程块的尺寸有多大，这用于计算所需的块数。

％ 设置t = tic（）;x = parallel.gpu.gpuarray.linspace（xlim（1），xlim（2），gridsize）;y = parallel.gpu.gpuarray.linspace（ylim（1），ylim（2），gridsize）;[xGrid，ygrid] = meshgrid（x，y）;%加载内核kernel = parallel.gpu.cudakernel（“processMandelbrotElement.ptx”,......'processmandelbrotelement.cu');％确保我们有足够的块来覆盖整个阵列numElements = numel（xgrid）;kernel.threadblocksize = [kernel.maxthreadsperblock，1,1];kernel.gridsize = [ceil（numElements / kernel.maxthreadsperblock），1];%调用内核count=parallel.gpu.GPUArray.zero（大小（xGrid））；count=feval（内核、count、xGrid、yGrid、maxIterations、numElements）；%展示gpuCUDAKernelTime=toc（t）；图像SC（x，y，计数）轴形象标题（Sprintf（“%1.2秒（GPU CUDAKernel）=%1.1秒更快”,......gpucudakertime，cpuTime/gpucudakertime））；