nBodyCS<I>学习笔记之计算着色器

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;

}