GPU上的Mandelbrot Set

发布的罗兰美国舒尔，2011年7月18日

27次观看(过去30天)|0喜欢|8的评论

今天我欢迎客座博主回来本Tordoff谁上次在这里写过关于如何做的博客雕刻恐龙．这些天本在做并行计算工具箱特别关注gpu。

内容

关于这个演示

Benoit Mandelbrot于去年秋天去世，尽管他发明了“分形”这个术语，并且是分形几何学的先驱之一，但对大多数人来说，他留下的永恒遗产无疑是以他的名字命名的分形。很少有数学家或数学算法能像曼德布罗特在80年代末和90年代初设计的那样，装饰了那么多青少年卧室的墙壁。

我还怀疑，在不同的平台上，没有哪一种分形语言能像它一样被如此多的语言所实现。我的第一个版本是在一台老旧的ZX Spectrum上的32x24块状版本，有七种绚丽的颜色。现在我们有了更多的图形功能，也有了更多的颜色。事实上，正如我们将看到的，像Mandelbrot集这样的算法非常适合运行在GPU(图形处理单元)上，而GPU通常闲置在您的PC中。

作为起点，我将使用一个Mandelbrot集的版本，它松散地取自Cleve Moler的集合MATLAB实验电子书。然后，我们将看看并行计算工具箱™提供的三种不同的方式来加速使用GPU:

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

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

设置

下面的值指定了Mandelbrot集的一个高度放大的部分，该部分位于主心脏线和它左边的p/q灯泡之间的山谷。

下面的代码通过在这些极限之间创建一个1000x1000的实部(X)和虚部(Y)网格，形成了每个计算的起点集。对于这个特定的位置，我碰巧知道500次迭代就足够画出一幅漂亮的图了。

maxIterations = 500;gridSize = 1000;Xlim = [-0.748766713922161， -0.748766707771757];Ylim = [0.123640844894862, 0.123640851045266];

MATLAB中的Mandelbrot集

下面是使用标准MATLAB命令在CPU上运行的Mandelbrot Set的实现。这个计算是向量化的，因此每个位置都是一次更新的。

%设置T = tic();x = linspace(xlim(1)， xlim(2)， gridSize);y = linspace(ylim(1)， ylim(2)， gridSize);[xGrid,yGrid] = meshgrid(x, y);z0 = xGrid + 1i*yGrid;Count = 0 (size(z0));%计算Z = z0;为n = 0:maxIterations z = z.*z + z0;Inside = abs(z)<=2;计数=计数+内部;结束Count = log(Count +1);%显示cpuTime = toc(t);集(gcf,“位置”， [200 200 600 600]);Imagesc (x, y, count);轴图像Colormap ([jet();flipud(jet());0 0 0]);标题(sprintf ('%1.2fsecs(不含GPU)'， cpuTime));

使用parallel.gpu.GPUArray -快16倍

parallel.gpu.GPUArray提供了许多可用于创建数据数组的函数的GPU版本，包括linspace，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 = complex( xGrid, yGrid ); count = parallel.gpu.GPUArray.zeros( size(z0) );%计算Z = z0;为n = 0:maxIterations z = z.*z + z0;Inside = abs(z)<=2;计数=计数+内部;结束Count = log(Count +1);%显示naiveGPUTime = toc(t);Imagesc (x, y，计数)轴图像标题(sprintf ('%1.2fsecs(原生GPU) = %1.1fx更快'，.．.naiveGPUTime, cpuTime/naiveGPUTime))

使用元素操作- 164倍快

注意到算法在输入的每个元素上都是平等的，我们可以将代码放在一个helper函数中，并使用它来调用它arrayfun．对于GPU数组输入，arrayfun使用的函数被编译成原生GPU代码。在本例中，我们将循环放入processMandelbrotElement.m如下所示:

函数count = processMandelbrotElement(x0,y0,maxIterations) z0 = complex(x0,y0);Z = z0;Count = 1;while (count <= maxIterations)…& &((真正的(z) * (z) +图像放大(z) *图像放大(z)) < = 4)数=计数+ 1;Z = Z * Z + z0;结束计数= log(count);

请注意，由于此函数只处理单个元素，因此引入了早期中止。对于Mandelbrot Set的大多数视图，大量的元素很早就停止了，这可以节省大量的处理。的为循环也被替换为而循环，因为它们通常更有效率。这个函数没有提到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);Imagesc (x, y，计数)轴图像标题(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代表“并行线程执行语言”)。

由于c++中缺少复数，CUDA/ c++代码比我们目前看到的MATLAB版本要复杂一些。但是，算法的本质没有改变:

__device__ unsigned int doIterations(double const realPart0, double const imagPart0, unsigned int const max){//初始化:z = z0 double realPart = realPart0;double imagPart = imagPart0;Unsigned int count = 0;//循环直到转义while ((count <= max) && ((realPart*realPart + imagPart*imagPart) <= 4.0)) {++count;//更新:z = z*z + z0;double const oldRealPart = realPart;realPart = realPart*realPart - imagPart*imagPart + 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);内核。ThreadBlockSize = [kernel.MaxThreadsPerBlock,1,1];内核。GridSize = [ceil(numElements/kernel.MaxThreadsPerBlock)，1];调用内核count = parallel.gpu.GPUArray。0 (size(xGrid));count = feval(内核，计数，xGrid, yGrid, maxIterations, numElements);%显示gpuCUDAKernelTime = toc(t);Imagesc (x, y，计数)轴图像标题(sprintf ('%1.2fsecs (GPU CUDAKernel) = %1.1fx更快'，.．.gpuCUDAKernelTime, cpuTime/gpuCUDAKernelTime));