nBodyCS<I>學習筆記之計算著色器

Nvidia-SDK（1）

版權聲明：本文為博主原創文章，未經博主允許不得轉載。

DirectX一直是Windows上圖形和遊戲開發的核心技術。DirectX提供了一種在顯卡上運行的程序——著色器（Shader）。從DirectX11開始，DirectX增加了一種計算著色器（Compute Shader），它是專門為與圖形無關的通用計算設計的。因此DirectX就變成了一個通用GPU計算的平臺。鑒於GPU擁有極其強大的並行運算能力，學習使用DirectCompute是很有意義的。基本上，DirectCompute需要通過計算著色器5.0（Compute Shader）編程模型（即CS 5.0）才能完全實現。然而CS 5.0需要DirectX 11硬件才能支持，本文默認機器支持DirectX11硬件的。

1.計算著色器創建步驟
1.1初始化設備和上下文

1.2從HLSL文件加載著色器程序並編譯

1.3為著色器創建並初始化資源（如緩沖區）

1.4設定著色器狀態，並執行

1.5取回運算結果

詳細的創建過程已經在高斯模糊一篇詳細介紹，本文大概簡單闡述下：

成功運行後，這段代碼將產生一個設備指針，一個上下文指針還有一個特性等級的Flag。

註意：為簡單起見以上代碼省略了許多變量聲明的代碼。完整示例代碼需要補上這些代碼。這裏的代碼片段僅用來展示程序中發生的事情。

選擇要用的顯卡

用IDXGIFactory對象即可枚舉系統中安裝的顯卡，如下面代碼所示。首先創建一個IDXGIFactory對象，然後調用EnumAdapters並傳入一個代表正在枚舉顯卡的整數。如果不存在，它會返回DXGI_ERROR_NOT_FOUND。

接下來，在調用D3DCreateDevice創建設備的時候從第一個參數傳入想用的顯卡適配器指針，並且將驅動類型設為D3D_DRIVER_TYPE_UNKNOWN。詳細信息請參見D3D11文檔中D3DCreateDevice函數的幫助。

運行計算著色器

著色器（Shader）是在顯卡上運行的程序，它並不同於CPU上執行的程序，可以用HLSL來編寫（見後文）。

DirectCompute程序中的計算著色器是通過Dispatch函數執行的：

以上語句分派了16x16個線程組。

註意，著色器的輸入通常考慮成“狀態”。就是說你應當在分派著色器程序之前設定狀態，而一旦分派了，“狀態”決定輸入變量的值。所以著色器分派代碼通常應該像這樣：

以上所有常量緩沖（constant buffer），緩沖等可以在著色器程序中看到東東都是在分派線程之前通過調用CSSet…()設定的

與CPU進行同步

請註意上面所有的調用都是異步的。CPU方面總是會立即返回然後才具體執行。如果有必要，其後調用的緩沖區“映射”操作（詳見下文緩沖區部分）時CPU的調用線程才會停下來等待所有異步操作的完成。

事件：基本剖析和同步操作

DirectCompute提供一種基於“查詢”的事件機制API。你可以創建、插入並等待特定狀態的查詢來判斷著色器（或其他異步調用）具體在何時執行。下面的例子創建了一個查詢，然後通過等待查詢來確保運行到某一點時所有該執行的操作都已經執行了，再分派著色器，最後等待另一個查詢並確認著色器程序已執行完畢。

創建查詢對象：

然後在一系列調用中插入“籬笆”，再等待之。如果查詢的信息不存在，GetData()將返回S_FALSE。

最後用這條語句釋放查詢對象：

請小心創建和釋放查詢對象以免弄出太多的查詢來（特別是你處理一幀畫面的時候）。

點擊此處查看MSDN關於查詢的幫助

DirectCompute中的資源

譯註：資源是指可以被GPU或CPU訪問的數據，是著色器的輸入與輸出。包括緩沖區和紋理等類型。

DirectX中資源是按照以下步驟創建出來的：

1.首先創建一個資源描述器，用來描述所要創建的資源。資源描述器是一種內含許多Flag和所需資源信息的結構體。

2.調用某種Create系方法，傳入描述器作為參數並創建資源。

CPU與GPU之間的通訊

gD3DContext->CopyResouce()函數可以用來讀取或復制資源。這裏復制是指兩個資源之間的復制。如果要在CPU和GPU之間（譯註：就是在內存和顯存之間）進行復制的話，先要創建一個CPU這邊的“中轉”資源。中轉資源可以映射到CPU的內存指針上，這樣就可以從中轉資源中讀取數據或者復制數據。之後解除中轉資源的映射，再用CopyResource()方法進行與GPU之間的復制。

CPU與GPU之間緩沖區復制的性能

CUDA-C語言（CUDA是nVidia的GPU通用計算平臺）可以分配定址（pinned）宿主指針和寫入聯合（write combined）宿主指針，通過它們可以進行性能最佳的GPU數據復制。而在DirectCompute中，緩沖區的“usage”屬性決定了內存分配的類型和訪問時的性能。

· D3D11_USAGE_STAGING 這種usage的資源是系統內存，可以直接由GPU進行讀寫。但是他們僅能用作復制操作（CopyResource(),CopySubresourceRegion()）的源或目標，而不能直接在著色器中使用。

· 如果資源創建的時候指定了D3D11_CPU_ACCESS_WRITE flag那麽從CPU到GPU復制的性能最佳。

· 如果用了D3D11_CPU_ACCESS_READ該資源將是一個由CPU緩存的資源，性能較低（但是支持取回操作）

· 如果同時指定，READ比WRITE優先。

· D3D11_USAGE_DYNAMIC （僅能用於緩沖區型資源，不能用於紋理資源）用於快速的CPU->GPU內存傳輸。這種資源不但可以作為復制和源和目標，還可以作為紋理（用D3D的術語說，叫做著色器資源視圖ShaderResourceView）在著色器中讀取。但是著色器不能寫入這種資源。這些資源的版本由驅動程序來控制，每次你用DISCARD flag映射內存的時候，如果這塊內存還在被GPU所使用，驅動程序就會產生一塊新的內存來，而不會等GPU的操作結束。它的意義在於提供一種流的方式將數據輸送到GPU。

二．示例代碼《nBodyCS》
//-----------------------------------------------------------------------------

// Reset the body system to its initial configuration

//-----------------------------------------------------------------------------

HRESULT NBodySystemCS::resetBodies(BodyDataconfigData)

{

HRESULThr =S_OK;

m_numBodies= configData.nBodies;

//for compute shader on CS_4_0, we can only have a single UAV per shader, so wehave to store particle

//position and velocity in the same array: all positions followed by all velocities

D3DXVECTOR4*particleArray =newD3DXVECTOR4[m_numBodies* 3];

for(unsignedinti=0; i < m_numBodies; i++){

particleArray[i] =D3DXVECTOR4(configData.position[i*3 + 0],

configData.position[i*3 +1],

configData.position[i*3 +2],

1.0);

particleArray[i +m_numBodies]=particleArray[i];

particleArray[i + 2 *m_numBodies]=D3DXVECTOR4(configData.velocity[i*3 +0],

configData.velocity[i*3 + 1],

configData.velocity[i*3 +2],

1.0);

}

//------------------------------------------------------------------------------------------------------

結構化緩沖區和亂序訪問視圖

ComputeShader的一個很重要的特性是結構化緩沖區和亂序訪問視圖。結構化緩沖區（structuredbuffer）在計算著色器中可以像數組一樣訪問。任意線程可以讀寫任意位置（即並行程序的散發scatter和收集gather動作）。亂序訪問視圖（unordered access view，UAV）是一種將調用方創建的資源綁定到著色器中的機制，並且允許……亂序訪問。

聲明結構化緩沖區

我們可以用D3D11_RESOURCE_MISC_BUFFER_STRUCTURED來創建結構化緩沖區。下面指定的綁定flag表示允許著色器亂序訪問。下邊采用的默認usage表示它可以被GPU進行讀寫，但需要復制到中轉資源當中才能被CPU讀寫。

//------------------------------------------------------------------------------------------------------

D3D11_SUBRESOURCE_DATAinitData = {particleArray,0, 0 };

// 創建結構化緩沖區

D3D11_BUFFER_DESCsbDesc;

sbDesc.BindFlags =D3D11_BIND_UNORDERED_ACCESS|D3D11_BIND_SHADER_RESOURCE;

sbDesc.CPUAccessFlags = 0;

sbDesc.MiscFlags =D3D11_RESOURCE_MISC_BUFFER_STRUCTURED;

sbDesc.StructureByteStride =sizeof(D3DXVECTOR4);

sbDesc.ByteWidth =sizeof(D3DXVECTOR4) *m_numBodies* 3;

sbDesc.Usage =D3D11_USAGE_DEFAULT;

V_RETURN(m_pd3dDevice->CreateBuffer(&sbDesc, &initData,&m_pStructuredBuffer) );

//create the Shader Resource View (SRV) for the structured buffer

D3D11_SHADER_RESOURCE_VIEW_DESCsbSRVDesc;

sbSRVDesc.Buffer.ElementOffset = 0;

sbSRVDesc.Buffer.ElementWidth =sizeof(D3DXVECTOR4);

sbSRVDesc.Buffer.FirstElement = 0;

sbSRVDesc.Buffer.NumElements =m_numBodies* 3;

sbSRVDesc.Format =DXGI_FORMAT_UNKNOWN;

sbSRVDesc.ViewDimension =D3D11_SRV_DIMENSION_BUFFER;

V_RETURN(m_pd3dDevice->CreateShaderResourceView(m_pStructuredBuffer, &sbSRVDesc,&m_pStructuredBufferSRV) );

聲明亂序訪問視圖

下面我們聲明一個亂序訪問視圖。註意需要給他一個結構化緩沖區的指針

// 創建一個亂序訪問視圖，指向結構化緩沖區

D3D11_UNORDERED_ACCESS_VIEW_DESCsbUAVDesc;

sbUAVDesc.Buffer.FirstElement = 0;

sbUAVDesc.Buffer.Flags = 0;

sbUAVDesc.Buffer.NumElements =m_numBodies* 3;

sbUAVDesc.Format =DXGI_FORMAT_UNKNOWN;

sbUAVDesc.ViewDimension =D3D11_UAV_DIMENSION_BUFFER;

V_RETURN(m_pd3dDevice->CreateUnorderedAccessView(m_pStructuredBuffer, &sbUAVDesc,&m_pStructuredBufferUAV) );

delete[] particleArray;

returnhr;

}

之後，在分派著色器線程之前，我們需要激活著色器使用的結構化緩沖：

分派線程之後，如果使用CS 4.x硬件，一定要將其解除綁定。因為CS4.x每條渲染流水線僅支持綁定一個UAV。

三．DirectCompute中的常量緩沖

常量緩沖（constantbuffer）是一組計算著色器運行時不能更改的數據。用作圖形程序是，常量緩沖可以是視角矩陣或顏色常量。在通用計算程序中，常量緩沖可以存放諸如信號過濾的權重和圖像處理的說明等數據。

如果要使用常量緩沖：

· 創建緩沖區資源

· 用內存映射的方式初始化數據（也可以用效果接口）

· 用CSSetConstantBuffers設定常量緩沖的值

下面代碼創建了三個常量緩沖。註意常量緩沖的尺寸，這裏我們知道在HLSL中它是一個四元矢量。

// 創建常量緩沖

D3D11_BUFFER_DESC cbDesc;

cbDesc.Usage =D3D11_USAGE_DYNAMIC;

cbDesc.BindFlags=D3D11_BIND_CONSTANT_BUFFER;

// CPU 可寫, 這樣我們可以每幀更新數據

cbDesc.CPUAccessFlags=D3D11_CPU_ACCESS_WRITE;

cbDesc.MiscFlags= 0;

cbDesc.ByteWidth=sizeof(CB_DRAW);

V_RETURN( pd3dDevice->CreateBuffer( &cbDesc,NULL, &m_pcbDraw) );

cbDesc.ByteWidth=sizeof(CB_UPDATE);

V_RETURN( pd3dDevice->CreateBuffer( &cbDesc,NULL, &m_pcbUpdate) );

cbDesc.Usage =D3D11_USAGE_IMMUTABLE;

cbDesc.BindFlags=D3D11_BIND_CONSTANT_BUFFER;

cbDesc.CPUAccessFlags= 0;

cbDesc.ByteWidth=sizeof(CB_IMMUTABLE);

D3D11_SUBRESOURCE_DATA initData= { &cbImmutable, 0, 0 };

V_RETURN( pd3dDevice->CreateBuffer(&cbDesc,&initData, &m_pcbImmutable) );

接下來用內存映射的方式將數據發送到常量緩沖。通常程序員會在CPU程序裏定義和HLSL一樣的結構體，因此會用sizeof取得尺寸，然後將緩沖區的指針映射到結構體來填充數據

// 必須用 D3D11_MAP_WRITE_DISCARD

D3D11_MAPPED_SUBRESOURCEmappedResource;

V( m_pd3dImmediateContext->Map( m_pcbUpdate,0, D3D11_MAP_WRITE_DISCARD, 0, &mappedResource ) );

CB_UPDATE* pcbUpdate= (CB_UPDATE*)mappedResource.pData;

pcbUpdate->g_timestep=dt;

pcbUpdate->g_softeningSquared = 0.01f;

pcbUpdate->g_numParticles=m_numBodies;

pcbUpdate->g_readOffset =m_readBuffer*m_numBodies;

pcbUpdate->g_writeOffset = (1 -m_readBuffer)*m_numBodies;

m_pd3dImmediateContext->Unmap(m_pcbUpdate,0);

註意計算著色器的輸入變量（在這裏就是常量緩沖）是當成“狀態”變量的，因此在分派計算著色器之前需要用CSSetShader()函數設置狀態。這樣計算著色器執行的時候就能訪問到這些變量

// 在計算著色器中激活

m_pd3dImmediateContext->CSSetShader(m_pCSUpdatePositionAndVelocity,NULL, 0 );

m_pd3dImmediateContext->CSSetConstantBuffers( 0, 1, &m_pcbUpdate );

// Run the CS

m_pd3dImmediateContext->Dispatch(m_numBodies/ 256, 1, 1 );

最後當計算著色器運行的時候，m_pCSUpdatePositionAndVelocity所指向的著色器就可以訪問這個m_pcbUpdate常量緩沖

4.計算著色器（CS）HLSL編程

運行在顯卡上的計算著色器是用HLSL（High Level Shader Language 高級著色器語言）寫成的。在我們的例子中它是以文本形式存在，並且在運行時動態編譯的。計算著色器是一種單一程序被許多線程並行執行的程序。這些線程分成多個“線程組”，在線程組內的線程之間可以共享數據或互相同步。

GPU硬件架構

GPU硬件結構主要由以下幾個關鍵模塊組成：內存（全局的，常量的，共享的）；流處理器簇（SM）；流處理器(SP)。如圖所示：

SP: 最基本的處理單元，streamingprocessor 最後具體的指令和任務都是在sp上處理的。GPU進行並行計算，也就是很多個sp同時做處理

SM:多個sp加上其他的一些資源組成一個sm, streaming multiprocessor. 其他資源也就是存儲資源，共享內存，寄儲器等。

WARP:GPU執行程序時的調度單位，目前cuda的warp的大小為32，同在一個warp的線程，以不同數據資源執行相同的指令。

GPU實際上是一個SM的陣列，每個SM包含N個核。一個GPU設備中包含一個或多個SM。SM內部組成結構圖如下所示：

4.1線程網格

一個線程網格是由若幹線程塊組成的，每個線程塊是二維的，擁有X軸Y軸。此時我們最多能開啟Y*X*T個線程。thread-->block-->grid：在利用cuda進行編程時，一個grid分為多個block，而一個block分為多個thread。其中任務劃分到是否影響最後的執行效果。劃分的依據是任務特性和GPU本身的硬件特性。一個sm只會執行一個block裏的warp，當該block裏warp執行完才會執行其他block裏的warp。進行劃分時，最好保證每個block裏的warp比較合理，那樣可以一個sm可以交替執行裏面的warp，從而提高效率，此外，在分配block時，要根據GPU的sm個數，分配出合理的block數，讓GPU的sm都利用起來，提利用率。分配時，也要考慮到同一個線程block的資源問題，不要出現對應的資源不夠。

假設我們在看一張高清的圖片，這張圖片的分辨率為1920 x1080。通常線程塊中的線程數量最好是一個線程束大小的整數倍，即32的整數倍。本例中我們在線程塊上開啟192個線程。每個線程塊192個線程，很容易計算出一行圖形需要10個線程塊（如圖4.1）。在這裏選在192這個是因為X軸方向處理的數據大小1920是它的整數倍，192又是線程束大小的整數倍。GPU上的一個線程束的大小是32（英偉達公司保留著對這個參數修改的權利），他們提供一個固有變量-WRAPSIZE，我們可以通過這個變量來獲取硬件支持的線程束的大小。

圖4.1按行分布的線程塊

在X軸方向的頂部我們可以得到線程的索引，在Y軸方向我們可以得到行號。由於每一行只處理一行像素，每一行有10個線程塊，因此我們需要1080行來處理整張圖片，一共1080*10=10800個線程塊。在費米架構的硬件上，一個SM可以處理8個線程塊，所以從應用層角度來說一共需要1350個（總共10800個線程塊 / 每個SM能調度的8個線程塊）SM來完全實現並行。但當前費米架構的硬件只有16個SM可供使用（GTx580）即每個SM將被分配675個線程塊進行處理。

上述例子很簡單，數據分布整齊容易理解。但是往往我們的數據可能不是一維的，這時我們可以使用二維模塊或者三維矩陣來存儲數據。

4.2網格（Grid）、線程塊（Block）和線程（Thread）的組織關系

CUDA的軟件架構由網格（Grid）、線程塊（Block）和線程（Thread）組成，相當於把GPU上的計算單元分為若幹（2~3）個網格，每個網格內包含若幹（65535）個線程塊，每個線程塊包含若幹（512）個線程，三者的關系如下圖：

Thread，block，grid是CUDA編程上的概念，為了方便程序員軟件設計，組織線程。

· thread：一個CUDA的並行程序會被以許多個threads來執行。

· block：數個threads會被群組成一個block，同一個block中的threads可以同步，也可以通過shared memory通信。

· grid：多個blocks則會再構成grid。

網格（Grid）、線程塊（Block）和線程（Thread）的最大數量

CUDA中可以創建的網格數量跟GPU的計算能力有關，可創建的Grid、Block和Thread的最大數量參看以下表格：

在單一維度上，程序的執行可以由多達3*65535*512=100661760（一億）個線程並行執行，這對在CPU上創建並行線程來說是不可想象的。

線程索引的計算公式

一個Grid可以包含多個Blocks，Blocks的組織方式可以是一維的，二維或者三維的。block包含多個Threads，這些Threads的組織方式也可以是一維，二維或者三維的。

CUDA中每一個線程都有一個唯一的標識ID—ThreadIdx，這個ID隨著Grid和Block的劃分方式的不同而變化，這裏給出Grid和Block不同劃分方式下線程索引ID的計算公式。

1、 grid劃分成1維，block劃分為1維

int threadId = blockIdx.x *blockDim.x + threadIdx.x;

2、 grid劃分成1維，block劃分為2維

int threadId = blockIdx.x * blockDim.x * blockDim.y+ threadIdx.y* blockDim.x + threadIdx.x;

3、 grid劃分成1維，block劃分為3維

int threadId = blockIdx.x * blockDim.x * blockDim.y * blockDim.z
+ threadIdx.z * blockDim.y * blockDim.x
+ threadIdx.y * blockDim.x + threadIdx.x;

4、 grid劃分成2維，block劃分為1維

int blockId = blockIdx.y * gridDim.x + blockIdx.x;
int threadId = blockId * blockDim.x + threadIdx.x;
5、 grid劃分成2維，block劃分為2維

int blockId = blockIdx.x + blockIdx.y * gridDim.x;
int threadId = blockId * (blockDim.x * blockDim.y)
+ (threadIdx.y * blockDim.x) + threadIdx.x;
6、 grid劃分成2維，block劃分為3維

int blockId = blockIdx.x + blockIdx.y * gridDim.x;
int threadId = blockId * (blockDim.x * blockDim.y * blockDim.z)
+ (threadIdx.z * (blockDim.x * blockDim.y))
+ (threadIdx.y * blockDim.x) + threadIdx.x;
7、 grid劃分成3維，block劃分為1維

int blockId = blockIdx.x + blockIdx.y * gridDim.x
+gridDim.x * gridDim.y * blockIdx.z;
int threadId = blockId * blockDim.x + threadIdx.x;


8、 grid劃分成3維，block劃分為2維

int blockId = blockIdx.x + blockIdx.y * gridDim.x
+gridDim.x * gridDim.y * blockIdx.z;
int threadId = blockId * (blockDim.x * blockDim.y)
+ (threadIdx.y * blockDim.x) + threadIdx.x;


9、 grid劃分成3維，block劃分為3維

int blockId = blockIdx.x + blockIdx.y * gridDim.x
+gridDim.x * gridDim.y * blockIdx.z;
int threadId = blockId * (blockDim.x * blockDim.y * blockDim.z)
+ (threadIdx.z * (blockDim.x * blockDim.y))
+ (threadIdx.y * blockDim.x) + threadIdx.x;
4.3 nBodyCS.hlsl

m_pd3dImmediateContext->Dispatch( m_numBodies / 256,1, 1 );

Dispatch(m_numBodies/ 256, 1, 1)創建了（m_numBodies/ 256）x 1 x 1個線程組。而每一個線程組中的線程是在著色器代碼中用這個語法來指定[numthreads(BLOCK_SIZE,1,1)]

/* 這表示線程組中的線程數，本例中是BLOCK_SIZE x1x1 = 256個線程

[numthreads(BLOCK_SIZE,1,1)]

void NBodyUpdate(uint threadId : SV_GroupIndex,

uint3 groupId : SV_GroupID,

uint3 globalThreadId :SV_DispatchThreadID)

{

float4pos = particles[g_readOffset + globalThreadId.x];

float4 vel = particles[2 * g_numParticles + globalThreadId.x];

//compute acceleration

float3 accel = computeBodyAccel(pos, threadId, groupId);

//Leapfrog-Verlet integration of velocity and position

vel.xyz+= accel * g_timestep;

pos.xyz+= vel * g_timestep;

particles[g_writeOffset+ globalThreadId.x] = pos;

particles[2* g_numParticles + globalThreadId.x] = vel;

}

// Computes the total acceleration on the body with position myPos

// caused by the gravitational attraction of all other bodies in

// the simulation

float3 computeBodyAccel(float4 bodyPos, uint threadId, uint blockId)

{

float3 acceleration = {0.0f, 0.0f, 0.0f};

uint p = BLOCK_SIZE;

uint n = g_numParticles;

uint numTiles = n / p;

for (uint tile = 0; tile < numTiles; tile++)

{

// 取所有線程組中相同索引的threadid 索引threadid 【0,255】範圍

sharedPos[threadId] = particles[g_readOffset + tile * p + threadId];

// 阻止組中所有線程的執行，直到所有組共享訪問完成，組中的所有線程都

//已達到此調用。

GroupMemoryBarrierWithGroupSync();

//計算 bodyPos同相同threakid不同的線程組的粒子坐標的重力引力影響值

acceleration = gravitation(bodyPos, acceleration);

GroupMemoryBarrierWithGroupSync();

}

return acceleration;

}

nBodyCS<I>學習筆記之計算著色器