""" ============================================================================= CS349E · 前置讲 P2 共享模块：simt_lab.py —— 一台"用 numpy 模拟的 GPU" ============================================================================= 【这个模块要解决什么问题】 没有 GPU，怎么学 GPU 编程？我们做一台"纸上的 GPU"：你像写 CUDA 内核一样写一个 Python 函数， 模拟器负责"启动"它（按网格 grid 和线程块 block 的形状），并且**数出**真实 GPU 上决定性能的那些事件： - 全局内存（HBM）访问：每个 warp 的一条 load/store 指令碰到了几个 32 字节的 sector？ （合并访存 coalescing 的好坏就看这个，P2 §5） - 共享内存访问：一个 warp 的 32 个线程落在 32 个 bank 上，有没有冲突（bank conflict）？（P2 §6） - 分支发散（divergence）：一个 warp 里有的线程走 if、有的走 else，两条路都要执行（P2 §4） - 同步：`blk.sync()` 相当于 CUDA 的 __syncthreads()。如果你**忘了写**，模拟器会报告数据竞争（race）。 - 占用率（occupancy）、波次（wave）、Amdahl 定律、内核启动开销的小公式（P2 §7、§8）。 【模拟器怎样"执行"一个内核（重要，和真实 GPU 的差别）】 真实 GPU 上，一个线程块里的每个线程各自执行同一段代码，彼此以不确定的顺序推进。 我们的模拟器把**一个线程块的所有线程放在 numpy 数组里，一条语句一次算完全部线程**： blk.tid → 形状 (线程数,) 的数组：0, 1, 2, ... x = blk.load(A, i) → 所有线程同时从 A 读各自的 A[i] 这相当于"所有线程严格同步地一起执行每条语句"，所以算出来的**结果**总是确定的； 但真实 GPU 上线程并不严格同步，所以模拟器另外记录"谁在哪个同步阶段写过共享内存"， 在你读到别的 warp 写、但中间没有 sync 的数据时报告竞争——真实 GPU 上那样写的程序结果是不确定的。 这种"整块向量化"的写法其实很接近 Triton、numpy 和 TPU（Pallas）的编程模型（P2 §9）。 【写内核的约定】 def my_kernel(blk, A, B, n): # 第一个参数总是 Block 对象 i = blk.bx * blk.bdx + blk.tx # 全局线程下标（数组） m = i < n # 边界检查：一个布尔数组（mask） a = blk.load(A, i, mask=m) # 只有 mask 为 True 的线程真的读 blk.store(B, i, a * 2, mask=m) stats = launch(my_kernel, grid=(n_blocks,), block=(256,), args=(A, B, n)) 【运行】python3 simt_lab.py # 自检：向量加法的结果和计数器是否符合手算 【依赖】numpy。 """ import math from dataclasses import dataclass, field import numpy as np # --------------------------------------------------------------------------- # 1. GPU 参数（H100 SXM，计算能力 9.0） # --------------------------------------------------------------------------- @dataclass(frozen=True) class GPUConfig: """ 一块 GPU 的结构参数。默认值是 H100 SXM（计算能力 9.0）。 【出处】线程/warp/块/寄存器/共享内存的上限：CUDA C++ 编程指南的"Technical Specifications per Compute Capability"表 （9.0 列：每 SM 最多 32 个块、64 个 warp、2048 个线程，64K 个 32 位寄存器，每线程最多 255 个寄存器， 每 SM 最多 228 KB 共享内存，每块最多 227 KB，另有 1 KB 保留给系统，32 个 bank）； SM 数 132、加速频率 1.98 GHz：How to Scale Your Model 第 12 章；HBM 带宽 3.35 TB/s：NVIDIA 数据手册。 【假设值】mem_latency_s（HBM 访问延迟约 500 ns，Gimlet 博客引用的测量）、launch_overhead_s（每次内核启动约 4 µs， 是"几微秒"这个常见量级里的一个代表值，真实值取决于驱动、CPU 和是否用 CUDA Graphs）、 reg_alloc_unit 和 smem_alloc_unit（占用率计算器里的分配粒度）。 """ name: str = "H100-SXM（模拟）" n_sm: int = 132 warp_size: int = 32 schedulers_per_sm: int = 4 fp32_lanes_per_sm: int = 128 max_threads_per_sm: int = 2048 max_warps_per_sm: int = 64 max_blocks_per_sm: int = 32 max_threads_per_block: int = 1024 regs_per_sm: int = 65536 max_regs_per_thread: int = 255 reg_alloc_unit: int = 256 smem_per_sm: int = 228 * 1024 max_smem_per_block: int = 227 * 1024 smem_reserved_per_block: int = 1024 smem_alloc_unit: int = 128 n_banks: int = 32 bank_bytes: int = 4 sector_bytes: int = 32 clock_hz: float = 1.98e9 mem_bw: float = 3.35e12 mem_latency_s: float = 500e-9 launch_overhead_s: float = 4e-6 H100 = GPUConfig() # --------------------------------------------------------------------------- # 2. 统计 # --------------------------------------------------------------------------- @dataclass class KernelStats: """ 一次内核启动的全部计数器。 【字段】 blocks, threads ：启动了多少个块、多少个线程 warp_instr ：发射了多少条 warp 指令（只要 warp 里有一个线程活跃，这条指令就要发射一次） active_lane_slots ：这些指令里一共有多少个"活跃线程 × 指令"——除以 (warp_instr × 32) 就是 SIMT 效率 gld_instr / gst_instr ：全局内存 load / store 的 warp 指令数 gld_sectors / gst_sectors：它们碰到的 32 字节 sector 数（≈ 真正从 HBM 搬的数据，假设没有 cache 复用） g_bytes_requested ：线程们真正想要的字节数（活跃线程数 × 元素大小） smem_instr ：共享内存 load/store 的 warp 指令数 smem_wavefronts ：这些指令实际花的"轮数"——没有冲突时 1 条指令 1 轮，k 路冲突要 k 轮 barriers ：sync() 的次数（按块计） branches / divergent ：遇到分支的 warp 数 / 其中发生发散的 warp 数 races ：检测到的数据竞争描述（最多记前 10 条） race_count ：竞争的总次数 max_smem_per_block ：每块用了多少共享内存 """ blocks: int = 0 threads: int = 0 warp_instr: int = 0 active_lane_slots: int = 0 gld_instr: int = 0 gst_instr: int = 0 gld_sectors: int = 0 gst_sectors: int = 0 g_bytes_requested: int = 0 smem_instr: int = 0 smem_wavefronts: int = 0 barriers: int = 0 branches: int = 0 divergent: int = 0 races: list = field(default_factory=list) race_count: int = 0 max_smem_per_block: int = 0 def dram_bytes(self, sector_bytes=32): """估计的 HBM 流量 = (load + store 的 sector 数) × 32 字节（假设没有 L1/L2 复用）。""" return (self.gld_sectors + self.gst_sectors) * sector_bytes def coalescing_efficiency(self, sector_bytes=32): """合并效率 = 线程想要的字节 / 实际搬运的字节（1.0 表示一个字节都没浪费）。没有访存时返回 1。""" moved = self.dram_bytes(sector_bytes) return self.g_bytes_requested / moved if moved else 1.0 def simt_efficiency(self, warp_size=32): """SIMT 效率 = 平均每条 warp 指令有多少比例的线程在干活。没有指令时返回 1。""" return self.active_lane_slots / (self.warp_instr * warp_size) if self.warp_instr else 1.0 def bank_conflict_factor(self): """平均每条共享内存指令要几轮（1.0 = 没有冲突，32 = 最坏）。没有共享内存访问时返回 1。""" return self.smem_wavefronts / self.smem_instr if self.smem_instr else 1.0 # --------------------------------------------------------------------------- # 3. 共享内存数组 # --------------------------------------------------------------------------- class SharedArray: """ 一个线程块的共享内存数组（只支持 4 字节的类型：float32 / int32 / uint32）。 【额外记录的东西（用来检测竞争）】 writer_warp / write_epoch：每个元素最近一次被哪个 warp、在第几个同步阶段写过 reader_warp / read_epoch ：每个元素最近一次被哪个 warp、在第几个同步阶段读过 "同步阶段"（epoch）= 两次 blk.sync() 之间的那段代码。同一阶段里，一个 warp 读到别的 warp 写的值（写后读）， 或者一个 warp 覆盖了别的 warp 刚读过的值（读后写），在真实 GPU 上结果都是不确定的。 """ def __init__(self, n, dtype): dtype = np.dtype(dtype) if dtype.itemsize != 4: raise ValueError("模拟器的共享内存只支持 4 字节的类型（float32/int32/uint32）") self.data = np.zeros(n, dtype=dtype) self.writer_warp = np.full(n, -1, dtype=np.int64) self.write_epoch = np.full(n, -1, dtype=np.int64) self.reader_warp = np.full(n, -1, dtype=np.int64) self.read_epoch = np.full(n, -1, dtype=np.int64) @property def nbytes(self): return self.data.nbytes # --------------------------------------------------------------------------- # 4. 线程块：内核里能用的全部"内置变量"和"指令" # --------------------------------------------------------------------------- def _dim3(d): """把 int 或 1–3 个元素的元组统一成 (x, y, z)。""" if isinstance(d, int): return (d, 1, 1) d = tuple(d) return d + (1,) * (3 - len(d)) class Block: """ 一个线程块的执行上下文。内核函数的第一个参数就是它。 【内置变量（对应 CUDA）】 tid ：块内线性线程号 0..n-1（数组） ← threadIdx 展平 tx, ty, tz ：块内三维线程号（数组） ← threadIdx.x/y/z bx, by, bz ：块号（整数） ← blockIdx.x/y/z bdx, bdy, bdz ：块的形状（整数） ← blockDim.x/y/z gdx, gdy, gdz ：网格的形状（整数） ← gridDim.x/y/z warp, lane ：每个线程属于块内第几个 warp、是 warp 里第几个（数组） 【指令】load / store（全局内存）、shared / sload / sstore（共享内存）、sync、branch、compute、shfl_down。 【Python 小知识】这里的"数组"都是 numpy 数组，长度 = 块内线程数；对它们做 + - * // % < 等运算， 就等于"每个线程各算一次"。 """ def __init__(self, gpu, stats, block_idx, block_dim, grid_dim): self.gpu, self.stats = gpu, stats self.bx, self.by, self.bz = block_idx self.bdx, self.bdy, self.bdz = block_dim self.gdx, self.gdy, self.gdz = grid_dim n = self.bdx * self.bdy * self.bdz self.n_threads = n self.tid = np.arange(n, dtype=np.int64) self.tx = self.tid % self.bdx self.ty = (self.tid // self.bdx) % self.bdy self.tz = self.tid // (self.bdx * self.bdy) self.warp = self.tid // gpu.warp_size self.lane = self.tid % gpu.warp_size self.n_warps = math.ceil(n / gpu.warp_size) self.epoch = 0 self.smem_bytes = 0 # ---------------- 辅助 ---------------- def _mask(self, mask): """把 mask 变成长度为 n_threads 的布尔数组；None 表示全部线程活跃。""" if mask is None: return np.ones(self.n_threads, dtype=bool) m = np.asarray(mask, dtype=bool) return np.broadcast_to(m, (self.n_threads,)).copy() def _idx(self, idx): """把下标变成长度为 n_threads 的 int64 数组（允许传一个标量，表示所有线程用同一个下标）。""" return np.broadcast_to(np.asarray(idx, dtype=np.int64), (self.n_threads,)) def _issue(self, act, n=1): """ 记录 n 条 warp 指令：只要一个 warp 里有活跃线程，它就要发射这条指令。 【返回】这条指令涉及的 warp 数 """ per_warp = np.bincount(self.warp[act], minlength=self.n_warps) warps = int(np.count_nonzero(per_warp)) self.stats.warp_instr += warps * n self.stats.active_lane_slots += int(act.sum()) * n return warps def _race(self, msg): self.stats.race_count += 1 if len(self.stats.races) < 10: self.stats.races.append(f"块 ({self.bx},{self.by}) 第 {self.epoch} 个同步阶段：{msg}") # ---------------- 计算与分支 ---------------- def compute(self, n=1, mask=None): """ 记录 n 条普通运算指令（加法、乘法……）。真正的数值运算你直接用 numpy 写；这里只负责计数。 【为什么需要】模拟器看不见你写的 numpy 表达式，所以要你告诉它"这里做了几条运算指令"， 它才能统计 SIMT 效率和估计指令发射时间。不调用也不影响结果的正确性。 """ self._issue(self._mask(mask), n) def branch(self, cond, mask=None): """ 一个 if 分支：返回"走 if 这条路的线程"的 mask（= cond 且原本活跃）。 【发散】一个 warp 里如果有的线程 cond 为真、有的为假，硬件只能先执行 if 这条路（让走 else 的线程闲着）， 再执行 else 这条路——两条路都要花时间。模拟器记录发生发散的 warp 数。 【用法】 t = blk.branch(blk.tid % 2 == 0) # 走 if 的线程 ...if 路径的语句用 mask=t... ...else 路径的语句用 mask=~t... """ act = self._mask(mask) c = np.asarray(cond, dtype=bool) & act n_act = np.bincount(self.warp[act], minlength=self.n_warps) n_true = np.bincount(self.warp[c], minlength=self.n_warps) has = n_act > 0 self.stats.branches += int(has.sum()) self.stats.divergent += int(((n_true > 0) & (n_true < n_act)).sum()) self._issue(act, 1) return c def shfl_down(self, val, offset, mask=None): """ warp 内的数据交换（CUDA 的 __shfl_down_sync）：lane i 拿到 lane i + offset 的值； i + offset 超出 warp 的线程拿回自己的值。不经过共享内存，也不需要 sync。 【用途】warp 内的归约：offset = 16, 8, 4, 2, 1 五步就能把 32 个数加到 lane 0（03 脚本演示 4）。 """ act = self._mask(mask) self._issue(act, 1) val = np.asarray(val) src = self.tid + offset ok = (self.lane + offset < self.gpu.warp_size) & (src < self.n_threads) out = val.copy() out[ok] = val[src[ok]] return out # ---------------- 全局内存 ---------------- def _global_access(self, arr, idx, act, is_store): """统计一条全局内存 warp 指令：碰到的 sector 数、请求的字节数。返回检查过的下标。""" if not arr.flags["C_CONTIGUOUS"]: raise ValueError("全局内存数组必须是内存连续的（C-contiguous）；切片或转置后请先 np.ascontiguousarray") if not act.any(): return idx flat_size = arr.size bad = act & ((idx < 0) | (idx >= flat_size)) if bad.any(): t = int(np.flatnonzero(bad)[0]) raise IndexError(f"线程 {t}（块 {self.bx},{self.by}）越界访问下标 {int(idx[t])}，数组大小 {flat_size}：" "真实 GPU 上这是一个可能悄悄读到垃圾、也可能直接崩溃的 bug——记得加边界检查 mask") addr = idx[act] * arr.itemsize sectors = addr // self.gpu.sector_bytes key = self.warp[act] * (1 << 40) + sectors n_sectors = len(np.unique(key)) warps = self._issue(act, 1) if is_store: self.stats.gst_instr += warps self.stats.gst_sectors += n_sectors else: self.stats.gld_instr += warps self.stats.gld_sectors += n_sectors self.stats.g_bytes_requested += int(act.sum()) * arr.itemsize return idx def load(self, arr, idx, mask=None): """ 从全局内存数组 arr（numpy 数组，按展平后的一维下标访问）读：每个活跃线程读 arr.flat[idx]。 【参数】arr：numpy 数组（视为从 0 号地址开始、对齐到 sector 的一块显存）；idx：每个线程的下标； mask：哪些线程真的执行（None = 全部） 【返回】长度为 n_threads 的数组；不活跃线程的位置填 0 【边界】活跃线程的下标越界时抛出 IndexError（提醒你加边界检查）。 """ act = self._mask(mask) idx = self._idx(idx) self._global_access(arr, idx, act, is_store=False) out = np.zeros(self.n_threads, dtype=arr.dtype) out[act] = arr.ravel()[idx[act]] return out def store(self, arr, idx, val, mask=None): """ 写全局内存：每个活跃线程把 val 写到 arr.flat[idx]。 【边界】两个活跃线程写同一个位置时，真实 GPU 上谁最后写进去是不确定的；模拟器记为一次竞争。 """ act = self._mask(mask) idx = self._idx(idx) self._global_access(arr, idx, act, is_store=True) val = np.broadcast_to(np.asarray(val), (self.n_threads,)) tgt = idx[act] if len(np.unique(tgt)) < len(tgt): self._race("多个线程在同一条 store 里写同一个全局内存位置") flat = arr.reshape(-1) flat[tgt] = val[act] # ---------------- 共享内存 ---------------- def shared(self, n, dtype=np.float32): """ 申请一块 n 个元素的共享内存（相当于 CUDA 的 __shared__ float buf[n]）。 【边界】一个块的共享内存总量超过 GPU 的每块上限（H100：227 KB）时抛出 MemoryError。 """ sh = SharedArray(n, dtype) self.smem_bytes += sh.nbytes if self.smem_bytes > self.gpu.max_smem_per_block: raise MemoryError(f"每块共享内存 {self.smem_bytes} 字节，超过上限 {self.gpu.max_smem_per_block}") self.stats.max_smem_per_block = max(self.stats.max_smem_per_block, self.smem_bytes) return sh def _shared_access(self, sh, idx, act): """统计一条共享内存 warp 指令的 bank 冲突：同一个 bank 上不同的 4 字节字有几个，就要几轮。""" if not act.any(): return n = len(sh.data) bad = act & ((idx < 0) | (idx >= n)) if bad.any(): t = int(np.flatnonzero(bad)[0]) raise IndexError(f"线程 {t} 越界访问共享内存下标 {int(idx[t])}，大小 {n}") word = idx[act] # 4 字节的类型：一个元素 = 一个字 bank = word % self.gpu.n_banks w = self.warp[act] key = (w * self.gpu.n_banks + bank) * (1 << 20) + word uniq = np.unique(key) # 同一个字被多个线程读，只算一次（广播） wb, counts = np.unique(uniq >> 20, return_counts=True) per_warp_max = np.zeros(self.n_warps, dtype=np.int64) np.maximum.at(per_warp_max, wb // self.gpu.n_banks, counts) warps = self._issue(act, 1) self.stats.smem_instr += warps self.stats.smem_wavefronts += int(per_warp_max.sum()) def sload(self, sh, idx, mask=None): """ 读共享内存。同时检查"写后读"竞争：本同步阶段里别的 warp 写过、而你没有 sync 就读。 【返回】长度为 n_threads 的数组；不活跃线程填 0 """ act = self._mask(mask) idx = self._idx(idx) self._shared_access(sh, idx, act) i = idx[act] w = self.warp[act] raw = (sh.write_epoch[i] == self.epoch) & (sh.writer_warp[i] != w) if raw.any(): self._race(f"warp {int(w[raw][0])} 读共享内存 [{int(i[raw][0])}]，它在本阶段刚被 warp " f"{int(sh.writer_warp[i[raw][0]])} 写过，中间没有 sync（写后读）") sh.reader_warp[i] = w sh.read_epoch[i] = self.epoch out = np.zeros(self.n_threads, dtype=sh.data.dtype) out[act] = sh.data[i] return out def sstore(self, sh, idx, val, mask=None): """ 写共享内存。同时检查"读后写"竞争：本同步阶段里别的 warp 刚读过这个位置，你就覆盖了它。 """ act = self._mask(mask) idx = self._idx(idx) self._shared_access(sh, idx, act) i = idx[act] w = self.warp[act] war = (sh.read_epoch[i] == self.epoch) & (sh.reader_warp[i] != w) & (sh.reader_warp[i] >= 0) if war.any(): self._race(f"warp {int(w[war][0])} 覆盖共享内存 [{int(i[war][0])}]，它在本阶段刚被 warp " f"{int(sh.reader_warp[i[war][0]])} 读过，中间没有 sync（读后写）") val = np.broadcast_to(np.asarray(val), (self.n_threads,)) sh.data[i] = val[act] sh.writer_warp[i] = w sh.write_epoch[i] = self.epoch def sync(self): """块内同步（CUDA 的 __syncthreads()）：之前的共享内存读写对之后的所有线程都可见。""" self.epoch += 1 self.stats.barriers += 1 self._issue(np.ones(self.n_threads, dtype=bool), 1) # --------------------------------------------------------------------------- # 5. 启动内核 # --------------------------------------------------------------------------- def launch(kernel, grid, block, args=(), gpu=H100): """ 启动一个内核：对网格里的每个块，构造一个 Block，调用 kernel(blk, *args)。 【参数】kernel：内核函数；grid、block：int 或最多 3 个元素的元组；args：传给内核的其他参数；gpu：GPUConfig 【返回】KernelStats 【边界】块内线程数超过上限（1024）时抛出 ValueError——真实 CUDA 会启动失败。 块是按顺序一个一个执行的（真实 GPU 上是并行、顺序不确定）；正确的内核不应该依赖块之间的执行顺序。 """ g, b = _dim3(grid), _dim3(block) n = b[0] * b[1] * b[2] if n > gpu.max_threads_per_block: raise ValueError(f"每块 {n} 个线程，超过上限 {gpu.max_threads_per_block}") stats = KernelStats() for bz in range(g[2]): for by in range(g[1]): for bx in range(g[0]): blk = Block(gpu, stats, (bx, by, bz), b, g) kernel(blk, *args) stats.blocks += 1 stats.threads += n return stats def cdiv(a, b): """向上取整的除法：cdiv(1000, 256) = 4。算网格大小时用：块数 = ⌈元素数 / 每块线程数⌉。""" return -(-a // b) # --------------------------------------------------------------------------- # 6. 粗略的时间估计 # --------------------------------------------------------------------------- def estimate_time(stats, gpu=H100): """ 用计数器粗估这个内核在真 GPU 上至少要多久（roofline 风格的下界 + 启动开销）。 【模型】 t_mem = HBM 流量 / 带宽 （搬数据） t_issue = (warp 指令数 + 共享内存冲突多出来的轮数) / (SM 数 × 每 SM 调度器数 × 频率) （每个调度器每周期发射一条 warp 指令） t = max(t_mem, t_issue) + 启动开销 【局限】不考虑 L1/L2 cache 的复用（HBM 流量可能被高估）、不考虑延迟是否被隐藏（需要足够的占用率，P2 §7）、 不考虑 Tensor Core。只用来比较同一个问题的不同写法谁更省，不是真实时间。 """ t_mem = stats.dram_bytes(gpu.sector_bytes) / gpu.mem_bw extra = stats.smem_wavefronts - stats.smem_instr t_issue = (stats.warp_instr + extra) / (gpu.n_sm * gpu.schedulers_per_sm * gpu.clock_hz) return {"t_mem": t_mem, "t_issue": t_issue, "t": max(t_mem, t_issue) + gpu.launch_overhead_s, "bound": "memory" if t_mem >= t_issue else "issue"} # --------------------------------------------------------------------------- # 7. 占用率、波次、延迟隐藏、Amdahl # --------------------------------------------------------------------------- def occupancy(threads_per_block, regs_per_thread, smem_per_block, gpu=H100): """ 占用率计算器：一个 SM 上最多能同时驻留几个块、几个 warp？被哪种资源卡住？ 【四种限制】 1. 线程（warp）数：每 SM 最多 64 个 warp 2. 块数：每 SM 最多 32 个块 3. 寄存器：每个 warp 占 ⌈每线程寄存器 × 32 / 256⌉ × 256 个寄存器，每 SM 共 65,536 个 4. 共享内存：每块占 ⌈(用量 + 1 KB 系统保留) / 128⌉ × 128 字节，每 SM 共 228 KB 【返回】字典：blocks_per_sm、active_warps、occupancy（= active_warps / 64）、limiter（卡住的资源名）、 以及各项限制下的块数 【边界】每线程寄存器超过 255 或每块共享内存超过 227 KB 时，内核无法启动，返回 blocks_per_sm = 0。 """ warps_per_block = math.ceil(threads_per_block / gpu.warp_size) if regs_per_thread > gpu.max_regs_per_thread or smem_per_block > gpu.max_smem_per_block \ or threads_per_block > gpu.max_threads_per_block: return {"blocks_per_sm": 0, "active_warps": 0, "occupancy": 0.0, "limiter": "无法启动", "limits": {}} regs_per_warp = math.ceil(regs_per_thread * gpu.warp_size / gpu.reg_alloc_unit) * gpu.reg_alloc_unit lim = { "warps": gpu.max_warps_per_sm // warps_per_block, "blocks": gpu.max_blocks_per_sm, "registers": (gpu.regs_per_sm // regs_per_warp) // warps_per_block if regs_per_warp else gpu.max_blocks_per_sm, } if smem_per_block > 0: per_block = math.ceil((smem_per_block + gpu.smem_reserved_per_block) / gpu.smem_alloc_unit) * gpu.smem_alloc_unit lim["shared_mem"] = gpu.smem_per_sm // per_block else: lim["shared_mem"] = gpu.max_blocks_per_sm blocks = min(lim.values()) limiter = min(lim, key=lambda k: lim[k]) active = blocks * warps_per_block return {"blocks_per_sm": blocks, "active_warps": active, "occupancy": active / gpu.max_warps_per_sm, "limiter": limiter, "limits": lim} def waves(n_blocks, blocks_per_sm, gpu=H100): """ 波次（wave）：所有 SM 一次最多同时跑 n_sm × blocks_per_sm 个块；块更多时要分几"波"。 【返回】(波数, 最后一波的利用率, 整体利用率) 整体利用率 = 块数 / (波数 × 每波容量)。例如 133 个块、每 SM 1 个：要 2 波，第二波只有 1 个块在跑。 【边界】blocks_per_sm = 0 时返回 (0, 0, 0)。 """ cap = gpu.n_sm * blocks_per_sm if cap == 0 or n_blocks == 0: return 0, 0.0, 0.0 w = cdiv(n_blocks, cap) last = n_blocks - (w - 1) * cap return w, last / cap, n_blocks / (w * cap) def latency_hiding(active_warps, bytes_per_warp_request, requests_in_flight_per_warp=1, gpu=H100): """ Little 定律在一个 SM 上：在途的字节够不够填满这个 SM 分到的那份带宽？ 【模型】需要的在途字节 = (HBM 带宽 / SM 数) × 访存延迟； 实际在途字节 = 活跃 warp 数 × 每个 warp 同时发出的请求数 × 每个请求的字节数。 能达到的带宽比例 ≈ min(1, 实际 / 需要)。 【返回】(能达到的带宽比例, 需要的在途字节, 实际在途字节) 【局限】假设每个 warp 一直在发请求、没有别的瓶颈；只用来理解"为什么要足够多的 warp 或者每个 warp 多发几个请求"。 """ need = gpu.mem_bw / gpu.n_sm * gpu.mem_latency_s have = active_warps * requests_in_flight_per_warp * bytes_per_warp_request return min(1.0, have / need), need, have def amdahl(serial_frac, n): """Amdahl 定律：加速比 = 1 / (s + (1 − s)/n)。问题规模固定，加更多处理器。""" return 1.0 / (serial_frac + (1.0 - serial_frac) / n) def gustafson(serial_frac, n): """Gustafson 定律：加速比 = s + (1 − s)·n。处理器越多，问题也一起变大（串行部分的时间不变）。""" return serial_frac + (1.0 - serial_frac) * n # --------------------------------------------------------------------------- # 自检 # --------------------------------------------------------------------------- def _vector_add(blk, a, b, c, n): 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) if __name__ == "__main__": n = 1000 a = np.arange(n, dtype=np.float32) b = np.ones(n, dtype=np.float32) c = np.zeros(n, dtype=np.float32) st = launch(_vector_add, grid=(cdiv(n, 256),), block=(256,), args=(a, b, c, n)) print("结果正确：", bool(np.array_equal(c, a + b))) # 1000 个 float32 = 4000 字节 = 125 个 sector；读 a、读 b、写 c 各 125 个 print("sector 数（应为 125 + 125 = 250 次读，125 次写）：", st.gld_sectors, st.gst_sectors) print("合并效率（应为 1.0）：", st.coalescing_efficiency()) print("块数 / 线程数：", st.blocks, st.threads, "（最后一块有 24 个线程越界，被 mask 掉）") print("SIMT 效率：", round(st.simt_efficiency(), 4)) occ = occupancy(256, 32, 0) print("256 线程/块、32 寄存器/线程的占用率：", occ["occupancy"], occ["limiter"])