文章目录
- reduce baseline实现
- reduce_v0普通并行版本
- reduce_v1消除warp divergence版本
- reduce_v1取余替换为位与版本
- reduce_v2解决bank_conflict版本
- reduce_v3让空闲线程干活版本
- reduce_v4展开最后一维减少同步
- reduce_v5完全展开
- reduce_v6在合适选择block数量
- reduce_v7用shuffle做warp reduce
- 参考链接
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内的idunsigned 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内的idint 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拎出来单独作reduceif (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 warpif (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会自动填充0sum += __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规约求和(已完结~)