GPU编程-Thread Hierarchy（3）

1. 如果处理的数据是二维的或者三维的，应该怎么办呢？

针对的，我们可以按照二维或者三维的方式，组织线程。老规矩，先代码、后解释

// Kernel definition
__global__ void MatAdd(float A[N][N], float B[N][N],
float C[N][N])
{
int i = threadIdx.x;
int j = threadIdx.y;
C[i][j] = A[i][j] + B[i][j];
}
int main()
{
...
// Kernel invocation with one block of N * N * 1 threads
int numBlocks = 1;
dim3 threadsPerBlock(N, N);
MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
...
}

线程可以一维、二维或者三维的方式，组织成Block，在上述代码中，我们指定有一个Block，这个Block按照NxN的二维结构进行组织。如果N就是矩阵相应的维度，那么上述代码块完成的功能就是矩阵对应元素相加。

2.受GPU资源的限制，每一个Block所含线程个数有限（一般情况下，最多为1024个），如果矩阵的维度超过了线程个数上限，是不是就计算不了大型矩阵的对应元素相加了呢？

答案是否定的。如果将Block看做一个基本组成单元，Block又可以按照一维、二维或者三维的形式组织成grid。Blcok、grid、thread的关系如下图所示

如果矩阵的维度超过了Block能够包含线程的上限，我们可采取以下方式应对（先代码，后解释）

// Kernel definition
__global__ void MatAdd(float A[N][N], float B[N][N],
float C[N][N])
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
if (i < N && j < N)
C[i][j] = A[i][j] + B[i][j];
}
int main()
{
...
// Kernel invocation
dim3 threadsPerBlock(16, 16);
dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
...
}

在上述代码中，N代表矩阵的维度，每一个Block按照16x16的二维结构组织，这样每一个Block只能够处理大型矩阵一个很小的patch。一般情况下，grid所有的thread是自然是顺序排列的（此时的Block索引可以理解为一种二级索引，一级索引指的是直接索引thread）。上述代码就是先将大型矩阵分解为Block，然后由Block里的线程完成具体的矩阵对应元素相加操作。

“The number of thread blocks in a grid is usually dictated by the size of the data being processed or the number of processors in the system, which it can greatly exceed.”

3. Block是并行执行的，假如所需Block数量超出GPU所能提供的Block的限制，会出现什么情况呢？

如上图所示，grid内的Block根据GPU的具体情况，选择顺序执行或者并行执行。

总结：线程的组织方式既能够匹配GPU硬件又能够处理大量数据，是一种很巧妙的安排。