问题介绍
通俗的来说,reduce就是要对一个数组求 sum,min,max,avg 等等。reduce又被叫作规约,意思就是递归约减,最后获得的输出相比于输入一般维度上会递减。
比如 nvidia 博客上这个 reduce sum 问题,一个长度为 8 的数组求和之后得到的输出只有一个数,从 1 维数组变成一个标量。本文就以 reduce sum 为例来记录 reduce 优化。
硬件环境
nvidia a100-pcie-40gb , 峰值带宽在 1555 gb/s , cuda版本为11.8.
构建baseline
在问题介绍一节中的 reduce 求和图实际上就指出了 baseline 的执行方式,我们将以树形图的方式去执行数据累加,最终得到总和。但由于gpu没有针对 global memory 的同步操作,所以博客指出我们可以通过将计算分成多个阶段的方式来避免 global memrory 的操作。如下图所示:
接着 nvidia 博客给出了 baseline 算法的实现:
在这里插入图片描述
这里的 g_idata 表示的是输入数据的指针,而 g_odata 则表示输出数据的指针。然后首先把 global memory 数据 load 到 shared memory 中,接着在 shared memory 中对数据进行 reduce sum 操作,最后将 reduce sum 的结果写会 global memory 中。
但接下来的这页 ppt 指出了 baseine 实现的低效之处:
这里指出了2个问题,一个是warp divergent,另一个是取模这个操作很昂贵。这里的warp divergent 指的是对于启动 baseline kernel 的一个 block 的 warp 来说,它所有的 thread 执行的指令都是一样的,而 baseline kernel 里面存在 if 分支语句,一个 warp 的32个 thread 都会执行存在的所有分支,但只会保留满足条件的分支产生的结果。
我们可以在第8页ppt里面看到,对于每一次迭代都会有两个分支,分别是有竖直的黑色箭头指向的小方块(有效计算的线程)以及其它没有箭头指向的方块,所以每一轮迭代实际上都有大量线程是空闲的,无法最大程度的利用gpu硬件。
从这个ppt我们可以计算出,对于一个 block 来说要完成reduce sum,一共有8次迭代,并且每次迭代都会产生warp divergent。
接下来我们先把 baseline 的代码抄一下,然后我们设定好一个 gridsize 和 blocksize 启动 kernel 测试下性能。在ppt的代码基础上,我么补充一下内存申请以及启动 kernel 的代码。
#include #include #include #define n 32*1024*1024#define block_size 256__global__ void reduce_v0(float *g_idata,float *g_odata){ __shared__ float sdata[block_size]; // 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] = g_idata[i]; __syncthreads(); // 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(); } // write result for this block to global mem if (tid == 0) g_odata[blockidx.x] = sdata[0];}int main() { float *input_host = (float*)malloc(n*sizeof(float)); float *input_device; cudamalloc((void **)&input_device, n*sizeof(float)); for (int i = 0; i < n; i++) input_host[i] = 2.0; cudamemcpy(input_device, input_host, n*sizeof(float), cudamemcpyhosttodevice); int32_t block_num = (n + block_size - 1) / block_size; float *output_host = (float*)malloc((n / block_size) * sizeof(float)); float *output_device; cudamalloc((void **)&output_device, (n / block_size) * sizeof(float)); dim3 grid(n / block_size, 1); dim3 block(block_size, 1); reduce_v0<>(input_device, output_device); cudamemcpy(output_device, output_host, block_num * sizeof(float), cudamemcpydevicetohost); return 0;}
我们这里设定输入数据的长度是 32*1024*1024 个float数 (也就是ppt中的4m数据),然后每个 block 的线程数我们设定为 256 (block_size = 256, 也就是一个 block 有 8 个 warp)并且每个 block 要计算的元素个数也是 256 个,然后当一个 block 的计算完成后对应一个输出元素。
所以对于输出数据来说,它的长度是输入数据长度处以 256 。我们代入 kernel 加载需要的 gridsize(n / 256) 和 blocksize(256) 再来理解一下 baseline 的 reduce kernel。
首先,第 tid 号线程会把 global memroy 的第 i 号数据取出来,然后塞到 shared memroy 中。接下来针对已经存储到 shared memroy 中的 256 个元素展开多轮迭代,迭代的过程如 ppt 的第8页所示。完成迭代过程之后,这个 block 负责的256个元素的和都放到了 shared memrory 的0号位置,我们只需要将这个元素写回global memory就做完了。
接下来我们使用 nvcc -o bin/reduce_v0 reduce_v0_baseline.cu 编译一下这个源文件,并且使用nsight compute去profile一下。
性能和带宽的测试情况如下:
优化手段 耗时(us) 带宽利用率 加速比
reduce_baseline 990.66us 39.57% ~
优化手段1: 交错寻址(interleaved addressing)
接下来直接nvidia的ppt给出了优化手段1:
这里是直接针对 baseline 中的 warp divergent 问题进行优化,通过调整baseline中的分支判断代码使得更多的线程可以走到同一个分支里面,降低迭代过程中的线程资源浪费。具体做法就是把 if (tid % (2*s) == 0) 替换成 strided index的方式也就是int index = 2 * s * tid,然后判断 index 是否在当前的 block 内。
虽然这份优化后的代码没有完全消除if语句,但是我们可以来计算一下这个版本的代码在8次迭代中产生 warp divergent 的次数。对于第一次迭代,0-3号warp的index都是满足=blockdim.x的,也就是说这次迭代根本不会出现warp divergent的问题,因为每个warp的32个线程执行的都是相同的分支。
接下来对于第二代迭代,0,1两个warp是满足=blockdim.x,依然不会出现warp divergent,以此类推直到第4次迭代时0号warp的前16个线程和后16线程会进入不同的分支,会产生一次warp divergent,接下来的迭代都分别会产生一次warp divergent。
但从整体上看,这个版本的代码相比于baseline的代码产生的warp divergent次数会少得多。
我们继续抄一下这个代码然后进行profile一下。
#define n 32*1024*1024#define block_size 256__global__ void reduce_v1(float *g_idata,float *g_odata){ __shared__ float sdata[block_size]; // 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] = g_idata[i]; __syncthreads(); // 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]; // } int index = 2 * s * tid; if (index 0; s >>= 1) { if (tid 0; s >>= 1) { if (tid < s){ sdata[tid] += sdata[tid + s]; } __syncthreads(); } // write result for this block to global mem if (tid == 0) g_odata[blockidx.x] = sdata[0];}
性能和带宽的测试情况如下:
优化手段 耗时(us) 带宽利用率 加速比
reduce_baseline 990.66us 39.57% ~
reduce_v1_interleaved_addressing 479.58us 81.74% 2.06
reduce_v2_bank_conflict_free 462.02us 84.81% 2.144
reduce_v3_idle_threads_free 244.16us 83.16% 4.057
优化手段4: 展开最后一个warp
首先来看 ppt 的第20页:
这里的意思是,对于 reduce_v3_idle_threads_free 这个 kernel 来说,它的带宽相比于理论带宽还差得比较远,因为 reduce 操作并不是算术密集型的算子。
因此,一个可能的瓶颈是指令的开销。这里说的指令不是加载,存储或者给计算核心用的辅算术指令。换句话说,这里的指令就是指地址算术指令和循环的开销。
接下来 ppt 21指出了减少指令开销的优化方法:
这里的意思是当reduce_v3_idle_threads_free kernel里面的s<=32时,此时的block中只有一个warp0在干活时,但线程还在进行同步操作。这一条语句造成了极大的指令浪费。
由于一个warp的32个线程都是在同一个simd单元上,天然保持了同步的状态,所以当s32; s >>= 1) { if (tid < s){ sdata[tid] += sdata[tid + s]; } __syncthreads(); } // write result for this block to global mem if (tid = 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 __global__ void reduce_v5(float *g_idata,float *g_odata){ __shared__ float sdata[block_size]; // 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] = g_idata[i] + g_idata[i + blockdim.x]; __syncthreads(); // do reduction in shared mem if(blocksize>=512){ if(tid=256){ if(tid=128){ if(tid<64){ sdata[tid]+=sdata[tid+64]; } __syncthreads(); } // write result for this block to global mem if(tid= 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 __global__ void reduce_v6(float *g_idata,float *g_odata){ __shared__ float sdata[block_size]; // each thread loads one element from global to shared mem unsigned int tid = threadidx.x; unsigned int i = blockidx.x*(blockdim.x * num_per_thread) + threadidx.x; sdata[tid] = 0; #pragma unroll for(int iter=0; iter=512){ if(tid=256){ if(tid=128){ if(tid<64){ sdata[tid]+=sdata[tid+64]; } __syncthreads(); } // write result for this block to global mem if(tid<32)warpreduce(sdata,tid); if (tid == 0) g_odata[blockidx.x] = sdata[0];}int main() { float *input_host = (float*)malloc(n*sizeof(float)); float *input_device; cudamalloc((void **)&input_device, n*sizeof(float)); for (int i = 0; i < n; i++) input_host[i] = 2.0; cudamemcpy(input_device, input_host, n*sizeof(float), cudamemcpyhosttodevice); const int block_num = 1024; const int num_per_block = n / block_num; const int num_per_thread = num_per_block / block_size; float *output_host = (float*)malloc((block_num) * sizeof(float)); float *output_device; cudamalloc((void **)&output_device, (block_num) * sizeof(float)); dim3 grid(block_num, 1); dim3 block(block_size, 1); reduce_v6<>(input_device, output_device); cudamemcpy(output_device, output_host, block_num * sizeof(float), cudamemcpydevicetohost); return 0;}
profile结果:
性能和带宽的测试情况如下:
优化手段 耗时(us) 带宽利用率 加速比
reduce_baseline 990.66us 39.57% ~
reduce_v1_interleaved_addressing 479.58us 81.74% 2.06
reduce_v2_bank_conflict_free 462.02us 84.81% 2.144
reduce_v3_idle_threads_free 244.16us 83.16% 4.057
reduce_v4_unroll_last_warp 167.10us 54.10% 5.928
reduce_v5_completely_unroll 158.78us 56.94% 6.239
reduce_v6_multi_add 105.47us 85.75% 9.392
在把block_num从65536调整到1024之后,无论是性能还是带宽都达到了最强,相比于最初的baseline加速了9.4倍。
总结
我这里的测试结果和nvidia ppt里提供的结果有一些出入,nvidia ppt的34页展示的结果是对于每一种优化相比于前一种无论是性能还是带宽都是稳步提升的。但我这里的测试结果不完全是这样,对于 reduce_v4_unroll_last_warp 和 reduce_v5_completely_unroll 这两个优化,虽然耗时近一步减少但是带宽却降低了,我也还没想清楚原因。
并且最终的kernel带宽利用率为 73 / 86.4 = 84.5% ,和我在a100上的reduce_v6_multi_add kernel的测试结果基本相当。
“Home+智能联接”三层体验升级:开辟家庭多元化业务场景
七夕情人节不知道送什么?送这几款蓝牙耳机对方肯定喜欢!
OPPO刷脸,华为用“芯”,OPPO干的过华为吗?
轴承座发热原因及处理方法
BT工作原理
NIVDIA的reduce优化笔记
飞兆半导体开发出PowerTrench MOSFET器件FDMB2307NZ
Q3俄罗斯手机市场线上销量华为夺冠
基于51单片机的酒精浓度监测仪
R型电源变压器有载调压是什么?
荣耀9怎么样?荣耀9与一加5谁才是颜值担当?荣耀9与一加5区别对比评测
遇到信号满格却无法打电话时,该如何解决
Team Owl Works携全树脂桌面3D打印机全新亮相
硅胶洗墙灯的特性介绍
土壤重金属检测仪的功能有哪些
如何将交流电源变成直流电源?
上海关区集成电路进口达4036.2亿元 为最大类进口产品
小米PocophoneF1Lite跑分曝光 搭载骁龙660或为千元级水准
什么是生物多样性
万用表检测多种类型电容器的方法盘点