【CUDA 】第3章 CUDA执行模型——3.5循环展开（3）

待解决的问题：

第三章 CUDA执行模型

3.5 循环展开

3.5.2 展开 线程的归约（最后一个warp内32个线程展开）

__syncthreads()是用于块内同步的，在归约核函数中，用来保证线程进入下一轮之前，每一轮的所有线程都把局部结果写入全局内存。

当只剩下一个线程束时（线程<=32），因为线程束执行是SIMT单指令多线程，每条指令之后有隐式的warp内同步过程。归约循环的最后6个迭代（最后一个warp内）可以用语句展开。

if(tid < 32){volatile int *vmem = idata;vmem[tid] += vmem[tid + 32];vmem[tid] += vmem[tid + 16];vmem[tid] += vmem[tid + 8];vmem[tid] += vmem[tid + 4];vmem[tid] += vmem[tid + 2];vmem[tid] += vmem[tid + 1];
}

这个线程束的展开避免了执行循环控制和线程同步逻辑。

volatile关键字：表明编译器对访问该变量的代码就不再进行优化，告诉编译器每次赋值的时候把vmem[tid]的值存回全局内存中。如果省略volatile关键字，就不能正常工作，因为编译器或缓存会对全局/共享内存优化读写。
如果全局/共享内存中的变量有volatile修饰，编译器假定它的值可以被其它线程任何时间修改使用，所以volatile修饰的变量强制读/写内存，而不是读写缓存或寄存器。

代码与reduceUnrolling8相比的区别：在写回最后的结果之前多了上面展开的代码，对最后一个warp块内进行手动线程规约，对应上面的就地归约中，限制条件也需要改成stride>32。

最后一个warp手动规约+一个线程块手动处理8个数据块

//展开线程规约4：手动展开八个块block的处理+最后一个块内手动规约
__global__ void reduceUnrollWarps8(int *g_idata, int *g_odata, unsigned int n){//设置线程idunsigned int tid = threadIdx.x;unsigned int idx = blockIdx.x * blockDim.x * 8 + threadIdx.x;//全局数据指针转换为当前块内数据指针int *idata = g_idata + blockIdx.x * blockDim.x * 8;//展开八个数据块，八个数据块加到第一个位置if(idx + 7 * blockDim.x < n){int a1 = g_idata[idx];int a2 = g_idata[idx + blockDim.x];int a3 = g_idata[idx + blockDim.x * 2];int a4 = g_idata[idx + blockDim.x * 3];int a5 = g_idata[idx + blockDim.x * 4];int a6 = g_idata[idx + blockDim.x * 5];int a7 = g_idata[idx + blockDim.x * 6];int a8 = g_idata[idx + blockDim.x * 7];g_idata[idx] = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8;}__syncthreads();//就地归约,一个数据块内for(int stride = blockDim.x/2; stride > 32; stride >>= 1){if(tid < stride){idata[tid] += idata[tid + stride];}__syncthreads();}//【新增！！！】展开warp，当归约到最后一个warp时，手动展开if(tid < 32){volatile int *vmem = idata;//volatile关键字说明编译器不对该变量优化vmem[tid] += vmem[tid + 32];vmem[tid] += vmem[tid + 16];vmem[tid] += vmem[tid + 8];vmem[tid] += vmem[tid + 4];vmem[tid] += vmem[tid + 2];vmem[tid] += vmem[tid + 1];}//把这个块的结果写回全局内容if(tid == 0) g_odata[blockIdx.x] = idata[0];
}

对应的main函数代码也需要修改：

//kernel7:reduceUnrollWarps8 展开规约5：一个线程块处理八个数据块cudaMemcpy(d_idata, h_idata, bytes, cudaMemcpyHostToDevice);cudaDeviceSynchronize();iStart = seconds();//因为现在每个线程块处理两个数据块，需要调整内核的执行配置，网格大小减为原来的一半reduceUnrollWarps8<<<grid.x/8, block>>>(d_idata, d_odata, size);//这里写错成了bytescudaDeviceSynchronize();iElaps = seconds() - iStart;cudaMemcpy(h_odata, d_odata, grid.x/8*sizeof(int), cudaMemcpyDeviceToHost);//需要除2gpu_sum = 0;for(int i = 0; i<grid.x/8; i++){//需要除gpu_sum += h_odata[i];}printf("gpu reduceUnrollWarps8  elapsed %f ms gpu_sum: %d <<<grid %d block %d>>>\n", iElaps, gpu_sum, grid.x/8, block.x);//需要除2

执行结果：

./3-3reduceInteger2 starting reduction atdevice 0: NVIDIA GeForce RTX 3090 with array size 16777216 grid 32768 block 512
cpu reduce         elapsed 0.041231 ms cpu_sum: 2139353471
gpu Warmup          elapsed 0.031330 ms gpu_sum: 0 <<<grid 32768 block 512>>>
gpu Neighbored          elapsed 0.000483 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu NeighboredLess      elapsed 0.000276 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu reduceInterleaved   elapsed 0.000239 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu reduceUnrolling2    elapsed 0.000156 ms gpu_sum: 2139353471 <<<grid 16384 block 512>>>
gpu reduceUnrolling4    elapsed 0.000125 ms gpu_sum: 2139353471 <<<grid 8192 block 512>>>
gpu reduceUnrolling8    elapsed 0.000118 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>
gpu reduceUnrollWarps8  elapsed 0.000119 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>

存在的问题：书上说比reduceUnrolling8快，但是实际上并没有，和书上结论不一致

reduceUnrollWarps8更快的原因是由于__syncthreads()的同步，更少线程束阻塞。
说明：__syncthreads()能减少核函数中的阻塞。

3.5.3 （交错归约）完全展开 的归约

如果编译时已知一个循环中的迭代次数，可以把循环完全展开。
因为Fermi或Kep1er架构中，每个块的最大线程数是1024，并且归约核函数中循环迭代次数是基于一个块block的维度。

把原来交错归约的步长自动迭代for循环具体手写。

    //【新增!!!】把就地归约步长迭代，完全展开if(blockDim.x >= 1024 && tid < 512) idata[tid] += idata[tid + 512];//线程块大小>1024时，前512个线程把1024个元素归约为512个元素__syncthreads();if(blockDim.x >= 512 && tid < 256) idata[tid] += idata[tid + 256];//线程块大小>512时，前256个线程把512个元素归约为256个元素__syncthreads();if(blockDim.x >= 256 && tid < 128) idata[tid] += idata[tid + 128];//线程块大小>256时，前128个线程把256个元素归约为128个元素__syncthreads();if(blockDim.x >= 128 && tid < 64) idata[tid] += idata[tid + 64];//线程块大小>128时，前64个线程把128个元素归约为64个元素__syncthreads();//tid此时已经<64了，<32的情况下面处理

手动展开交错归约+最后一个warp手动规约+一个线程块手动处理8个数据块

//完全展开的规约5：手动展开八个块block的处理+最后一个块内手动规约
__global__ void reduceCompleteUnrollWrap8(int *g_idata, int *g_odata, unsigned int n){//设置线程idunsigned int tid = threadIdx.x;unsigned int idx = blockIdx.x * blockDim.x * 8 + threadIdx.x;//全局数据指针转换为当前块内数据指针int *idata = g_idata + blockIdx.x * blockDim.x * 8;//展开八个数据块，八个数据块加到第一个位置if(idx + blockDim.x * 7 < n){int a1 = g_idata[idx];int a2 = g_idata[idx + blockDim.x];int a3 = g_idata[idx + blockDim.x * 2];int a4 = g_idata[idx + blockDim.x * 3];int a5 = g_idata[idx + blockDim.x * 4];int a6 = g_idata[idx + blockDim.x * 5];int a7 = g_idata[idx + blockDim.x * 6];int a8 = g_idata[idx + blockDim.x * 7];g_idata[idx] = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8;}__syncthreads();//【新增!!!】把就地归约步长迭代，完全展开if(blockDim.x >= 1024 && tid < 512) idata[tid] += idata[tid + 512];//线程块大小>1024时，前512个线程把1024个元素归约为512个元素__syncthreads();if(blockDim.x >= 512 && tid < 256) idata[tid] += idata[tid + 256];//线程块大小>512时，前256个线程把512个元素归约为256个元素__syncthreads();if(blockDim.x >= 256 && tid < 128) idata[tid] += idata[tid + 128];//线程块大小>256时，前128个线程把256个元素归约为128个元素__syncthreads();if(blockDim.x >= 128 && tid < 64) idata[tid] += idata[tid + 64];//线程块大小>128时，前64个线程把128个元素归约为64个元素__syncthreads();//tid此时已经<64了，<32的情况下面处理//展开warp，当归约到最后一个warp时，手动展开if(tid < 32){volatile int *vmem = idata;vmem[tid] += vmem[tid + 32];vmem[tid] += vmem[tid + 16];vmem[tid] += vmem[tid + 8];vmem[tid] += vmem[tid + 4];vmem[tid] += vmem[tid + 2];vmem[tid] += vmem[tid + 1];}//把这个块的结果写回全局内容if(tid == 0) g_odata[blockIdx.x] = idata[0];
}

执行结果：

./3-3reduceInteger2 starting reduction atdevice 0: NVIDIA GeForce RTX 3090 with array size 16777216 grid 32768 block 512
cpu reduce         elapsed 0.039004 ms cpu_sum: 2139353471
gpu Warmup          elapsed 0.038595 ms gpu_sum: 0 <<<grid 32768 block 512>>>
gpu Neighbored          elapsed 0.001485 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu NeighboredLess      elapsed 0.001385 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu reduceInterleaved   elapsed 0.001355 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu reduceUnrolling2    elapsed 0.000496 ms gpu_sum: 2139353471 <<<grid 16384 block 512>>>
gpu reduceUnrolling4    elapsed 0.000456 ms gpu_sum: 2139353471 <<<grid 8192 block 512>>>
gpu reduceUnrolling8    elapsed 0.000908 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>
gpu reduceUnrollWarps8  elapsed 0.000190 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>
gpu reduceCompleteUnrollWrap8  elapsed 0.000469 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>

书上说比reduceUnrollWarps8又有了小加速，实际运行并没有

3.5.4 模板函数 的归约

在手动展开交错归约的基础上，用模板函数进一步减少分支消耗，CUDA支持模板参数，指定块大小作为模板函数的参数。

与上一个代码的区别：
①参数模板iBlockSize代替原来的块大小blockDim.x。检查块大小的if语句在编译的时候会被评估，如果该条件为false，编译时就删除，循环更有效率。
例如：块大小永远是256，下面语句永远是false，直接不进入该归约步骤，编译器自动从执行内核中移除它。
②调用的时候用switch-case语句，编译器可以为不同大小的代码块自动优化代码。

iBlockSize >= 1024 && tid < 512

关键代码修改：

    //【新增改变!!!】把就地归约步长迭代，完全展开+iBlockSize代替原来if(iBlockSize >= 1024 && tid < 512) idata[tid] += idata[tid + 512];//线程块大小>1024时，前512个线程把1024个元素归约为512个元素__syncthreads();if(iBlockSize >= 512 && tid < 256) idata[tid] += idata[tid + 256];//线程块大小>512时，前256个线程把512个元素归约为256个元素__syncthreads();if(iBlockSize >= 256 && tid < 128) idata[tid] += idata[tid + 128];//线程块大小>256时，前128个线程把256个元素归约为128个元素__syncthreads();if(iBlockSize >= 128 && tid < 64) idata[tid] += idata[tid + 64];//线程块大小>128时，前64个线程把128个元素归约为64个元素__syncthreads();//tid此时已经<64了，<32的情况下面处理

模板函数+手动展开交错归约+最后一个warp手动规约+一个线程块手动处理8个数据块

//完全展开的规约6：模板参数替换块大小+手动展开八个块block的处理+最后一个块内手动规约
template <unsigned int iBlockSize>
__global__ void reduceCompleteUnroll(int *g_idata, int *g_odata, unsigned int n){//设置线程idunsigned int tid = threadIdx.x;unsigned int idx = blockIdx.x * blockDim.x * 8 + threadIdx.x;//全局数据指针转换为当前块内数据指针int *idata = g_idata + blockIdx.x * blockDim.x * 8;//展开八个数据块，八个数据块加到第一个位置if(idx + blockDim.x * 7 < n){int a1 = g_idata[idx];int a2 = g_idata[idx + blockDim.x];int a3 = g_idata[idx + blockDim.x * 2];int a4 = g_idata[idx + blockDim.x * 3];int a5 = g_idata[idx + blockDim.x * 4];int a6 = g_idata[idx + blockDim.x * 5];int a7 = g_idata[idx + blockDim.x * 6];int a8 = g_idata[idx + blockDim.x * 7];g_idata[idx] = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8;}__syncthreads();//【新增改变!!!】把就地归约步长迭代，完全展开+iBlockSize代替原来if(iBlockSize >= 1024 && tid < 512) idata[tid] += idata[tid + 512];//线程块大小>1024时，前512个线程把1024个元素归约为512个元素__syncthreads();if(iBlockSize >= 512 && tid < 256) idata[tid] += idata[tid + 256];//线程块大小>512时，前256个线程把512个元素归约为256个元素__syncthreads();if(iBlockSize >= 256 && tid < 128) idata[tid] += idata[tid + 128];//线程块大小>256时，前128个线程把256个元素归约为128个元素__syncthreads();if(iBlockSize >= 128 && tid < 64) idata[tid] += idata[tid + 64];//线程块大小>128时，前64个线程把128个元素归约为64个元素__syncthreads();//tid此时已经<64了，<32的情况下面处理//展开warp，当归约到最后一个warp时，手动展开if(tid < 32){volatile int *vmem = idata;vmem[tid] += vmem[tid + 32];vmem[tid] += vmem[tid + 16];vmem[tid] += vmem[tid + 8];vmem[tid] += vmem[tid + 4];vmem[tid] += vmem[tid + 2];vmem[tid] += vmem[tid + 1];}//把这个块的结果写回全局内容if(tid == 0) g_odata[blockIdx.x] = idata[0];
}

对应main函数调用也修改：

    //kernel9:reduceCompleteUnroll 展开规约5：一个线程块处理八个数据块//调用方式修改为在switch case结构cudaMemcpy(d_idata, h_idata, bytes, cudaMemcpyHostToDevice);cudaDeviceSynchronize();iStart = seconds();//改变【新增！！！】允许编译器为特定线程块大小自动优化代码switch(blocksize){case 1024:reduceCompleteUnroll<1024><<<grid.x/8, block>>>(d_idata, d_odata, size);break;case 512:reduceCompleteUnroll<512><<<grid.x/8, block>>>(d_idata, d_odata, size);break;case 256:reduceCompleteUnroll<256><<<grid.x/8, block>>>(d_idata, d_odata, size);break;case 128:reduceCompleteUnroll<128><<<grid.x/8, block>>>(d_idata, d_odata, size);break;case 64:reduceCompleteUnroll<64><<<grid.x/8, block>>>(d_idata, d_odata, size);break;}//因为现在每个线程块处理两个数据块，需要调整内核的执行配置，网格大小减为原来的一半//reduceCompleteUnroll<<<grid.x/8, block>>>(d_idata, d_odata, size);//这里写错成了bytescudaDeviceSynchronize();iElaps = seconds() - iStart;cudaMemcpy(h_odata, d_odata, grid.x/8*sizeof(int), cudaMemcpyDeviceToHost);//需要除2gpu_sum = 0;for(int i = 0; i<grid.x/8; i++){//需要除gpu_sum += h_odata[i];}printf("gpu reduceCompleteUnroll  elapsed %f ms gpu_sum: %d <<<grid %d block %d>>>\n", iElaps, gpu_sum, grid.x/8, block.x);//需要除2

执行结果：

./3-3reduceInteger2 starting reduction atdevice 0: NVIDIA GeForce RTX 3090 with array size 16777216 grid 32768 block 512
cpu reduce         elapsed 0.040380 ms cpu_sum: 2139353471
gpu Warmup          elapsed 0.053873 ms gpu_sum: 0 <<<grid 32768 block 512>>>
gpu Neighbored          elapsed 0.000478 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu NeighboredLess      elapsed 0.000273 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu reduceInterleaved   elapsed 0.000235 ms gpu_sum: 2139353471 <<<grid 32768 block 512>>>
gpu reduceUnrolling2    elapsed 0.000157 ms gpu_sum: 2139353471 <<<grid 16384 block 512>>>
gpu reduceUnrolling4    elapsed 0.000124 ms gpu_sum: 2139353471 <<<grid 8192 block 512>>>
gpu reduceUnrolling8    elapsed 0.000117 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>
gpu reduceUnrollWarps8  elapsed 0.000120 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>
gpu reduceCompleteUnrollWrap8  elapsed 0.000119 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>
gpu reduceCompleteUnroll  elapsed 0.000119 ms gpu_sum: 2139353471 <<<grid 4096 block 512>>>

最大的性能增益是reduceUnrolling8核函数获得的，每个线程规约前处理8个数据块，有8个独立的内存访问，更好让内存带宽饱和，并且还能隐藏加载/存储延迟。