""" ============================================================================= CS349E · 前置讲 P2 · 02_memory_coalescing.py —— 合并访存与矩阵转置 ============================================================================= 【这个脚本要回答的问题】 GPU 的显存（HBM）带宽很大，但它是以 32 字节的 **sector** 为单位搬数据的。一个 warp 的 32 个线程 同时发出一条 load 指令时，硬件把这 32 个地址落在的 sector 合并起来一起搬——这叫**合并访存**（coalescing）。 32 个线程读 32 个相邻的 float32（128 字节）只要 4 个 sector；如果它们的地址分散，就可能要 32 个 sector， 搬 1024 字节只用了 128 字节。本脚本： 演示 1：一个 warp 的各种访问模式各碰到几个 sector（连续、错位、跨步、广播、随机、更宽的类型、向量化读取）。 演示 2：矩阵转置——朴素写法（读合并、写不合并）、用共享内存分块（读写都合并，但有 bank 冲突）、 加一列填充消除 bank 冲突、以及忘了 sync 的版本。对照 Mark Harris 在 NVIDIA 博客上的实测。 【怎么读输出】全部是模拟（真实计算 + 事件计数），每次运行完全相同。 "估计的 HBM 流量"假设每个 sector 都来自 HBM（不考虑 L1/L2 cache 的复用），会高估不合并写法的代价（见观察）。 【运行】python3 02_memory_coalescing.py （约 2 秒） python3 02_memory_coalescing.py --quick （更快，矩阵更小） 【依赖】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, estimate_time # noqa: E402 QUICK = "--quick" in sys.argv # --------------------------------------------------------------------------- # 演示 1：一个 warp 的访问模式 # --------------------------------------------------------------------------- def pattern_kernel(blk, A, fn): """ 每个线程按 fn(全局线程号) 给出的下标从 A 读一个元素（只统计读，不写结果）。 【参数】A：被读的数组；fn：把全局线程号数组映射成下标数组的函数 【Python 小知识】把函数当参数传进来（fn），就能用同一个内核试很多种访问模式。 """ i = blk.bx * blk.bdx + blk.tx blk.load(A, fn(i)) def demo1_patterns(n_blocks=64, threads=256): print_header("演示 1：一个 warp 读 32 个数，碰到几个 32 字节的 sector？") n = n_blocks * threads big = 64 * n + 64 A32 = np.zeros(big, dtype=np.float32) A64 = np.zeros(big, dtype=np.float64) A16B = np.zeros(4 * big, dtype=np.float32).view("V16") # 每个元素 16 字节（相当于一次读 4 个 float32） rng = np.random.default_rng(0) perm = rng.integers(0, big, size=n) cases = [ ("连续：线程 i 读第 i 个 float32", A32, lambda i: i), ("错位 1 个元素：读第 i+1 个", A32, lambda i: i + 1), ("跨 2 个元素", A32, lambda i: 2 * i), ("跨 4 个元素", A32, lambda i: 4 * i), ("跨 8 个元素（32 B）", A32, lambda i: 8 * i), ("跨 32 个元素（128 B）", A32, lambda i: 32 * i), ("广播：所有线程读同一个数", A32, lambda i: np.zeros_like(i)), ("随机下标", A32, lambda i: perm[i]), ("连续读 float64（8 B）", A64, lambda i: i), ("连续读 16 B（向量化读取，如 float4）", A16B, lambda i: i), ] rows = [] for name, arr, fn in cases: st = launch(pattern_kernel, grid=(n_blocks,), block=(threads,), args=(arr, fn)) per_instr = st.gld_sectors / st.gld_instr req = st.g_bytes_requested / st.gld_instr rows.append([name, f"{per_instr:.2f}", f"{req:.0f} B", f"{per_instr * 32:.0f} B", f"{st.coalescing_efficiency() * 100:.1f}%"]) print_table(["访问模式（一个 warp 32 个线程）", "每条指令的 sector 数", "线程要的字节", "实际搬的字节", "合并效率"], rows) print_observations([ "连续访问：32 个 float32 正好落在 4 个 sector（128 字节）里，一个字节都不浪费。这是 GPU 上最重要的性能规则。", "错位一个元素：首地址没有对齐到 sector，128 字节横跨 5 个 sector，效率 80%。", "跨 8 个元素（32 字节）以后，每个线程独占一个 sector：搬 1024 字节只用 128 字节，效率掉到 12.5%——" "和 P1 §2.5 的'跨度达到一整行，cache line 只用了一点点'是同一件事，只是单位从 cache line 换成了 sector。", "广播（所有线程读同一个地址）只要 1 个 sector，'合并效率'超过 100% 只是因为 32 个线程要的是同一份 4 字节；" "随机下标几乎每个线程一个 sector。", "更宽的类型（float64、16 字节向量读取）连续访问时效率仍是 100%，而且一条指令搬的字节更多（256 B、512 B）——" "同样多的 warp 能有更多数据在途，这正是 P1 §4.4 说的'向量化读取'帮助隐藏延迟的原因（04 脚本）。", ]) # --------------------------------------------------------------------------- # 演示 2：矩阵转置 # --------------------------------------------------------------------------- TILE, BLOCK_ROWS = 32, 8 # 和 Harris 的博客一样：块 32×8 个线程，每块负责 32×32 的一块，每个线程处理 4 个元素 def copy_kernel(blk, A, B, N): """基线：B = A（读写都按行连续，完全合并）。性能目标：转置应该尽量接近它。""" x = blk.bx * TILE + blk.tx y = blk.by * TILE + blk.ty for j in range(0, TILE, BLOCK_ROWS): v = blk.load(A, (y + j) * N + x) blk.store(B, (y + j) * N + x, v) def transpose_naive(blk, A, B, N): """ 朴素转置 B[x][y] = A[y][x]：读 A 时相邻线程读同一行的相邻元素（合并）； 写 B 时相邻线程写同一**列**的相邻元素——地址相隔 N 个 float（不合并）。 【边界】为了简洁，要求 N 是 32 的倍数（不做边界检查）。 """ x = blk.bx * TILE + blk.tx y = blk.by * TILE + blk.ty for j in range(0, TILE, BLOCK_ROWS): v = blk.load(A, (y + j) * N + x) blk.store(B, x * N + (y + j), v) def transpose_tiled(blk, A, B, N, pad=0, do_sync=True): """ 用共享内存分块转置（Harris 的 transposeCoalesced / transposeNoBankConflicts）。 【思路】 1. 整块 32×32 先按行读进共享内存 tile（读全局内存：合并）； 2. sync：确保整块都写进来了（下一步每个线程要读的是**别的线程**写的元素）； 3. 按列从 tile 里读出来，按行写到 B 的转置位置（写全局内存：合并）。 "转置"这一步在共享内存里完成，全局内存的读写都是连续的。 【bank 冲突】tile 每行 32 个 float 时，第 3 步按列读：同一个 warp 的 32 个线程读 tile[0][c]、tile[1][c]…， 地址相隔 32 个字（4 字节），全部落在同一个 bank——32 路冲突，要 32 轮。 每行填充 1 个元素（pad=1，宽 33），相隔 33 个字，32 个线程正好落在 32 个不同的 bank 上。 【参数】pad：每行填充几个元素；do_sync：是否在第 2 步同步（False 用来演示竞争） """ tw = TILE + pad tile = blk.shared(TILE * tw) x = blk.bx * TILE + blk.tx y = blk.by * TILE + blk.ty for j in range(0, TILE, BLOCK_ROWS): v = blk.load(A, (y + j) * N + x) blk.sstore(tile, (blk.ty + j) * tw + blk.tx, v) if do_sync: blk.sync() x = blk.by * TILE + blk.tx # 转置后的块的位置 y = blk.bx * TILE + blk.ty for j in range(0, TILE, BLOCK_ROWS): v = blk.sload(tile, blk.tx * tw + blk.ty + j) blk.store(B, (y + j) * N + x, v) def demo2_transpose(N=256 if QUICK else 512): print_header(f"演示 2：{N}×{N} float32 矩阵转置（块 32×8 线程，每块处理 32×32）") A = np.arange(N * N, dtype=np.float32).reshape(N, N) grid = (N // TILE, N // TILE) kernels = [ ("copy（基线，不是转置）", copy_kernel, {}, A), ("朴素转置", transpose_naive, {}, A.T), ("共享内存分块", transpose_tiled, {"pad": 0}, A.T), ("共享内存分块 + 填充一列", transpose_tiled, {"pad": 1}, A.T), ("分块 + 填充，但忘了 sync", transpose_tiled, {"pad": 1, "do_sync": False}, A.T), ] rows = [] t_copy = None for name, fn, kw, expect in kernels: B = np.zeros_like(A) def kern(blk, A_, B_, N_, fn=fn, kw=kw): fn(blk, A_, B_, N_, **kw) st = launch(kern, grid=grid, block=(TILE, BLOCK_ROWS), args=(A, B, N)) est = estimate_time(st) t_copy = t_copy or est["t_mem"] rows.append([name, f"{st.gld_sectors:,}", f"{st.gst_sectors:,}", f"{st.coalescing_efficiency() * 100:.0f}%", f"{st.bank_conflict_factor():.1f}", st.race_count, f"{est['t_mem'] / t_copy:.2f}×", "对" if np.array_equal(B, expect) else "错"]) print_table(["内核", "读 sector", "写 sector", "合并效率", "共享内存每条指令几轮", "检测到的竞争", "估计 HBM 时间（相对 copy）", "结果"], rows) print("\n对照：Mark Harris《An Efficient Matrix Transpose in CUDA C/C++》（NVIDIA 博客，1024×1024，有效带宽 GB/s）") print_table(["内核", "Tesla M2050", "Tesla K20c"], [["copy", "105.2", "136.0"], ["copySharedMem", "104.6", "152.3"], ["transposeNaive", "18.8", "55.3"], ["transposeCoalesced（分块）", "51.3", "97.6"], ["transposeNoBankConflicts（分块 + 填充）", "99.5", "144.3"]]) print_observations([ "朴素转置的读和 copy 一样，但写的时候相邻线程写相隔一整行的地址：每个线程独占一个 sector，写的流量是 copy 的 8 倍。", "分块后，全局内存的读写都合并了，HBM 流量回到和 copy 一样——'转置'这件不规整的事被挪到了共享内存里做。", "但不填充时，从 tile 按列读是 32 路 bank 冲突：读和写平均下来每条共享内存指令要 16.5 轮。填充一列后回到 1 轮。" "这和 P1 §2.5 的 cache 冲突缺失是同一个道理、同一个解法：2 的幂次方的跨度，填一个元素错开。", "Harris 的实测印证了这个顺序：朴素转置只有 copy 的 18%（M2050）到 41%（K20c）；分块后提到 49%–72%；" "消除 bank 冲突后达到 copy 的约 95%。", "模拟器估计朴素转置的 HBM 时间是 copy 的约 4.5 倍；Harris 实测在 M2050 上慢 5.6 倍、在 K20c 上只慢 2.5 倍。" "真实 GPU 的 L2 cache 会把相邻 warp 写同一个 sector 的部分合并起来（新一代做得更好），模拟器没有建模这一点。" "计数器是用来比较写法、找出问题的，不是用来预测时间的。", "忘了 sync 的版本，模拟器算出的结果是对的（它让所有线程严格同步地执行每条语句），但检测到了大量'写后读'竞争：" "真实 GPU 上，一个 warp 可能在别的 warp 还没把数据写进 tile 时就去读，结果是错的、而且每次运行可能不一样。", ]) if __name__ == "__main__": demo1_patterns() demo2_transpose()