锚文本对网站济南市建设工程交易网

锚文本对网站,济南市建设工程交易网,天津做网站报价,优质企业网站推广1.Ruduce/规约 定义 对整个张量进行一个操作#xff0c;得到一个标量结果。这个操作可以是max,min,summax,min,summax,min,sum等。 2.原理 我们用单线程的思路来实现的话#xff0c;就是遍历整个张量#xff0c;然后来做reducereducereduce的操作即可。 for(int i 0; i…1.Ruduce/规约 定义对整个张量进行一个操作得到一个标量结果。这个操作可以是m a x , m i n , s u m max,min,summax,min,sum等。2.原理我们用单线程的思路来实现的话就是遍历整个张量然后来做r e d u c e reducereduce的操作即可。for(inti0;in;i){ansreduce_op(ans,a[i]);}但我们这里用c u d a cudacuda了想要加速就要挖掘可并行的地方。注意到r e d u c e reducereduce操作是有结合律的可以分块计算每一块的结果然后再合并。那么每一块的计算彼此之间是没有数据依赖的是可并行的。需要同步的地方只有需要等到每一块都算完才能开始合并。我们这里为了简单起见实际上实现的是每个b l o c k blockblock计算出规约结果即可接下来把每个b l o c k blockblock的结果再规约到一个标量的过程是类似的。3.实现以下实现主要借鉴视频【CUDA】Reduce规约求和已完结~【CUDA】硬件与工具探索合集更新ing以下的性能分析均在R T X 4090 RTX4090RTX4090上实现3.0 暴力/朴素规约算法说是暴力但其实理解起来也没那么简单。需要理解基础的规约算法还有线程模型基础的规约算法思路是每一个分块内我们每一次都让问题规模折半让当前的结果数组中元素两两进行合并。假设初始长度是8 88那么第一次合并后结果位于[ 0 , 2 , 4 , 6 ] [0,2,4,6][0,2,4,6]再合并一次结果位于[ 0 , 4 ] [0,4][0,4]以此类推如下图这其实类似一个位运算t r i c k tricktrick给你一个i n t 64 int64int64型变量如何在6 66次位运算操作内计算出二进制位上1 11的个数的奇偶。答案是每次都s h i f t shiftshift然后异或这和规约算法是等价的这个过程我们的并行度在于每一次问题规模的缩小进行的操作都是可并行的读写的都是不同位置。需要同步的地方只有每一层需要全算完了才能开始下一层。这个算法的实现上可以每一轮都遍历当前有效位置让相邻的有效位置进行规约操作想要遍历有效位置每轮的步长需要倍增。每个分块我们安排等同于块长的线程数第i ii个线程就处理第i ii个位置的操作所以第一轮需要操作的线程就是0 , 2 , 4 , 6 0,2,4,60,2,4,6第二轮就是0 , 4 0,40,4以此类推当然这是第一份代码所以我们需要搭建整体框架包括d e v i c e devicedevice侧的暴力计算资源申请释放调用核函数检查结果正确性。#includestdio.h#defineBLOCK_SIZE256boolcheck(int*a,int*b,intN){for(inti0;iN;i){if(a[i]!b[i]){returnfalse;}}returntrue;}__global__voidreduce(int*d_input,int*d_output){//计算当前块的起始地址int*startd_inputblockIdx.x*blockDim.x;for(inti1;iblockDim.x;i*2){if(threadIdx.x%(i*2)0){start[threadIdx.x]start[threadIdx.xi];}//每层需要同步 全算完了才能进入下一层__syncthreads();}//全部计算完了0号线程保存的就是这个块的结果拷贝到结果数组if(threadIdx.x0){d_output[blockIdx.x]start[0];}}voidsolve(){// 申请host device侧空间intN32*1024*1024;intblock_numN/BLOCK_SIZE;int*input(int*)malloc(N*sizeof(int));int*d_input;cudaMalloc((int**)d_input,N*sizeof(int));int*output(int*)malloc(block_num*sizeof(int));int*d_output;cudaMalloc((int**)d_output,block_num*sizeof(int));int*result(int*)malloc(block_num*sizeof(int));// 初始化值srand(666);for(inti0;iN;i){input[i]rand()%1000;}// cpu暴力计算for(inti0;iblock_num;i){intsum0;for(intj0;jBLOCK_SIZE;j){suminput[i*BLOCK_SIZEj];}result[i]sum;}cudaMemcpy(d_input,input,N*sizeof(int),cudaMemcpyHostToDevice);dim3Grid(block_num);dim3Block(BLOCK_SIZE);reduceGrid,Block(d_input,d_output);cudaMemcpy(output,d_output,block_num*sizeof(int),cudaMemcpyDeviceToHost);if(!check(output,result,block_num)){for(inti0;iblock_num;i){if(i0output[i]!result[i]){printf(blockidx:%d,output:%d,ans:%d\n,i,output[i],result[i]);}}}else{printf(test passed!\n);}cudaFree(d_output);cudaFree(d_input);free(input);free(output);free(result);}intmain(){solve();return0;}n c u ncuncu分析性能看D u r a t i o n DurationDuration字段377 u s 377 us377us以此作为基准衡量后面的优化效果3.1 使用共享内存/shared_memory一个比较显然的优化是上一份代码中我们操作的都是全局内存但每个块都只用自己块内的数据那么可以把每个块的数据拷贝到自己的共享内存中访问速度更快。这个简单的优化其实包含一个深刻的思想局部性和缓存。在G P U GPUGPU编程诞生前人们很早就意识到程序的空间局部性因此发明了层次存储结构和多级缓存把频繁访问的数据放在访问速度更快的缓存中只不过操作缓存的接口并没有留给编程开发者而是留给系统自行调度开发者想要利用这个局部性只能通过调整数组存储顺序或者在遍历时进行分块来实现//调整存储顺序for(inti0;in;i){for(intj0;jm;j){//ans a[j][i];ansa[i][j];}}//分块//for(int i 0; i n; i ){////}for(inti0;in/B;iB){for(intj0;jB;j){}}c u d a cudacuda的共享内存其实就是一个高速缓存并且把操作缓存的权限给了开发者于是我们如果知道每个线程块内频繁使用的数据是哪些就可以直接拷贝进去而不需要等着编译器或者操作系统来给我们优化#includestdio.h#defineBLOCK_SIZE256// 使用共享内存__global__voidreduce1(int*d_input,int*d_output){// 计算全局内存下标拷贝到共享内存int*startd_inputblockIdx.x*blockDim.x;__shared__intsdata[BLOCK_SIZE];sdata[threadIdx.x]start[threadIdx.x];//拷贝这里也需要同步__syncthreads();for(inti1;iblockDim.x;i*2){if(threadIdx.x%(i*2)0){sdata[threadIdx.x]sdata[threadIdx.xi];}__syncthreads();}if(threadIdx.x0){d_output[blockIdx.x]sdata[0];}}boolcheck(int*a,int*b,intN){for(inti0;iN;i){if(a[i]!b[i]){returnfalse;}}returntrue;}voidsolve(){// 申请host device侧空间intN32*1024*1024;intblock_numN/BLOCK_SIZE;int*input(int*)malloc(N*sizeof(int));int*d_input;cudaMalloc((int**)d_input,N*sizeof(int));int*output(int*)malloc(block_num*sizeof(int));int*d_output;cudaMalloc((int**)d_output,block_num*sizeof(int));int*result(int*)malloc(block_num*sizeof(int));// 初始化值srand(666);for(inti0;iN;i){input[i]rand()%1000;}// cpu暴力计算for(inti0;iblock_num;i){intsum0;for(intj0;jBLOCK_SIZE;j){suminput[i*BLOCK_SIZEj];}result[i]sum;}cudaMemcpy(d_input,input,N*sizeof(int),cudaMemcpyHostToDevice);dim3Grid(block_num);dim3Block(BLOCK_SIZE);reduce1Grid,Block(d_input,d_output);cudaMemcpy(output,d_output,block_num*sizeof(int),cudaMemcpyDeviceToHost);if(!check(output,result,block_num)){for(inti0;iblock_num;i){if(i0output[i]!result[i]){printf(blockidx:%d,output:%d,ans:%d\n,i,output[i],result[i]);}}}else{printf(test passed!\n);}cudaFree(d_output);cudaFree(d_input);free(input);free(output);free(result);}intmain(){solve();return0;}338 u s 338 us338us有一定效果不是很明显主要是因为我么访存次数不算太多共享内存的读取加速没有掩盖掉从全局内存拷贝到贡献内存的开销。如果是g e m m gemmgemm之类访存很多的算子加速效果才会比较明显3.2 消除 线程束分化/warp_diverse进一步的优化需要我们理解c u d a cudacuda的执行过程。c u d a cudacuda的多线程执行过程实际上是每次取出一个b l o c k blockblock内的32 3232个线程组成一个线程束w a r p warpwarp然后统一执行这里一个问题一个w a r p warpwarp内的线程如果执行的操作不一样怎么办比如有个分支指令一部分线程进入了A AA分支一部分进入了B BB分支答案是这两部分之间会退化成串行的如下图来看一下我们目前的实现没有没线程束分化我们前面说了我们是让第i ii个线程负责第i ii个位置那么在第一轮中0 , 2 , . . 30 0,2,..300,2,..30这些线程是有运算的其余是没有运算的走的是两个不同分支第二轮中0 , 4 , . . . , 28 0,4,...,280,4,...,28这些线程是有运算的其余是不进行运算的也是走了两个不同分支。可以看到是有线程束分化的。for(inti1;iblockDim.x;i*2){if(threadIdx.x%(i*2)0){sdata[threadIdx.x]sdata[threadIdx.xi];}__syncthreads();}我们想避免这一点去掉分支指令是不太可能的但是可以让一个w a r p warpwarp内的分支结果全都相同。仔细来看产生分化的地方如果一共有64 6464个元素第一轮执行运算的线程是0 , 2 , 4.. , 62 0,2,4..,620,2,4..,62这里只用了32 3232个线程我们可以让前32 3232个线程0 , 1 , . . , 31 0,1,..,310,1,..,31来执行这些运算操作然后后32 3232个线程不执行运算操作这样两个w a r p warpwarp内的线程就都无分化了这里在实现上线程i ii负责的位置就不再是i ii了而是需要计算的我们计算出这个下标i d x idxidx看他是否超出当前线程块的大小如果不超说明这个线程需要进行运算操作超出则说明不需要这样我们每轮执行运算操作线程都是从0 00号线程开始的一段就能消除线程束分化了。如下图红色和橙色分别表示有无计算操作的线程#includestdio.h#defineBLOCK_SIZE256// 使用共享内存// 消除wrap分化__global__voidreduce2(int*d_input,int*d_output){// 计算全局内存下标拷贝到共享内存int*startd_inputblockIdx.x*blockDim.x;__shared__intsdata[BLOCK_SIZE];sdata[threadIdx.x]start[threadIdx.x];__syncthreads();for(inti1;iblockDim.x;i*2){intidxthreadIdx.x*i*2;if(idxBLOCK_SIZE){sdata[idx]sdata[idxi];}__syncthreads();}if(threadIdx.x0){d_output[blockIdx.x]sdata[0];}}boolcheck(int*a,int*b,intN){for(inti0;iN;i){if(a[i]!b[i]){returnfalse;}}returntrue;}voidsolve(){// 申请host device侧空间intN32*1024*1024;intblock_numN/BLOCK_SIZE;int*input(int*)malloc(N*sizeof(int));int*d_input;cudaMalloc((int**)d_input,N*sizeof(int));int*output(int*)malloc(block_num*sizeof(int));int*d_output;cudaMalloc((int**)d_output,block_num*sizeof(int));int*result(int*)malloc(block_num*sizeof(int));// 初始化值srand(666);for(inti0;iN;i){input[i]rand()%1000;}// cpu暴力计算for(inti0;iblock_num;i){intsum0;for(intj0;jBLOCK_SIZE;j){suminput[i*BLOCK_SIZEj];}result[i]sum;}cudaMemcpy(d_input,input,N*sizeof(int),cudaMemcpyHostToDevice);dim3Grid(block_num);dim3Block(BLOCK_SIZE);reduce2Grid,Block(d_input,d_output);cudaMemcpy(output,d_output,block_num*sizeof(int),cudaMemcpyDeviceToHost);if(!check(output,result,block_num)){for(inti0;iblock_num;i){if(i0output[i]!result[i]){printf(blockidx:%d,output:%d,ans:%d\n,i,output[i],result[i]);}}}else{printf(test passed!\n);}cudaFree(d_output);cudaFree(d_input);free(input);free(output);free(result);}intmain(){solve();return0;}235 u s 235us235us这个效果还是比较明显的3.3 消除bank conflict更进一步的优化需要我们了解c u d a cudacuda的内存架构为了加速内存的读取底层不是只有一个内存可读而是把内存拆成32 3232份每一份的读写时独立的可并行这里设计成32 3232也是为了配合前面提到的w a r p warpwarp中有32 3232个线程理想情况下w a r p warpwarp中的每个线程正好都去读不同的b a n k bankbank内存带宽提升32 3232倍。比如下图b a n k 0 bank0bank0存的就是内存上0 , 32 , 64 0,32,640,32,64这些位置的值这里可能出现的问题是可能多个线程会去同时读相同的b a n k bankbank如果有n nn个线程去读同一b a n k bankbank里的不同数据称为n nn路b a n k bankbank冲突那么在这个b a n k bankbank的读就又变成串行的了。注意如果多个线程读b a n k bankbank里的同一数据不会触发冲突而是会触发广播b a n k bankbank会把这个值广播出去给所有需要的线程。如下图第二种情况在b a n k 0 bank0bank0两个线程来读不同的数据就是一个2 22路b a n k bankbank冲突而最后两种情况多个线程读同一个b a n k bankbank但是读的都是同一数据不会冲突而是广播回到我们的程序有b a n k bankbank冲突吗显然是有的首先我们只用关注一个w a r p warpwarp内的线程因为他们只有是同时执行可能触发冲突然后如下图0 , 8 , 16 , 24 0,8,16,240,8,16,24线程一块取加法第一个操作数的值时会同时访问0 , 32 , 64 , 96 0,32,64,960,32,64,96的位置触发一个至少四路的冲突如何解决呢可以调整每个线程负责的的元素位置我们让每个线程负责的元素不再是相邻的两个而是左右半区各一个如下图这样线程0 , 8 , 16 , 24 0,8,16,240,8,16,24访问的内存就是0 , 8 , 16 , 24 0,8,16,240,8,16,24了他们之间显然是不存在冲突的。并且之前线程编号和负责的内存位置不是直接对应的有点别扭现在就简单了线程i ii负责的就是i , i s t r i d e i,istridei,istride这两个位置。s t r i d e stridestride是当前有效元素个数的一半每轮折半。这样实现起来也更清晰。#includestdio.h#defineBLOCK_SIZE256// 使用共享内存// 消除wrap分化// 消除bank conflict__global__voidreduce3(int*d_input,int*d_output){int*startd_inputblockIdx.x*blockDim.x;__shared__intsdata[BLOCK_SIZE];sdata[threadIdx.x]start[threadIdx.x];__syncthreads();// if (blockIdx.x 1 threadIdx.x 0)// {// for (int i 0; i BLOCK_SIZE; i)// {// printf((%d)%d , i, sdata[i]);// }// printf(\n);// }// 消除wrap分化for(intiBLOCK_SIZE/2;i1;i/2){if(threadIdx.xi){sdata[threadIdx.x]sdata[threadIdx.xi];}__syncthreads();// if (blockIdx.x 1 threadIdx.x 0)// {// for (int i 0; i BLOCK_SIZE; i)// {// printf((%d)%d , i, sdata[i]);// }// printf(\n);// }}if(threadIdx.x0){d_output[blockIdx.x]sdata[0];}}boolcheck(int*a,int*b,intN){for(inti0;iN;i){if(a[i]!b[i]){returnfalse;}}returntrue;}voidsolve(){// 申请host device侧空间intN32*1024*1024;intblock_numN/BLOCK_SIZE;int*input(int*)malloc(N*sizeof(int));int*d_input;cudaMalloc((int**)d_input,N*sizeof(int));int*output(int*)malloc(block_num*sizeof(int));int*d_output;cudaMalloc((int**)d_output,block_num*sizeof(int));int*result(int*)malloc(block_num*sizeof(int));// 初始化值srand(666);for(inti0;iN;i){input[i]rand()%1000;}// cpu暴力计算for(inti0;iblock_num;i){intsum0;for(intj0;jBLOCK_SIZE;j){suminput[i*BLOCK_SIZEj];}result[i]sum;}cudaMemcpy(d_input,input,N*sizeof(int),cudaMemcpyHostToDevice);dim3Grid(block_num);dim3Block(BLOCK_SIZE);reduce3Grid,Block(d_input,d_output);cudaMemcpy(output,d_output,block_num*sizeof(int),cudaMemcpyDeviceToHost);if(!check(output,result,block_num)){for(inti0;iblock_num;i){if(i0output[i]!result[i]){printf(blockidx:%d,output:%d,ans:%d\n,i,output[i],result[i]);}}}else{printf(test passed!\n);}cudaFree(d_output);cudaFree(d_input);free(input);free(output);free(result);}intmain(){solve();return0;}221 u s 221us221us略有提高内存吞吐率和S M SMSM吞吐率都提高了5 55个点还是有点效果的。3.4 提升线程工作量线程块个数减半到这一步小优化已经提不了太多了需要从每个线程的操作着手。注意到按我们目前的思路第一轮只有一半的线程有计算任务第二轮只有四分之一的线程有计算任务以此类推有计算任务的线程越来越少剩下的线程都闲着那么一个优化思路就是来一波降本增效把不干活的人开了。我们可以让线程数减半每个线程的工作量翻倍这样虽然还是有线程闲着但是闲着的线程数减少了。这里有两个思路一个是线程块个数减半每个线程块大小不变另一个是线程块个数不变每个线程块大小减半。我们先实现第一个实现起来其实也比较简单从全局内存搬运到共享内存时我们就让每个线程搬两个假设有个位于线程块i ii,块内编号j jj的线程就让他搬运块i ii的j jj位置和块i 1 i1i1的j jj位置。然后在调用核函数时把线程块个数减半其余都不用改动。#includestdio.h#defineBLOCK_SIZE256// 增加每个线程的工作量减少block数__global__voidreduce4(int*d_input,int*d_output){int*startd_inputblockIdx.x*blockDim.x*2;__shared__intsdata[BLOCK_SIZE];sdata[threadIdx.x]start[threadIdx.x]start[threadIdx.xblockDim.x];__syncthreads();// 消除wrap分化for(intiBLOCK_SIZE/2;i1;i1){if(threadIdx.xi){sdata[threadIdx.x]sdata[threadIdx.xi];}__syncthreads();}if(threadIdx.x0){d_output[blockIdx.x]sdata[0];}}boolcheck(int*a,int*b,intN){for(inti0;iN;i){if(a[i]!b[i]){returnfalse;}}returntrue;}voidsolve(){// 申请host device侧空间intN32*1024*1024;intblock_numN/(BLOCK_SIZE*2);int*input(int*)malloc(N*sizeof(int));int*d_input;cudaMalloc((int**)d_input,N*sizeof(int));int*output(int*)malloc(block_num*sizeof(int));int*d_output;cudaMalloc((int**)d_output,block_num*sizeof(int));int*result(int*)malloc(block_num*sizeof(int));// 初始化值srand(666);for(inti0;iN;i){input[i]rand()%1000;}// cpu暴力计算for(inti0;iblock_num;i){intsum0;for(intj0;jBLOCK_SIZE*2;j){suminput[i*BLOCK_SIZE*2j];}result[i]sum;}cudaMemcpy(d_input,input,N*sizeof(int),cudaMemcpyHostToDevice);dim3Grid(block_num);dim3Block(BLOCK_SIZE);reduce4Grid,Block(d_input,d_output);cudaMemcpy(output,d_output,block_num*sizeof(int),cudaMemcpyDeviceToHost);if(!check(output,result,block_num)){for(inti0;iblock_num;i){if(i0output[i]!result[i]){printf(blockidx:%d,output:%d,ans:%d\n,i,output[i],result[i]);}}}else{printf(test passed!\n);}cudaFree(d_output);cudaFree(d_input);free(input);free(output);free(result);}intmain(){solve();return0;}146 u s 146us146us几乎减半这还是很有效的3.5 增加线程工作量线程块大小减半类似地我们也可以让线程块大小减半注意计算起始地址偏移量时用的是线程块大小乘二是因为这里的线程块大小实际上是线程块所对应的数据块大小的一半。后面同理搬运两个数据另一个数据的下标加上线程块大小。#includestdio.h#defineBLOCK_SIZE256// block数不变 减少block里的线程数__global__voidreduce5(int*d_input,int*d_output){int*startd_inputblockIdx.x*blockDim.x*2;__shared__intsdata[BLOCK_SIZE/2];sdata[threadIdx.x]start[threadIdx.x]start[threadIdx.xblockDim.x];__syncthreads();// 消除wrap分化for(intiBLOCK_SIZE/4;i1;i/2){if(threadIdx.xi){sdata[threadIdx.x]sdata[threadIdx.xi];}__syncthreads();}if(threadIdx.x0){d_output[blockIdx.x]sdata[0];}}boolcheck(int*a,int*b,intN){for(inti0;iN;i){if(a[i]!b[i]){returnfalse;}}returntrue;}voidsolve(){// 申请host device侧空间intN32*1024*1024;intblock_numN/BLOCK_SIZE;int*input(int*)malloc(N*sizeof(int));int*d_input;cudaMalloc((int**)d_input,N*sizeof(int));int*output(int*)malloc(block_num*sizeof(int));int*d_output;cudaMalloc((int**)d_output,block_num*sizeof(int));int*result(int*)malloc(block_num*sizeof(int));// 初始化值srand(666);for(inti0;iN;i){input[i]rand()%1000;}// cpu暴力计算for(inti0;iblock_num;i){intsum0;for(intj0;jBLOCK_SIZE;j){suminput[i*BLOCK_SIZEj];}result[i]sum;}cudaMemcpy(d_input,input,N*sizeof(int),cudaMemcpyHostToDevice);dim3Grid(block_num);dim3Block(BLOCK_SIZE/2);reduce5Grid,Block(d_input,d_output);cudaMemcpy(output,d_output,block_num*sizeof(int),cudaMemcpyDeviceToHost);if(!check(output,result,block_num)){for(inti0;iblock_num;i){if(i32output[i]!result[i]){printf(blockidx:%d,output:%d,ans:%d\n,i,output[i],result[i]);}}}else{printf(test passed!\n);}cudaFree(d_output);cudaFree(d_input);free(input);free(output);free(result);}intmain(){solve();return0;}146 u s 146us146us表现差不多3.6 不足一个warp时循环展开到这里其实已经没什么优化的了基本到极限了注意上一节的性能指标内存吞吐量已经到了极限了。剩下的优化如果不能改进内存读取的效率都没什么加速。当然理论上这个优化还是有意义的我们不妨实现一下。线程同步也是一个耗时的开销我们来考虑优化这里。注意到我们计算每一层由于计算任务可能被分成了多个w a r p warpwarp都要等到所有w a r p warpwarp都计算完才能计算下一层也就是需要一个同步。但最后有效位置不足一个w a r p warpwarp也就是不到32 3232个时只有第一个w a r p warpwarp在计算其它w a r p warpwarp根本不干活同步等待他们是无意义的。甚至去创建这些w a r p warpwarp都是无意义的。所以我们可以在不足32 3232个元素时不再进行循环而是把循环展开手动累加这样省去同步和循环的开销。也就是前32 3232个线程执行w a r p R e d u c e warpReducewarpReduce函数。这里初看可能有点迷惑因为实际上还是有数据依赖的每一层规约要全部执行玩才能进行下一层为什么我们这里没有同步也不会出错呢因为这些线程始终在一个w a r p warpwarp里每次都只会同步的执行一个指令也就是会先同步的执行a[tid] a[tid 32]再同步执行a[tid] a[tid 16]以此类推。w a r p warpwarp的s i m t simtsimt机制自动帮我们实现了同步。另一个可能迷惑的点是根据之前我们的思路第一轮有32 3232个线程有计算任务下一轮就只有16 1616个了以此类推每轮有计算任务的线程数减半但这里我们让前32 3232个线程都执行了完整的6 66轮规约计算这对吗实际上对结果无影响我们只是让一些原本无计算任务的线程也一起执行计算了比如第二轮16 − 31 16-3116−31这些线程原本不用求和了但现在我们也让他们进行计算了。这其实没有影响因为从结果的角度讲[ 16 , 31 ] [16,31][16,31]这些位置只有第一轮规约才会读取后面都不用读取了那它们的值变化了也无所谓了从性能的角度讲就算不给他们安排计算任务也要跟着w a r p warpwarp里的[ 0 , 15 ] [0,15][0,15]线程一起计算等待他们计算完成那就算让[ 15 , 31 ] [15,31][15,31]参与计算也不会导致额外的时间开销。这里还有一个陷阱w a r p R e d u c e warpReducewarpReduce函数内我们给传入的数组指针这个参数设置成了v o l a t i l e volatilevolatile也就是不可变这是告诉编译器不要对这个内存的操作进行优化。这是因为这个函数内对a [ t i d ] a[tid]a[tid]多次读取如果不加关键字编译器会优化成把a [ t i d ] a[tid]a[tid]这个位置的值取到寄存器里然后在寄存器内进行累加再写回内存这看似是对的但我们这里是多线程且有数据依赖的每一轮规约后a [ t i d ] a[tid]a[tid]的值都是会变的我们累计是要用的是最新的值而寄存器里的值相当于是一个缓存没有同步到最新的值结果就错了。加上不可变关键字强制每次都去共享内存重新读取才是对的。#includestdio.h#defineBLOCK_SIZE256// 有任务的线程不足一个warp时省去同步__device__voidwarpReduce(volatileint*a,inttid){a[tid]a[tid32];a[tid]a[tid16];a[tid]a[tid8];a[tid]a[tid4];a[tid]a[tid2];a[tid]a[tid1];}__global__voidreduce6(int*d_input,int*d_output){int*startd_inputblockIdx.x*blockDim.x*2;volatile__shared__intsdata[BLOCK_SIZE];sdata[threadIdx.x]start[threadIdx.x]start[threadIdx.xblockDim.x];__syncthreads();// 消除wrap分化for(intiBLOCK_SIZE/2;i32;i/2){if(threadIdx.xi){sdata[threadIdx.x]sdata[threadIdx.xi];}__syncthreads();}if(threadIdx.x32){warpReduce(sdata,threadIdx.x);}if(threadIdx.x0){d_output[blockIdx.x]sdata[0];}}boolcheck(int*a,int*b,intN){for(inti0;iN;i){if(a[i]!b[i]){returnfalse;}}returntrue;}voidsolve(){// 申请host device侧空间intN32*1024*1024;intblock_numN/(BLOCK_SIZE*2);int*input(int*)malloc(N*sizeof(int));int*d_input;cudaMalloc((int**)d_input,N*sizeof(int));int*output(int*)malloc(block_num*sizeof(int));int*d_output;cudaMalloc((int**)d_output,block_num*sizeof(int));int*result(int*)malloc(block_num*sizeof(int));// 初始化值srand(666);for(inti0;iN;i){input[i]rand()%1000;}// cpu暴力计算for(inti0;iblock_num;i){intsum0;for(intj0;jBLOCK_SIZE*2;j){suminput[i*BLOCK_SIZE*2j];}result[i]sum;}cudaMemcpy(d_input,input,N*sizeof(int),cudaMemcpyHostToDevice);dim3Grid(block_num);dim3Block(BLOCK_SIZE);reduce6Grid,Block(d_input,d_output);cudaMemcpy(output,d_output,block_num*sizeof(int),cudaMemcpyDeviceToHost);if(!check(output,result,block_num)){for(inti0;iblock_num;i){if(i0output[i]!result[i]){printf(blockidx:%d,output:%d,ans:%d\n,i,output[i],result[i]);}}}else{printf(test passed!\n);}cudaFree(d_output);cudaFree(d_input);free(input);free(output);free(result);}intmain(){solve();return0;}145 u s 145us145us几乎无优化因为内存带宽已经快用满了瓶颈在内存。不过这其实减小了计算量因此注意到计算吞吐率降低了。3.7 全部循环展开根据3.6 3.63.6的思路我们其实可以把全部循环都展开这样所有的同步都可以省去了循环的开销也可以省去。#includestdio.h#defineBLOCK_SIZE256__device__voidwarpReduce(volatileint*a,inttid){a[tid]a[tid32];a[tid]a[tid16];a[tid]a[tid8];a[tid]a[tid4];a[tid]a[tid2];a[tid]a[tid1];}__global__voidreduce7(int*d_input,int*d_output){int*startd_inputblockIdx.x*blockDim.x*2;volatile__shared__intsdata[BLOCK_SIZE];sdata[threadIdx.x]start[threadIdx.x]start[threadIdx.xblockDim.x];__syncthreads();// 消除wrap分化if(threadIdx.x128){sdata[threadIdx.x]sdata[threadIdx.x128];}__syncthreads();if(threadIdx.x64){sdata[threadIdx.x]sdata[threadIdx.x64];}__syncthreads();if(threadIdx.x32){warpReduce(sdata,threadIdx.x);}if(threadIdx.x0){d_output[blockIdx.x]sdata[0];}}boolcheck(int*a,int*b,intN){for(inti0;iN;i){if(a[i]!b[i]){returnfalse;}}returntrue;}voidsolve(){// 申请host device侧空间intN32*1024*1024;intblock_numN/(BLOCK_SIZE*2);int*input(int*)malloc(N*sizeof(int));int*d_input;cudaMalloc((int**)d_input,N*sizeof(int));int*output(int*)malloc(block_num*sizeof(int));int*d_output;cudaMalloc((int**)d_output,block_num*sizeof(int));int*result(int*)malloc(block_num*sizeof(int));// 初始化值srand(666);for(inti0;iN;i){input[i]rand()%1000;}// cpu暴力计算for(inti0;iblock_num;i){intsum0;for(intj0;jBLOCK_SIZE*2;j){suminput[i*BLOCK_SIZE*2j];}result[i]sum;}cudaMemcpy(d_input,input,N*sizeof(int),cudaMemcpyHostToDevice);dim3Grid(block_num);dim3Block(BLOCK_SIZE);reduce7Grid,Block(d_input,d_output);cudaMemcpy(output,d_output,block_num*sizeof(int),cudaMemcpyDeviceToHost);if(!check(output,result,block_num)){for(inti0;iblock_num;i){if(i0output[i]!result[i]){printf(blockidx:%d,output:%d,ans:%d\n,i,output[i],result[i]);}}}else{printf(test passed!\n);}cudaFree(d_output);cudaFree(d_input);free(input);free(output);free(result);}intmain(){solve();return0;}145 u s 145us145us没什么优化但是计算吞吐量有一定的降低计算任务确实变小了。这里应该是内存b o u n d boundbound了最后这几个优化确实是优化了计算量和同步时间的看起来没有效果主要是时间都在访存上了。4.番外 h100上测试4.0 朴素版本比4090 40904090还略慢神秘h 100 h100h100上n c u ncuncu要r o o t rootroot用不了只能n s y s nsysnsys凑合着用也有个核函数计时只是单位变成n s 1 e − 9 s ns1e-9sns1e−9s而n c u ncuncu是u s 1 e − 6 s us1e-6sus1e−6s换算一下即可相当于479 u s 479us479us比4090 40904090慢百分之三十左右。4.1 共享内存还是比4090慢429 u s 429us429us但是优化后加速比比4090 40904090要大怀疑是h 100 h100h100的H B M HBMHBM高速内存带来的加速因为我们换到共享内存上相当于是在优化访存H B M HBMHBM内存带宽大优化一点访存就能提速很多4.2 消除线程束分化260 u s 260us260us已经快追上4090了4.3 消除bank冲突199 u s 199us199us终于反超40904.4 增加线程工作量减少线程块个数109 u s 109us109us快了一倍比4090的140 u s 140us140us快很多了4.5 增加线程工作量减小线程块大小104 u s 104us104us还要更快一点后续优化以这个版本为基础4.6 循环展开最后一个warp87.4 u s 87.4us87.4us证明这个优化是有用的我们前面的4090没效果完全是因为内存b o u n d boundbound了4.7 全部循环展开87.1 u s 87.1us87.1us还是有提升可能是因为我们的线程块大小是256 256256线程数超过32 3232的规约轮数只有两轮128 , 64 128,64128,64而且这两轮的计算是在不同w a r p warpwarp的w a r p warpwarp间需要同步同步的时间不能省只剩去了f o r forfor循环的开销