""" ============================================================================= CS349E · 前置讲 P2 · 01_cuda_programming_model.py —— CUDA 编程模型：kernel、thread、block、grid ============================================================================= 【这个脚本要回答的问题】 GPU 上的程序（内核，kernel）长什么样？成千上万个线程怎样知道"自己该处理哪个数据"？ 本脚本用 simt_lab 模拟器把 CUDA 编程模型的几个核心概念走一遍： 演示 1：向量加法——CUDA C++ 原文 vs 模拟器写法；块大小怎么选；为什么最好是 32 的倍数。 演示 2：忘了边界检查会怎样。 演示 3：二维网格——给矩阵的每一行加一个偏置；线程怎样按二维编号、又怎样被打包成 warp。 演示 4：grid-stride 循环——用固定数量的线程处理任意大小的数组。 【怎么读输出】全部是模拟（真实计算 + 事件计数），每次运行完全相同。 【运行】python3 01_cuda_programming_model.py （约 1 秒） 【依赖】numpy；本目录的 simt_lab；../lecture01 的 console。 """ import os import sys import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, os.path.join(HERE, "..", "lecture01")) from console import print_header, print_observations, print_table # noqa: E402 from simt_lab import launch, cdiv # noqa: E402 CUDA_VECTOR_ADD = r""" // CUDA C++：每个线程负责 c[i] = a[i] + b[i] 中的一个 i __global__ void vector_add(const float* a, const float* b, float* c, int n) { int i = blockIdx.x * blockDim.x + threadIdx.x; // 我是全局第几个线程？ if (i < n) { // 边界检查：最后一个块可能有多余的线程 c[i] = a[i] + b[i]; } } // 在 CPU（host）上启动：⌈n / 256⌉ 个块，每块 256 个线程 int threads = 256; int blocks = (n + threads - 1) / threads; vector_add<<>>(d_a, d_b, d_c, n); """ def vector_add(blk, a, b, c, n): """ 模拟器版的向量加法。和上面的 CUDA 代码逐行对应： i = blockIdx.x * blockDim.x + threadIdx.x → i = blk.bx * blk.bdx + blk.tx （i 是一个数组：每个线程一个值） if (i < n) → m = i < n，之后的读写都带 mask=m c[i] = a[i] + b[i] → load、load、加法、store 【参数】blk：模拟器传进来的 Block；a、b、c：numpy 数组（"显存"）；n：元素个数 【边界】i ≥ n 的线程被 mask 掉，不读不写。 """ i = blk.bx * blk.bdx + blk.tx m = i < n x = blk.load(a, i, mask=m) y = blk.load(b, i, mask=m) blk.compute(1, mask=m) blk.store(c, i, x + y, mask=m) def vector_add_no_check(blk, a, b, c, n): """故意漏掉边界检查的版本（演示 2）。""" i = blk.bx * blk.bdx + blk.tx x = blk.load(a, i) y = blk.load(b, i) blk.store(c, i, x + y) # --------------------------------------------------------------------------- # 演示 1：向量加法与块大小 # --------------------------------------------------------------------------- def demo1_vector_add(n=10_000): print_header("演示 1：向量加法——CUDA 原文与模拟器写法；块大小怎么选") print(CUDA_VECTOR_ADD) a = np.arange(n, dtype=np.float32) b = np.full(n, 0.5, dtype=np.float32) rows = [] for bs in (32, 100, 128, 256, 1000, 1024): c = np.zeros(n, dtype=np.float32) st = launch(vector_add, grid=(cdiv(n, bs),), block=(bs,), args=(a, b, c, n)) warps_per_block = cdiv(bs, 32) rows.append([bs, st.blocks, st.threads, st.blocks * warps_per_block, st.threads - n, f"{st.simt_efficiency() * 100:.1f}%", f"{st.coalescing_efficiency() * 100:.0f}%", "对" if np.array_equal(c, a + b) else "错"]) print_table(["每块线程数", "块数", "线程总数", "warp 总数", "多余线程", "SIMT 效率", "合并效率", "结果"], rows) print_observations([ "每个线程只算一个元素：'我是第几个线程'= 块号 × 块大小 + 块内线程号。这一行几乎出现在每一个 CUDA 内核里。", "块数用向上取整的除法 ⌈n / 块大小⌉ 算出来，所以最后一个块通常有多余的线程；if (i < n) 让它们什么都不做。", "块大小 100 时每块有 4 个 warp，但最后一个 warp 只有 4 个线程在干活（100 = 3×32 + 4）：硬件按 warp 调度，" "不满的 warp 照样占一个调度名额。块大小应取 32 的倍数——CUDA 编程指南也这样建议。", "块大小 100 还有一个意外的代价：合并效率只有 86%。第 k 个块从第 100k 个元素（第 400k 字节）开始，" "大多数块的起点没有对齐到 32 字节的 sector，一个 warp 的访问会跨进相邻的 sector，多搬了数据（02 脚本详细讲）。", "块大小 128、256、1024（以及起点恰好对齐的 1000：4000 字节是 32 的倍数）的计数几乎一样，合并效率都是 100%：块大小主要影响的是" "占用率和资源分配（04 脚本），不是访存模式。", ]) # --------------------------------------------------------------------------- # 演示 2：忘了边界检查 # --------------------------------------------------------------------------- def demo2_bounds(n=1000): print_header("演示 2：忘了边界检查会怎样？") a = np.arange(n, dtype=np.float32) b = np.ones(n, dtype=np.float32) c = np.zeros(n, dtype=np.float32) try: launch(vector_add_no_check, grid=(cdiv(n, 256),), block=(256,), args=(a, b, c, n)) print(" （没有报错——说明 n 恰好是 256 的倍数）") except IndexError as e: print(f" 模拟器报错：{e}") print_observations([ "n = 1000 不是 256 的倍数，最后一个块的第 232–255 号线程会去读 a[1000]…a[1023]——数组外面的内存。", "模拟器在这里直接报错；真实 GPU 上越界读可能读到别的数据（结果悄悄出错），也可能触发非法地址错误让整个程序崩溃，" "越界写更危险。这是 CUDA 新手最常见的 bug 之一，排查工具是 NVIDIA 的 compute-sanitizer。", ]) # --------------------------------------------------------------------------- # 演示 3：二维网格 # --------------------------------------------------------------------------- def add_row_bias(blk, X, bias, M, N): """ X[r, c] += bias[r]，X 是 M × N、按行存储。用二维的块和二维的网格。 【下标】行 r = blockIdx.y × blockDim.y + threadIdx.y；列 c = blockIdx.x × blockDim.x + threadIdx.x。 x 方向对应列：这样 threadIdx.x 相邻的线程访问同一行里相邻的元素（地址连续，合并访存）。 【边界】行、列都要检查：(r < M) & (c < N)。 """ r = blk.by * blk.bdy + blk.ty c = blk.bx * blk.bdx + blk.tx m = (r < M) & (c < N) x = blk.load(X, r * N + c, mask=m) bb = blk.load(bias, r, mask=m) blk.compute(1, mask=m) blk.store(X, r * N + c, x + bb, mask=m) def add_row_bias_swapped(blk, X, bias, M, N): """同样的运算，但让 x 方向对应**行**（相邻线程访问相隔一整行的元素）——一个常见的错误写法。""" r = blk.bx * blk.bdx + blk.tx c = blk.by * blk.bdy + blk.ty m = (r < M) & (c < N) x = blk.load(X, r * N + c, mask=m) bb = blk.load(bias, r, mask=m) blk.compute(1, mask=m) blk.store(X, r * N + c, x + bb, mask=m) def demo3_2d(M=100, N=70): print_header(f"演示 3：二维网格——给 {M}×{N} 矩阵的每一行加偏置（块 16×16）") X0 = np.arange(M * N, dtype=np.float32).reshape(M, N) bias = np.arange(M, dtype=np.float32) * 1000 expect = X0 + bias[:, None] rows = [] for fn, label in ((add_row_bias, "x 方向 = 列（正确的映射）"), (add_row_bias_swapped, "x 方向 = 行（错误的映射）")): X = X0.copy() if fn is add_row_bias: grid = (cdiv(N, 16), cdiv(M, 16)) else: grid = (cdiv(M, 16), cdiv(N, 16)) st = launch(fn, grid=grid, block=(16, 16), args=(X, bias, M, N)) rows.append([label, f"{grid[0]}×{grid[1]}", st.gld_sectors + st.gst_sectors, f"{st.coalescing_efficiency() * 100:.0f}%", "对" if np.array_equal(X, expect) else "错"]) print_table(["映射方式", "网格", "访存 sector 数", "合并效率", "结果"], rows) print("\n一个 16×16 的块怎样被切成 warp（warp = 线程号 tid = tx + 16·ty 连续的 32 个）：") tid = np.arange(256) tx, ty = tid % 16, tid // 16 for w in range(3): sel = tid // 32 == w print(f" warp {w}：ty ∈ {{{', '.join(str(v) for v in sorted(set(ty[sel])))}}}，tx = 0..15") print_observations([ "二维问题用二维的块和网格来写更自然：每个线程用 (blockIdx, threadIdx) 的 x、y 分量算出自己负责的行和列。", "线程按 x 最快变化的顺序编号，再每 32 个打包成一个 warp：16×16 的块里，一个 warp 是两行各 16 个线程。", "所以 x 方向应该对应内存里连续的那一维（按行存储时是列）。两种映射结果都对，但错误的映射让相邻线程访问相隔一整行的地址，" "要多搬好几倍的数据——这就是下一个脚本的主题：合并访存。", ]) # --------------------------------------------------------------------------- # 演示 4：grid-stride 循环 # --------------------------------------------------------------------------- def scale_grid_stride(blk, x, y, alpha, n): """ y[i] = alpha · x[i]，用 grid-stride 循环：每个线程从自己的全局编号出发，每次跳过"全网格的线程总数"。 for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += gridDim.x * blockDim.x) 【好处】网格大小和数据大小脱钩：可以按"刚好填满 GPU"来定网格大小，数组再大也不用改启动参数； 线程的初始化开销（算下标、读参数）被摊到多个元素上。 """ stride = blk.gdx * blk.bdx start = blk.bx * blk.bdx + blk.tx k = 0 while True: i = start + k * stride m = i < n if not m.any(): break v = blk.load(x, i, mask=m) blk.compute(1, mask=m) blk.store(y, i, alpha * v, mask=m) k += 1 def demo4_grid_stride(n=100_000): print_header(f"演示 4：grid-stride 循环——固定的网格处理 {n:,} 个元素") x = np.arange(n, dtype=np.float32) rows = [] for blocks, threads in ((4, 128), (16, 256), (cdiv(n, 256), 256)): y = np.zeros(n, dtype=np.float32) st = launch(scale_grid_stride, grid=(blocks,), block=(threads,), args=(x, y, np.float32(3.0), n)) total = blocks * threads rows.append([f"{blocks} × {threads}", total, f"{cdiv(n, total)}", st.gld_instr, f"{st.coalescing_efficiency() * 100:.0f}%", "对" if np.array_equal(y, 3 * x) else "错"]) print_table(["网格 × 块", "线程总数", "每线程处理几个元素", "load 指令（warp 级）", "合并效率", "结果"], rows) print_observations([ "三种网格大小都得到正确结果，load 指令总数也一样：工作量是由数据决定的，网格只决定'谁来做'。", "每一轮里，全网格的线程访问一段连续的地址（相邻线程相邻地址），所以合并效率仍是 100%。", "实践中常把网格大小设成'SM 数的几倍'，让每个 SM 都有活干，又不至于启动太多块（04 脚本的波次）。", ]) if __name__ == "__main__": demo1_vector_add() demo2_bounds() demo3_2d() demo4_grid_stride()