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

时间:2009-12-13   来源:博客园   网友评论:0   人气: 313 作者: 装配中的脑袋

近年来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,我们采用以下公式:

clip_image002

进行迭代,所得的所有复数的集合就叫曼德勃罗特集,其中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);
}


 

文章评论