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*128addKernel <<<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[N];    float b[N];    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;}

参考资料：

1. 《CUDA by Example: An Introduction to General-Purpose GPU Programming》 中文名《GPU高性能编程CUDA实战》
2. 深入理解CUDA点积运算 https://segmentfault.com/a/1190000007519944