1. 执行模型：warp、block、grid 与启动配置

CPU 程序员习惯“开 \(N\) 个线程，每个线程跑一段活”。GPU 的并行是三层嵌套：grid 里有很多 block，block 里有很多 warp，warp 由 32 个 lane 组成。写 kernel 时你只写“一个线程做什么”，但硬件实际按 warp 发指令。

1.1 三层结构

grid（一次 kernel 启动）
├── block (0,0)                     ← 一个 block 固定在一个 SM 上
│   ├── warp 0 : lane  0 … 31
│   ├── warp 1 : lane 32 … 63
│   └── …                           ← block 内可以 __syncthreads()
├── block (1,0)  → 另一个 SM      ┐
├── block (0,1)  → 再一个 SM      ├ 执行顺序不确定，不能互相等待
└── …                              ┘

三条规则就够了：

block 是同步的边界。block 内的线程可以 barrier、可以共享 shared memory；block 之间不能，只能靠全局内存 + 原子/fence（第 3 节），或者拆成两个 kernel。
block 是调度的单位。一个 block 一旦被放到某个 SM 上就待在那里跑完；block 到 SM 的映射由硬件决定，没有顺序保证。
warp 是执行的单位。32 个 lane 一起取指、一起执行；分支不同时才拆开。

写索引的时候只用四个内建变量：threadIdx、blockIdx、blockDim、gridDim。二维版本如下：

index2d.cu

__global__ void scale2d(float* out, const float* in, int rows, int cols, float k) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (col >= cols || row >= rows) return;      // 边界检查：最后一个 block 会越界
    out[row * cols + col] = in[row * cols + col] * k;
}

// host 侧
dim3 block(32, 8);                               // 256 线程
dim3 grid((cols + block.x - 1) / block.x,
          (rows + block.y - 1) / block.y);
scale2d<<<grid, block>>>(d_out, d_in, rows, cols, 2.0f);

注意两点。第一，blockDim 是运行时值，编译器不知道它是 32，所以索引乘加不会被常量折叠；追求极致时可以用模板把 block 大小变成编译期常量。第二，二维只是索引方便，硬件调度仍是一维的，block 内的线程编号按 x 最快变化，这会影响你后面设计合并访问的方式。

1.2 warp 与 SIMT：分支为什么会变贵

一个 warp 的 32 个 lane 共享同一个程序计数器。执行 if / else 时，如果 lane 分成两派，硬件会：

用掩码只让“条件为真”的 lane 执行 then 分支；
再用另一个掩码让“条件为假”的 lane 执行 else 分支；
在汇合点（immediate post-dominator）恢复全掩码。

也就是说，两条路径是串行执行的。下面这段代码里，每个 warp 正好一半 lane 走一条路：

divergence.cu

__global__ void divergence_demo(const float* x, float* y, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= n) return;
    float v = x[i];
    if ((threadIdx.x & 31) < 16) v = __expf(v);              // lane 0–15
    else                          v = __logf(fabsf(v) + 1.0f); // lane 16–31
    y[i] = v;
}

如果两条路径各有 10 条指令，这个 warp 实际发射 20 条，而理想情况是 10 条，利用率 50%。更糟的是循环：如果每个线程的循环次数不同，warp 会一直执行到最长的那个，短的 lane 全程空转。

反过来，块内一致（uniform）的分支是免费的：

if (blockIdx.x == 0) { … }        // 整个 block 都走同一边，没有掩码切换
if (threadIdx.x >= n) return;     // 只有最后一个 block 会部分返回，代价可以忽略

warp 级原语可以让你“看见”掩码，从而把分支变成数据操作：

bool positive = v > 0.0f;
unsigned mask = __ballot_sync(0xffffffffu, positive);  // 位 i = lane i 的条件
int npos = __popc(mask);                                // 这个 warp 里有多少个正数

Volta 之后的一个细节：每线程有了独立的程序计数器（independent thread scheduling），warp 内可以更细粒度地交错，但调度单位仍是 warp，__shfl_sync / __syncwarp 也仍以 warp 为单位。不要依赖“两个 lane 会自动在同一时刻汇合”，需要同步就写 __syncwarp(mask)。

1.3 启动配置

语法是 kernel<<<grid, block, sharedBytes, stream>>>(args...)。后两个参数可以省略。取值经验：

block：128–256 线程最稳。必须是 32 的倍数（否则尾 warp 有 lane 永远闲置）；上限 1024。
grid：ceil(N / block)。若要封顶，取 blocksPerSM × SM 数，循环体内用 grid-stride（1.4）。
动态 shared memory：第三个参数传入字节数，kernel 里用 extern __shared__ 取地址。超过 48 KB 时要在 host 侧显式提高上限。

launch.cu

__global__ void scale(float* out, const float* in, int n, float k) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) out[i] = in[i] * k;
}

// 普通启动：grid、block 两个参数
scale<<<(n + 255) / 256, 256>>>(d_out, d_in, n, 2.0f);

// 需要 96 KiB 动态 shared 的 kernel：先提高单 block 上限（H100 约 227 KB）
cudaFuncSetAttribute(myKernel, cudaFuncAttributeMaxDynamicSharedMemorySize, 96 << 10);
myKernel<<<grid, 256, 96 << 10>>>(args...);   // 第三个参数就是字节数

启动配置本身也有成本：每次 launch 约 5 µs。如果 kernel 只跑 2 µs，把 100 次 launch 合成 1 次、用 grid-stride 在设备内循环，收益会比任何微优化都大。

1.4 grid-stride 循环

固定 grid、让每个线程处理多个元素，是 GPU 上的标准写法：

gridstride.cu

__global__ void saxpy(size_t n, float a, const float* x, float* y) {
    size_t stride = (size_t)blockDim.x * gridDim.x;
    for (size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x; i < n; i += stride)
        y[i] = a * x[i] + y[i];
}

// host：grid 封顶在“能常驻的 block 数”，和 N 无关
int perSM = 0;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&perSM, saxpy, 256, 0);
int sms = 0;
cudaDeviceGetAttribute(&sms, cudaDevAttrMultiProcessorCount, 0);
int blocks = perSM * sms;
saxpy<<<blocks, 256>>>(n, 2.0f, d_x, d_y);

好处有三个：任意 \(N\) 都不会越界；launch 开销固定，不会因为 \(N\) 大就启动几百万个 block；索引计算可以复用，每个线程的第一次访存天然连续。代价是每个线程要做一次循环判断，以及尾部的负载不均——相比带宽，这些都不值一提。

注意 stride 用 blockDim.x * gridDim.x，不是 gridDim.x。这是最常见的下标错误之一。

1.5 占用率：有多少 warp 能同时在岗

占用率的定义很直接：

\[ \text{occupancy} \;=\; \frac{\text{每 SM 常驻的 warp 数}}{\text{每 SM 最大 warp 数}}, \]

H100/A100 的上限都是 64 个 warp（2048 线程）。常驻多少由四个资源里最紧的那个决定：

寄存器：每 SM 65536 个 32 位寄存器。每线程 32 个寄存器时，256 线程的 block 占 8192 个 → 8 个 block → 64 warp（100%）。每线程 64 个 → 4 个 block → 32 warp（50%）。
shared memory：每 block 48 KB、SM 有 228 KB → 最多 4 个 block → 50%。
block 数：每 SM 最多 32 个 block。block 只有 32 线程时，32 个 block 也只有 32 warp → 50%。
线程数：每 SM 最多 2048，block 1024 线程时只能放 2 个。

查询真实值不要手算，用运行时 API：

int perSM = 0;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&perSM, myKernel, 256, /*dynamicSMem=*/0);
// 实际占用率 = perSM * (256/32) / 64

经验

占用率不是越高越好。它只是“有多少 warp 可以用来藏延迟”。一个每线程读 16 B、ILP 很深的 kernel，40% 占用率也可能跑满带宽；反过来，100% 占用率但每个 warp 只读 4 B，依然喂不饱内存（回到 0.3 的账）。先把字节/指令比做上去，再看占用率。

控制寄存器用量可以用 __launch_bounds__(maxThreads, minBlocksPerSM) 提示编译器，或编译选项 -maxrregcount=N。前者更精细：它只约束这个 kernel，不会牵连别的。

1.6 和 CPU 的对照

CPU 概念	GPU 对应	差别
std::thread	一个 lane	GPU 线程极轻，创建/切换由硬件完成
SIMD intrinsics	warp 内 32 lane	自动，不需要手写 `_mm256_*`
线程数 8–128	线程数 \(10^4\)–\(10^6\)	并行度高 3–4 个数量级
抢占式调度	warp 轮换，block 不抢占	长 kernel 会占住 SM 直到跑完
pthread_barrier	`__syncthreads()`	只作用于 block 内
thread_local	寄存器 + 局部内存	寄存器溢出会落到“local memory”（其实是显存），非常慢

1.7 容易写错的地方

常见错误

blockDim 不是 32 的倍数，却用 0xffffffff 做 shuffle。最后一个 warp 的 mask 必须按实际活跃 lane 计算，否则是 UB。
用 int 索引大数组。\(2^{31}\) 个 float 只有 8 GiB，超过就溢出。size_t 或 long long。
假设 block 之间有执行顺序。“先算 block 0 再算 block 1”是错的；需要全局顺序就拆 kernel。
launch 之后立刻在 host 上读设备指针。kernel 是异步的，host 直接解引用 → 段错误或垃圾值。
不检查 launch 错误。参数错误在 launch 时就能报，越界要等同步；两者都要查（见 6.1 的 CK 宏）。