您的位置：UltraLAB圖形工作站方案網(wǎng)站 > GPU計(jì)算 > GPU開發(fā)工具 > 用GPU通用并行計(jì)算繪制曼德勃羅特集圖形上篇

用GPU通用并行計(jì)算繪制曼德勃羅特集圖形上篇

時(shí)間：2009-12-13 13:16:00 來源：UltraLAB圖形工作站方案網(wǎng)站 人氣：8336 作者：admin

近年來PC的計(jì)算能力發(fā)生了天翻地覆的變化。CPU逐漸趨向于多核發(fā)展，同時(shí)內(nèi)存帶寬和緩存不斷增加，如今的PC已經(jīng)成為小型的統(tǒng)一地址空間的并行計(jì)算機(jī)。然而我們的PC中還有一個(gè)設(shè)備可以提供比CPU更加強(qiáng)大的并行計(jì)算設(shè)備——顯卡，它在進(jìn)行充分并行的任務(wù)時(shí)可以提供高達(dá)數(shù)TFLOPS的峰值運(yùn)算能力，這幾乎是2000-2001年間國產(chǎn)超級(jí)計(jì)算機(jī)的運(yùn)算能力。在顯卡剛出現(xiàn)時(shí)，顯卡內(nèi)的模塊都是為特定的圖形任務(wù)而設(shè)計(jì)的。比如會(huì)有光照和坐標(biāo)變換單元以及光柵單元等。隨著顯卡和圖形技術(shù)的發(fā)展，在DirectX 8出現(xiàn)之后顯卡開始允許在這些特殊單元處理的基礎(chǔ)上增加編程控制。比如光照和坐標(biāo)變換階段可以用頂點(diǎn)著色器（Vertex Shader）進(jìn)行控制，而光柵后的階段則可以使用像素著色器（Pixel Shader）進(jìn)行控制。這些“著色器（Shader）”其實(shí)就是運(yùn)行在顯卡上的程序。在DirectX 10階段之后，整個(gè)渲染階段可編程控制階段大大增強(qiáng)，整個(gè)渲染流水線中各個(gè)階段都可以通過著色器語言進(jìn)行控制。因此DirectX 10后的顯卡，都將執(zhí)行著色器的運(yùn)算單元從渲染過程中獨(dú)立出來，成為一個(gè)統(tǒng)一著色單元，這些統(tǒng)一著色單元是可以大規(guī)模并行執(zhí)行的運(yùn)算單元。人們逐漸發(fā)現(xiàn)這些統(tǒng)一著色單元實(shí)際上可以進(jìn)行與圖形技術(shù)無關(guān)的通用技術(shù)。顯卡的性能幾乎按照三倍于摩爾定律的速度發(fā)展，GPU用于通用計(jì)算的潛力也非常巨大。比如本文編寫時(shí)最強(qiáng)的GPU芯片——Ati Radeon HD5870可以提供2.72TFLOPS的單精度峰值浮點(diǎn)運(yùn)算性能和500GFLOPS左右的雙精度峰值浮點(diǎn)運(yùn)算性能。而時(shí)下最先進(jìn)的酷睿i7四核處理器只能提供60-80GFLOPS的峰值浮點(diǎn)運(yùn)算能力。

GPU芯片制造商nVidia和Ati是GPU通用計(jì)算的先行者，他們分別提供了CUDA和Ati Stream技術(shù)用于GPU程序的開發(fā)。隨后由蘋果等廠商領(lǐng)導(dǎo)并獲得諸多廠商支持的OpenCL出臺(tái)，成為GPU通用計(jì)算的統(tǒng)一標(biāo)準(zhǔn)。而微軟也看到了這個(gè)潛力，在DirectX 11中提供了不與具體渲染流程綁定的計(jì)算著色器（Compute Shader）。雖然還叫“著色器”，但是Compute Shader已經(jīng)是完全通用的編程技術(shù)。編程語言方面CUDA采用的編程語言有擴(kuò)展的C、C++和Fortran等，OpenCL采用的是擴(kuò)展的標(biāo)準(zhǔn)C。而DirectX Compute Shader是和C類似的HLSL語言。從開發(fā)難度來講DirectX要比CUDA和OpenCL難用。HLSL沒有指針，而且需要?jiǎng)討B(tài)編譯和執(zhí)行。但是目前Compute Shader 5.0支持雙精度浮點(diǎn)數(shù)、完善的原子操作和三維線程組支持。它還支持32KB的線程組共享內(nèi)存。因此采用Compute Shader進(jìn)行通用計(jì)算仍然是一個(gè)值得研究的領(lǐng)域。我也是DirectX的初學(xué)者，對(duì)它的圖形技術(shù)一竅不通，完全是本著通用計(jì)算的目的來學(xué)習(xí)HLSL。下面就動(dòng)手來實(shí)現(xiàn)一個(gè)非常適合并行計(jì)算的程序——曼德勃羅特集。

曼德勃羅特（Mandelbrot）集是一種復(fù)平面上的點(diǎn)集。對(duì)任意復(fù)數(shù)c，我們采用以下公式：

進(jìn)行迭代，所得的所有復(fù)數(shù)的集合就叫曼德勃羅特集，其中z從0開始。這個(gè)公式如此的簡(jiǎn)單，但效果卻非同一般。我們會(huì)發(fā)現(xiàn)許多點(diǎn)在迭代中處于亞穩(wěn)定的狀態(tài)，它們會(huì)在迭代數(shù)次之后逃逸（發(fā)散）。而這個(gè)迭代次數(shù)對(duì)每一點(diǎn)c而言都是不同的。如果我們畫出復(fù)平面中各個(gè)點(diǎn)的迭代次數(shù)，就會(huì)發(fā)現(xiàn)它能組成難以置信的圖案。

首先我們能夠發(fā)現(xiàn)每一個(gè)復(fù)數(shù)c都能夠獨(dú)立地采用這個(gè)這個(gè)公式進(jìn)行運(yùn)算，它無需與其他的點(diǎn)進(jìn)行任何交互。因此該算法默認(rèn)就是完全并行的。我們只需要按照每一個(gè)輸入進(jìn)行劃分就可以輕易地寫出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);
}

#p#page_title#e#

如你所見，HLSL和C語言十分相似，只是有一些特別的擴(kuò)展。一開始的cbuffer稱為Constant Buffer，是顯卡里所有線程都能訪問的資源。我們調(diào)用Compute Shader程序時(shí)傳的參數(shù)就都要放入Constant Buffer當(dāng)中。下面就是HLSL的函數(shù)。HLSL沒有棧，因此函數(shù)不能遞歸調(diào)用。所有函數(shù)的參數(shù)和本地變量都儲(chǔ)存在寄存器中。共有4096個(gè)寄存器可以用于這個(gè)目的。

ComposeColor是我為了畫出好看的顏色而編寫的，可以自由發(fā)揮，下面跳過它看Calculate，這就是我們的主角。DirectX中每一個(gè)著色器程序都綁定到一個(gè)特定的渲染流程上，因此著色器程序通常都有一個(gè)作用目標(biāo)。例如像素著色器針對(duì)每個(gè)像素執(zhí)行，頂點(diǎn)著色器針對(duì)每個(gè)頂點(diǎn)執(zhí)行。而計(jì)算著色器不與這些圖形概念綁定，計(jì)算著色器針對(duì)的是“線程”。每一個(gè)線程執(zhí)行一次。在這段代碼里[numthreads(16, 16, 1)]是對(duì)線程進(jìn)行分組的attribute。Compute Shader 5.0可以將線程分為三維的線程組。這個(gè)attribute就用來描述每一組的線程數(shù)。我們這個(gè)應(yīng)用實(shí)際不需要進(jìn)行分組，但是分組可以增加處理數(shù)據(jù)的規(guī)模。比如采用1024個(gè)線程的線程組，就能處理多達(dá)67108864個(gè)輸入數(shù)據(jù)。請(qǐng)注意，計(jì)算著色器總是大量線程并行執(zhí)行的狀態(tài)，編寫它的時(shí)候一定要變換思想，不要采用順序編程的思路。

我們的Calculate函數(shù)輸入?yún)?shù)DTid是一個(gè)三維向量，它有一個(gè)顯式用途SV_DispatchThreadID，這表示線程的分布ID（DispatchID），也就是線程組的ID與組內(nèi)線程ID的乘積。我們用它來表示復(fù)平面的輸入點(diǎn)坐標(biāo)。首先我們把輸入點(diǎn)換算成從-2i – 2i和-2 - 2這個(gè)區(qū)間的點(diǎn)坐標(biāo)。這就是我們公式里的c。然后就簡(jiǎn)單了，按照公式迭代就行了。迭代終止條件有兩個(gè)：一是迭代變量z的模大于2（表示這個(gè)點(diǎn)發(fā)散了），二是達(dá)到了最大的迭代次數(shù)MAX_ITER。每個(gè)點(diǎn)迭代次數(shù)具體是多少是無法知道的，不過它不會(huì)大于MAX_ITER。這里進(jìn)行復(fù)數(shù)運(yùn)算時(shí)所用的類型是float2，它是一個(gè)有兩個(gè)float分量的向量。

最后我們要把迭代次數(shù)輸出回CPU，這里使用的是一個(gè)RWStructuredBuffer，它是DirectX 11新增的可寫緩沖區(qū)類型，是通過亂序訪問緩沖區(qū)視圖（unordered access view）建立的。其聲明中的寄存器u0指向的就是亂序訪問視圖。使用這個(gè)緩沖區(qū)我們就可以讓各個(gè)線程以順序無關(guān)的方式寫入自己的結(jié)果。比如我們先利用DTid的分量計(jì)算出該像素在輸出緩沖區(qū)中的位置，然后用一條賦值語句寫入該緩沖。

比起編寫這個(gè)HLSL，讓他運(yùn)行起來可要麻煩多了。這真是個(gè)不幸的事實(shí)。我們需要在CPU端創(chuàng)建Direct X的設(shè)備，逐個(gè)創(chuàng)建各個(gè)緩沖區(qū)、緩沖區(qū)視圖并且填充緩沖區(qū)。還需要?jiǎng)?chuàng)建拷貝結(jié)果的輔助緩沖。我完全是按照Direct X SDK中的例子編寫的，否則還真沒法輕松的寫出來。大家也只好先看SDK的例子學(xué)習(xí)吧……也許使用OpenCL不需要這么多麻煩的輔助工作。

先讓我們看看結(jié)果吧！首先是從-2i ~ 2i以及-2 ~ 2區(qū)間的圖形：

黑色區(qū)域要么是逃逸的點(diǎn)，要么是迭代次數(shù)超過了最大限度的點(diǎn)（穩(wěn)定的點(diǎn)）。剩下的亞穩(wěn)定點(diǎn)組成了神奇的美麗圖形。當(dāng)我們重新計(jì)算局部區(qū)域的時(shí)候，會(huì)驚喜地發(fā)現(xiàn)它有極其復(fù)雜的微觀結(jié)構(gòu)