先说清楚要做什么

用 Mojo 把矩阵乘法从朴素三重循环加速到 SIMD 加并行，是官方展示语言能力的第一案例；但网上流传的复现大多只有”我跑出了 xx 倍”的结论，没有可复现的实验设计。这篇反其道而行：倍数是结果，实验纪律才是主题。同一份代码走四步演进——朴素 ijk、循环重排、行级并行、SIMD 向量化——每一步都能在你自己的机器上复跑、验证、计时。语法细节不展开，缺口查《Mojo 语法要点》与《Mojo 基础知识参考字典》。”矩阵乘法”这个约束会让每个设计决策都变重：数据怎么进内存、循环按什么顺序走、每一步改动的收益怎么归因，全都绕不开。

实验纪律：三条反模式先立起来

性能实验的结论可信度取决于纪律，代码之前先立规矩。第一条反模式：结果没人用，计算被整体删除。编译器看到矩阵乘出来没人读，有权把整个循环优化掉——计时器测的是空气。解法是算完后打印 checksum，强迫计算真实发生。第二条：单次计时。操作系统调度、频率迁移、缓存冷热都会让单次耗时抖动，预热一次后重复计时取最小值，才能让版本间的差异大于噪声。第三条：跨语言倍数当结论。官方系列博客《How Mojo gets a 35,000x speedup over Python》（2023 年公开口径）里的倍数，基线是纯 Python 朴素循环——换成 NumPy 基线，倍数立刻缩百倍。报倍数必须同时报基线、环境与矩阵规模，本篇的口径是”同语言内四步演进的相对变化”，不背官方倍数。

基础设施：自己搭 Matrix，不用 Tensor

官方论坛在 2025 年 5 月宣布淘汰 Tensor API，25.4 前后的旧教程大量使用 Tensor，照抄会撞上弃用警告甚至编译失败。本篇直接用 UnsafePointer 搭一个最小 Matrix——既绕开 API 迁移，又顺手把 Mojo 的内存模型看明白：

@fieldwise_init
struct Matrix:
    var rows: Int
    var cols: Int
    var data: UnsafePointer[Float64]   # 行主序裸指针，绕开正被淘汰的 Tensor API

    fn __init__(out self, rows: Int, cols: Int):
        self.rows = rows
        self.cols = cols
        self.data = UnsafePointer[Float64].alloc(rows * cols)

    fn get(self, r: Int, c: Int) -> Float64:
        return self.data[r * self.cols + c]

    fn set(mut self, r: Int, c: Int, v: Float64):
        self.data[r * self.cols + c] = v

两个刻意为之的细节。取值用普通方法 get/set 而不是下标协议：下标协议的方法命名在 25.x 间有变动，普通方法跨版本更稳。内存只分配不释放：这是实验代码，跑完即弃；生产代码必须补齐生命周期管理——更要紧的是，持裸指针的 struct 默认拷贝是浅拷贝，两份”Matrix”指向同一块内存，谁先释放另一方就悬空，这是 UnsafePointer 最典型的坑，也解释了为什么名字里带 Unsafe。

fn fill(mut m: Matrix):
    for r in range(m.rows):
        for c in range(m.cols):
            m.set(r, c, Float64((r + c) % 13))   # 确定性数据，结果可复现

fn zero(mut m: Matrix):
    for r in range(m.rows):
        for c in range(m.cols):
            m.set(r, c, 0.0)

fn checksum(m: Matrix) -> Float64:
    var s = 0.0
    for r in range(m.rows):
        for c in range(m.cols):
            s += m.get(r, c)
    return s

四步演进

V0：朴素 ijk 三重循环

fn matmul_naive(a: Matrix, b: Matrix, mut out: Matrix):
    for i in range(a.rows):
        for j in range(b.cols):
            var acc = 0.0
            for k in range(a.cols):
                acc += a.get(i, k) * b.get(k, j)
            out.set(i, j, acc)

这是基线，也是正确性对照的参照物：任何加速版必须和它逐元素一致。a、b 按默认约定只读借用传入，零拷贝；out 标 mut 就地写。

V1：循环重排 ikj，让内存访问连续

fn matmul_ikj(a: Matrix, b: Matrix, mut out: Matrix):
    for i in range(a.rows):
        for k in range(a.cols):
            var aik = a.get(i, k)                 # 提到外层，避免反复下标
            for j in range(b.cols):
                out.set(i, j, out.get(i, j) + aik * b.get(k, j))

数学完全等价，改的是访存模式：行主序存储下，内层 j 沿内存连续推进，b 的一行、out 的一行都能整段吃进缓存行；ijk 版本内层跳着访问 b 的不同行，缓存命中率是天然的短板。这一步展示的是”算法没变，数据摆放方式改变性能”——系统语言的第一课。

V2：行级并行，把核用满

from algorithm import parallelize

fn matmul_parallel(a: Matrix, b: Matrix, mut out: Matrix):
    @parameter
    fn row_work(r: Int):
        for k in range(a.cols):
            var aik = a.get(r, k)
            for j in range(b.cols):
                out.set(r, j, out.get(r, j) + aik * b.get(k, j))

    parallelize[row_work](out.rows)   # 行号分发给各线程，核心数自动定

并行是安全的，因为划分是安全的：每行只被一个线程写，线程间无共享写、无竞争——并行的正确性必须由数据划分保证，不靠运气。这也是矩阵乘法适合并行演示的原因：行与行天然独立。V1 的 ikj 顺序直接搬进并行版，两步优化是叠加关系。

V3：SIMD 向量化，一次算四个

alias W = 4   # Float64 每 4 个一组；CPU 支持 AVX-512 时可试 8，以官方文档为准

fn matmul_simd(a: Matrix, b: Matrix, mut out: Matrix):
    for i in range(a.rows):
        for k in range(a.cols):
            var aik = SIMD[DType.float64, W](a.get(i, k))       # 标量广播成 W 通道
            for j in range(0, b.cols - W + 1, W):               # 主块：整组处理
                var bv = b.data.load[W](k * b.cols + j)         # 一次读 W 个连续元素
                var ov = out.data.load[W](i * b.cols + j)
                ov += aik * bv
                out.data.store[W](i * b.cols + j, ov)           # 一次写回 W 个
            for j in range(b.cols - (b.cols % W), b.cols):      # 尾部：标量收尾
                out.set(i, j, out.get(i, j) + a.get(i, k) * b.get(k, j))

这是 Mojo 与众不同的地方：向量化不是调库，是把”宽度 1 的标量”改写成”宽度 W 的 SIMD”，load/store 按组读写。尾数不足 W 的部分标量收尾——b.cols 不整除 W 时正确性不受影响，正确性验收会兜住它。load/store 的参数拼写随版本微调，以官方文档为准。

计时与验证：可复现的最小骨架

from time import perf_counter_ns

def main():
    var N = 256
    var a = Matrix(N, N)
    var b = Matrix(N, N)
    fill(a)
    fill(b)

    1. 正确性：加速版必须与朴素版逐元素一致
    var expect = Matrix(N, N)
    zero(expect)
    matmul_naive(a, b, expect)
    print("checksum(expect):", checksum(expect))

    1. 计时：预热一次，5 次取最小
    var out = Matrix(N, N)
    var best = -1
    for _ in range(5):
        zero(out)
        var t0 = perf_counter_ns()
        matmul_simd(a, b, out)      # 这里以 V3 为例；四个版本各跑一组互相对比
        var dt = perf_counter_ns() - t0
        if best < 0 or dt < best:
            best = dt
    print("best ms:", Float64(best) / 1e6)
    print("checksum:", checksum(out))   # 防止编译器删除"没人用"的计算

    1. 误差核对
    var max_err = 0.0
    for r in range(N):
        for c in range(N):
            var e = abs(out.get(r, c) - expect.get(r, c))
            if e > max_err:
                max_err = e
    print("max_abs_err:", max_err)

跑之前记录环境表，这是”可复现”的最后一环：CPU 型号、物理核数、操作系统、pixi run mojo --version 的输出。没有这张表，别人无法判断你的趋势是否成立。

结果怎么解读：定性预期，自己复跑

本篇不给编造的实测数字——没在你的机器上跑过的数字没有意义。四步改动的预期方向如下，具体倍数请按上面的骨架自行复跑：

两条诚实边界。第一，N=256 的 Float64 矩阵不足以把缓存体系吃透，更大规模该上分块（tiling）——官方 matmul 系列的下一步就是它，生产环境直接用成熟 BLAS 库，手写循环是学习手段不是替代方案。第二，本篇全部数字口径是”待复跑”，公开口径只有一个：官方博客系列给出的 35,000x（纯 Python 基线、特定硬件），别把它当自己机器的预期。

总结

这次实验教会的东西：性能优化的第一性原理是访存模式（V1 的重排不写一行”快代码”就能改变性能）；并行的正确性来自数据划分而不是小心翼翼；SIMD 在 Mojo 里是类型参数而不是外部库；计时必须防死代码消除、防单次抖动。反过来说，这套流程什么时候不值得用：矩阵规模小到毫秒级以下、瓶颈根本不在计算、或者生产环境里一个 BLAS 调用就能解决——那就直接调库。下一步要把学习过程排成顺序，读《Mojo 学习路线》；开山篇的定位判断（谁该现在跟进）在《Mojo 观察报告》。

版本核实（以官方文档为准）：代码按 25.4（2025-06-18 发布）时点语法书写——参数约定 read/mut/owned/out、@fieldwise_init、UnsafePointer 的 alloc/load/store、SIMD 类型参数化、algorithm 模块的 parallelize、time 模块的 perf_counter_ns 均为该时点已发布事实；Tensor API 正在被官方淘汰（官方论坛 2025-05 声明）故未使用；parallelize 的模块归属、load/store 参数拼写在版本间可能调整，以 docs.modular.com 当前文档为准；官方 35,000x 为 2023 年公开口径（Modular 博客系列，纯 Python 基线），本篇未做实测，所有耗时结论需读者自行复跑。