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

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

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

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

首先我們能夠發(fā)現(xiàn)每一個復(fù)數(shù)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，是顯卡里所有線程都能訪問的資源。我們調(diào)用Compute Shader程序時傳的參數(shù)就都要放入Constant Buffer當中。下面就是HLSL的函數(shù)。HLSL沒有棧，因此函數(shù)不能遞歸調(diào)用。所有函數(shù)的參數(shù)和本地變量都儲存在寄存器中。共有4096個寄存器可以用于這個目的。

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

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

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

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

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

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

（點擊看大圖）。

而最大迭代次數(shù)也會對圖形產(chǎn)生很大的影響。下面是同一局部分別迭代最大256次-32768次的不同圖形

（點擊看大圖）

怎么樣，非常神奇吧。好像看到了宇宙誕生的奧秘

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