从算法骨架开始，而不是从线程开始

大多数 GPU 算法可由少数模式组合：

先识别骨架，通常能直接复用 CUB/Thrust/cuBLAS 等成熟实现。自己写 kernel 的主要理由应是融合、特殊布局、特殊数值语义或经过测量的性能缺口。

Map：带宽基线

__global__ void relu(const float* x, float* y, int n) {
  for (int i = blockIdx.x * blockDim.x + threadIdx.x;
       i < n; i += blockDim.x * gridDim.x) {
    y[i] = fmaxf(x[i], 0.0f);
  }
}

Map 通常算术强度低，适合作为内存带宽基线。优化方向：连续访问、减少中间数组、把多次 map 融成一次。但融合也会增加寄存器和编译复杂度，不能无限进行。

Reduce：从串行到树形

求和满足结合律的数学形式：

$$((((x_0+x_1)+x_2)+x_3)+\cdots)$$

可改成深度$O(\log N)$ 的树：

x0 x1 x2 x3 x4 x5 x6 x7
 \ /   \ /   \ /   \ /
 s0    s1    s2    s3
   \  /        \  /
    t0          t1
       \       /
          sum

浮点加法并不严格结合，所以不同树形会得到细微不同结果。这是数值语义变化，不一定是 bug。

Block reduce：shuffle + shared

以下版本要求 blockDim.x 是 32 的整数倍且不超过 1024：

__device__ float warp_sum(float v) {
  constexpr unsigned full = 0xffffffffu;
  for (int offset = warpSize / 2; offset > 0; offset >>= 1)
    v += __shfl_down_sync(full, v, offset);
  return v;
}

__global__ void reduce_sum(const float* x, float* partial, int n) {
  __shared__ float warp_sums[32];
  int tid = threadIdx.x;
  int lane = tid & (warpSize - 1);
  int warp = tid / warpSize;

  float sum = 0.0f;
  for (int i = blockIdx.x * blockDim.x * 2 + tid;
       i < n; i += gridDim.x * blockDim.x * 2) {
    sum += x[i];
    if (i + blockDim.x < n) sum += x[i + blockDim.x];
  }

  sum = warp_sum(sum);
  if (lane == 0) warp_sums[warp] = sum;
  __syncthreads();

  int warp_count = (blockDim.x + warpSize - 1) / warpSize;
  float block_sum = (tid < warp_count) ? warp_sums[lane] : 0.0f;
  if (warp == 0) {
    block_sum = warp_sum(block_sum);
    if (lane == 0) partial[blockIdx.x] = block_sum;
  }
}

设计层次：

每线程：grid-stride 累加多个元素到 register
每 warp：shuffle 在 register 间归约
每 block：每 warp 写一个 shared partial，首 warp 再归约
全 grid：第二次 kernel 归约 partial，或调用 CUB

这比“每个输入直接 atomicAdd(global_sum)”减少几个数量级的全局原子争用。

生产代码优先评估 CUB：cub::DeviceReduce、cub::BlockReduce、cub::WarpReduce。CUB 的层次恰好对应 device/block/warp，且能随架构选择策略。

Scan：并行算法的瑞士军刀

Exclusive scan：

input:          [3, 1, 7, 0, 4, 1, 6, 3]
exclusive sum:  [0, 3, 4,11,11,15,16,22]
inclusive sum:  [3, 4,11,11,15,16,22,25]

Hillis–Steele 每轮向前读距离$2^d$ 的元素，深度$O(\log N)$，但总工作$O(N\log N)$。Blelloch scan 使用 upsweep 构建归约树、把根置单位元，再 downsweep，工作量$O(N)$：

upsweep:   叶子 → 分段和 → 总和
downsweep: 总和位置置 0，父前缀分发到左右子树

完整 device scan 通常分三层：

1. 每 block 扫描自己的 tile，并输出 block total；
2. 扫描所有 block totals；
3. 把对应 block prefix 加回每个 tile。

典型派生：

• Stream compaction：predicate → 0/1 flags → exclusive scan → 有效元素 scatter；
• Radix sort：按某一位分桶，scan 得到每桶写入位置；
• CSR 构造：每行非零数 scan 得到 row offsets；
• 队列分配：每线程请求量 scan 得到无冲突区间。

除非是学习或要与其他操作融合，直接使用 cub::DeviceScan 或 thrust::inclusive_scan/exclusive_scan。

Histogram：解决热点，而不只是使用 atomic

最朴素：

atomicAdd(&bins[input[i]], 1);

当 key 分布偏斜，许多线程争用同一地址。优化层次：

global histogram
    ▲ 每 block 合并一次
block-private histogram（shared）
    ▲ 每 warp 合并一次
warp-private / lane-private partials
    ▲
input items

选择取决于 bin 数、数据分布和 shared memory 容量：

• bin 少：每 block/warp 私有副本可行；
• bin 多：完整副本太贵，可按 tile 处理 key range；
• key 高度偏斜：warp 内先用 ballot/match 聚合同 key，再由 leader atomic；
• 稀疏 key：排序 + reduce-by-key 可能更适合。

优化 histogram 必须用多种分布测试：均匀、单热点、Zipf、全相同。只测均匀随机会掩盖争用。

Stream compaction

目标：从输入中保留满足 predicate 的元素，且保持顺序。

input:   [8, 1, 6, 2, 9]
keep:    [1, 0, 1, 0, 1]
scan:    [0, 1, 1, 2, 2]
scatter: out[0]=8, out[1]=6, out[2]=9

代码骨架：

flag[i] = pred(input[i]);
exclusive_scan(flag, position);
if (flag[i]) output[position[i]] = input[i];

若输出顺序无要求，可用 block-local queues + global atomic reservation，减少全局 scan；但这改变了语义。先写清“稳定/不稳定”的合同。

Transpose：重排访问，而不是增加计算

朴素转置无法同时让 global load 与 store 都合并：

out[x * height + y] = in[y * width + x];

Shared tile 将问题拆成两次合并访问，中间在片上转置：

Global input      Shared tile       Global output
连续读一行  ───► 写 [row][col]
                  读 [col][row] ───► 连续写一行
                  第二维 +1 padding 避免 bank conflict

它展示了 shared memory 的本质：不是“缓存一切”，而是在片上重排，让片外两端都高效。

Stencil：Halo 与时间分块

一维 3-point stencil：

$$y_i = a x_{i-1} + b x_i + c x_{i+1}$$

每 block 把中心区域加左右 halo 放入 shared：

global: ... [halo L | block interior ........ | halo R] ...
shared:      [   L   | reusable interior       |   R   ]

注意：

• 边界条件（zero/clamp/wrap）必须明确；
• halo 加载要保持尽量合并；
• 邻域小且缓存命中好时，shared 版本未必更快；
• 多时间步 temporal blocking 增加复用，但扩大 halo、shared/register 压力和边界复杂度。

Gather 与 Scatter

out[i] = values[index[i]];       // gather：读地址不规则，写连续
out[index[i]] = values[i];       // scatter：读连续，写地址不规则

Gather 的主要问题是 cache locality 和内存事务；scatter 还可能有写冲突。策略：

• 对索引排序/分桶，改善局部性；
• 压缩索引宽度；
• 将重复 index 先局部归约；
• scatter 冲突时明确使用 atomic 或唯一所有权；
• 若重排成本超过收益，不要重排。

Sort、Top-k 与图算法

图和稀疏问题常受不规则访存与负载不均限制，峰值 FLOP/s 不是有用目标。应观察 edge throughput、memory transactions、L2 hit、atomic contention 和尾部 block。

算法设计的五个量

1. Work：总操作数是否接近最优复杂度？
2. Depth：最长依赖链多长？能并行到什么程度？
3. Traffic：各层级搬运多少字节？
4. Contention：多少线程竞争同一资源/地址？
5. Balance：warp/block 间工作是否均匀？

例如 scan 的工作高效版本可能指令更多、同步结构更复杂，但对大规模输入更可扩展；一个原子版本在小输入上可能更快。正确结论总是带输入规模和分布。

练习顺序

• Vector map：建立带宽与计时基线。
• Reduce：naive atomic → shared tree → shuffle → CUB。
• Scan：先实现单 block，再实现分层 device scan，对照 CUB。
• Histogram：测试四种分布，做 block-private 与 warp aggregation。
• Transpose：naive → coalesced shared → padded shared。
• Compact：predicate + scan + scatter，并验证稳定性。
• Stencil：global-only 与 shared halo 两版对照缓存指标。

官方参考

• Cooperative Groups
• CUB Developer Overview
• CUB API
• Thrust API
• CUDA Samples