用GPU通用并行计算绘制曼德勃罗特集图形上篇（转）

http://www.cnblogs.com/Ninputer/archive/2009/11/24/1609364.html

近年来PC的计算能力发生了天翻地覆的变化。CPU逐渐趋向于多核发展，同时内存带宽和缓存不断增加，如今的PC已经成为小型的统一地址空间的并行计算机。然而我们的PC中还有一个设备可以提供比CPU更加强大的并行计算设备——显卡，它在进行充分并行的任务时可以提供高达数TFLOPS的峰值运算能力，这几乎是2000-2001年间国产超级计算机的运算能力。在显卡刚出现时，显卡内的模块都是为特定的图形任务而设计的。比如会有光照和坐标变换单元以及光栅单元等。随着显卡和图形技术的发展，在DirectX 8出现之后显卡开始允许在这些特殊单元处理的基础上增加编程控制。比如光照和坐标变换阶段可以用顶点着色器（Vertex Shader）进行控制，而光栅后的阶段则可以使用像素着色器（Pixel Shader）进行控制。这些“着色器（Shader）”其实就是运行在显卡上的程序。在DirectX 10阶段之后，整个渲染阶段可编程控制阶段大大增强，整个渲染流水线中各个阶段都可以通过着色器语言进行控制。因此DirectX 10后的显卡，都将执行着色器的运算单元从渲染过程中独立出来，成为一个统一着色单元，这些统一着色单元是可以大规模并行执行的运算单元。人们逐渐发现这些统一着色单元实际上可以进行与图形技术无关的通用技术。显卡的性能几乎按照三倍于摩尔定律的速度发展，GPU用于通用计算的潜力也非常巨大。比如本文编写时最强的GPU芯片——Ati Radeon HD5870可以提供2.72TFLOPS的单精度峰值浮点运算性能和500GFLOPS左右的双精度峰值浮点运算性能。而时下最先进的酷睿i7四核处理器只能提供60-80GFLOPS的峰值浮点运算能力。

GPU芯片制造商nVidia和Ati是GPU通用计算的先行者，他们分别提供了 CUDA和Ati Stream技术用于GPU程序的开发。随后由苹果等厂商领导并获得诸多厂商支持的OpenCL出台，成为GPU通用计算的统一标准。而微软也看到了这个潜力，在DirectX 11中提供了不与具体渲染流程绑定的计算着色器（Compute Shader）。虽然还叫“着色器”，但是Compute Shader已经是完全通用的编程技术。编程语言方面CUDA采用的编程语言有扩展的C、C++和Fortran等，OpenCL采用的是扩展的标准C。而DirectX Compute Shader是和C类似的HLSL语言。从开发难度来讲DirectX要比CUDA和OpenCL难用。HLSL没有指针，而且需要动态编译和执行。但是目前Compute Shader 5.0支持双精度浮点数、完善的原子操作和三维线程组支持。它还支持32KB的线程组共享内存。因此采用Compute Shader进行通用计算仍然是一个值得研究的领域。我也是DirectX的初学者，对它的图形技术一窍不通，完全是本着通用计算的目的来学习HLSL。下面就动手来实现一个非常适合并行计算的程序——曼德勃罗特集。

曼德勃罗特（Mandelbrot）集是一种复平面上的点集。对任意复数c，我们采用以下公式：

进行迭代，所得的所有复数的集合就叫曼德勃罗特集，其中z从0开始。这个公式如此的简单，但效果却非同一般。我们会发现许多点在迭代中处于亚稳定的状态，它们会在迭代数次之后逃逸（发散）。而这个迭代次数对每一点c而言都是不同的。如果我们画出复平面中各个点的迭代次数，就会发现它能组成难以置信的图案。

首先我们能够发现每一个复数c都能够独立地采用这个这个公式进行运算，它无需与其他的点进行任何交互。因此该算法默认就是完全并行的。我们只需要按照每一个输入进行划分就可以轻易地写出HLSL的代码。

我们先来看看代码，然后再来分析：

#define MAX_ITER 32768

//--------------------------------------------------------------------------------------
// Constant Buffers
//--------------------------------------------------------------------------------------
cbuffer CB : register( b0 )
{
    unsigned int c_Stride;
    unsigned int c_Width;
    unsigned int c_Height;
    float c_RealMin;
    float c_ImagMin;
    float c_ScaleReal;
    float c_ScaleImag;
};


RWStructuredBuffer<unsigned int> Data : register( u0 );

inline uint ComposeColor(uint index)
{    
    index = index * clamp(0, 1, MAX_ITER - index);
    uint red, green, blue;

    float phase = index * 3.0f / MAX_ITER;
    red = (uint)(max(0.0f, phase - 2.0f) * 255.0f);
    green = (uint)(smoothstep(0.0f, 1.0f, phase - 1.3f) * 255.0f);
    blue = (uint)(max(0.0f, 1.0f - abs(phase - 1.0f)) * 255.0f);
    
    return 0xff000000 | (red << 16) | (green << 8) | blue;
}

[numthreads(16, 16, 1)]
void Calculate(uint3 DTid : SV_DispatchThreadID)
{
    float2 c;
    c.x = c_RealMin + (DTid.x * c_ScaleReal);
    c.y = c_ImagMin + ((c_Width - DTid.y) * c_ScaleImag);
    
    float2 z;
    z.x = 0.0f;
    z.y = 0.0f;
    
    float temp, lengthSqr;
    uint count = 0;
    
    do
    {
        temp = z.x * z.x - z.y * z.y + c.x;
        z.y = 2 * z.x * z.y + c.y;
        z.x = temp;
        
        lengthSqr = z.x * z.x + z.y * z.y;
        count++;
    }
    while ((lengthSqr < 4.0f) && (count < MAX_ITER));    
    
    //write to result
    uint currentIndex = DTid.x + DTid.y * c_Width;
    
    Data[currentIndex] = ComposeColor(log(count) / log(MAX_ITER) * MAX_ITER);
}

如你所见，HLSL和C语言十分相似，只是有一些特别的扩展。一开始的cbuffer称为Constant Buffer，是显卡里所有线程都能访问的资源。我们调用Compute Shader程序时传的参数就都要放入Constant Buffer当中。下面就是HLSL的函数。HLSL没有栈，因此函数不能递归调用。所有函数的参数和本地变量都储存在寄存器中。共有4096个寄存器可以用于这个目的。

ComposeColor是我为了画出好看的颜色而编写的，可以自由发挥，下面跳过它看Calculate，这就是我们的主角。DirectX中每一个着色器程序都绑定到一个特定的渲染流程上，因此着色器程序通常都有一个作用目标。例如像素着色器针对每个像素执行，顶点着色器针对每个顶点执行。而计算着色器不与这些图形概念绑定，计算着色器针对的是“线程”。每一个线程执行一次。在这段代码里[numthreads(16, 16, 1)]是对线程进行分组的attribute。Compute Shader 5.0可以将线程分为三维的线程组。这个attribute就用来描述每一组的线程数。我们这个应用实际不需要进行分组，但是分组可以增加处理数据的规模。比如采用1024个线程的线程组，就能处理多达67108864个输入数据。请注意，计算着色器总是大量线程并行执行的状态，编写它的时候一定要变换思想，不要采用顺序编程的思路。

我们的Calculate函数输入参数DTid是一个三维向量，它有一个显式用途SV_DispatchThreadID，这表示线程的分布 ID（DispatchID），也就是线程组的ID与组内线程ID的乘积。我们用它来表示复平面的输入点坐标。首先我们把输入点换算成从-2i – 2i和-2 - 2这个区间的点坐标。这就是我们公式里的c。然后就简单了，按照公式迭代就行了。迭代终止条件有两个：一是迭代变量z的模大于2（表示这个点发散了），二是达到了最大的迭代次数MAX_ITER。每个点迭代次数具体是多少是无法知道的，不过它不会大于MAX_ITER。这里进行复数运算时所用的类型是 float2，它是一个有两个float分量的向量。

最后我们要把迭代次数输出回CPU，这里使用的是一个RWStructuredBuffer，它是DirectX 11新增的可写缓冲区类型，是通过乱序访问缓冲区视图（unordered access view）建立的。其声明中的寄存器u0指向的就是乱序访问视图。使用这个缓冲区我们就可以让各个线程以顺序无关的方式写入自己的结果。比如我们先利用 DTid的分量计算出该像素在输出缓冲区中的位置，然后用一条赋值语句写入该缓冲。

比起编写这个HLSL，让他运行起来可要麻烦多了。这真是个不幸的事实。我们需要在CPU端创建Direct X的设备，逐个创建各个缓冲区、缓冲区视图并且填充缓冲区。还需要创建拷贝结果的辅助缓冲。我完全是按照Direct X SDK中的例子编写的，否则还真没法轻松的写出来。大家也只好先看SDK的例子学习吧……也许使用OpenCL不需要这么多麻烦的辅助工作。

先让我们看看结果吧！首先是从-2i ~ 2i以及-2 ~ 2区间的图形：