cuda逐步优化实现reduce sum 操作

归约是一种常见的数据并行原语，它将数组中的元素通过某种二元操作（如加法）合并成一个单一的值。通过逐步展示不同的CUDA实现版本，来演示重要的优化策略。
由于规约的算术操作很简单，对算力要求不高，因此我们逐步优化目标是尽可能达到最高的带宽利用率，基本想法是：

• 树状归约方法：在每个线程块内使用基于树的方法进行局部归约，然后需要处理如何跨线程块通信部分结果。

• 全局同步问题：CUDA没有全局同步机制，因为这样做在硬件上成本高昂，并且会限制程序运行的线程块数量，影响整体效率。

• 内核分解：通过分解计算为多个内核调用来避免全局同步，内核启动点作为全局同步点，具有较低的硬件和软件开销。

• 优化目标：对于归约操作，由于其算术强度很低（每个加载的元素仅有一次浮点操作），优化目标是达到峰值带宽。

基础实现

__global__ void reduceSum(int *g_idata, int* g_odata)
{
    extern __shared__ int sdata[];
    uint tid = threadIdx.x;
    uint i = blockIdx.x*blockDim.x+threadIdx.x;

    sdata[tid] = g_idata[i];
    // printf("blockIdx=%d,sdata[%d]=%d ",blockIdx.x,tid,sdata[tid]);

    __syncthreads();

    for(uint s=1; s<blockDim.x; s*=2){
        if (tid %(2*s) == 0){
            sdata[tid] += sdata[tid+s];
        }
        __syncthreads();
    }
    if(tid ==0) {
        g_odata[blockIdx.x] = sdata[0];
        // atomicAdd(g_odata, sdata[0]);
    }
}

Warp thread divergent

在 CUDA 编程中，高度发散的 warps 和使用 %（取模）运算符都会对性能产生负面影响。

高度发散的 Warps Warp 是 CUDA 中的一个基本执行单元。一个 Warp 包含 32 个线程，这些线程在同一个流多处理器（SM）中并行执行相同的指令。
如果一个 Warp 中的所有线程都执行相同的指令，则 Warp 是一致的，性能最好。
Warp 发散 发生在同一个 Warp 中的线程执行不同的指令路径时。通常是因为条件分支语句（如 if-else）导致不同线程走不同的代码路径。

• 当 Warp 发散时，CUDA 硬件必须序列化不同的执行路径。这意味着，虽然所有线程在逻辑上是并行的，但实际上它们不得不逐路径地执行不同的指令，这大大降低了并行效率。
• 举例来说，如果一个 Warp 中一半的线程执行一个路径，另一半执行另一个路径，那么两个路径将被顺序执行，每个路径只利用了一半的线程，效率降低。

% 运算符很慢

• % 运算符在很多硬件架构上实现起来比较复杂和耗时，因为它通常需要进行除法运算，而除法比加法、减法和乘法慢很多。
• 在 CUDA 编程中，特别是对于 GPU 的流多处理器（SM）来说，整数除法和取模操作更为耗时，因为这些操作需要更多的时钟周期来完成。

解决方案

• 减少 Warp 发散
  • 最小化条件分支：尽量减少 if-else 语句的使用，特别是在 Warp 内部。
  • 数据重构：尝试重构数据，使得同一个 Warp 中的线程能够执行相同的指令。
  • 避免复杂的条件判断：如果条件判断无法避免，尝试使用其它算法或数据结构来最小化发散。
• 优化取模操作
  • 使用位操作：如果取模的数是 2 的幂，可以使用位操作来代替 %。例如，x % 4 可以替换为 x & 3。
  • 查找表：对于小范围的取模操作，可以使用查找表来替代计算。
  • 简化算法：如果可能，重构算法以减少或避免取模操作。

__global__ void reduceSum1(int *g_idata, int* g_odata)
{
    extern __shared__ int sdata[];
    uint tid = threadIdx.x;
    uint i = blockIdx.x*blockDim.x+threadIdx.x;

    sdata[tid] = g_idata[i];
    // printf("blockIdx=%d,sdata[%d]=%d ",blockIdx.x,tid,sdata[tid]);

    __syncthreads();

    for(uint s=1; s<blockDim.x; s*=2){
        int index = 2*s*tid;
        
        if (index< blockDim.x){
            sdata[tid] += sdata[tid+s];
        }
        __syncthreads();
    }
    if(tid ==0) {
        // g_odata[blockIdx.x] = sdata[0];
        atomicAdd(g_odata, sdata[0]);
    }
}

Bank conflict

在归约过程中，我们确保访问模式没有bank冲突。在 CUDA 设备上，默认情况下每个内存Bank有 32 个 int 元素，因此访问模式 sdata[tid] 和 sdata[tid + s] 在没有填充的情况下通常不会导致bank冲突，特别是当 blockDim.x 是 2 的幂时。

// solve bank conflict, use sequence addressing
__global__ void reduceSum2(int *g_idata, int* g_odata)
{
    extern __shared__ int sdata[];
    uint tid = threadIdx.x;
    uint i = blockIdx.x*blockDim.x+threadIdx.x;

    sdata[tid] = g_idata[i];
    // printf("blockIdx=%d,sdata[%d]=%d ",blockIdx.x,tid,sdata[tid]);

    __syncthreads();

    for(uint s=blockDim.x/2; s>0; s>>=1){
        if (tid< s){
            sdata[tid] += sdata[tid+s];
        }
        __syncthreads();
    }
    if(tid ==0) {
        // g_odata[blockIdx.x] = sdata[0];
        atomicAdd(g_odata, sdata[0]);
    }
}

Unroll last warp

指令开销（Instruction Overhead） 是指那些与核心计算无直接关系的辅助指令的开销。这些辅助指令包括地址计算、循环控制、条件判断等。尽管这些指令在 CPU 编程中可能看起来很轻量，但在 GPU 上，由于大量线程的并行执行，甚至少量的开销也会累积成显著的性能瓶颈。

如何缓解这些开销

1. 减少地址计算
   尽量减少复杂的索引计算，将索引计算移出循环体或频繁执行的代码块中。
   使用共享内存来存储中间结果，减少全局内存访问的复杂索引计算。
2. 优化循环
   尽量减少循环的层数和每次迭代的复杂度。
   使用 unrolling 技术手动展开循环，减少循环控制指令的开销。
3. 减少分支和条件判断
   避免 Warp 分歧，尽量减少条件分支。
   使用更简单的逻辑和数据结构来避免复杂的条件判断。

// unroll last warp, minimize loop and other condition code, and address arithmetric instruction 
__global__ void reduceSum3(int *g_idata, int* g_odata)
{
    extern __shared__ int sdata[];
    uint tid = threadIdx.x;
    uint i = blockIdx.x*blockDim.x+threadIdx.x;

    sdata[tid] = g_idata[i];
    // printf("blockIdx=%d,sdata[%d]=%d ",blockIdx.x,tid,sdata[tid]);

    __syncthreads();

    for(uint s=blockDim.x/2; s>32; s>>=1){
        if (tid< s){
            sdata[tid] += sdata[tid+s];
        }
        __syncthreads();
    }

    //last 32 thread all in the same warp, so do not need to use syncthread in a single warp
    // this saves useless work in all warps, not just the last one
    // without unrolling, all warps execute every iteration of the for loop and if statement
    if (tid<32) {
    	//use volatile to avoid compiler optimization, keep 
    	//fetch value from the memory(not the register) everytime
        volatile int* vdmem = sdata;
        vdmem[tid] += vdmem[tid+32];
        vdmem[tid] += vdmem[tid+16];
        vdmem[tid] += vdmem[tid+8];
        vdmem[tid] += vdmem[tid+4];
        vdmem[tid] += vdmem[tid+2];
        vdmem[tid] += vdmem[tid+1];
    };

    if(tid ==0) {
        // g_odata[blockIdx.x] = sdata[0];
        atomicAdd(g_odata, sdata[0]);
    }
}

Unroll all loop

template<unsigned int blockSize>
__device__ void warp_reduce(volatile int* sdata, int tid)
{
    sdata[tid] += sdata[tid+32];
    sdata[tid] += sdata[tid+16];
    sdata[tid] += sdata[tid+8];
    sdata[tid] += sdata[tid+4];
    sdata[tid] += sdata[tid+2];
    sdata[tid] += sdata[tid+1];
}

//unroll all loop code, use template to keep kernel be generic
template<unsigned int blockSize>
__global__ void reduceSum4(int *g_idata, int* g_odata)
{
    extern __shared__ int sdata[];
    uint tid = threadIdx.x;
    uint i = blockIdx.x*blockDim.x+threadIdx.x;
    uint gridSize = 2*blockSize*gridDim.x;
    sdata[tid] = 0;
    uint n = 4000000;
    while(i<n){
        sdata[tid] +=g_idata[i] + g_idata[i+blockSize];
        i +=gridSize;
    }
    __syncthreads();

    //blocksize condition judgement will be evaluated at compile time
    if (blockSize >= 1024) { if (tid < 512) { sdata[tid] += sdata[tid + 512]; } __syncthreads(); }
    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(); }

    //last 32 thread all in the same warp, so do not need to use syncthread in a single warp
    // this saves useless work in all warps, not just the last one
    // without unrolling, all warps execute every iteration of the for loop and if statement
    if (tid<32) warp_reduce<blockSize>(sdata, tid);

    if(tid ==0) {
        // g_odata[blockIdx.x] = sdata[0];
        atomicAdd(g_odata, sdata[0]);
    }
}

• 线程步长和循环遍历：
  每个线程在初始化时会根据 blockIdx.x 和 threadIdx.x 计算其全局索引 i。随后，通过 while (i < n) 循环确保每个线程在遍历其负责的所有数据。i += gridSize 会使线程跳到其下一个负责的数据块。

• 合并内存访问：
  当每个线程遍历数据时，gridSize 作为步长使得每个线程块中的线程能够均匀分布在全局数据范围内，确保内存访问的合并。

• 减少内存访问次数：每个线程在while循环中尽可能多地从全局内存加载数据到共享内存，减少了对全局内存的访问次数。全局内存访问通常比共享内存访问要慢得多，因此减少全局内存访问可以提高性能（如果输入gridDim是一维的，while 循环只能执行一次，并不能累积数据）。

• 提高内存访问的效率：通过while循环，每个线程可以加载更多的数据，这样可以通过增加每个线程的工作量来提高内存访问的效率。在CUDA编程中，内存访问的效率（如内存吞吐量）对于性能至关重要。

• 减少内核启动开销：通过让每个线程做更多的工作，可以减少需要启动的内核数量。内核启动有一定的开销，减少内核数量可以减少这种开销。

• 算法级联（Algorithm Cascading）：这种方法实际上是算法级联的一个应用，即将顺序算法和并行算法结合起来。每个线程加载多个元素，然后在共享内存中进行树状归约。这种方法可以减少递归内核调用的层级，从而减少内核启动的开销。

• 保持内存访问的连贯性：在while循环中，通过使用gridSize作为步长，可以保持内存访问的连贯性（coalescing），这有助于进一步提高内存访问的性能。

• 减少线程空闲时间：在前面的kernel中，有一半的线程在第一次循环迭代时是空闲的。而在这个kernel中，通过while循环确保了所有线程都在忙碌地执行工作，从而更充分地利用了GPU的计算资源。

Average elasped time:(0.000252) second, N size:(4000000), bandwidth:(63.516257 GB/s)
Average elasped time:(0.000200) second, N size:(4000000), bandwidth:(80.102532 GB/s)
Average elasped time:(0.000197) second, N size:(4000000), bandwidth:(81.168833 GB/s)
Average elasped time:(0.000161) second, N size:(4000000), bandwidth:(99.423344 GB/s)
Average elasped time:(0.000172) second, N size:(4000000), bandwidth:(92.885010 GB/s)