GPU编程自学5 —— 线程协作

深度学习的兴起，使得多线程以及GPU编程逐渐成为算法工程师无法规避的问题。这里主要记录自己的GPU自学历程。

目录

《GPU编程自学1 —— 引言》

《GPU编程自学2 —— CUDA环境配置》

《GPU编程自学3 —— CUDA程序初探》

《GPU编程自学4 —— CUDA核函数运行参数》

《GPU编程自学5 —— 线程协作》

《GPU编程自学6 —— 函数与变量类型限定符》

《GPU编程自学7 —— 常量内存与事件》

《GPU编程自学8 —— 纹理内存》

《GPU编程自学9 —— 原子操作》

《GPU编程自学10 —— 流并行》

五、 线程协作

5.1 并行程序块的分解

首先回顾我们之前实现的矢量相加程序：

// 核函数：每个线程负责一个分量的加法
__global__ void addKernel(int *c, const int *a, const int *b)
{
int i = threadIdx.x; // 获取线程ID
c[i] = a[i] + b[i];
}

// 运行核函数，运行设置为1个block，每个block中size个线程
addKernel << <1, size >> >(dev_c, dev_a, dev_b);

通过前面小节，我们知道一个Block中的可开辟的线程数量是有限的(不超过1024)。

如果矢量特别长，上面的操作是会出现问题的。于是我们可以采用多个线程块(Block)来解决线程不足的问题。 假如我们设定每个线程块包含128个线程，则需要的线程块的数量为 size / 128。 为了避免不能整除带来的问题，我们可以稍微多开一点 (size + 127) / 128，但需要增加判断条件来避免越界。

// 核函数：每个线程负责一个分量的加法
__global__ void addKernel(int *c, const int *a, const int *b, const int size)
{
int i = blockIdx.x * blockDim.x + threadIdx.x; // 获取线程ID
if (i < size)
c[i] = a[i] + b[i];
}

// 运行核函数，运行设置为多个block，每个block中128个线程
addKernel <<<(size + 127) / 128, 128 >>>(dev_c, dev_a, dev_b, size);

通过前面小节，我们同时也知道一个Grid中可开辟的Block数量也是有限的。

如果数据量大于 Block_num * Thread_num，那么我们就无法为每个分量单独分配一个线程了。 不过，一个简单的解决办法就是在核函数中增加循环。

// 核函数：每个线程负责多个分量的加法
__global__ void addKernel(int *c, const int *a, const int *b, const int size)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
while (i < size)
{
c[i] = a[i] + b[i];
// 偏移分量等于一个Grid中包含的线程数量
i += blockDim.x * gridDim.x;
}
}

// 运行核函数，运行设置为1个Grid包含128个block，每个block包含128个线程
// 其中已经假设 size > 128*128
addKernel <<<128, 128 >>>(dev_c, dev_a, dev_b, size);

5.2 共享内存与同步

上面提到线程块的分解似乎是为了增加可用的线程数量，但这种说法并不靠谱，因为这完全可以由CUDA在幕后全权控制。 其实，分解线程块的重要原因是因为内存。

在“4.3 内存结构”中我们已经知道，同一个Block中的线程可以访问一块共享内存。由于共享内存缓冲区驻留在物理GPU上，而不是GPU之外的系统内存上，因此访问共享内存的延迟要远远低于访问普通缓冲区的延迟。

不同Block之间存在隔离，如果我们需要不同线程之间进行通信，那么还需要考虑线程同步的问题。比如线程A将某个数值写入内存，然后线程B会对该数值进行一些操作，那么显然必须等A完成之后B才可以操作，如果没有同步，程序将会因进入“竞态条件”而产生意想不到的错误。

接下来我们通过一个矢量点积的例子来说明上述问题。

矢量点积的定义如下：

\((x_1,x_2,x_3,x_4)\cdot (y_1,y_2,y_3,y_4)=x_1y_1+x_2y_2+x_3y_3+x_4y_4\)

由上面的定义来看，点积的实现可以分为两步：

（1）计算每个分量的乘积，并暂存该结果；

（2）将所有临时结果加和。

我们先来简单实现下第一步：

// 定义我们的矢量长度
const int N = 64 * 256;

// 定义每个Block中包含的Thread数量
const int threadsPerBlock = 256;

// 定义每个Grid中包含的Block数量, 这里32 < 64， 是为了模拟线程数量不足的情况
const int blocksPerGrid = 32;

__global__ void dot( float *a, float *b, float *c )
{
// 声明共享内存用于存储临时乘积结果，内存大小为1个Block中的线程数量
// PS. 每个Block都相当于有一份程序副本，因此相当于每个Block都有这样的一份共享内存
__shared__ float cache[threadsPerBlock];

// 线程索引
int tid = threadIdx.x + blockIdx.x * blockDim.x;

// 一个Block中的线程索引
int cacheIndex = threadIdx.x;

// 计算分量乘积，同时处理线程不足的问题
float   temp = 0;
while (tid < N) {
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}

// 存储临时乘积结果
cache[cacheIndex] = temp;
}

执行完上面的部分，我们剩下的就是把cache中的值相加求和。 但是，我们必须要保证所有乘积都已经计算完成，才能去计算求和。 命令如下：

// 对线程块中的所有线程进行同步
// 线程块中的所有线程都执行完前面的代码后才会继续往后执行
__syncthreads();

求和最直接的方式就是循环累加，此时复杂度与数组长度成正比。不过我们可以再用一种更加高效的方法，其复杂度与数组长度的log成正比：将值两两合并，于是数据量减小一半，再重复两两合并直至全部计算完成。代码如下：

// 合并算法要求长度为2的指数倍
int i = blockDim.x/2;
while (i != 0)
{
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}

// 最后将一个Block的求和结果进行保存
if (cacheIndex == 0)
c[blockIdx.x] = cache[0];

下面给出完整的代码(简单起见，不再做错误检查)：

#include <iostream>

#include "cuda_runtime.h"

//定义矢量长度
const int N = 64 * 256;

// 定义每个Block中包含的Thread数量
const int threadsPerBlock = 256;

// 定义每个Grid中包含的Block数量, 这里32 < 64， 是为了模拟线程数量不足的情况
const int blocksPerGrid = 32;

// 核函数：矢量点积
__global__ void dot(float* a, float* b, float* c)
{
// 声明共享内存用于存储临时乘积结果，内存大小为1个Block中的线程数量
// PS. 每个Block都相当于有一份程序副本，因此相当于每个Block都有这样的一份共享内存
__shared__ float cache[threadsPerBlock];

// 线程索引
int tid = blockIdx.x * blockDim.x + threadIdx.x;

// 一个Block中的线程索引
int cacheIndex = threadIdx.x;

// 计算分量乘积，同时处理线程不足的问题
float temp = 0.0f;
while (tid < N)
{
temp += a[tid] * b[tid];
tid  += gridDim.x * blockDim.x;
}

// 存储临时乘积结果
cache[cacheIndex] = temp;

// 对线程块中的所有线程进行同步
// 线程块中的所有线程都执行完前面的代码后才会继续往后执行
__syncthreads();

// 合并算法要求长度为2的指数倍
int i = threadsPerBlock / 2;
while (i != 0)
{
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}

if (cacheIndex == 0)
c[blockIdx.x] = cache[0];
}

int main()
{
// 在主机端创建数组
float a
;
float b
;
float c[threadsPerBlock];
for (size_t i = 0; i < N; i++)
{
a[i] = 1.f;
b[i] = 1.f;
}

// 申请GPU内存
float* dev_a = nullptr;
float* dev_b = nullptr;
float* dev_c = nullptr;
cudaMalloc((void**)&dev_a, N * sizeof(float));
cudaMalloc((void**)&dev_b, N * sizeof(float));
cudaMalloc((void**)&dev_c, blocksPerGrid * sizeof(float));

//将数据从主机copy进GPU
cudaMemcpy(dev_a, a, N * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, N * sizeof(float), cudaMemcpyHostToDevice);

//进行点积计算
dot<<<32, 256>>>(dev_a, dev_b, dev_c);

//将计算结果copy回主机
cudaMemcpy(c, dev_c, blocksPerGrid * sizeof(float), cudaMemcpyDeviceToHost);

//将每个block的结果进行累加
for (size_t i = 1; i < blocksPerGrid; i++)
c[0] += c[i];

// 输出结果
std::cout << "The ground truth is 16384, our answer is " << c[0] << std::endl;

//释放内存
cudaFree(dev_a); cudaFree(dev_b); cudaFree(dev_c);

system("pause");
return 0;
}

参考资料：

《CUDA by Example: An Introduction to General-Purpose GPU Programming》 中文名《GPU高性能编程CUDA实战》

深入理解CUDA点积运算 https://segmentfault.com/a/1190000007519944