规约reduce算法

内存类型	速度	容量	可见范围
全局内存	~1 TB/s (HBM)	大 (GB级)	所有线程
共享内存	~200 TB/s (SRAM)	小	同一个 Block
寄存器	最快 (0延迟)	极少	单个线程

• 规约算法的核心思路：把数据从全局内存搬到共享内存和寄存器，在最快的层级上完成计算并写回

SIMT 和 Warp Divergence

GPU使用SIMT，对应的一个warp的3个线程在同一时刻执行同一条指令
如果写了分支代码，warp内有些线程走if，有些走else/此时硬件只能先执行if，再执行else分支。实际上代码背串行化了。
因此再规约算法中要小心处理 if (threadIdx.x < stride) 这类条件，一旦stride小于32，分支会在同一warp内产生divergence

共享内存和访存冲突

共享内存被划分为若干Bank，通常为32个。每个bank在同一时钟周期内只能服务于一个线程。如果同一个Warp中的多个线程访问同一个Bank中的不同地址，就会发生Bank conflict，导致访问串行化
- 无冲突：所有线程访问不同bank，或访问同一bank的同一地址
- 2路冲突：两个线程访问同一bank的不同地址，访存时间翻倍
- 最坏情况：32个线程挤在同一bank，串行化32倍
在规约算法中内存访问模式的微小差异会对性能造成巨大影响

规约是什么？

规约就是“把一堆东西合并成一个东西”，形式化的说，给定集合$A = {a_1,a_2,…,a_n}$，和一个二元操作⊕，规约计算的是：
$$
/text{result} = a_1 ⊕ a_2 ⊕ … ⊕ a_n
$$
规约的特点是：输入高维，输出低维

特性	逐元素操作	规约操作
输出形状	与输入相同	维度降低
数据依赖	无	前后依赖
并行策略	直接映射	需要多轮合并
典型例子	add sigmoid gelu	sum max mean softmax

常见的规约算子有：

算子	操作	可能出现的位置
Sum	a + b	损失函数求和、 BatchNorm均值
Max / Min	max(a,b)	MaxPooling、Softmax的稳定化
Argmax	比较值，保留索引	分类模型最后的预测
Mean	sum再除以N	LayerNorm/ BatchNorm
Variance	先sum再求方差	LayerNorm/BatchNorm
Dot Product	逐元素乘再sum	矩阵乘法的核心
其中Sum是所有规约算子的基础，满足结合律和交换律，单位元是0，硬件实现高效，存在于几乎所有深度学习框架

算子优化(向量求和为例)

def sequential_sum(arr):
    result = 0
    for x in arr:
        result += x
    return result

以上为最基础的求和方法，是线性链式结构，每一步依赖前一步的结果。为了让算法更快，首先改变计算的拓扑结构

树形规约

既然串行是链式的，我们将其修改为属性，基本思想是凉凉相加，具体来说：
每一轮相距步长的元素相加，步长逐轮翻倍

__global__ void reduce_vo(const float* in, float* out, int N) {
    __shared__ float sdata[1024];

    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + tid;

    sdata[tid] = (idx < N) ? in[idx] : 0.0f;
    __syncthreads();

    // 步长翻倍
    for(int s = 1; s < blockDim.x; s *= 2) {
        if (tid % (2 * s) == 0) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    // 每个block的结果用原子加法写回全局内存
    if (tid == 0) atomicAdd(out, sdata[0]);
}

shared

为什么这里使用共享内存呢？因为相比全局内存拿取写回，共享内存会给同一个block里面的线程提供一个空间，物理位置位于SM内部，延迟极低，带宽极高。一个Block内的所有线程几乎都可以零成本的读写这一块共享空间，这让数据交换快了很多
声明 __shard__ float sdata[1024],意味着每个block都会有一份sdata，大小为1024，如果我的block启动1024个线程，每个线程有一个专属的sdata[tid]的各自来存放从全局内存搬运过来的数。

int tid = threadIdx.x;
int idx = blockIdx.x * blockDim.x + tid;
sdata[tid] = (idx < N) ? in[idx] : 0.0f;
__syncthreads();

这一步表示每个线程计算自己的idx，然后把自己负责的元素从in数组拿进来，放进共享内存对应的小格子里面。
如果数组长度不是block大小整数倍，最后几个线程会越界，通过(idx < N)?in[idx]:0.0f来保护
最终通过 _syncthreads()确保所有线程把自己的数据放妥当，等到所有线程就绪后进行下一步

循环部分

现在共享内存摆好每个线程带进来的数，接下来的任务是在block内部把这些数相加，浓缩为一个和。
最基础的算法是树状折半规约，这里采用循环实现

for (int s = 1;s<blockDim.x;s*=2){
    if (tid %(2*s) == 0){
        sdata[tid] += sdata[tid+s];
    }
    __syncthreads();
}

最终输出

if (tid == 0) atomicAdd(out, sdata[0]);

最终通过原子加法将局部和累加到全局输出out中
因为有很多block同时做这件事，使用atomicAdd保证不会发生读写冲突

继续优化！

从逻辑上，规约算法已经完整了，但是从性能上，实际还有很大差距：

• 取模运算
  在GPU上取模非常慢，2*s在编译器中没有办法高效优化为移位和加法，绝大部分在排队等待除法结果，计算资源被浪费
• Bank conflict
  共享内存虽然快，但是容易发生访存冲突，它被分为32个bank，回到我们的循环：第一轮 s=1，干活的是偶数线程 0,2,4,…，它们访问的地址是 (0,1), (2,3), (4,5)，相邻干活线程之间隔了一个线程，但访问的 bank 还勉强算分散。 可随着 s 增大，干活线程越来越稀疏，到 s=16 时，参与运算的线程 0 和线程 32 访问的地址分别落在 bank 0 和 bank 32——注意 bank 数是 32，bank 32 在硬件上就是 bank 0 的另一行，于是 bank 0 被同时访问，冲突发生。 再往后，s=32 时干活线程的间隔更大，几乎所有访问都挤压在极少数 bank 上，冲突严重到让共享内存的带宽优势荡然无存。

优化取模

__global__ void reduce_v1(const float* in, float* out, int N) {
    __shared__ float sdata[1024];
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    sdata[tid] = (idx < N) ? in[idx] : 0.0f;
    __syncthreads();

    for(int s=1; s < blockDim.x; s *= 2) {
        int index = 2 * s * tid;
        if (index < blockDim.x) {
            sdata[index] += sdata[index + s];
        }
        __syncthreads();
    }

    if (tid == 0) atomicAdd(out, sdata[0]);
}

使用乘加代替取模

让每个线程直接算出自己要服务的下标index。这个计算包含一次整数乘法和一次加法，整数乘法在GPU上的吞吐量远高于除法和取模。编译器看到这种固定模式的乘法和移位，往往能够进一步优化指令组合

优化bank conflict

__global__ void reduce_v2(const float* in, float* out, int N) {
    __shared__ float smem[1024];
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + tid;

    smem[tid] = (idx < N) ? in[idx] : 0.0f;
    __syncthreads();

    for (int s = BLOCK >> 1; s > 0; s >>= 1) {
        if (tid < s) smem[tid] += smem[tid + s];
        __syncthreads();
    }

    if (tid == 0) atomicAdd(out, smem[0]);
}

循环方式被彻底修改了：

for (int s = BLOCK >> 1; s > 0; s >>= 1) {
    if (tid < s) smem[tid] += smem[tid + s];
    __syncthreads();
}

我们的步长不再是从小到大，原来是1，2，4，6，8。 而我们现在的计算方式是 s = BLOCKK / 2 比如 BLOCK = 1024，那么s一开始就是512.
-第一轮 0-511干活，剩下的闲着
-第二轮 0-255干活…
-最后剩下线程0
并且工作的线程是连续的，但是每个线程访存的地址发生了变化是连续的smem[tid]和另一段连续的smem[tid+s]。
整个规约过程中，每一轮都减半，确保bank冲突从头到尾都是零

优化！优化！

__global__ void reduce_v3(const float* in, float* out, int N) {
    __shared__ float smem[1024];
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockDim.x * 2) + threadIdx.x;

    float val1 = (idx < N) ? in[idx] : 0.0f;
    float val2 = (idx + blockDim.x < N) ? in[idx + blockDim.x] : 0.0f;
    smem[tid] = val1 + val2;
    __syncthreads();

    for (int s = BLOCK >> 1; s > 0; s >>= 1) {
        if (tid < s) smem[tid] += smem[tid + s];
        __syncthreads();
    }

    if (tid == 0) atomicAdd(out, smem[0]);
}

针对数据加载模块优化，从全局内存一次搬运一个元素优化为搬运多个元素，例如两个，并直接在本地寄存器加起来

int idx = blockIdx.x * (blockDim.x * 2) + threadIdx.x;

带来了两种收益：
1. 减少了atomicAdd的竞争
2. 提升计算密度。寄存器的带宽元朝共享内存，并且没有bank conflict，线程在等待全局内存时就可以做加法

• 为什么不一次读4个，8个？
  当然可以，而且许多高性能库确实是这么做的。但是增加单线程负载会降低活跃线程块数，降低硬件占用率，因此在”每线程工作量“和”并行度“之间找到一个平衡点

展开最后一个Warp

优化完数据加载，再来优化掉最后一次加载。
如果最后一次加载不满足32个线程，直接打表计算（将所有倾向列出来）

__device__ void warpReduce_v4(volatile float* smem, int tid) {
    smem[tid] += smem[tid + 32];
    smem[tid] += smem[tid + 16];
    smem[tid] += smem[tid + 8];
    smem[tid] += smem[tid + 4];
    smem[tid] += smem[tid + 2];
    smem[tid] += smem[tid + 1];
}

__global__ void reduce_v4(const float* in, float* out, int N) {
    __shared__ float smem[1024];
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockDim.x * 2) + threadIdx.x;

    float val1 = (idx < N) ? in[idx] : 0.0f;
    float val2 = (idx + blockDim.x < N) ? in[idx + blockDim.x] : 0.0f;
    smem[tid] = val1 + val2;
    __syncthreads();

    // 规约到只剩 64 个数时停手
    for (int s = BLOCK >> 1; s > 32; s >>= 1) {
        if (tid < s) smem[tid] += smem[tid + s];
        __syncthreads();
    }

    // 最后 64 → 1 的规约，完全在 warp 内完成
    if (tid < 32) warpReduce_v4(smem, tid);

    if (tid == 0) atomicAdd(out, smem[0]);
}

提前停手

当规约到s=64时直接停止循环，此时共享内存的钱64个严肃？包含了64个和，生效的可以手动实现

if (tid < 32) warpReduce_v4(smem, tid);

这里只让前32个线程进入函数，其他warp此时已经解放，不需要继续做同步

__device__ void warpReduce_v4(volatile float* smem, int tid) {
    smem[tid] += smem[tid + 32];   // 64 → 32
    smem[tid] += smem[tid + 16];   // 32 → 16
    smem[tid] += smem[tid + 8];    // 16 → 8
    smem[tid] += smem[tid + 4];    // 8  → 4
    smem[tid] += smem[tid + 2];    // 4  → 2
    smem[tid] += smem[tid + 1];    // 2  → 1
}

每一步步幅和之前保持一致，但是去掉了同步

• 为什么可以这么干？
  Warp内的32条线程在硬件上是锁步执行的，它们执行的是同一条命令。只要32个线程一起进入函数，它们内部的每一步加法对外都是天然同步的，不存在写后读的风险，也不需要显式同步

Volatile的作用

__device__ void warpReduce_v4(volatile float* smem, int tid)

注意到smem的参数类型是volatile float*？
- 硬件保证不了编译器会老老实实地每次把读写都落到共享内存上。编译器会做非常激进地优化：把频繁访问的共享内存缓存在寄存器里。

smem[tid] += smem[tid + 32];
smem[tid] += smem[tid + 16];

比如以上这段代码，smem[tid]在第一句之后已经被加载到寄存器了，第二局再次使用，可能直接使用寄存器中的旧值，而不是回到共享内存。但是smem[tid]在第一句中被写回了新值，而且读smem[tid+16]的线程可能在同一个warp内的另一个位置，他需要看到的是写进共享内存的最新结果。
一旦编译器smem[tid]缓存到寄存器，就会发生：
- 线程0执行完第一句，smem[0]新值被写回共享内存
- 线程0执行第二句时，编译器可能直接使用寄存器中smem[0]值加上smem[16]，共享内存的最新smem[0]被无视
- 最终smem[0]得到一个错误的部分和
volatile 告诉编译器每次用到这个变量时，必须从共享内存里老老实实读写。加上关键字之后，编译器不再对smem指向的地址做任何寄存器和缓存优化。每一步累加都从共享内存取最新的写入，保证计算完全正确。

• 总结：
  在CUDA编程中，volatile核心作用时强制编译器每次从显存中读取变量，以避免因为优化导致的数据不一致问题
• 注意：
  volatile并不等于同步或者原子操作。如果多个线程并发写，那么仍然会产生数据竞争。并不能保证原子性和顺序

试试模板化

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

template <unsigned int blockSize>
__global__ void reduce_v5_impl(const float* in, float* out, int N) {
    __shared__ float smem[1024];
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockSize * 2) + tid;

    float val1 = (idx < N) ? in[idx] : 0.0f;
    float val2 = (idx + blockSize < N) ? in[idx + blockSize] : 0.0f;
    smem[tid] = val1 + val2;
    __syncthreads();

    if (blockSize >= 512) { if (tid < 256) smem[tid] += smem[tid + 256]; __syncthreads(); }
    if (blockSize >= 256) { if (tid < 128) smem[tid] += smem[tid + 128]; __syncthreads(); }
    if (blockSize >= 128) { if (tid < 64) smem[tid] += smem[tid + 64]; __syncthreads(); }

    if (tid < 32) warpReduce_v5<blockSize>(smem, tid);

    if (tid == 0) atomicAdd(out, smem[0]);
}

• 模板+编译期常量展开
  之前的版本都是通过for循环实现的，每次迭代都有读取，判断，右移，跳转等额外开销。而这个版本的想法是：既然BLOCKSIZE在编译时就是常数，干脆将所有可能的不正常直接展开成为条件判断

  if (blockSize >= 512) { ... }
  if (blockSize >= 256) { ... }
  if (blockSize >= 128) { ... }

  比如调用 reduce_v5_impl<256>时，最终只有等价的两行：

  if (tid < 128) smem[tid] += smem[tid + 128]; __syncthreads();
  if (tid < 64)  smem[tid] += smem[tid + 64];  __syncthreads();

  此时循环就彻底消失了，没有循环变量、条件跳转、没有任何控制指令，只剩下加法运行算。
  模板把可配置性与性能这两个看似矛盾的事情联系在了一起，在为任何blockSize生成专属的优化版本，但是对于不同blockSize，我们同样可以进行模板化

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

  多了一个if判断，但是整个后边的逻辑自动适配blockSize，不再需要人工调整
• 为何使用模板：
  如果blockSize时普通函数参数，在运行时仍然要进行一次判断；但是作为模板参数，它属于编译期常数，编译器会为每一个不同的blockSize生成专属内核函数版本

Grid Loop优化

虽然每个线程都会加载元素并做加法，但是实际上越到规约后期。大量的线程是在干等。如果数组过大，需要启动机器多线程块，此时线程启动的开销无法被忽视，并且原子加法也会因为线程太多而产生激烈的竞争

Grid Loop是做什么的?

• 与其让大量的线程只算一次，不如只开几十个线程块，让每个线程不停计算

  template <unsigned int blockSize>
  __device__ void warpReduce_v6(volatile float* smem, int tid) {
      if (blockSize >= 64) smem[tid] += smem[tid + 32];
      if (blockSize >= 32) smem[tid] += smem[tid + 16];
      if (blockSize >= 16) smem[tid] += smem[tid + 8];
      if (blockSize >= 8) smem[tid] += smem[tid + 4];
      if (blockSize >= 4) smem[tid] += smem[tid + 2];
      if (blockSize >= 2) smem[tid] += smem[tid + 1];
  }
  
  
  template <unsigned int blockSize>
  __global__ void reduce_v6_impl(const float* in, float* out, int N) {
      __shared__ float smem[blockSize];
      int tid = threadIdx.x;
      int idx = blockIdx.x * (blockSize * 2) + tid;
      int gridSize = blockSize * 2 * gridDim.x;
  
      float sum = 0.0f;
  
      while (idx < N) {
          sum += in[idx];
          if (idx + blockSize < N) sum += in[idx + blockSize];
          idx += gridSize;
      }
  
      smem[tid] = sum;
      __syncthreads();
  
      // 树形规约（模板展开）
      if (blockSize >= 512) { if (tid < 256) smem[tid] += smem[tid + 256]; __syncthreads(); }
      if (blockSize >= 256) { if (tid < 128) smem[tid] += smem[tid + 128]; __syncthreads(); }
      if (blockSize >= 128) { if (tid < 64)  smem[tid] += smem[tid + 64];  __syncthreads(); }
  
      if (tid < 32) warpReduce_v6<blockSize>(smem, tid);
  
      if (tid == 0) atomicAdd(out, smem[0]);
  }

  在这个版本的循环中，sum是一个存储在寄存器中的变量，原本每一次加法都要放在共享内存中，虽然它也很快，但是和寄存器差一个数量级，在这个版本中只有最后一步规约才启动共享内存。
  在这个版本中，不需要根据N的大小计算开多少Block，而是可以根据硬件SM的数量来固定Grid大小。
  此外我们还可以把向量化也加进来一起试一下：

  template <unsigned int blockSize, unsigned int elementsPerThread = 8>
  __global__ void reduce_v7(const float* __restrict__ in, float* __restrict__ out, int N) {
      __shared__ float smem[blockSize];
      int tid = threadIdx.x;
  
      const float4* in4 = reinterpret_cast<const float4*>(in);
      int n_vec = N / 4;
  
      int idx = blockIdx.x * blockSize + tid;
      int grid_stride = gridDim.x * blockSize;
  
      float sum = 0.0f;
  
      for (int i = idx; i < n_vec; i += grid_stride) {
          float4 val = in4[i];
          sum += val.x + val.y + val.z + val.w;
      }
  
      // 处理剩余的末尾元素 (如果 N 不是 4 的倍数)
      if (idx == 0) {
          for (int i = n_vec * 4; i < N; i++) atomicAdd(out, in[i]);
      }
  
      smem[tid] = sum;
      __syncthreads();
  
      if (blockSize >= 512) { if (tid < 256) smem[tid] += smem[tid + 256]; __syncthreads(); }
      if (blockSize >= 256) { if (tid < 128) smem[tid] += smem[tid + 128]; __syncthreads(); }
      if (blockSize >= 128) { if (tid < 64)  smem[tid] += smem[tid + 64];  __syncthreads(); }
      if (tid < 32) warpReduce_v5<blockSize>(smem, tid);
      if (tid == 0) atomicAdd(out, smem[0]);
  }

  __resctrict__表示函数的两个指针指向的内存区域不会重叠，此时编译器就会进行相应的优化，例如寄存器服用，指令重排等。

配置Grid和Block

在前三个版本中，每个线程都对应一个数据，计算公式为：
Gridsize = [N / Blocksize]

(N + BLOCK - 1) / BLOCK

在中间三个版本中，一个线程对应两个元素，只是因为我们对加载做了优化，每次把相邻的一个元素读进来

(N + 2 * BLOCK - 1) / (2 * BLOCK)

最终两个版本的Grid大小并不取决于数据总量，而是GPU硬件规模，将数据与算子解耦