""" ============================================================================= CS349E · 第 10 讲共享模块：kernel_lab.py —— 在"纸上的 GPU"上一步步写矩阵乘法内核 ============================================================================= 【这个模块要解决什么问题】 第 10 讲要从最朴素的矩阵乘法内核出发，一步步做到接近厂商库（cuBLAS）的写法： 合并访存 → 共享内存分块 → 一维 / 二维寄存器分块 → 向量化访存 → （warp 分块、Tensor Core）。 我们没有 GPU，所以分两步： 1. **正确性 + 事件计数**：把每个版本的内核按 CUDA 的写法、一句一句写成 Python，放在前置讲 P2 的 numpy GPU 模拟器（simt_lab）上真的执行，和 numpy 的结果比较；同时数出每个版本发了多少条 全局内存 / 共享内存 / FFMA 指令、碰到多少个 cache line 和 sector、共享内存有多少 bank 冲突。 2. **性能估计**：把计数代入一个"多层屋顶"模型（FP32 运算、指令发射、L1、共享内存、显存）， 看每一步优化**解除了哪个瓶颈、又碰到了哪个新瓶颈**，再和公开的实测（Boehm 的博客，RTX A6000）对照。 【本模块在 simt_lab 上加了什么】 - VBlock：向量访存（float4 = 128 位的 LDG.128 / LDS.128 / STS.128）；共享内存 bank 冲突按"相位"计 （32 位访问：一个 warp 一个相位；64 位：两个半 warp；128 位：四个四分之一 warp，见 smem_wavefronts）。 - 全局内存访问额外记录"碰到几个 128 字节的 cache line"（L1 的 wavefront 数，见 l1_wavefronts）。 - 六个矩阵乘法内核 k1–k6（对应 Boehm 博客的 kernel 1–6）。 - 性能模型 predict()、A6000 的参数、Tensor Core 的 mma.m16n8k16 片段布局与 ldmatrix 的 bank 冲突、 一个"Triton 风格"的块级编程小解释器（mini-Triton）。 【运行】python3 kernel_lab.py # 自检：六个内核的结果正确、计数和手算一致 【依赖】numpy；从 ../prereq2 导入 simt_lab（只读）。 """ import os import sys from dataclasses import dataclass, replace import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, os.path.join(HERE, "..", "prereq2")) from simt_lab import Block, KernelStats, cdiv, occupancy, H100 as H100_SIMT # noqa: E402 # --------------------------------------------------------------------------- # 1. 两块 GPU 的结构参数 # --------------------------------------------------------------------------- # RTX A6000（GA102，Ampere，计算能力 8.6）：84 个 SM，每 SM 128 个 FP32 单元，加速频率 1.80 GHz # → FP32 峰值 84 × 128 × 2 × 1.8e9 = 38.7 TFLOPS；GDDR6 768 GB/s；每 SM 最多 1536 个线程（48 个 warp）、 # 16 个块；共享内存每 SM 最多 100 KB、每块最多 99 KB（CUDA 编程指南计算能力 8.6 的那一列）。 # Boehm 的博客说 kernel 3 的占用率是"48 个 warp 里的 32 个"，和这里的 48 一致。 A6000 = replace(H100_SIMT, name="RTX A6000（模拟）", n_sm=84, fp32_lanes_per_sm=128, max_threads_per_sm=1536, max_warps_per_sm=48, max_blocks_per_sm=16, smem_per_sm=100 * 1024, max_smem_per_block=99 * 1024, clock_hz=1.80e9, mem_bw=768e9) H100 = H100_SIMT # simt_lab 里的 H100 SXM：132 SM、1.98 GHz、3.35 TB/s # A100 SXM（GA100，计算能力 8.0）：108 SM，每 SM 只有 64 个 FP32 单元（GA102 是 128），1.41 GHz → FP32 19.5 TFLOPS； # 每 SM 最多 2048 个线程、32 个块、164 KB 共享内存（每块最多 163 KB）；HBM2 1.555 TB/s（40 GB 版）。 A100 = replace(H100_SIMT, name="A100（模拟）", n_sm=108, fp32_lanes_per_sm=64, smem_per_sm=164 * 1024, max_smem_per_block=163 * 1024, clock_hz=1.41e9, mem_bw=1.555e12) LINE_BYTES = 128 # L1 / L2 的 cache line 大小（4 个 32 字节的 sector） # --------------------------------------------------------------------------- # 2. 计数器与带向量访存的 Block # --------------------------------------------------------------------------- @dataclass class KStats(KernelStats): """ 在 simt_lab.KernelStats 上多记几项。 【新增字段】 ffma_instr ：FFMA（FP32 乘加）的 warp 指令数 fma ：线程级的乘加次数（= 活跃线程数 × 每条指令的乘加数） gld_lines / gst_lines：全局 load / store 碰到的 128 字节 cache line 数（按"每条 warp 指令里不同的 line"计） lds_instr / sts_instr：共享内存 load / store 的 warp 指令数（smem_instr 是两者之和） lds_wavefronts / sts_wavefronts：共享内存 load / store 实际花的轮数（含 bank 冲突） ideal_smem_wavefronts：同样的访问如果没有 bank 冲突需要几轮（用来算"冲突倍数"） 【Python 小知识】dataclass 可以继承：子类的字段接在父类后面；父类字段都有默认值，所以子类的也必须有。 """ ffma_instr: int = 0 fma: int = 0 gld_lines: int = 0 gst_lines: int = 0 lds_instr: int = 0 sts_instr: int = 0 lds_wavefronts: int = 0 sts_wavefronts: int = 0 ideal_smem_wavefronts: int = 0 def smem_wavefronts(word, warp, lane, width, n_banks=32, warp_size=32): """ 一条共享内存 warp 指令要几轮（wavefront）才能完成：bank 冲突的"分相位"模型。 【模型】共享内存有 32 个 bank，每个 bank 每个周期送出一个 4 字节的字（word）；一个周期最多送 128 字节。 - 32 位访问（LDS.32）：整个 warp 是一个相位； - 64 位访问（LDS.64）：分成两个相位，每个相位 16 个线程； - 128 位访问（LDS.128，float4）：分成四个相位，每个相位 8 个线程。 每个相位里：同一个 bank 上有 k 个**不同的**字，就要 k 轮（同一个字被多个线程读只算一次 = 广播）； 相位之间不合并。这是 NVIDIA 论坛和 profiler 文档里常见的描述，NVIDIA 没有给出完整的官方规格。 【推论】LDS.128 至少要 4 轮（每个线程拿 16 字节，32 个线程共 512 字节，每轮最多 128 字节）—— 向量化**不会**减少共享内存的轮数，只减少指令数（03 脚本演示 3）。 【参数】word：每个活跃线程访问的**起始**字号（数组）；warp、lane：这些线程的 warp 号和 lane 号； width：每个线程连续访问几个字（1、2、4） 【返回】(实际轮数, 无冲突时的轮数) 【边界】没有活跃线程时返回 (0, 0)。 """ if len(word) == 0: return 0, 0 per_phase = warp_size // width phase = lane // per_phase words = (word[:, None] + np.arange(width)[None, :]).ravel() grp = np.repeat(warp * 8 + phase, width) # 每个 (warp, 相位) 一组；相位最多 4 个 key = np.unique(grp * (1 << 32) + words) # 去掉重复的字（广播） g = key >> 32 bank = (key & ((1 << 32) - 1)) % n_banks gb, counts = np.unique(g * n_banks + bank, return_counts=True) grp_of = gb // n_banks ug, inv = np.unique(grp_of, return_inverse=True) per_group = np.zeros(len(ug), dtype=np.int64) np.maximum.at(per_group, inv, counts) # 无冲突时：每个相位要搬的不同字数 / 32，向上取整（至少 1） _, n_words = np.unique(g, return_counts=True) ideal = np.maximum(1, -(-n_words // n_banks)) return int(per_group.sum()), int(ideal.sum()) def l1_wavefronts(byte_addr, warp, line_bytes=LINE_BYTES): """ 一条全局内存 warp 指令碰到几个不同的 128 字节 cache line（按 warp 分别数，再求和）。 【为什么要数 line】性能模型假设 L1 每个周期处理一个 cache line（一个 wavefront）： 合并得好的 warp 访问（32 个连续的 float = 128 字节）只碰 1 个 line； 朴素内核里 32 个线程各读一行 A，碰 32 个 line，要 32 个周期。 """ if len(byte_addr) == 0: return 0 return len(np.unique(warp * (1 << 40) + byte_addr // line_bytes)) class VBlock(Block): """ simt_lab.Block 加上：FFMA 计数、cache line 计数、向量访存（width = 1、2、4 个元素）。 【新指令】 ffma(n, mask) ：n 条 FFMA warp 指令（只计数；数值请自己用 numpy 算） load_vec(arr, idx, w, mask) ：每个线程从全局内存读 w 个连续元素 arr[idx : idx+w]（一条 LDG.32/64/128） store_vec(arr, idx, val, w) ：写 w 个连续元素（val 形状 [线程数, w]） sload_vec / sstore_vec ：共享内存的向量读写（LDS.128 / STS.128 等） 【边界】向量访问要求 idx 是 w 的倍数（float4 要 16 字节对齐，否则真实 GPU 报"misaligned address"）—— 这里检查并抛出 ValueError。 """ # ---------- 计数 ---------- def ffma(self, n=1, mask=None): act = self._mask(mask) warps = self._issue(act, n) self.stats.ffma_instr += warps * n self.stats.fma += int(act.sum()) * n # ---------- 全局内存 ---------- def _gcount_lines(self, arr, idx, act, w, is_store): if not act.any(): return addr = (idx[act][:, None] + np.arange(w)[None, :]).ravel() * arr.itemsize wv = np.repeat(self.warp[act], w) lines = l1_wavefronts(addr, wv) if is_store: self.stats.gst_lines += lines else: self.stats.gld_lines += lines def load(self, arr, idx, mask=None): act = self._mask(mask) self._gcount_lines(arr, self._idx(idx), act, 1, False) return super().load(arr, idx, mask) def store(self, arr, idx, val, mask=None): act = self._mask(mask) self._gcount_lines(arr, self._idx(idx), act, 1, True) return super().store(arr, idx, val, mask) def _check_vec(self, idx, act, w, n): if w not in (1, 2, 4): raise ValueError("向量宽度只能是 1、2、4") if np.any(idx[act] % w != 0): raise ValueError(f"向量访问没有对齐：下标必须是 {w} 的倍数（真实 GPU 会报 misaligned address）") bad = act & ((idx < 0) | (idx + w > n)) if bad.any(): t = int(np.flatnonzero(bad)[0]) raise IndexError(f"线程 {t} 越界访问下标 {int(idx[t])}..{int(idx[t]) + w - 1}，数组大小 {n}") def load_vec(self, arr, idx, w, mask=None): """ 每个活跃线程读 arr.flat[idx : idx + w]：一条 warp 指令。 【返回】形状 [线程数, w] 的数组；不活跃的线程填 0 """ act = self._mask(mask) idx = self._idx(idx) self._check_vec(idx, act, w, arr.size) out = np.zeros((self.n_threads, w), dtype=arr.dtype) if not act.any(): return out flat = arr.reshape(-1) addr = (idx[act][:, None] + np.arange(w)[None, :]).ravel() * arr.itemsize wv = np.repeat(self.warp[act], w) warps = self._issue(act, 1) self.stats.gld_instr += warps self.stats.gld_sectors += len(np.unique(wv * (1 << 40) + addr // self.gpu.sector_bytes)) self.stats.gld_lines += l1_wavefronts(addr, wv) self.stats.g_bytes_requested += int(act.sum()) * w * arr.itemsize out[act] = flat[idx[act][:, None] + np.arange(w)[None, :]] return out def store_vec(self, arr, idx, val, w, mask=None): """每个活跃线程把 val[t, :w] 写到 arr.flat[idx[t] : idx[t] + w]。""" act = self._mask(mask) idx = self._idx(idx) self._check_vec(idx, act, w, arr.size) if not act.any(): return flat = arr.reshape(-1) tgt = (idx[act][:, None] + np.arange(w)[None, :]) addr = tgt.ravel() * arr.itemsize wv = np.repeat(self.warp[act], w) warps = self._issue(act, 1) self.stats.gst_instr += warps self.stats.gst_sectors += len(np.unique(wv * (1 << 40) + addr // self.gpu.sector_bytes)) self.stats.gst_lines += l1_wavefronts(addr, wv) self.stats.g_bytes_requested += int(act.sum()) * w * arr.itemsize if len(np.unique(tgt.ravel())) < tgt.size: self._race("多个线程在同一条向量 store 里写同一个全局内存位置") flat[tgt.ravel()] = np.asarray(val)[act].reshape(-1) # ---------- 共享内存 ---------- def _smem_count(self, sh, idx, act, w, is_store): n = len(sh.data) bad = act & ((idx < 0) | (idx + w > n)) if bad.any(): t = int(np.flatnonzero(bad)[0]) raise IndexError(f"线程 {t} 越界访问共享内存下标 {int(idx[t])}，大小 {n}") wf, ideal = smem_wavefronts(idx[act], self.warp[act], self.lane[act], w) warps = self._issue(act, 1) self.stats.smem_instr += warps self.stats.smem_wavefronts += wf self.stats.ideal_smem_wavefronts += ideal if is_store: self.stats.sts_instr += warps self.stats.sts_wavefronts += wf else: self.stats.lds_instr += warps self.stats.lds_wavefronts += wf def sload(self, sh, idx, mask=None): act = self._mask(mask) idx2 = self._idx(idx) if act.any(): self._smem_count(sh, idx2, act, 1, False) return self._sread(sh, idx2, act, 1)[:, 0] def sstore(self, sh, idx, val, mask=None): act = self._mask(mask) idx2 = self._idx(idx) if act.any(): self._smem_count(sh, idx2, act, 1, True) v = np.broadcast_to(np.asarray(val), (self.n_threads,))[:, None] self._swrite(sh, idx2, act, v, 1) def sload_vec(self, sh, idx, w, mask=None): """共享内存的向量读（LDS.64 / LDS.128）。返回 [线程数, w]。""" act = self._mask(mask) idx2 = self._idx(idx) self._check_vec(idx2, act, w, len(sh.data)) if act.any(): self._smem_count(sh, idx2, act, w, False) return self._sread(sh, idx2, act, w) def sstore_vec(self, sh, idx, val, w, mask=None): """共享内存的向量写（STS.128 等）。val 形状 [线程数, w]。""" act = self._mask(mask) idx2 = self._idx(idx) self._check_vec(idx2, act, w, len(sh.data)) if act.any(): self._smem_count(sh, idx2, act, w, True) self._swrite(sh, idx2, act, np.asarray(val).reshape(self.n_threads, w), w) def _sread(self, sh, idx, act, w): """读数据，并做和 simt_lab 相同的"写后读"竞争检查（逐个元素）。""" out = np.zeros((self.n_threads, w), dtype=sh.data.dtype) if not act.any(): return out i = (idx[act][:, None] + np.arange(w)[None, :]).ravel() wv = np.repeat(self.warp[act], w) raw = (sh.write_epoch[i] == self.epoch) & (sh.writer_warp[i] != wv) if raw.any(): self._race(f"warp {int(wv[raw][0])} 读共享内存 [{int(i[raw][0])}]，它在本阶段刚被 warp " f"{int(sh.writer_warp[i[raw][0]])} 写过，中间没有 sync（写后读）") sh.reader_warp[i] = wv sh.read_epoch[i] = self.epoch out[act] = sh.data[i].reshape(-1, w) return out def _swrite(self, sh, idx, act, val, w): """写数据，并做"读后写"竞争检查。""" if not act.any(): return i = (idx[act][:, None] + np.arange(w)[None, :]).ravel() wv = np.repeat(self.warp[act], w) war = (sh.read_epoch[i] == self.epoch) & (sh.reader_warp[i] != wv) & (sh.reader_warp[i] >= 0) if war.any(): self._race(f"warp {int(wv[war][0])} 覆盖共享内存 [{int(i[war][0])}]，它在本阶段刚被 warp " f"{int(sh.reader_warp[i[war][0]])} 读过，中间没有 sync（读后写）") sh.data[i] = np.asarray(val)[act].reshape(-1) sh.writer_warp[i] = wv sh.write_epoch[i] = self.epoch def launch(kernel, grid, block, args=(), gpu=A6000): """ 和 simt_lab.launch 相同，只是用 VBlock 和 KStats。 【返回】KStats 【边界】块内线程数超过 1024 时抛出 ValueError。块按顺序一个个执行（真实 GPU 上并行、顺序不定）。 """ def d3(d): d = (d,) if isinstance(d, int) else tuple(d) return d + (1,) * (3 - len(d)) g, b = d3(grid), d3(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 = KStats() for bz in range(g[2]): for by in range(g[1]): for bx in range(g[0]): blk = VBlock(gpu, stats, (bx, by, bz), b, g) kernel(blk, *args) stats.blocks += 1 stats.threads += n return stats # --------------------------------------------------------------------------- # 3. 六个矩阵乘法内核：C[M,N] = A[M,K] · B[K,N]，行主序（row-major），FP32 # --------------------------------------------------------------------------- # 写法尽量逐句对应 Boehm 的 CUDA 代码（github.com/siboehm/SGEMM_CUDA），只去掉了 alpha、beta # （C = αAB + βC 的那两个系数，不影响内存访问的模式；去掉后少读一遍 C）。 # 每个内核的 CUDA 版本写在讲义里；这里的 Python 版本里，"一个 numpy 数组" = "块里每个线程各一个值"。 def k1_naive(blk, A, B, C, M, N, K): """ kernel 1：朴素版。每个线程算 C 的一个元素：一行 A 和一列 B 的点积。 【线程映射】块 32×32；x = 行 = blockIdx.x·32 + threadIdx.x，y = 列 = blockIdx.y·32 + threadIdx.y。 一个 warp 是 threadIdx.x = 0..31（threadIdx.y 相同）→ 32 个线程读 32 个**不同的行**的 A： A[x·K + i] 的地址相差 K·4 字节，碰 32 个 cache line（不合并）；B[i·N + y] 32 个线程是同一个地址（广播）。 【边界】有边界检查（x < M 且 y < N），任何形状都能算。 """ x = blk.bx * blk.bdx + blk.tx y = blk.by * blk.bdy + blk.ty m = (x < M) & (y < N) acc = np.zeros(blk.n_threads, dtype=np.float32) for i in range(K): a = blk.load(A, x * K + i, mask=m) b = blk.load(B, i * N + y, mask=m) blk.ffma(1, mask=m) acc += a * b blk.store(C, x * N + y, acc, mask=m) def k2_coalesced(blk, A, B, C, M, N, K): """ kernel 2：合并访存。只改了"线程 → C 元素"的映射：块改成一维 1024 个线程， 行 = blockIdx.x·32 + threadIdx.x / 32，列 = blockIdx.y·32 + threadIdx.x % 32。 一个 warp（threadIdx.x 连续的 32 个）现在是**同一行、连续 32 列**： A[行·K + i] 32 个线程同一个地址（广播，1 个 line）；B[i·N + 列] 32 个连续的 float = 128 字节（1 个 line）。 """ row = blk.bx * 32 + blk.tid // 32 col = blk.by * 32 + blk.tid % 32 m = (row < M) & (col < N) acc = np.zeros(blk.n_threads, dtype=np.float32) for i in range(K): a = blk.load(A, row * K + i, mask=m) b = blk.load(B, i * N + col, mask=m) blk.ffma(1, mask=m) acc += a * b blk.store(C, row * N + col, acc, mask=m) def k3_smem(blk, A, B, C, M, N, K, BS=32, drop_sync=None): """ kernel 3：共享内存分块（BS×BS = 32×32）。每个块先把 A 的 32×32 块和 B 的 32×32 块搬进共享内存， 再让 1024 个线程从共享内存里读着做 32 步乘加。 【复用】从显存（全局内存）读进来的每个数被块里 32 个线程用到：全局访问少了 32 倍。 但每次乘加仍要从共享内存读 2 个数——瓶颈挪到了共享内存（Boehm 的 profiler："Stall MIO Throttle"）。 【两次 sync】第一次：保证块装满了再读；第二次：保证大家读完了，下一轮才能覆盖。 drop_sync="first" / "second" 去掉其中一次（02 脚本演示 3：模拟器会报告"写后读" / "读后写"竞争）。 【边界】要求 M、N、K 都是 BS 的倍数（Boehm 的版本也这样假设；不整除需要补边界检查，见讲义 §5）。 """ tr, tc = blk.tid // BS, blk.tid % BS As = blk.shared(BS * BS) Bs = blk.shared(BS * BS) acc = np.zeros(blk.n_threads, dtype=np.float32) r0, c0 = blk.bx * BS, blk.by * BS for bk in range(0, K, BS): a = blk.load(A, (r0 + tr) * K + bk + tc) b = blk.load(B, (bk + tr) * N + c0 + tc) blk.sstore(As, tr * BS + tc, a) blk.sstore(Bs, tr * BS + tc, b) if drop_sync != "first": blk.sync() for d in range(BS): av = blk.sload(As, tr * BS + d) bv = blk.sload(Bs, d * BS + tc) blk.ffma(1) acc += av * bv if drop_sync != "second": blk.sync() blk.store(C, (r0 + tr) * N + c0 + tc, acc) def k4_blocktile_1d(blk, A, B, C, M, N, K, BM=64, BN=64, BK=8, TM=8): """ kernel 4：一维寄存器分块。块算 BM×BN = 64×64 的 C，512 个线程，每个线程算**一列里连续 TM = 8 个**元素。 【为什么快】内层循环每读一个 Bs 的数（tmpB），和 8 个 As 的数各做一次乘加： 每次乘加的共享内存读从 2 次降到 (1 + 8)/8 = 1.125 次（Boehm："9K/8 SMEM accesses per result"）。 【线程映射】threadCol = tid % BN、threadRow = tid / BN（0..7）；搬 A 时 innerRowA = tid / BK、innerColA = tid % BK （64 行 × 8 列正好 512 个线程各搬一个）；搬 B 时 innerRowB = tid / BN、innerColB = tid % BN（8 行 × 64 列）。 【边界】M、N、K 是 BM、BN、BK 的倍数；线程数 = BM·BN/TM。 """ tcol, trow = blk.tid % BN, blk.tid // BN iAc, iAr = blk.tid % BK, blk.tid // BK iBc, iBr = blk.tid % BN, blk.tid // BN As = blk.shared(BM * BK) Bs = blk.shared(BK * BN) res = np.zeros((blk.n_threads, TM), dtype=np.float32) r0, c0 = blk.bx * BM, blk.by * BN for bk in range(0, K, BK): blk.sstore(As, iAr * BK + iAc, blk.load(A, (r0 + iAr) * K + bk + iAc)) blk.sstore(Bs, iBr * BN + iBc, blk.load(B, (bk + iBr) * N + c0 + iBc)) blk.sync() for d in range(BK): tb = blk.sload(Bs, d * BN + tcol) for r in range(TM): av = blk.sload(As, (trow * TM + r) * BK + d) blk.ffma(1) res[:, r] += av * tb blk.sync() for r in range(TM): blk.store(C, (r0 + trow * TM + r) * N + c0 + tcol, res[:, r]) def k5_blocktile_2d(blk, A, B, C, M, N, K, BM=128, BN=128, BK=8, TM=8, TN=8, vec_b=False): """ kernel 5：二维寄存器分块。块算 128×128，256 个线程，每个线程算 TM×TN = 8×8 = 64 个元素。 【外积】内层每一步从共享内存读 TM 个 A（regM）和 TN 个 B（regN），做 TM·TN 次乘加： 每次乘加的共享内存读 = (TM + TN)/(TM·TN) = 16/64 = 0.25 次（Boehm："K/4 SMEM accesses per result"）。 【代价】每个线程 64 个累加器 + 16 个 regM/regN 放在寄存器里 → 每线程寄存器多、占用率下降（讲义 §6.3）。 【bank 冲突】一个 warp = threadRow 的 2 个值 × threadCol 的 16 个值： As[(threadRow·TM + i)·BK + d]：两个 threadRow 的地址差 64 个字 → 同一个 bank → 2 路冲突； Bs[d·BN + threadCol·TN + i]：16 个 threadCol 步长 8 个字 → 只落在 4 个 bank 上 → 4 路冲突。 模拟器会数出来（03 脚本）。 【vec_b】True 时把 regN 的 8 个连续数改成两条 LDS.128 读（编译器在能证明对齐时可能自动这样做； 我们看不到 nvcc 生成的 SASS，03 脚本把两种情况都算出来）。As 的 8 个数不连续（步长 BK），不能这样合并。 【边界】M、N、K 是 BM、BN、BK 的倍数；线程数 = (BM/TM)·(BN/TN)。 """ nt = (BM // TM) * (BN // TN) tcol, trow = blk.tid % (BN // TN), blk.tid // (BN // TN) iAc, iAr, strideA = blk.tid % BK, blk.tid // BK, nt // BK iBc, iBr, strideB = blk.tid % BN, blk.tid // BN, nt // BN As = blk.shared(BM * BK) Bs = blk.shared(BK * BN) res = np.zeros((blk.n_threads, TM, TN), dtype=np.float32) r0, c0 = blk.bx * BM, blk.by * BN for bk in range(0, K, BK): for off in range(0, BM, strideA): blk.sstore(As, (iAr + off) * BK + iAc, blk.load(A, (r0 + iAr + off) * K + bk + iAc)) for off in range(0, BK, strideB): blk.sstore(Bs, (iBr + off) * BN + iBc, blk.load(B, (bk + iBr + off) * N + c0 + iBc)) blk.sync() for d in range(BK): regM = np.stack([blk.sload(As, (trow * TM + i) * BK + d) for i in range(TM)], axis=1) if vec_b: regN = np.concatenate([blk.sload_vec(Bs, d * BN + tcol * TN + q, 4) for q in range(0, TN, 4)], axis=1) else: regN = np.stack([blk.sload(Bs, d * BN + tcol * TN + j) for j in range(TN)], axis=1) blk.ffma(TM * TN) res += regM[:, :, None] * regN[:, None, :] blk.sync() for i in range(TM): for j in range(TN): blk.store(C, (r0 + trow * TM + i) * N + c0 + tcol * TN + j, res[:, i, j]) def k6_vectorized(blk, A, B, C, M, N, K, BM=128, BN=128, BK=8, TM=8, TN=8): """ kernel 6：向量化访存。和 kernel 5 一样的分块，但 1. 从显存读 A、B 用 float4（LDG.128，一条指令搬 16 字节）； 2. A 的块在共享内存里**转置**存放（As[k][m]），这样 regM 的 8 个数是连续的，可以用两条 LDS.128 读； 3. regN 也用两条 LDS.128；写 C 用 float4（STG.128）。 【效果】指令数大幅下降（每 8 个数 2 条指令而不是 8 条）；共享内存的**轮数**并不减少（LDS.128 至少 4 轮）， 而 Bs 的读仍有 2 路冲突（Boehm："remaining problems: shared-memory bank conflicts"）。 【边界】同 kernel 5；另外要求 N、K 是 4 的倍数（float4 对齐）。 """ nt = (BM // TM) * (BN // TN) tcol, trow = blk.tid % (BN // TN), blk.tid // (BN // TN) iAr, iAc = blk.tid // (BK // 4), blk.tid % (BK // 4) iBr, iBc = blk.tid // (BN // 4), blk.tid % (BN // 4) rowsA = nt // (BK // 4) # 一轮能搬几行 A（256 / 2 = 128 = BM，一轮搬完） rowsB = nt // (BN // 4) # 一轮能搬几行 B（256 / 32 = 8 = BK） As = blk.shared(BK * BM) Bs = blk.shared(BK * BN) res = np.zeros((blk.n_threads, TM, TN), dtype=np.float32) r0, c0 = blk.bx * BM, blk.by * BN for bk in range(0, K, BK): for off in range(0, BM, rowsA): t = blk.load_vec(A, (r0 + iAr + off) * K + bk + iAc * 4, 4) for q in range(4): # 转置：4 条标量 STS blk.sstore(As, (iAc * 4 + q) * BM + iAr + off, t[:, q]) for off in range(0, BK, rowsB): t = blk.load_vec(B, (bk + iBr + off) * N + c0 + iBc * 4, 4) blk.sstore_vec(Bs, (iBr + off) * BN + iBc * 4, t, 4) blk.sync() for d in range(BK): regM = np.concatenate([blk.sload_vec(As, d * BM + trow * TM + q, 4) for q in range(0, TM, 4)], axis=1) regN = np.concatenate([blk.sload_vec(Bs, d * BN + tcol * TN + q, 4) for q in range(0, TN, 4)], axis=1) blk.ffma(TM * TN) res += regM[:, :, None] * regN[:, None, :] blk.sync() for i in range(TM): for q in range(0, TN, 4): blk.store_vec(C, (r0 + trow * TM + i) * N + c0 + tcol * TN + q, res[:, i, q:q + 4], 4) @dataclass(frozen=True) class KernelSpec: """ 一个内核版本的"启动配置"。 【字段】name：名字；fn：内核函数；block：每块线程数；tile：(BM, BN) 每块算的 C 的大小； smem：每块共享内存字节数；regs：每线程寄存器数的**假设值**（我们看不到编译器的分配， 按"累加器 + 操作数 + 地址"估计，讲义 §6.3）；multiple：要求 M、N、K 是它的倍数； boehm：Boehm 博客报告的 GFLOP/s（RTX A6000，4092² 的 FP32 矩阵乘法） """ name: str fn: object block: tuple tile: tuple smem: int regs: int multiple: int boehm: float KERNELS = [ KernelSpec("k1 朴素", k1_naive, (32, 32), (32, 32), 0, 32, 1, 309.0), KernelSpec("k2 合并访存", k2_coalesced, (1024,), (32, 32), 0, 32, 1, 1986.5), KernelSpec("k3 共享内存分块", k3_smem, (1024,), (32, 32), 2 * 32 * 32 * 4, 40, 32, 2980.3), KernelSpec("k4 一维寄存器分块", k4_blocktile_1d, (512,), (64, 64), (64 * 8 + 8 * 64) * 4, 48, 64, 8474.7), KernelSpec("k5 二维寄存器分块", k5_blocktile_2d, (256,), (128, 128), 2 * 128 * 8 * 4, 128, 128, 15971.7), KernelSpec("k6 向量化访存", k6_vectorized, (256,), (128, 128), 2 * 128 * 8 * 4, 128, 128, 18237.3), ] BOEHM_LATER = {"k9 自动调优（BK=16 等）": 19721.0, "k10 warp 分块": 21779.3, "cuBLAS": 23249.6} def run_kernel(spec, A, B, gpu=A6000, **kw): """ 在模拟器上跑一个内核版本，返回 (C, stats)。 【参数】spec：KernelSpec；A：[M, K]、B：[K, N] 的 float32 数组；kw：传给内核的额外参数 【边界】形状不是 spec.multiple 的倍数时抛出 ValueError（只有 k1、k2 有边界检查）。 """ M, K = A.shape K2, N = B.shape assert K == K2 if any(v % spec.multiple for v in (M, N, K)): raise ValueError(f"{spec.name} 要求 M、N、K 是 {spec.multiple} 的倍数") A = np.ascontiguousarray(A, dtype=np.float32) B = np.ascontiguousarray(B, dtype=np.float32) C = np.zeros((M, N), dtype=np.float32) BM, BN = spec.tile grid = (cdiv(M, BM), cdiv(N, BN)) fn = (lambda blk, *a: spec.fn(blk, *a, **kw)) if kw else spec.fn st = launch(fn, grid, spec.block, args=(A, B, C, M, N, K), gpu=gpu) return C, st # --------------------------------------------------------------------------- # 4. 从小规模的计数外推到大矩阵 # --------------------------------------------------------------------------- COUNT_FIELDS = ("warp_instr", "ffma_instr", "gld_instr", "gst_instr", "gld_lines", "gst_lines", "gld_sectors", "gst_sectors", "lds_instr", "sts_instr", "lds_wavefronts", "sts_wavefronts", "ideal_smem_wavefronts", "barriers", "fma") def per_fma(st): """把计数器都除以乘加次数：这些比值对方阵的大小 n 几乎不变（03 脚本用两个 n 验证）。""" return {f: getattr(st, f) / st.fma for f in COUNT_FIELDS} def scale_counts(st, n_big): """ 把小规模（n_small³）跑出来的计数按乘加次数等比例放大到 n_big³。 【为什么可以】主循环里的每类指令都和乘加次数成正比（每次乘加读几个数、每 8 次乘加一条 LDS……）； 只有"写 C"这类每个输出一次的指令和 n² 成正比，在大矩阵里占比 1/n，可以忽略。 【返回】字典：每个字段放大后的值 """ fma_big = float(n_big) ** 3 return {f: v * fma_big for f, v in per_fma(st).items()} # --------------------------------------------------------------------------- # 5. 性能模型：多层屋顶 # --------------------------------------------------------------------------- @dataclass(frozen=True) class PerfModel: """ 一个 SM 级的"多层屋顶"性能模型（每个资源每周期能处理多少，乘以 SM 数和频率）。 【字段（每个 SM、每个周期）】 fma_per_clk ：FP32 乘加（GA102 / Hopper 都是 128） issue_per_clk ：warp 指令发射（4 个调度器，每个每周期 1 条） l1_lines_per_clk：L1 每周期处理的 cache line 数（假设 1：一个 wavefront 处理一个 128 字节的 line） smem_wf_per_clk ：共享内存每周期 1 轮（128 字节） 【全局】dram_bw（字节/秒） 【假设与局限】 - 各资源完美重叠，时间 = 最慢的那个资源（roofline 的思想，第 1 讲）； - 不模拟延迟和依赖：如果在途的请求不够多（占用率低、没有双缓冲），真实时间会更长（讲义 §7）； - 指令数只数了访存、FFMA 和 sync，没数地址计算和循环控制（真实 SASS 里还有不少）； - "L1 每周期一个 line"是根据 Nsight Compute 文档里 wavefront 概念做的简化假设，不是 NVIDIA 的公开规格。 """ name: str n_sm: int clock_hz: float dram_bw: float fma_per_clk: int = 128 issue_per_clk: int = 4 l1_lines_per_clk: int = 1 smem_wf_per_clk: int = 1 def peak_flops(self): return 2 * self.n_sm * self.fma_per_clk * self.clock_hz PERF_A6000 = PerfModel("RTX A6000", 84, 1.80e9, 768e9) PERF_H100 = PerfModel("H100 SXM（FP32 CUDA core）", 132, 1.98e9, 3.35e12) PERF_A100 = PerfModel("A100（FP32 CUDA core）", 108, 1.41e9, 1.555e12, fma_per_clk=64) def tile_order(tiles_m, tiles_n, group=None): """ 程序实例（线程块）的启动顺序 → 它负责 C 的哪个块 (块行, 块列)。 【参数】group=None：行主序（先走完一行块，再下一行）；group=G：Triton 教程的"分组顺序"—— 先在 G 个块行里竖着走一列，再走下一列（G 行一组，组内列主序）。 【返回】(pm, pn)：两个整数数组，第 i 个元素是第 i 个启动的块的块行、块列 【为什么】同时在跑的那一波块，碰到的不同 A 块行、B 块列越少，L2 的复用就越好（06 脚本演示 2）。 """ pid = np.arange(tiles_m * tiles_n) if group is None: return pid // tiles_n, pid % tiles_n per_group = group * tiles_n gid = pid // per_group first_m = gid * group gsize = np.minimum(tiles_m - first_m, group) pm = first_m + (pid % per_group) % gsize pn = (pid % per_group) // gsize return pm, pn def dram_bytes_tiled(n, bm, bn, blocks_per_sm, n_sm, elem=4, group=None): """ 显存流量的"波次 + K 同步推进"模型：同一波里同时在跑的块按 K 方向齐头并进， 同一时刻用到的 A、B 小块在 L2 里只从显存读一次。 【公式】一波 W = n_sm × blocks_per_sm 个块：每波读 A 的 (碰到的块行数 × bm) 行、B 的 (碰到的块列数 × bn) 列， 各 K = n 个数；最后再加上写 C 一遍。块的顺序由 tile_order(group) 决定。 【边界】假设 L2 能装下一波在一个 K 步里的工作集（几 MB 以内都成立）；不考虑跨波的复用（偏保守）。 """ tiles_m, tiles_n = n // bm, n // bn pm, pn = tile_order(tiles_m, tiles_n, group) total = tiles_m * tiles_n W = max(1, n_sm * blocks_per_sm) byts = 0 for start in range(0, total, W): rows = len(np.unique(pm[start:start + W])) cols = len(np.unique(pn[start:start + W])) byts += (rows * bm + cols * bn) * n * elem return byts + n * n * elem def predict(spec, counts_big, n, perf=PERF_A6000, gpu=A6000): """ 用放大后的计数估计 n×n×n 矩阵乘法的时间和 GFLOP/s。 【返回】字典：每个资源的时间（t_fma、t_issue、t_l1、t_smem、t_dram）、t（最大值）、bound（是哪个）、 gflops、occ（占用率计算器的结果） """ occ = occupancy(spec.block[0] * (spec.block[1] if len(spec.block) > 1 else 1), spec.regs, spec.smem, gpu) cyc = perf.n_sm * perf.clock_hz t = { "t_fma": counts_big["fma"] / (perf.fma_per_clk * cyc), "t_issue": counts_big["warp_instr"] / (perf.issue_per_clk * cyc), "t_l1": (counts_big["gld_lines"] + counts_big["gst_lines"]) / (perf.l1_lines_per_clk * cyc), "t_smem": (counts_big["lds_wavefronts"] + counts_big["sts_wavefronts"]) / (perf.smem_wf_per_clk * cyc), "t_dram": dram_bytes_tiled(n, spec.tile[0], spec.tile[1], max(1, occ["blocks_per_sm"]), perf.n_sm) / perf.dram_bw, } bound = max(t, key=t.get) tt = t[bound] return {**t, "t": tt, "bound": bound[2:], "gflops": 2 * n ** 3 / tt / 1e9, "occ": occ} # --------------------------------------------------------------------------- # 6. Tensor Core：mma.m16n8k16 的片段布局与 ldmatrix 的 bank 冲突 # --------------------------------------------------------------------------- def mma_frag_coords(): """ PTX 指令 mma.sync.aligned.m16n8k16（FP16 输入、FP32 累加）里，每个 lane 持有矩阵的哪些元素。 【布局（PTX ISA "Matrix Fragments for mma.m16n8k16"）】groupID = lane >> 2，t = lane % 4： A（16×16）每个 lane 8 个元素 a0..a7： a0,a1：(groupID, 2t), (groupID, 2t+1) a2,a3：行 + 8 a4,a5：(groupID, 2t+8), (groupID, 2t+9) a6,a7：行 + 8、列 + 8 B（16×8）每个 lane 4 个元素 b0..b3：b0,b1：(2t, groupID), (2t+1, groupID)；b2,b3：行 + 8 C / D（16×8）每个 lane 4 个元素 c0..c3：c0,c1：(groupID, 2t), (groupID, 2t+1)；c2,c3：行 + 8 【返回】字典 {"A": [32, 8, 2], "B": [32, 4, 2], "C": [32, 4, 2]} 的 (行, 列) 坐标数组 【直觉】一条 mma 指令由一个 warp 的 32 个线程**合作**完成：16×16×8 = 2048 次乘加， 每个线程只拿着输入、输出的一小片；数据怎样分布在线程之间是指令规定死的。 """ lane = np.arange(32) g, t = lane >> 2, lane % 4 A = np.zeros((32, 8, 2), dtype=np.int64) for i in range(8): A[:, i, 0] = g + 8 * ((i >> 1) & 1) A[:, i, 1] = 2 * t + (i & 1) + 8 * (i >> 2) Bc = np.zeros((32, 4, 2), dtype=np.int64) for i in range(4): Bc[:, i, 0] = 2 * t + (i & 1) + 8 * (i >> 1) Bc[:, i, 1] = g Cc = np.zeros((32, 4, 2), dtype=np.int64) for i in range(4): Cc[:, i, 0] = g + 8 * (i >> 1) Cc[:, i, 1] = 2 * t + (i & 1) return {"A": A, "B": Bc, "C": Cc} def warp_mma(a_frag, b_frag, c_frag): """ 模拟一条 mma.m16n8k16：输入是 32 个 lane 各自的片段，输出 D 的片段。D = A·B + C（FP32 累加）。 【参数】a_frag：[32, 8]（FP16 值）；b_frag：[32, 4]；c_frag：[32, 4]（FP32） 【返回】d_frag：[32, 4]（FP32） 【思路】按 mma_frag_coords 把片段拼回 16×16 的 A、16×8 的 B、16×8 的 C，算 A·B + C，再拆回各 lane。 真实硬件里没有"拼回"这一步——Tensor Core 直接从这些寄存器取数；拼回只是为了让我们能用 numpy 算。 """ co = mma_frag_coords() A = np.zeros((16, 16), dtype=np.float32) Bm = np.zeros((16, 8), dtype=np.float32) Cm = np.zeros((16, 8), dtype=np.float32) A[co["A"][..., 0], co["A"][..., 1]] = a_frag.astype(np.float32) Bm[co["B"][..., 0], co["B"][..., 1]] = b_frag.astype(np.float32) Cm[co["C"][..., 0], co["C"][..., 1]] = c_frag D = A @ Bm + Cm return D[co["C"][..., 0], co["C"][..., 1]] def swizzle_offset(r, c, row_elems, elem_bytes=2, swizzle=True): """ 共享内存里一个 [行, row_elems] 的 FP16 数据块，逻辑位置 (r, c) 的物理字节地址。 【不 swizzle】行主序：(r·row_elems + c)·2。 【XOR swizzle】把地址看成一串 128 字节的"line"，每个 line 有 8 个 16 字节的"chunk"； 物理 chunk = 逻辑 chunk XOR (line 号 % 8)。同一列的连续 8 行于是落在 8 个不同的 chunk 上 （CuTe 里写作 Swizzle<3,4,3>；CUTLASS、Triton、FlashAttention 都用这类布局）。 【边界】只是重新排列，不浪费空间（和"每行补几个元素"的 padding 办法不同）。 """ o = (np.asarray(r) * row_elems + np.asarray(c)) * elem_bytes if not swizzle: return o line = o // 128 chunk = (o % 128) // 16 return line * 128 + ((chunk ^ (line % 8)) * 16) + o % 16 def ldmatrix_wavefronts(row0, col0, row_elems, swizzle): """ ldmatrix.x4：一个 warp 读 4 个 8×8 的 FP16 矩阵（一个 16×16 的块），每个 lane 提供一行的地址（16 字节）。 每个 8×8 矩阵是一个相位（8 行 × 16 字节 = 128 字节）；数每个相位的 bank 冲突。 【参数】row0, col0：16×16 块的左上角在共享内存数据块里的逻辑坐标；row_elems：每行几个 FP16 【返回】(实际轮数, 理想轮数 = 4) """ total = 0 for mat in range(4): rows = row0 + (mat % 2) * 8 + np.arange(8) col = col0 + (mat // 2) * 8 start = swizzle_offset(rows, col, row_elems, swizzle=swizzle) // 4 # 起始字号 words = (start[:, None] + np.arange(4)[None, :]).ravel() bank = words % 32 cnt = {} for w, b in zip(words.tolist(), bank.tolist()): cnt.setdefault(b, set()).add(w) total += max(len(s) for s in cnt.values()) return total, 4 # --------------------------------------------------------------------------- # 7. mini-Triton：块级编程模型的一个小解释器 # --------------------------------------------------------------------------- class TL: """ 一个"程序实例"（Triton 里叫 program，对应 CUDA 的一个线程块）看到的 tl 命名空间。 【和 CUDA 的区别】CUDA 里你写**一个线程**做什么；Triton 里你写**一个块**做什么： offs = pid * BLOCK + tl.arange(0, BLOCK) # 一个块负责的下标（一个向量） x = tl.load(x_ptr + offs, mask=offs < n) # 一次读一整块 块内怎样分给线程、怎样合并访存、要不要经过共享内存、用不用 Tensor Core、怎样流水——都由编译器决定。 这里的"指针"直接用 (numpy 数组, 下标数组) 表示：load(arr, offs, mask) 读 arr.flat[offs]。 【计数】每个 load / store 记录读写的字节数（只算 mask 为真的元素），供 06 脚本统计 IO。 """ def __init__(self, pid, counters): self._pid = pid self.c = counters def program_id(self, axis): return self._pid[axis] @staticmethod def arange(start, end): return np.arange(start, end) @staticmethod def zeros(shape, dtype=np.float32): return np.zeros(shape, dtype=dtype) def load(self, arr, offs, mask=None, other=0.0): offs = np.asarray(offs) flat = arr.reshape(-1) if mask is None: mask = np.ones(offs.shape, dtype=bool) mask = np.broadcast_to(mask, offs.shape) out = np.full(offs.shape, other, dtype=arr.dtype) out[mask] = flat[offs[mask]] self.c["load_bytes"] += int(mask.sum()) * arr.itemsize self.c["loads"] += 1 return out def store(self, arr, offs, val, mask=None): offs = np.asarray(offs) flat = arr.reshape(-1) if mask is None: mask = np.ones(offs.shape, dtype=bool) mask = np.broadcast_to(mask, offs.shape) flat[offs[mask]] = np.broadcast_to(val, offs.shape)[mask] self.c["store_bytes"] += int(mask.sum()) * arr.itemsize self.c["stores"] += 1 def dot(self, a, b, acc=None): """块级矩阵乘法：FP32 累加（Triton 的 tl.dot 在 GPU 上会变成 Tensor Core 指令）。""" r = a.astype(np.float32) @ b.astype(np.float32) self.c["dot_flops"] += 2 * a.shape[0] * a.shape[1] * b.shape[1] return r if acc is None else acc + r def triton_launch(kernel, grid, *args, **meta): """ 按网格逐个执行程序实例：kernel(tl, *args, **meta)。 【参数】grid：整数或元组（程序实例的个数）；meta：编译期常量（BLOCK_SIZE 等，Triton 里写成 tl.constexpr） 【返回】计数字典：programs、loads、stores、load_bytes、store_bytes、dot_flops、order（程序实例的执行顺序） """ grid = (grid,) if isinstance(grid, int) else tuple(grid) c = {"programs": 0, "loads": 0, "stores": 0, "load_bytes": 0, "store_bytes": 0, "dot_flops": 0} for idx in np.ndindex(*grid[::-1]): pid = idx[::-1] # 第 0 维变化最快（和 CUDA 的 blockIdx.x 一样） kernel(TL(pid, c), *args, **meta) c["programs"] += 1 return c # --------------------------------------------------------------------------- # 自检 # --------------------------------------------------------------------------- def _selftest(): rng = np.random.default_rng(0) n = 128 A = rng.standard_normal((n, n)).astype(np.float32) B = rng.standard_normal((n, n)).astype(np.float32) ref = A.astype(np.float64) @ B.astype(np.float64) for spec in KERNELS: C, st = run_kernel(spec, A, B) err = float(np.abs(C - ref).max()) print(f" {spec.name:<10} 最大误差 {err:.2e} 竞争 {st.race_count} FMA {st.fma}") assert err < 1e-3 and st.race_count == 0 and st.fma == n ** 3 # 手算：k1 的一个 warp 每步读 A 碰 32 个 line、读 B 碰 1 个 line _, st1 = run_kernel(KERNELS[0], A, B) assert st1.gld_lines == (n * n // 32) * n * (32 + 1) # 手算：k2 的一个 warp 每步读 A 1 个 line（广播）、读 B 1 个 line _, st2 = run_kernel(KERNELS[1], A, B) assert st2.gld_lines == (n * n // 32) * n * 2 # 分相位的 bank 模型：32 个线程读连续 32 个字 → 1 轮；LDS.128 读连续 128 个字 → 4 轮 lane = np.arange(32) assert smem_wavefronts(lane, np.zeros(32, int), lane, 1) == (1, 1) assert smem_wavefronts(lane * 4, np.zeros(32, int), lane, 4) == (4, 4) # 步长 32 个字：32 路冲突 assert smem_wavefronts(lane * 32, np.zeros(32, int), lane, 1)[0] == 32 # mma 片段：随机 A、B、C，拆成片段 → warp_mma → 拼回，等于 A·B + C co = mma_frag_coords() Am = rng.standard_normal((16, 16)).astype(np.float16) Bm = rng.standard_normal((16, 8)).astype(np.float16) Cm = rng.standard_normal((16, 8)).astype(np.float32) d = warp_mma(Am[co["A"][..., 0], co["A"][..., 1]], Bm[co["B"][..., 0], co["B"][..., 1]], Cm[co["C"][..., 0], co["C"][..., 1]]) D = np.zeros((16, 8), np.float32) D[co["C"][..., 0], co["C"][..., 1]] = d assert np.allclose(D, Am.astype(np.float32) @ Bm.astype(np.float32) + Cm, atol=1e-5) # 每个元素恰好被一个 lane 持有一次 assert len({(int(r), int(c)) for r, c in co["A"].reshape(-1, 2)}) == 256 # ldmatrix：每行 64 个 FP16（128 字节）时 8 路冲突，swizzle 后无冲突 assert ldmatrix_wavefronts(0, 0, 64, False) == (32, 4) assert ldmatrix_wavefronts(0, 0, 64, True) == (4, 4) print("kernel_lab 自检通过。") if __name__ == "__main__": print("== kernel_lab 自检 ==") _selftest()