CUDA算子：reduce优化

最新推荐文章于 2024-09-24 16:30:29 发布

2_33_

最新推荐文章于 2024-09-24 16:30:29 发布

阅读量1.1k

点赞数 21

文章标签：人工智能大数据

本文链接：https://blog.csdn.net/qq_56047026/article/details/139053428

版权

如何对GPU中的reduce算法进行优化。官方地址

https://developer.download.nvidia.cn/assets/cuda/files/reduction.pdf

reduce算法的本质： $x = x_{0}\bigotimes x_{1}\bigotimes x_{2}...\bigotimes x_{n-1}\bigotimes x_{n}$

并行算法设计

GPU中，reduce采用了树形的计算方式

从上到下，将数据不断累加，直到得出最后的结果，25。GPU没有针对global数据的同步操作，只能针对block的数据进行同步。所以将reduce分为两个阶段，示意图如下：

假设给定一个长度为N的数组，需要计算该数组的所有元素之和。首先需要将数组分为m个小份。在第一阶段，开启m个block计算出m个小份的reduce值。在第二阶段，使用一个block将m个小份再次进行reduce，得到最终结果。由于第二段本质上可以调用第一阶段的kernel，所以本文只探索第一阶段的优化技巧。

kernel接口

__global__void reduce(T *input, T* output)

input：输入的数组，一个长度为N的数组；output：输出的数组，第一阶段的结果，长度为M的数组。

在开始CUDA编程前，设置三个参数：

1.BlockNum：开启的block数量，即M，代表需要将数组切成几份

2.Thread_per_block:每个block中开启的线程束，一般：128，256，512，1024

3.Num_per_block:每个block需要进行reduce操作的长度

其中：BlockNum*Num_per_block = N

**reduce baseline算法：if (tid%(2*s)==0)**

三个步骤：

1.将数据load至shared memory中

2.在shared memory中对数据进行reduce操作

3.将最后的结果写回global memory中。

在第一个步骤中，让Num_per_block和Thread_per_block一致，每个block设置256个线程，一个block负责256个数据的reduce工作。

假设需要处理32M（ $2^{5}*2^{20}$ ）的数据，则有128K（ $2^{7}*2^{10}$ ）个block。

tid：线程号；i：原始数组中的索引号。第tid号线程将第i号的数据从global中取出，放到shared memory的第tid元素中。比如在第0号block中，0号线程将0号元素取出，放到shared memory的第0号位置。

从硬件角度分析，为了执行代码，GPU需要分配存储资源和计算资源。存储资源包括在global memory中分配的一块32M*sizeof（float）的空间以及在shared memory中分配的256*sizeof（float）的空间。shared memory存在bank冲突的问题。计算资源根据thread数量确定，一个block中分配256个thread线程，32个线程一组，绑定在一个SIMD（单指令多数据流）单元。（实际的硬件资源分配不是这样的，因为一个SM的计算资源有限，不可能给每一个block都分配这么多的SIMD）。

第一阶段：tid号线程将i号数据从全局显存中取出，再放进共享显存中，中间走一遍寄存器再到共享显存中。

#include <bits/stdc++.h>
#include <cuda.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <time.h>
#include <sys/time.h>

#define THREAD_PER_BLOCK 256

__global__ void reduce0(float *d_in,float *d_out){
    __shared__ float sdata[THREAD_PER_BLOCK];

    // each thread loads one element from global to shared mem
    unsigned int tid = threadIdx.x;//线程号
    unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;//在原始数组中的索引号
    sdata[tid] = d_in[i];//第tid号线程将第i号的数据从global中取出，放到shared memory的第tid元素中
    __syncthreads();//确保线程块中的所有线程都已到达该调用点

第二阶段，block中需要计算的256个元素全部被存储在了shared memory中，此时需要进行reduce操作。多轮迭代，第一轮迭代中，如果tid%2==0,则第tid号线程将shared memory中第tid号位置的值和第tid+1号的值进行相加，放在tid位置；第二轮迭代，如果tid%4==0，则将tid和tid+2号的值相加，放在第tid号位置，不断迭代，所有元素都被累加到第0号位置。

判断tid%2*s == 0再进行下一步，产生了分支

    // do reduction in shared mem
    for(unsigned int s=1; s < blockDim.x; s *= 2) {
        if (tid % (2*s) == 0) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

在第三个阶段，block负责的256个元素值和都放置在shared memory的0号位置，此时，只需要将0号位置的元素写回即可。

    // write result for this block to global mem
    if (tid == 0) d_out[blockIdx.x] = sdata[0];

实验结果

nsight-sys打开分析系统

**优化技巧1：解决warp divergence 问题线程束分散 if (2stid < blockDim.x)**

对于一个block而言，所有thread都是执行同一条命令，如果存在if-else这样的分支情况的话，thread会执行所有的分支。只是不满足条件的分支，所产生的结果不会被记录下来。每一轮迭代都会产生两个分支，严重影响效率

解决方式

尽可能地让所有线程走到同一个分支里面

只要2*s*tid<blockDim.x，都进行计算

    // do reduction in shared mem
    for(unsigned int s=1; s < blockDim.x; s *= 2) {
        int index = 2 * s * tid;
        if (index < blockDim.x) {
            sdata[index] += sdata[index + s];
        }
        __syncthreads();
    }

假定block中存在256个thread，即拥有256/32=8个warp。当进行第一次迭代时，0-3号warp的index<blockDim.x，4-7号warp的index>=blockDim.x。对于每个warp而言，都只是进入到一个分支内，不存在warp divergence的情况。当进行第二次迭代时，0、1号两个warp进入计算分支。当进行第3次迭代时，只有0号warp进入计算分支。进行第4此迭代时，只有0号warp的前16个线程进入分支。此时开始产生warp divergence。通过这种方式，消除了前3次迭代的warp。

优化技巧2：解决bank冲突反循环 int s = blockDim.x/2,

reduce1最大问题是bank冲突。目光聚焦在for循环，并且只聚焦在0号warp。第一次迭代，0号线程需要去load shared memory中的0号地址以及1号地址的数，然后写回到0号地址。此时，这个warp中的16号线程，需要去load shared memory中的32号地址和33号地址。可以发现，0号地址跟32号地址产生了两路bank冲突。

在第二次迭代中，0号线程需要load shared memory中的0号地址和2号地址。这个warp中的8号线程需要load shared memory中的32号地址和34号地址，16号线程需要load shared memory中的64号地址和68号地址，24号线程需要load shared memory中的32号地址以及34号地址，16号线程需要load shared memory中的64号地址和68号地址，24号线程需要load shared memory中的96号地址和100号地址。又因为0、32、64、96号地址对应着同一个银行，所以此时产生了4路的银行冲突。现在，可以继续算下去，8路bank冲突，16路冲突。由于bank冲突，所以reduce1性能受限。下图说明了在加载第一个数据时所产生的bank冲突。

解决方式

在reduce中，解决bank冲突的方式就是把for循环逆着来。原来stride从0到256，现在stride从128到0。其伪代码如下：

for(unsigned int s = blockDim.x/2; s > 0; s>>=1){
    if (tid < s){
        sdata[tid] += sdata[tid + s];
    }
    _syncthreads();
}

0号线程需要加载共享内存的0号元素和128号元素。1号线程需要加载共享内存中的1号元素和129号元素。在这一迭代中，在读取第一个数时，warp中的32个线程刚好加载一行shared memory数据。再分析第2轮迭代，0号线程加载0号元素和64号元素，1号元素加载1号元素和65号元素。每次加载共享内存的一行。分析第3轮迭代，0号线程加载0号元素和32号元素。一个warp load shared memory的一行。没有bank冲突，问题：到了第4轮迭代，0号线程加载0号元素和16号元素。16号线程啥也不干，s=16,16-31啥也不干，跳过去了。

优化技巧3：解决idle线程

reduce2最大的问题是线程的浪费。可以看到启动了256个线程，但是在第1轮迭代时只有128个线程在干活，第二轮迭代只有64个线程在干活，每次干活的线程都会减少一半。第一轮迭代示意图如下，只有前128个线程在load数据。后128个线程啥也不干。

解决方式：

每一次迭代都有一半的线程不干活。而且，128-255号线程没有贡献，让他们做一次加法。除了去全局内存中取数外，再做一次加法。为了实现这个，block数量减少，Num_per_block增加一倍。原来一个block需要管256个数，现在需要管512个数

__global__ void reduce3(float *d_in, float *d_out) {
    __shared__ float sdata[THREAD_PER_BLOCK];
    
    //ench thread loads one element from global memory to shared mem
    unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x;//当前的block和thread的索引，每个block有两个warp
    unsigned int tid = threadIdx.x;
    sdata[tid] = d_in[i] + d_in[i+blockDim.x];
    __syncthreads();

}
    //do reduction in shared mem
    for(unsigned int s = blockDim.x/2; s>0; s>>=1) {
        if (tid < s){
            sdata[tid] += sdata[tid+s];
        }
        __syncthreads();
    }
    //write result for this block to global mem
    if(tid == 0)d_out[blockIdx.x] = sdata[tid];
}

通过这种方式，将一些idle的线程利用起来了

优化技巧4：展开最后减少同步

现有问题

当进行到最后几轮迭代时，此时的block中只有warp0在干活时，线程还在进行同步操作。这一条语句造成了极大浪费。

解决方式

由于一个warp中的32个线程，其实是在一个SIMD单元上，这32个线程每次都是执行同一条命令，保持了同步状态。当s=32时，即只有一个SIMD单元在工作时，完全可以将__syncthreads()这条同步代码去掉。所以我们将最后一维展开以减少同步。伪代码如下：

__device__ void warpReduce(volatile float* cache,int tid){
    cache[tid]+=cache[tid+32];
    cache[tid]+=cache[tid+16];
    cache[tid]+=cache[tid+8];
    cache[tid]+=cache[tid+4];
    cache[tid]+=cache[tid+2];
    cache[tid]+=cache[tid+1];
}

__global__ void reduce3(float *d_in, float *d_out) {
    __shared__ float sdata[THREAD_PER_BLOCK];
    
    //ench thread loads one element from global memory to shared mem
    unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x;//当前的block和thread的索引，每个block有两个warp
    unsigned int tid = threadIdx.x;
    sdata[tid] = d_in[i] + d_in[i+blockDim.x];
    __syncthreads();

}
    //do reduction in shared mem
    for(unsigned int s = blockDim.x/2; s>0; s>>=1) {
        if (tid < s){
            sdata[tid] += sdata[tid+s];
        }
        __syncthreads();
    }
    //write result for this block to global mem
    if(tid < 32) warpReduce(sdata,tid);
    if(tid == 0)d_out[blockIdx.x] = sdata[tid];
}

优化技巧5：完全展开减少计算

对于for循环进行完全展开，主要减少for循环的开销

template <unsigned int blockSize>
__device__ void warpReduce(volatile float* cache, unsigned int tid){
    if (blockSize >= 64)cache[tid]+=cache[tid+32];
    if (blockSize >= 32)cache[tid]+=cache[tid+16];
    if (blockSize >= 16)cache[tid]+=cache[tid+8];
    if (blockSize >= 8)cache[tid]+=cache[tid+4];
    if (blockSize >= 4)cache[tid]+=cache[tid+2];
    if (blockSize >= 2)cache[tid]+=cache[tid+1];
}

template <unsigned int blockSize>
__global__ void reduce5(float *d_in,float *d_out){
    __shared__ float sdata[THREAD_PER_BLOCK];

    // each thread loads one element from global to shared mem
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    sdata[tid] = d_in[i] + d_in[i + blockDim.x];
    __syncthreads();

    // do reduction in shared mem
    if (blockSize >= 512) {
        if (tid < 256) { 
            sdata[tid] += sdata[tid + 256]; 
        } 
        __syncthreads(); 
    }
    if (blockSize >= 256) {
        if (tid < 128) { 
            sdata[tid] += sdata[tid + 128]; 
        } 
        __syncthreads(); 
    }
    if (blockSize >= 128) {
        if (tid < 64) { 
            sdata[tid] += sdata[tid + 64]; 
        } 
        __syncthreads(); 
    }
    if (tid < 32) warpReduce<blockSize>(sdata, tid);

    // write result for this block to global mem
    if (tid == 0) d_out[blockIdx.x] = sdata[0];
}

优化技巧6：合理设置block数量

调整block和thread的数值，之前一个block开启256个线程，这个block负责256个元素的reduce操作，可以让block多管点数，开启block的数量少些。

因此需要重新定义block的取值，block数量多，block可以快速进行切换，去掩盖访存的延时。如果一个block被分配更多的工作，可能会更好的覆盖延时，如果线程有更多的work，对编译器而言，就可能对相关指令进行重新排列，从而覆盖延时，证明block少一些更好。block设置2048时要比512，1024，4096更好

优化技巧7 使用shuffle指令

目前绝大多数的访存类算子，如softmax，batch_norm，reduce等，都是用shuffle实现。如何把shuffle指令用在reduce优化上。

shuffle指令是一组针对warp（基本的并行执行单元）的指令，重要特性是warp内的寄存器可以相互访问。没有shuffle时，各个线程在进行通信时只能通过共享内存来访问彼此的寄存器。采用shuffle后，warp内的线程可以直接对其他线程的寄存器进行访问，从而减小访存的延时。此外，可编程性提高了，某些场景下，不用共享内存了。