【CUDA】Reduce归约求和（下）

前言

学习 UP 主比飞鸟贵重的多_HKL 的【CUDA】Reduce规约求和（已完结~）视频，记录下个人学习笔记，仅供自己参考😄

refer1：【CUDA】Reduce规约求和（已完结~）

refer2：Optimizing Parallel Reduction in CUDA

refer3：https://github.com/BBuf/how-to-optim-algorithm-in-cuda/tree/master/reduce

refer4：深入浅出GPU优化系列：reduce优化

refer5：https://chatgpt.com/

1. 优化技巧4：展开最后一个warp减少同步

我们接着上篇文章来讲解 reduce 的优化，上篇文章最后我们提到 reduce_v3 中最后几轮迭代时对于空闲的线程还在进行 __syncthreads() 同步，这造成了浪费

为了减少线程同步带来的指令开销，NVIDIA 博客中给出的优化技巧 4 的代码如下：

__device__ void warpReduce(volatile float* sdata, unsigned 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];
}

__global__ void reduce_v4(float* g_idata, float* g_odata){
    extern __shared__ float sdata[];

    // 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
    for(unsigned int s = blockDim.x / 2; s > 32; s >>= 1){
        if(tid < s){
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    // write result for this block to global mem
    if(tid < 32){
        warpReduce(sdata, tid);
    }
    if(tid == 0){
        g_odata[blockIdx.x] = sdata[0];
    }
}

相比于 reduce_v3 的做法，reduce_v4 主要是将原来的：

// do reduction in shared mem
for(unsigned int s = blockDim.x / 2; s > 0; s >>= 1){
    if(tid < s){
        sdata[tid] += sdata[tid + s];
    }
    __syncthreads();
}

替换为了：

// do reduction in shared mem
for(unsigned int s = blockDim.x / 2; s > 32; s >>= 1){
    if(tid < s){
        sdata[tid] += sdata[tid + s];
    }
    __syncthreads();
}

// write result for this block to global mem
if(tid < 32){
    warpReduce(sdata, tid);
}

其他部分的代码保持不变，下面我们来具体分析下 NVIDIA 的优化思路：

1. 缩短归约循环范围

原来的归约循环一直执行到 s > 0，每一次迭代都调用 __syncthreads() 进行全局同步
在 reduce_v4 中，将归约循环的条件修改为 s > 32，即只有还有超过一个 warp 的活跃线程时进行同步归约

2. 利用 warp 内同步特性

当 s <= 32 时，活跃的线程全部位于同一 warp 中，而同一 warp 内的线程是同步执行的（位于一个 SIMD 单元上），不再需要显式调用 __syncthreads()
这里可以避免在最后几轮归约中多余的同步操作，从而节省了大量的指令开销

3. 使用 unrolled 的 warpReduce 函数

在归约循环结束后，对于最后 32 个线程，使用一个专门的 warpReduce 函数来完成剩余的归约操作
该函数通过直接对共享内存中的值进行累加，实现了高效的无同步归约

总的来说，在 reduce_v4 中 NVIDIA 通过缩短归约循环范围并利用 warp 内同步，避免了最后几轮归约中空闲线程的同步开销，从而提升了整体性能

此外我们还发现 warpReduce 函数中对于 sdata 参数还加入了一个 volatile 修饰符，它的作用是什么呢？为什么要加入呢？🤔

在 warpReduce 函数中，volatile 修饰符的作用是告诉编译器，该内存区域可能会在不同线程间发生变化，因此不要对它的读写操作进行缓存或重排序优化

具体来说，volatile 修饰符包含两个功能：

保证内存可见性：在没有显式同步操作（__syncthreads()）的 warp 内部归约中，所有线程处于同步状态，但编译器仍可能优化内存访问。使用 volatile 可以确保每次对 sdata 的读写都直接从共享内存中完成，这样线程可以看到其他线程最新更新的值
避免编译器优化问题：在 warp 内进行连续的归约操作时，编译器可能会将数据保存在寄存器中，而不是每次都重新加载共享内存。volatile 强制每次都访问共享内存，从而避免由于优化导致的错误累加或数据不一致问题。

假设没有 volatile 修饰符，编译器为了优化可能会将共享内存变化缓存到寄存器中，例如在第一轮归约时，线程 0 将 sdata[0] 和 sdata[32] 相加，并将结果存入寄存器中，但这个更新可能没有立即写回共享内存，而在进行第二轮归约时，线程 0 再次读取 sdata[0] 时，由于编译器优化，重读共享内存时未能获取刚才更新后的值，而是读取的旧值，从而导致使用了错误的数值，具体大家可以参考：CUDA: In warp reduction and volatile keyword

博主测试了如果不加 volatile 修饰符时输出如下：

output_host[0] = 112.00
output_host[1] = 112.00
output_host[2] = 112.00
output_host[3] = 112.00
output_host[4] = 112.00

性能和带宽的测试情况如下：

优化手段	耗时(us)	Memory throughout(%)	DRAM throughout(%)	加速比
reduce_baseline	2490us	62.01%	15.44%	~
reduce_v1_interleaved_addressing	1800us	85.78%	21.35%	1.38
reduce_v2_bank_conflict_free	1730us	89.24%	22.01%	1.44
reduce_v3_idle_threads_free	896.19us	89.48%	42.89%	2.78
reduce_v4_unroll_last_warp	507.71us	76.08%	76.08%	4.90

2. 优化技巧5：完全展开循环

其实到了 reduce_v4 这一步，reduce 的效率已经足够高了，再进一步优化其实已经非常困难了。为了探索极致的性能表现，NVIDIA 博客中给出了对 for 循环进行完全展开的方案

代码如下：

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

template <unsigned int blockSize>
__global__ void reduce_v5(float* g_idata, float* g_odata){
    extern __shared__ float sdata[];

    // 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){
            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();
    }

    // write result for this block to global mem
    if(tid < 32){
        warpReduce<blockSize>(sdata, tid);
    }
    if(tid == 0){
        g_odata[blockIdx.x] = sdata[0];
    }
}

在 reduce_v5 中，NVIDIA 采用了模板和完全展开技术，将归约过程中的循环和条件判断在编译时就“摊平”，从而极大地减少了运行时的开销，详细来说：

1. 使用模板参数实现编译时常量

通过模板参数 <unsigned int blockSize> 使得 blockSize 成为编译时常量。这样所有关于 blockSize 的 if 条件（例如 if(blockSize >= 512)）都会在编译期间确定，从而使得代码可以完全展开，无需在运行时进行判断

2. 分阶段归约，逐步减少数据量

初始阶段，每个线程从全局内存加载两个元素并预归约，然后对于较大的 blockSize，先进行大步归约，当归约到 32 个线程以内时，调用 unrolled 的 warpReduce 函数完成最后的归约，而这时由于所有线程都在同一个 warp 内，因此不需要 __syncthreads()

3. 完全展开消除循环和同步开销

采用 if 语句和模板展开后，每一步归约操作都显式写出，而不是使用 for 循环。这样不仅减少了循环控制的开销，还避免了在归约末尾额外的分支和同步指令
最终，这种方式利用编译时展开将“硬编码”成一系列连续的加法操作，从而让编译器可以生成更加高效的汇编代码

性能和带宽的测试情况如下：

优化手段	耗时(us)	Memory throughout(%)	DRAM throughout(%)	加速比
reduce_baseline	2490us	62.01%	15.44%	~
reduce_v1_interleaved_addressing	1800us	85.78%	21.35%	1.38
reduce_v2_bank_conflict_free	1730us	89.24%	22.01%	1.44
reduce_v3_idle_threads_free	896.19us	89.48%	42.89%	2.78
reduce_v4_unroll_last_warp	507.71us	76.08%	76.08%	4.90
reduce_v5_completely_unroll	485.54us	79.48%	79.48%	5.13

从测试情况可以看到 reduce_v5 相比 reduce_v4 收益较少，这主要是因为 NVIDIA 博客中测试的设备比较老，而目前的 GPU 硬件架构在不断发展，另外 NVIDIA 在编译器上也做了较多的优化工作，所以这里测试出的加速比没有博客中的那么明显

3. 优化技巧6：调节GridSize和BlockSize

当走到这一步的时候，能调的东西基本上已经调完了，最后让我们把目光放在 block 和 thread 数量的设置上，之前在 reduce_v5 中我们默认一个 block 开启 256 个线程，负责 512 个元素的 reduce 操作，那我们可以让每个 block 处理的数量更多一点，这样开启的 block 数量就会少一些，以此对 block 设置进行调整，获得最优的 block 取值，这样或许能够带来一些性能上的收益

NVIDIA 博客中给出的优化技巧 6 的代码如下：

#include <cuda.h>
#include <cuda_runtime.h>
#include <time.h>
#include <stdio.h>

#define N 32*1024*1024
#define BLOCK_SIZE 256

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

template <unsigned int blockSize, int NUM_PER_THREAD>
__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 < NUM_PER_THREAD; iter++){
        sdata[tid] += g_idata[i + iter * blockSize];
    }
    __syncthreads();

    // do reduction in shared mem
    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();
    }

    // write result for this block to global mem
    if(tid < 32){
        warpReduce<blockSize>(sdata, tid);
    }
    if(tid == 0){
        g_odata[blockIdx.x] = sdata[0];
    }    
}

int main(){

    // input host
    float* input_host = (float*)malloc(N * sizeof(float));
    for(int i = 0; i < N; ++i){
        input_host[i] = 2.0;
    }

    // input device
    float* input_device;
    cudaMalloc((void**)&input_device, N * sizeof(float));

    // copy data from host to device
    cudaMemcpy(input_device, input_host, N * sizeof(float), cudaMemcpyHostToDevice);
    
    const int block_num = 1024;
    // 每个 block 要处理的数据量
    const int NUM_PER_BLOCK = N / block_num;
    // 每个 thread 要处理的数据量
    const int NUM_PER_THREAD = NUM_PER_BLOCK / BLOCK_SIZE;
    float* output_host = (float*)malloc(block_num * sizeof(float));

    // output device
    float* output_device;
    cudaMalloc((void**)&output_device, block_num * sizeof(float));

    // kernel launch
    dim3 grid(block_num, 1);
    dim3 block(BLOCK_SIZE, 1);
    reduce_v6<BLOCK_SIZE, NUM_PER_THREAD><<<grid, block>>>(input_device, output_device);

    // copy result from device to host
    cudaMemcpy(output_host, output_device, block_num * sizeof(float), cudaMemcpyDeviceToHost);

    // print some result
    for(int i = 0; i < 5; ++i){
        float data = output_host[i];
        printf("output_host[%d] = %.2f\n", i, data);
    }

    return 0;
}

在 reduce_v6 中，NVIDIA 的优化重点在于通过调整每个 block 处理的数据量来提高硬件资源的利用率，具体思路包括：

1. 增加每个 block 处理的数据量

原来版本（reduce_v5）：每个 block 处理 512 个元素（每个线程处理 2 个元素，加载两个数据相加）
当前优化（reduce_v6）：每个 block 处理的数据量更大，由 NUM_PER_BLOCK 定义，每个线程也不再只处理一个数据对，而是通过一个循环加载 NUM_PER_THREAD 个数据块，每次以 blockSize 为步长访问全局内存

2. 利用循环展开（Loop Unrolling）

使用 #pragma unroll 告诉编译器展开 for 循环，这样循环内的指令可以在编译时直接生成连续的加法指令，减少了循环控制开销，提高执行效率

3. 调整 block 和 thread 数量的权衡

reduce_v6 通过让每个线程加载并预处理多个数据元素，使用循环展开以及合理调整 block 和 thread 配置，在降低全局内存访问次数的同时减少了 kernel 启动开销，从而进一步提升了归约操作的整体性能。

性能和带宽的测试情况如下：

优化手段	耗时(us)	Memory throughout(%)	DRAM throughout(%)	加速比
reduce_baseline	2490us	62.01%	15.44%	~
reduce_v1_interleaved_addressing	1800us	85.78%	21.35%	1.38
reduce_v2_bank_conflict_free	1730us	89.24%	22.01%	1.44
reduce_v3_idle_threads_free	896.19us	89.48%	42.89%	2.78
reduce_v4_unroll_last_warp	507.71us	76.08%	76.08%	4.90
reduce_v5_completely_unroll	485.54us	79.48%	79.48%	5.13
reduce_v6_multi_add	400.06us	96.94%	96.94%	6.22

对于 block 的取值理论上来说取 SM 数量的倍数会比较合理，这里博主没有详细测试，表中的数据设置的 block 数量为 1024，可以看到此时带宽利用率得到了较大的提升

4. 优化技巧7：使用shuffle指令

其实对于 reduce 的优化 NVIDIA 的博客到这里就结束了，但是 NVIDIA 后面出了 shuffle 指令，对于 reduce 的优化有着非常好的效果。目前绝大多数访存类算子像 softmax、batchnorm、reduce 等都是用 shuffle 实现的，所以这里简单聊下怎么把 shuffle 指令用在 reduce 的优化上

Note：之前博主在学习 llama.cpp 中 RMSNorm 的实现时 warp 内的归约求和就是通过 shuffle 指令实现的，大家感兴趣的可以看看：RMSNorm算子的CUDA实现

shuffle 指令是一组针对 warp 的指令，它最重要的特性就是可以让 warp 内的寄存器相互访问，在没有 shuffle 指令的时候，各个线程在进行通信时只能通过 shared memory 来访问彼此的寄存器。而采用了 shuffle 指令之后，warp 内的线程可以直接对其他线程的寄存器进行访存，通过这种方式可以减少访存的延时

代码实现如下：

#include <cuda.h>
#include <cuda_runtime.h>
#include <time.h>
#include <stdio.h>

#define N 32*1024*1024
#define BLOCK_SIZE 256
#define WARP_SIZE 32

template <unsigned int blockSize>
__device__ __forceinline__ float warpReduceSum(float sum){
    if(blockSize >= 32) sum += __shfl_down_sync(0xffffffff, sum, 16);
    if(blockSize >= 16) sum += __shfl_down_sync(0xffffffff, sum, 8);
    if(blockSize >= 8) sum += __shfl_down_sync(0xffffffff, sum, 4);
    if(blockSize >= 4) sum += __shfl_down_sync(0xffffffff, sum, 2);
    if(blockSize >= 2) sum += __shfl_down_sync(0xffffffff, sum, 1);
    return sum;
}

template <unsigned int blockSize, int NUM_PER_THREAD>
__global__ void reduce_v7(float* g_idata, float* g_odata){
    float sum = 0;

    // each thread loads NUM_PER_THREAD element from global
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockSize * NUM_PER_THREAD) + threadIdx.x;
    #pragma unroll
    for(int iter = 0; iter < NUM_PER_THREAD; iter++){
        sum += g_idata[i + iter * blockSize];
    }
    
    // shared memory for partial sums (one per warp in the block)
    static __shared__ float warpLevelSum[WARP_SIZE];
    const int warpId = threadIdx.x / WARP_SIZE;
    const int laneId = threadIdx.x % WARP_SIZE;

    sum = warpReduceSum<blockSize>(sum);

    if(laneId == 0){
        warpLevelSum[warpId] = sum;
    }
    __syncthreads();

    // read from shared memory only if that warp existed
    if(threadIdx.x < blockDim.x / WARP_SIZE){
        sum = warpLevelSum[laneId];
    }else{
        sum = 0;
    }

    // final reduce using first warp
    if(warpId == 0){
        sum = warpReduceSum<blockSize / WARP_SIZE>(sum);
    }

    // write result for this block to global memory
    if(tid == 0){
        g_odata[blockIdx.x] = sum;
    }
}

int main(){

    // input host
    float* input_host = (float*)malloc(N * sizeof(float));
    for(int i = 0; i < N; ++i){
        input_host[i] = 2.0;
    }

    // input device
    float* input_device;
    cudaMalloc((void**)&input_device, N * sizeof(float));

    // copy data from host to device
    cudaMemcpy(input_device, input_host, N * sizeof(float), cudaMemcpyHostToDevice);
    
    const int block_num = 1024;
    // 每个 block 要处理的数据量
    const int NUM_PER_BLOCK = N / block_num;
    // 每个 thread 要处理的数据量
    const int NUM_PER_THREAD = NUM_PER_BLOCK / BLOCK_SIZE;
    float* output_host = (float*)malloc(block_num * sizeof(float));

    // output device
    float* output_device;
    cudaMalloc((void**)&output_device, block_num * sizeof(float));

    // kernel launch
    dim3 grid(block_num, 1);
    dim3 block(BLOCK_SIZE, 1);
    reduce_v7<BLOCK_SIZE, NUM_PER_THREAD><<<grid, block>>>(input_device, output_device);

    // copy result from device to host
    cudaMemcpy(output_host, output_device, block_num * sizeof(float), cudaMemcpyDeviceToHost);

    // print some result
    for(int i = 0; i < 5; ++i){
        float data = output_host[i];
        printf("output_host[%d] = %.2f\n", i, data);
    }

    return 0;
}

reduce_v7 的优化思路在于利用 shuffle 指令在 warp 内高效完成归约，减少共享内存的访问和同步开销，下面是具体流程的说明：

1. 数据加载与局部求和

每个线程负责加载多个数据元素：每个线程根据 blockIdx 和线程索引计算出起始全局内存地址，然后在一个展开的循环中累加 NUM_PER_THREAD 个数据：

unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * (blockSize * NUM_PER_THREAD) + threadIdx.x;
#pragma unroll
for(int iter = 0; iter < NUM_PER_THREAD; iter++){
    sum += g_idata[i + iter * blockSize];
}

这样，每个线程得到了自己局部的和，保存在变量 sum 中

2. 使用 shuffle 指令进行 warp 内归约

warpReduceSum 函数：该 device 函数利用 __shfl_down_sync 指令在 warp 内归约，具体实现我们后续再分析

3. 将 warp 级归约结果存入共享内存

写入共享内存：定义了一个大小为 WARP_SIZE(32) 的共享内存数组 warpLevelSum，每个 warp 内的第一个线程（即 laneId == 0）将其 warp 内的归约结果写入：

if(laneId == 0){
    warpLevelSum[warpId] = sum;
}

这里，warpId = threadIdx.x / WARP_SIZE 表示线程所属的 warp 编号

同步与读取：调用 __syncthreads() 保证所有 warp 内的部分结果都写入共享内存。接着，只有前 blockDim.x / WARP_SIZE 个线程（对应每个 warp 的归约结果）参与下一步：

if(threadIdx.x < blockDim.x / WARP_SIZE){
    sum = warpLevelSum[laneId];
} else {
    sum = 0;
}

4. 最终的 block 级归约

利用第一个 warp（warp0）进行归约：此时，归约结果已经在共享内存中存储为一个数组，每个元素代表一个 warp 的部分和。只有第一个 warp 中的线程会执行最终归约（if(warpId == 0)），再次调用 warpReduceSum：

if(warpId == 0){
    sum = warpReduceSum<blockSize / WARP_SIZE>(sum);
}

这一步骤将所有 warp 的部分和归并到一起，最终 threadIdx.x = 0 的线程得到整个 block 的归约结果

写回全局内存：最后，将整个 block 的归约结果写入全局内存：

if(tid == 0){
    g_odata[blockIdx.x] = sum;
}

reduce_v7 通过让每个线程预先累加多个数据，再利用 shuffle 指令在 warp 内高效归约、将部分结果写入共享内存后在第一个 warp 内完成最终归约，从而最大程度地减少了同步和内存访问开销，实现了极高的性能。

博主绘制了一张草图来帮助理解：

在这里插入图片描述

reduce_v7 核函数的整体思路还是比较清晰的，以一个 block 为例，reduce_v7 会先对 warp0~warp7 进行归约计算，并将其结果保存在共享内存数组 warpLevelSum 中，然后通过线程同步确保所有 warp 结果均已写入 warpLevelSum 中，接着将共享内存数组中的前 8 个数分别赋值给 warp0 的前 8 个线程的 sum 变量，而 warp0 的其它线程中的 sum 变量会赋值 0，最后对 warp0 的前 8 个线程做一次 warp 归约得到整个 block 的归约结果

下面我们来看下 warp 内归约函数 warpReduceSum 具体是如何做的

我们知道，一个 warp 中有 32 个线程，归约的目标是把这 32 个线程中的值累加到一个线程（通常是 landId = 0）中。__shfl_down_sync 指令允许同一 warp 内的线程直接交换寄存器中的数据，而无需借助共享内存或同步调用

具体来说，__shfl_down_sync(mask, value, delta) 的参数含义如下：

mask：表示参与交换的线程掩码，通常设置为 0xffffffff 表示整个 warp 的 32 个线程都参与
value：当前线程中要传递的数据
delta：表示向下（即向更高的 lane 编号方向）移动多少个线程的数据。例如 delta 为 16 时，每个线程会获得它自己所在 lane 号加 16 的线程中的值

在 warpReduceSum 中采用了依次使用 delta 为 16、8、4、2、1 的调用，形成了一棵二叉树结构的归约流程。举个例子，假设一个 warp 中线程的编号是 0~31，初始时每个线程保存一个数值，例如：

lane:   0  1  2  3  ... 15 | 16 17 18 19 ... 31
value:  a  b  c  d  ... p  | q  r  s  t  ... z

归约过程如下：

1. 第一步（delta = 16）

每个线程执行：

sum += __shfl_down_sync(0xffffffff, sum, 16);

线程 0~15：
线程 0：a = a + q
线程 1：b = b + r
…
线程 15：p = p + z
线程 16~31 的返回值未被用于累加，因为它们不需要再累加其他数据

结果：前 16 个线程保存了 16 个部分和

2. 第二步（delta = 8）

接下来，对于线程 0~7：

线程 0：之前的 a + q 再加上来自线程 8 的值
线程 1：同理，累加线程 9 的值
…

这一步将 16 个数归约为 8 个数

3. 第三步（delta = 4）

对线程 0~3，分别加上各自向下 4 个位置的数，归约到 4 个数

4. 第四步（delta = 2）

对线程 0~1，分别加上自己下面 2 个位置的数，归约到 2 个数

5. 第五步（delta = 1）

最后，线程 0 加上线程 1 的值，得到整个 warp 的累加和

__shfl_down_sync 允许在 warp 内直接从相对位置“下移”数据。通过连续使用 delta 值 16、8、4、2、1，我们形成了一棵二叉树归约结构，确保所有线程中的数据最终都被归约到线程 0 中，从而实现高效的 warp 内归约

那到这里博主其实有个困惑，在第二步（delta = 8）的时候线程 0 会累加线程 8 的值，线程 8 会累加线程 16 的值，那会不会出现这么一种情况，那就是线程 0 读取线程 8 寄存器的值累加的时候，恰好这时候线程 8 又在读取线程 16 寄存器的值进行累加，这个累加值会不会更新到线程 8 的寄存器中，导致线程 0 在取线程 8 寄存器值的时候其实是线程 8 在第二步中已经累加的结果，最终导致归约的结果出错🤔

其实在 CUDA 的 SIMT 模型中，warp 中所有线程是以锁步（lockstep）的方式执行指令的，__shfl_down_sync 指令就是在这种锁步机制下进行数据交换的，具体解释如下：

1. 锁步执行保证同时读取

当 warp 内的所有线程执行第二步的 __shfl_down_sync(sum, 8) 时，所有线程同时进入这一指令。此时，每个线程都会从目标线程（例如 landId = 8）的寄存器中读取值。

硬件保证这个读取操作是在所有线程已经完成第一步的更新后统一开始的，而不是某个线程先更新、另一个线程再读取。这意味着线程 0 读取线程 8 时，拿到的正是线程 8 在第一步结束时的值，而不会看到线程 8 在第二步中“更新一半”的状态

2. __shfl_down_sync 的工作机制

__shfl_down_sync(maks, value, delta) 指令在执行时，所有线程在同一时刻根据各自寄存器中的 value，读取各自 laneId + delta 线程的 value。

这个读取是原子进行的，所有线程都在同一个指令周期内完成读取和加法操作，因此不存在线程 8 读取线程 16 的值“抢先更新”导致线程 0 读到第二步结果的问题

性能和带宽的测试情况如下：

优化手段	耗时(us)	Memory throughout(%)	DRAM throughout(%)	加速比
reduce_baseline	2490us	62.01%	15.44%	~
reduce_v1_interleaved_addressing	1800us	85.78%	21.35%	1.38
reduce_v2_bank_conflict_free	1730us	89.24%	22.01%	1.44
reduce_v3_idle_threads_free	896.19us	89.48%	42.89%	2.78
reduce_v4_unroll_last_warp	507.71us	76.08%	76.08%	4.90
reduce_v5_completely_unroll	485.54us	79.48%	79.48%	5.13
reduce_v6_multi_add	400.06us	96.94%	96.94%	6.22
reduce_v7_shufl_down_sync	400.80us	96.88%	96.88%	6.21

可以看到基于 warp 原语 __shfl_down_sync 进行优化之后，耗时和带宽利用率和 reduce_v6 基本没有什么区别，这可能是 GPU 硬件架构迭代后内部更新的原因

不同优化技巧带来的性能表现如下：

在这里插入图片描述

通过上面一系列的优化技巧，我们对 reduce 进行了不断地优化，可以看到其实在 reduce_v6 的时候就已经很难再提升了，而这个数据跟 NVIDIA 博客中的数据还是有比较大的出入，主要原因可能是因为 GPU 硬件架构已经更新了好几代。总而言之，我们通过一系列的优化已经把 reduce 优化到一个非常好的程度了，还是非常 nice 的🤗

5. 拓展—CUDA工具链的使用

我们平时主要关注两个 CUDA 工具链的使用：NVIDIA Nsight Compute 和 NVIDIA Nsight Systems

NVIDIA Nsight Compute 和 NVIDIA Nsight Systems 都是 NVIDIA 提供的性能分析工具，但它们关注的层次和用途有所不同：

NVIDIA Nsight Compute

用途：专注于单个 CUDA 内核（kernel）的详细性能分析
功能：提供指令级别的计数、内存吞吐量、线程执行、寄存器使用、共享内存利用率等详细指标

NVIDIA Nsight Systems

用途：面向整个系统级别的性能分析，捕捉应用程序从 CPU 到 GPU 的整体行为
功能：提供时间轴视图，展示 CPU 线程、GPU 活动、驱动程序调用以及系统其他组件之间的交互与同步情况

Note：在本文的 reduce 核函数分析中，我们只需要用 Nsight Compute 即可

Nsight Compute 其实在安装完 CUDA 之后就已经安装了，具体安装位置在 /usr/local/cuda/bin，如下图所示：

在这里插入图片描述

Nsight Compute 提供两种界面：

1. ncu（命令行工具）：

用途：通过命令行运行 Nsight Compute，直接在终端中对 GPU 内核进行采样和分析
特点：适合自动化脚本、批量分析和生成报告文件（如 .ncu-rep），便于嵌入到开发和 CI/CD 流程中
使用方式：通过传递参数（例如 --launch-count）运行可执行文件，并将分析结果输出到终端或指定的报告文件中

2. ncu-ui（图形用户界面工具）

用途：提供直观的图形界面，用于加载和浏览由 ncu 生成的报告文件
特点：界面中会显示各种图表和 Breakdown 表格，如 GPU Throughput 图、Compute 和 Memory 的各项细分指标，便于交互地查看和分析性能数据
使用方式：启动 ncu-ui 后，通过菜单加载 .ncu-rep 报告文件，或直接配置目标程序进行交互式分析

关于 Nisight Compute 的使用大家可以参考：【模型分析】Nsight Compute使用入门，大家如果感兴趣的话也可以看看 UP 的视频：【CUDA调优指南】合并访存

另外一个工具 Nsight System 博主目前也暂未使用过，只是进行了安装，Nsight System 安装比较简单，大家可以参考：【cuda】Nsight System 下载，安装与使用，软件包外网访问下载可能比较慢，大家可以点击 here 下载（注意下载的 Nsight System 是 2025.1.1 版本，如果有其它版本需求可以自行下载）

关于服务器上的 Nsight System 安装大家可以参考：Nsight System的安装和使用

由于博主对这两个工具的使用也不熟练，自己也还在摸索中，所以这里只教大家如何利用 Nsight Compute 简单的查看 kernel 的耗时和带宽利用率😄

我们以 redece_v7 核函数为例，首先我们需要编写代码，新建一个 .cu 文件 reduce_v7_shfl_down_sync.cu，其内容如下：

#include <cuda.h>
#include <cuda_runtime.h>
#include <time.h>
#include <stdio.h>

#define N 32*1024*1024
#define BLOCK_SIZE 256
#define WARP_SIZE 32

template <unsigned int blockSize>
__device__ __forceinline__ float warpReduceSum(float sum){
    if(blockSize >= 32) sum += __shfl_down_sync(0xffffffff, sum, 16);
    if(blockSize >= 16) sum += __shfl_down_sync(0xffffffff, sum, 8);
    if(blockSize >= 8) sum += __shfl_down_sync(0xffffffff, sum, 4);
    if(blockSize >= 4) sum += __shfl_down_sync(0xffffffff, sum, 2);
    if(blockSize >= 2) sum += __shfl_down_sync(0xffffffff, sum, 1);
    return sum;
}

template <unsigned int blockSize, int NUM_PER_THREAD>
__global__ void reduce_v7(float* g_idata, float* g_odata){
    float sum = 0;

    // each thread loads NUM_PER_THREAD element from global
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockSize * NUM_PER_THREAD) + threadIdx.x;
    #pragma unroll
    for(int iter = 0; iter < NUM_PER_THREAD; iter++){
        sum += g_idata[i + iter * blockSize];
    }
    
    // shared memory for partial sums (one per warp in the block)
    static __shared__ float warpLevelSum[WARP_SIZE];
    const int warpId = threadIdx.x / WARP_SIZE;
    const int laneId = threadIdx.x % WARP_SIZE;

    sum = warpReduceSum<blockSize>(sum);

    if(laneId == 0){
        warpLevelSum[warpId] = sum;
    }
    __syncthreads();

    // read from shared memory only if that warp existed
    if(threadIdx.x < blockDim.x / WARP_SIZE){
        sum = warpLevelSum[laneId];
    }else{
        sum = 0;
    }

    // final reduce using first warp
    if(warpId == 0){
        sum = warpReduceSum<blockSize / WARP_SIZE>(sum);
    }

    // write result for this block to global memory
    if(tid == 0){
        g_odata[blockIdx.x] = sum;
    }
}

int main(){

    // input host
    float* input_host = (float*)malloc(N * sizeof(float));
    for(int i = 0; i < N; ++i){
        input_host[i] = 2.0;
    }

    // input device
    float* input_device;
    cudaMalloc((void**)&input_device, N * sizeof(float));

    // copy data from host to device
    cudaMemcpy(input_device, input_host, N * sizeof(float), cudaMemcpyHostToDevice);
    
    const int block_num = 1024;
    // 每个 block 要处理的数据量
    const int NUM_PER_BLOCK = N / block_num;
    // 每个 thread 要处理的数据量
    const int NUM_PER_THREAD = NUM_PER_BLOCK / BLOCK_SIZE;
    float* output_host = (float*)malloc(block_num * sizeof(float));

    // output device
    float* output_device;
    cudaMalloc((void**)&output_device, block_num * sizeof(float));

    // kernel launch
    dim3 grid(block_num, 1);
    dim3 block(BLOCK_SIZE, 1);
    reduce_v7<BLOCK_SIZE, NUM_PER_THREAD><<<grid, block>>>(input_device, output_device);

    // copy result from device to host
    cudaMemcpy(output_host, output_device, block_num * sizeof(float), cudaMemcpyDeviceToHost);

    // print some result
    for(int i = 0; i < 5; ++i){
        float data = output_host[i];
        printf("output_host[%d] = %.2f\n", i, data);
    }

    return 0;
}

接着我们利用 nvcc 工具生成一个可执行文件，指令如下：

mkdir bin
/usr/local/cuda/bin/nvcc -o bin/reduce_v7 reduce_v7_shfl_down_sync.cu

执行成功后我们会在 bin 文件夹下看到一个可执行文件 reduce_v7

接着我们可以利用 ncu 命令行工具分析核函数

sudo /usr/local/cuda/bin/ncu --launch-count 2000 ./bin/reduce_v7

Note：--launch-conut 参数会告诉 Nsight Compute 对目标程序中检测到的 kernel 调用重复采样 N 次，在这里是 2000 次

输出如下图所示：

在这里插入图片描述

可以看到执行指令后在终端输出了一系列的信息，我们依次来看下：

==PROF== Connected to process 805261 (/home/jarvis/Downloads/reduce/bin/reduce_v7)
==PROF== Profiling "reduce_v7" - 1 of 2000: 0%....50%....100% - 9 passes
==PROF== Disconnected from process 805261
[805261] reduce_v7@127.0.0.1
  void reduce_v7<(unsigned int)256, (int)128>(float *, float *), 2025-Mar-08 13:50:55, Context 1, Stream 7

首先显示的这部分基本信息记录了当前采样 kernel 的一些基本上下文信息，包括 kernel 名称、模板参数、执行时间戳、CUDA Contex、CUDA Stream 等等，这些信息可以帮助我们确认正在分析的是哪一个 kernel 以及运行时的上下文环境

接下来报告中主要有三个 Section，每个 Section 的作用如下：

1. GPU Speed Of Light Throughput

功能：

这部分反映了 kernel 在执行时各个硬件单元（例如 DRAM、SM、缓存）的吞吐能力和利用情况，基本上说明了内核的“极限速度情况”，因此对本文来说我们主要关注这个 Section 即可

主要指标：

Duration：表示 kernel 执行的时间（如 400.42 us），直接反映 kernel 的运行耗时
Memory [%] 和 DRAM Throughout：显示 kernel 在访问全局内存（DRAM）时所达到的利用率，当数值接近 100% 时，表明内存带宽已经接近瓶颈
L1/TEX Cache Throughput 和 L2 Cache Throughout：表示各级缓存的访问情况
Compute (SM) [%]：表示 SM 中的计算单元的利用率

作用：

通过这一节的数据，我们可以判断 kernel 是受内存带宽还是计算能力限制，进而为优化方向提供指令

2. Launch Statistics

功能：

这一部分展示了 kernel 的启动配置和资源分配情况，主要涉及 kernel 是如何在 GPU 上分布的，以及使用了多少硬件资源

主要指标：

Grid Size 和 Block Size：显示了启动时设置的 grid 和 block 的维度
Register Per Thread：每个线程使用的寄存器数，反映了寄存器资源的消耗
Shared Memory Configuration：包括静态和动态共享内存的分配情况
Threads 和 Waves Per SM：表示总线程数以及每个 SM 中的 warp 数量

作用：

这些信息帮助我们了解 kernel 的配置是否合理，是否存在资源使用过多或不足的问题，从而影响性能。

3. Occupancy

功能：

Occupancy 部分描述了 kernel 在 GPU 上的实际活跃程度，即 GPU 多大程度上被我们的 kernel 利用起来了

主要指标：

Theroretical Active Warps per SM：理论上每个 SM 能同时激活的 warp 数量
Achieved Occupancy：实际达到的占用率，通常以百分比表示
Achieved Active Warps Per SM：每个 SM 实际激活的 warp 数量
Block/Registers/Shared Memory Limit：表示由于硬件资源限制（如寄存器或共享内存）而可能限制 kernel 并行度的部分

作用：

Occupancy 数据可以帮助我们判断 kernel 的并行度是否达到预期。如果占用率低，可能意味着存在资源瓶颈（例如寄存器或共享内存分配不合理），这时可能需要优化 kernel 的资源使用或调整启动参数

这三个部分合起来为我们提供了一个从硬件资源利用、启动配置到实际执行效率的全面视角，从而帮助我们找出性能优化的重点。

在这些信息中我们关注的 reduce_v7 核函数的耗时和带宽利用率分别是 400.42 us 和 96.98 %

除了使用 ncu 直接在终端查看 kernel 性能分析外，我们还可以将其导出保存为 .ncu-rep 的报告，接着在 nuc-ui 图形化界面中查看，更加直观

具体流程如下：

首先我们需要利用 ncu 工具生成一份 .ncu-rep 报告，指令如下：

sudo /usr/local/cuda/bin/ncu --launch-count 2000 -o reduce_v7 ./bin/reduce_v7

执行后输出如下：

在这里插入图片描述

可以看到性能分析报告保存到了 reduce_v7.ncu-rep 中

接着我们启动 ncu-ui，在终端直接输入指令：

ncu-ui

启动完成后我们要将刚才导出的报告加载到界面中，依次点击 File->Open Files，接着选择我们刚才导出的 reduce_v7.ncu-rep 文件，最后点击 Open，接着我们就能在 UI 界面中看到相关的性能分析信息了，如下图所示：

在这里插入图片描述

可以看到 UI 界面展示的内容和终端输出的没什么差异，只是以界面的方式展示更加的直观

OK，以上就是使用 Nsight Compute 简单查看 kernel 耗时和带宽分析的内容了

结语

这篇文章我们继续学习了 reduce 的一些优化，包括展开最后一个 warp 减少同步带来的浪费、完全展开循环、调节 grid 和 block 参数以及使用 shuffle 指令，此外博主还简单介绍了 Nsight Compute 的使用。

OK，以上就是整篇文章的全部内容了

总的来说，整个 reduce 的优化过程还是挺有意思的，以解决问题的角度出发去优化 reduce sum，相比于直接枯燥的学习语法要有意思得多，大家感兴趣的可以看看 up 主的视频，还是很有收获的🤗