线程网络

cuda将线程抽象为grid、block、thread三个层次，构成两种视图：

视图1：一个device就是一个grid，grid的最小元素是block，一个grid由若干个block组成。
cuda通过一个dim3的变量描述一个grid里面的block的排列方式。一个grid可以是一维、二维、三维矩阵。

struct __device_builtin__ dim3
{
    unsigned int x, y, z;
#if defined(__cplusplus)
    __host__ __device__ dim3(unsigned int vx = 1, unsigned int vy = 1, unsigned int vz = 1) : x(vx), y(vy), z(vz) {}
    __host__ __device__ dim3(uint3 v) : x(v.x), y(v.y), z(v.z) {}
    __host__ __device__ operator uint3(void) { uint3 t; t.x = x; t.y = y; t.z = z; return t; }
#endif /* __cplusplus */
};

假设有一个二维的grid，描述这个grid的参数包括：

//cuda的内建变量
//当我们执行一个kernel函数时，cuda自动帮我们赋值了
gridDim.x   //x方向的block数量
gridDim.y   //y方向的block数量

blockIdx.x  //x方向的block索引
blockIdx.y  //y方向的block索引

其中，gridDim用于访问线性内存时，将二维索引转换为一维索引，即当我们知道thread在x方向上索引idx、y方向上的索引idy，通过gridDim.x*blockDim.x*idy+idx求出threadID；如果grid是三维的，则通过gridDim.x*blockDim.x*idy + gridDim.x*blockDim.x*idy + idx

视图2：一个block，block的最小元素是thread，一个block由若干个线程组成。同样的，cuda通过一个dim3变量描述一个block内部的thread排列形式。一个block可以是一维、二维、三维矩阵。
假设有一个二维的block，描述这个block的参数包括：

//cuda的内建变量
blockDim.x   //x方向的thread数量
blockDim.y   //y方向的thread数量

threadIdx.x  //x方向的thread索引
threadIdx.y  //y方向的thread索引

gridDim、blockDim
在访问线性内存时，需要将多维索引转换为一维索引，如对于一个block，二维索引(x,y)可以转换为(x+y*blockDim.x)，三维索引(x,y,z)可以转换为(x+y*blockDim.x+z*blockDim.x*blockDim.y)

举例

现在有一个grid，包含40*40个block，每个block包含16*16个thread

dim3 gridSize(40,40);
dim3 blockSize(16,16);

假设要完成对单精度浮点数据的处理

float* src;
float* dst;
int dataLen = gridSize.x*blockSize.x 
    * gridSize.y*blockSize.y * sizeof(float);
cudaMalloc(&src,dataLen);  //输入缓冲
cudaMalloc(&dst,dataLen);  //输出缓冲
myKernel<<<gridSize,blockSize>>>(
    float* src,float* dst);
...
cudaFree(src);
cudaFree(dst);  

其中myKernel定义如下

__global__ void 
    myKernel(float* src,float* dst)
{
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    int idy = blockIdx.y*blockDim.y + threadIdx.y;
    int thread_id = gridDim.x*blockDim.x*idy+x;  //二维索引转一维索引
    dst[idx] = src[thread_id]*0.5;   //通过一维索引访问线性内存
}

如果数据是一个二维数组

float src[640][640];
float dst[640][640];
__global__ void 
    myKernel(float* src,float* dst)
{
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    int idy = blockIdx.y*blockDim.y + threadIdx.y;
    dst[idy][idx] = src[idy][idx]*0.5;   //通过二维索引访问数组
}

当然，实际中不会将数据定义为数组，因为栈的资源是有限的，这里仅仅为了说明threadID的问题

线程分配

原则：

• 1、coalesced，确保对数据的合并访问，通过合理设置grid的x方向上的线程数，使每个线程中访问的数据大小是cache大小的整数倍，这样能够减少访存操作，优化访存效率（前提条件：数据本身是内存对齐的）。
  参考：《A developer’s guide to parallel …》Shane Cook一书的6.6和9.2.1
• 2、尽量保证block的数目是SM数目的整数倍，避免出现剩余线程块，剩余线程块会降低gpu的利用率
  参考：《A developer’s guide to parallel …》Shane Cook一书的5.6
• 3、合理设置block的thread数目，不能过少，也不能过多，从而使每个SM的利用率达到最大。当线程数满足一定条件，利用率可以达到100%
  参考：《A developer’s guide to parallel …》Shane Cook一书的5.5.2

cuda对每个SM能够容纳的最大block数目、最大thread数目有限制，具体参照cuda的compute capability列表。因此，我们需要根据实际gpu的compute capability，来分配线程数目。
Maximum number of resident blocks per multiprocessor - 每个SM能够容纳的最大block数目
Maximum number of resident threads per multiprocessor - 每个SM能够容纳的最大thread数目
注：这里指的是resident，即每一时刻SM中同时运行的thread或者block数目，并不是指我们能够分配的数目（上文所说的gridSize、blockSize），在计算能力3.x的设备上，我们甚至能够分配2147483648个block！但硬件中并不能同时运行这么多个block，每个SM能够容纳的最大block数目是有限的！

计算硬件利用率的方法： (每个block的线程数 × 每个SM能够容纳的最大block数目)/每个SM能够容纳的最大thread数目 × 100%

要确保利用率为100%，其必要条件是每个block的线程数 × 每个SM的block数目=每个SM能够容纳的最大thread数目

例如，GTX760，其compute capability为3.0，每个SM能够容纳的最大block数目为16，每个SM能够容纳的最大thread数目为2048。
令每个SM的block数目为16，要使利用率最大，则需要令每个block的线程数=2048/16=128，即blockSize.x*blockSize.y=128，可行的一个方案是blockSize（32，4）

上面说了，这是必要条件，因为每个SM并一定能够运行其最大的block数目，如上面例子的16，原因如下：

一旦我们的kernel中存在同步操作，则较早到达同步点的thread，需要等待其它还没运行到同步点的thread，只有所有thread都运行至同步点，程序才能继续向下执行。在这种情形下，执行一个block的时间（或者latency）是不确定的。而SM是以block为单位执行代码的，每个block只有当它执行完全部指令才会从SM撤走，但如果block在等待某一个线程（准确地讲应该是warp）执行完成，就会导致SM处于闲置状态，大大降低程序的性能。
由此可见，block包含的线程数越多，就会潜在增加等待执行较慢线程的可能性，从而导致SM阻塞。
参考：《A developer’s guide to parallel …》Shane Cook一书的5.5.2

为了实际的利用率达到最优，我们需要实验，不同的block数目、不同的block线程数量，哪一种组合使程序的耗时最低。