【项目复盘】CUDA编程——Reduce sum cuda核函数优化

reduce baseline实现

cuda核函数的实现基本流程

1、分配内存

cpu malloc申请内存

cudamalloc申请显存

初始化cpu数据并用cudaMemcpy拷贝到GPU，当然cudaMemcpy也可以将GPU数据拷贝回CPU，取决于拷贝类型：cudaMemcpyHostToDevice 或者 cudaMemcpyDeviceToHost

2、分配block和thread数量

通过dim3定义Grid和Block

dim3有三个成员变量x y z，证明可以定位三维的size，也证实后面的gridDim.x gridDim.y gridDim.z （x y z有多少个block），blockDim.x blockDim.y blockDim.z （x y z有多少个thread）

3、cuda核函数

global 关键字代表的就是核函数，主机上CPU调用 设备GPU上执行

host 由CPU调用且在CPU执行的函数

device 设备执行的程序，device前缀定义的函数只能在GPU上执行，所以device修饰的函数里面不能调用一般常见的函数

reduce_v0普通并行版本

baseline的实现只给了1个block和1个thread，其实相当于是串行程序，那v0版本就实现一次并行，并且在block层面利用上shared memory

首先核函数的设计有个规则：

1、设计算法按照Block设计，但是写程序按照thread写

2、索引最好为每一个block设计一个更进一步的索引

3、绘制清晰的block算法图很有必要

核函数内部index应该有两个，一个是内部block的id，也就是tid=threadIdx.x；再一个就是所有block范围内的全局id，也就是gtid=blockIdx.x * blockDimx.x + threadIdx.x

当然在v0版本里gtid=blockIdx.x * blockSize + threadIdx.x，这里blockSize是模版参数，这也是普遍的写法，因为用于静态shared memory的申请需要传入编译期常量指定大小（L10)

在核函数里__shared__ 就是分配了共享内存，另外注意一个点，共享内存的读写是否应该要加上__syncthreads，因为对于后面做运算是要同步block内所有线程都要完成赋值正常运行完的

这张图也代表reduce_v0 block内部的算法图，两两相加最后block内所有值的sum都在第0个线程上，但是要记住既然用到了shared memory就只是在block层面，要在所有block范围的话在最后还应该要对所有的block的0号元素做reduce，我的这个版本是在最后checkresult的时候再把d_out的所有元素全加起来（注意除每个blockIdx.x=0才会有值，其余都为0）

reduce_v1消除warp divergence版本

warp divergence是否存在？只有if语句而没有else语句

这个问题暂且先放着了，当前出现的warp分歧按照解释可以做优化，例如下图

这里从第一轮迭代中可以看到工作的线程（红色）是隔着的，这就会很容易形成在每个warp里会出现两种状态：工作和等待，那为何不让工作的都在同一个warp里呢？如下图：

注意：如何判断你写的核函数里哪个线程工作哪个线程不工作呢？比如下面的代码修改成让工作的线程在一块，等待的线程在一块。

for(int index = 1; index < blockDim.x; index *= 2) {
        // 算法思路和v0一致，仅仅是用位运算替代了v0 if语句中的除余操作
        // if ((tid & (2 * index - 1)) == 0){
        //     smem[tid] += smem[tid + index];
        // }

        int index_s = 2*index*tid;
        if(tid < blockDim.x / (2* index))
            smem[index_s] += smem[index_s + index];
        __syncthreads();
    }

一定要注意是什么产生运算，比如smem[index_s] += smem[index_s + index];前的if条件，那这里的满足条件的tid就是会工作的线程，否则就是不参与运算的等待线程，因此是不是可以体会到核函数要带着千军万马来抢一个大任务，当遇到share memory的时候要注意都停下来，休整好后再一同出发。因此这里的if条件每进行一次循环，工作的线程就会减半，并且都是前blockDim.x / (2* index)个线程。

reduce_v1取余替换为位与版本

nv的cuda example文档里写明了取余是一个很耗时的操作，其实不管是GPU还是其他如CPU、NPU都是这样的，因此就想到能否用按位与操作替换取余。

for(int index = 1; index < blockDim.x; index *= 2) {
        // 算法思路和v0一致，仅仅是用位运算替代了v0 if语句中的除余操作
        if ((tid & (2 * index - 1)) == 0){
            smem[tid] += smem[tid + index];
        }
        __syncthreads();
    }

reduce_v2解决bank_conflict版本

为什么会发生bank_conflict，是因为cuda为了更好的并行，共享内存被分割成32个等大小的内存块，内存块也就是bank。bank冲突的发生条件是一个warp里不同线程访问相同bank的不同字段，那这里cuda的设计初衷其实是你不同的线程应该访问不同的bank，例如你可以thread 0完全只访问bank0，thread 1完全只访问bank1等。以下解决bank冲突其实就是尽量让一个线程只在同一个bank中活动。

共享内存中，连续的4-bytes被分配到连续的32个bank中（每一个bank存放一个32-bits的数据）

上图就是bank conflict的示意图，一个warp里thread0 和thread16都会访问bank0 和bank1，因此会发生两路冲突，这里还要注意，有一种情况不会发生bank conflict。

上图如果不同线程访问的是一个bank中同一个字段，那么这个会发生broadcast，不会发生bank conflict。

上图通过调整shared memory每个工作线程做加法的数据，可以让warp里的每个线程都在一个bank中活动，因此解决了bank conflict。

代码也如下：

// 基于v1作出改进: 从之前的当前线程ID加2*线程ID位置然后不断加上*2位置上的数据，改成不断地对半相加，以消除bank conflict
    // 此时一个block对d_in这块数据的reduce sum结果保存在id为0的线程上面
    for (unsigned int index = blockDim.x / 2; index > 0; index >>= 1) {
        if (tid < index) {   // 这里其实还改变了工作线程的排布，让一个block里工作的线程都排在最前面
            smem[tid] += smem[tid + index];
        }
        __syncthreads();
    }

reduce_v3让空闲线程干活版本

在每次迭代你会发现最多只有一半的线程处于工作状态，难免会太浪费，因此想让“更多”的线程都参与进来。为何打引号？见下图

其实更多意思就是让一个线程能干多份活，这里最简单就是在shared memory赋值的同时做一次加法运算。例如plan A里smem[tid] += d_in[gtid] + d_in[gtid + blockSize];这里block数量减少了，thread数量不变

planB是保持block的数量，thread数量减少。

代码如下：

    __shared__ float smem[blockSize];
    // 泛指当前线程在其block内的id
    unsigned int tid = threadIdx.x;
    // 泛指当前线程在所有block范围内的全局id, *2代表当前block要处理2*blocksize的数据
    // ep. blocksize = 2, blockIdx.x = 1, when tid = 0, gtid = 4, gtid + blockSize = 6; when tid = 1, gtid = 5, gtid + blockSize = 7
    // ep. blocksize = 2, blockIdx.x = 0, when tid = 0, gtid = 0, gtid + blockSize = 2; when tid = 1, gtid = 1, gtid + blockSize = 3
    // so, we can understand L18, one thread handle data located in tid and tid + blockSize 
    unsigned int gtid = blockIdx.x * (blockSize * 2) + threadIdx.x;
    // load: 每个线程加载两个元素到shared mem对应位置
    smem[tid] = d_in[gtid] + d_in[gtid + blockSize];
    __syncthreads();

reduce_v4展开最后一维减少同步

众所周知shared memory的读写要千万注意sync，这是为了保证一个block里所有线程的数据都正确写入后再做后续的计算。那其实这个层面是为了同步不同warp的，在同一个warp是不需要sync，因为并行是以一个warp为单位，如下图，绿色框是指warp里所有线程都是在同一时间内执行相同的指令。所以很简单的想法是当迭代到工作线程tid<32时，就可以不做__syncthreads()了。

代码如下：

__device__ void WarpSharedMemReduce(volatile float* smem, int tid){
    // CUDA不保证所有的shared memory读操作都能在写操作之前完成，因此存在竞争关系，可能导致结果错误
    // 比如smem[tid] += smem[tid + 16] => smem[0] += smem[16], smem[16] += smem[32]
    // 此时L9中smem[16]的读和写到底谁在前谁在后，这是不确定的，所以在Volta架构后最后加入中间寄存器(L11)配合syncwarp和volatile(使得不会看见其他线程更新smem上的结果)保证读写依赖
    float x = smem[tid];
    if (blockDim.x >= 64) {
      x += smem[tid + 32]; __syncwarp();
      smem[tid] = x; __syncwarp();
    }
    x += smem[tid + 16]; __syncwarp();
    smem[tid] = x; __syncwarp();
    x += smem[tid + 8]; __syncwarp();
    smem[tid] = x; __syncwarp();
    x += smem[tid + 4]; __syncwarp();
    smem[tid] = x; __syncwarp();
    x += smem[tid + 2]; __syncwarp();
    smem[tid] = x; __syncwarp();
    x += smem[tid + 1]; __syncwarp();
    smem[tid] = x; __syncwarp();
}
// Note: using blockSize as a template arg can benefit from NVCC compiler optimization, 
// which is better than using blockDim.x that is known in runtime.
template<int blockSize>
__global__ void reduce_v4(float *d_in,float *d_out){
    __shared__ float smem[blockSize];
    // 泛指当前线程在其block内的id
    int tid = threadIdx.x;
    // 泛指当前线程在所有block范围内的全局id, *2代表当前block要处理2*blocksize的数据
    // ep. blocksize = 2, blockIdx.x = 1, when tid = 0, gtid = 4, gtid + blockSize = 6; when tid = 1, gtid = 5, gtid + blockSize = 7
    // ep. blocksize = 2, blockIdx.x = 0, when tid = 0, gtid = 0, gtid + blockSize = 2; when tid = 1, gtid = 1, gtid + blockSize = 3
    // so, we can understand L38, one thread handle data located in tid and tid + blockSize 
    int i = blockIdx.x * (blockSize * 2) + threadIdx.x;
    // load: 每个线程加载两个元素到shared mem对应位置
    smem[tid] = d_in[i] + d_in[i + blockSize];
    __syncthreads();

    // 基于v3改进：把最后一个warp抽离出来reduce，避免多做一次sync threads
    // 此时一个block对d_in这块数据的reduce sum结果保存在id为0的线程上面
    for (int s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) {
            smem[tid] += smem[tid + s];
        }
        __syncthreads();
    }

    // last warp拎出来单独作reduce
    if (tid < 32) {
        WarpSharedMemReduce(smem, tid);
    }
    // store: 哪里来回哪里去，把reduce结果写回显存
    // GridSize个block内部的reduce sum已得出，保存到d_out的每个索引位置
    if (tid == 0) {
        d_out[blockIdx.x] = smem[0];
    }
}

这里要注意的点是warpReduce里也要注意读写依赖。

reduce_v5完全展开

做到极致，这里的优化其实没有什么技巧，仅仅只是为了让指令运行的更少，因此做循环展开，减少for循环中的加法指令，代码如下：

template <int blockSize>
__device__ void BlockSharedMemReduce(float* smem) {

    //对v4 L45的for循环展开，以减去for循环中的加法指令，以及给编译器更多重排指令的空间
  if (blockSize >= 1024) {
    if (threadIdx.x < 512) {
      smem[threadIdx.x] += smem[threadIdx.x + 512];
    }
    __syncthreads();
  }
  if (blockSize >= 512) {
    if (threadIdx.x < 256) {
      smem[threadIdx.x] += smem[threadIdx.x + 256];
    }
    __syncthreads();
  }
  if (blockSize >= 256) {
    if (threadIdx.x < 128) {
      smem[threadIdx.x] += smem[threadIdx.x + 128];
    }
    __syncthreads();
  }
  if (blockSize >= 128) {
    if (threadIdx.x < 64) {
      smem[threadIdx.x] += smem[threadIdx.x + 64];
    }
    __syncthreads();
  }
  // the final warp
  if (threadIdx.x < 32) {
    volatile float* vshm = smem;
    if (blockDim.x >= 64) {
      vshm[threadIdx.x] += vshm[threadIdx.x + 32];
    }
    vshm[threadIdx.x] += vshm[threadIdx.x + 16];
    vshm[threadIdx.x] += vshm[threadIdx.x + 8];
    vshm[threadIdx.x] += vshm[threadIdx.x + 4];
    vshm[threadIdx.x] += vshm[threadIdx.x + 2]; 
    vshm[threadIdx.x] += vshm[threadIdx.x + 1];
  }
}

reduce_v6在合适选择block数量

这一块其实还是在压榨thread干活，即让一个thread在第一个迭代给shared memory赋值时，让他多加几个数。

代码如下：

// 基于v5的改进：不用显式指定一个线程处理2个元素，而是通过L58的for循环来自动确定每个线程处理的元素个数
    float sum = 0.0f;
    for (int32_t i = gtid; i < nums; i += total_thread_num) {
        sum += d_in[i];
    }
    smem[tid] = sum;
    __syncthreads();

reduce_v7用shuffle做warp reduce

NV出了Shuffle指令，对于reduce优化有着非常好的效果。目前绝大多数访存类算子，像是softmax，batch_norm，reduce等，都是用Shuffle实现。Shuffle指令是一组针对warp的指令。Shuffle指令最重要的特性就是warp内的寄存器可以相互访问。在没有shuffle指令的时候，各个线程在进行通信时只能通过shared memory来访问彼此的寄存器。而采用了shuffle指令之后，warp内的线程可以直接对其他线程的寄存器进行访存。通过这种方式可以减少访存的延时。
如下图，reduce sum的warp shuffle做相加如下操作：

代码如下：

template <int blockSize>
__device__ float WarpShuffle(float sum) {
    // __shfl_down_sync：前面的thread向后面的thread要数据
    // __shfl_up_sync: 后面的thread向前面的thread要数据
    // 1. 返回前面的thread向后面的thread要的数据，比如__shfl_down_sync(0xffffffff, sum, 16)那就是返回16号线程，17号线程的数据
    // 2. 使用warp shuffle指令的数据交换不会出现warp在shared memory上交换数据时的不一致现象，这一点是由GPU driver完成，故无需任何sync, 比如syncwarp
    // 3. 原先15-19行有5个if判断block size的大小，目前已经被移除，确认了一下__shfl_down_sync等warp shuffle指令可以handle一个block或一个warp的线程数量<32，不足32会自动填充0
    sum += __shfl_down_sync(0xffffffff, sum, 16); // 0-16, 1-17, 2-18, etc.
    sum += __shfl_down_sync(0xffffffff, sum, 8);// 0-8, 1-9, 2-10, etc.
    sum += __shfl_down_sync(0xffffffff, sum, 4);// 0-4, 1-5, 2-6, etc.
    sum += __shfl_down_sync(0xffffffff, sum, 2);// 0-2, 1-3, 4-6, 5-7, etc.
    sum += __shfl_down_sync(0xffffffff, sum, 1);// 0-1, 2-3, 4-5, etc.
    return sum;
}

要注意shuffle之后超出的部分可以不用考虑，直接忽视掉就行。

参考链接

深入浅出GPU优化系列：reduce优化
【CUDA】Reduce规约求和（已完结~）

原文地址：https://blog.csdn.net/npu_yuanfang/article/details/144332607

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！