前言

上回把三重循环换了下顺序，矩阵乘法从9秒干到1.5秒，当时对照的OpenBLAS成绩是0.05秒，中间还隔着三十倍，这三十倍总不能全赖人家多线程吧。翻了一堆资料把这笔账算明白了，大头其实藏在两个地方：CPU的指令级并行，和矩阵的分解策略。

CPU早就不是一次算一个数了

写代码的人很容易有个惯性思维：一行 a[i]*b[i] 就是做了一次乘法，CPU老老实实算完一个再算下一个。这是指令串行的想象，真实情况是现代CPU一条SIMD指令一次就能算一批：AVX2的寄存器256位宽，塞得下4个double，一条乘加指令（FMA）出发就是4乘4加；再加上主流CPU有两个FMA运算端口，一个时钟周期理论吞吐是16个浮点操作。标量代码再怎么优化，每个周期也就吃掉这一小口，理论峰值的大餐一口没动。

拿最熟悉的点积写个最小对比，标量串行版和AVX指令版：

#include <immintrin.h>

// 指令串行版: 一次乘一个, 乘和加分两步
double dot_scalar(const double *a, const double *b, int n) {
    double s = 0;
    for (int i = 0; i < n; i++)
        s += a[i] * b[i];
    return s;
}

// AVX2版: 一次装4个double, 一条FMA指令同时完成4乘4加
double dot_avx(const double *a, const double *b, int n) {
    __m256d s = _mm256_setzero_pd();
    for (int i = 0; i < n; i += 4)
        s = _mm256_fmadd_pd(_mm256_loadu_pd(a + i),
                            _mm256_loadu_pd(b + i), s);
    // 4条车道归约成1个数
    __m128d lo = _mm256_castpd256_pd128(s);
    __m128d hi = _mm256_extractf128_pd(s, 1);
    lo = _mm_add_pd(lo, hi);
    lo = _mm_add_pd(lo, _mm_permute_pd(lo, 1));
    return _mm_cvtsd_f64(lo);
}

编译时给足条件，gcc的自动向量化才有用武之地：

gcc -O2 -march=native dot.c -o dot        # AVX版直接生效
gcc -O2 -fno-tree-vectorize dot.c -o dot  # 强制标量串行做对照

一百万长度的点积，串行版约1.4毫秒，AVX版约0.4毫秒，三倍多。这就是"指令并行"四个字的直接含义：不是多开线程，而是单条指令内部就把并行吃干净了。矩阵乘法里最内层就是一个点积，把micro层面的这口理论峰值吃上，宏观成绩立刻改观。

寄存器分块：让数据坐上直达车

光有SIMD还不够，还得让喂给FMA的数据够快。寄存器离运算器最近，所以终极目标是让参与计算的小块常驻寄存器：把C切出4×4的小块，A的4行和B的4列分别一小条一小条轮着过，这个小C块从头到尾不落出寄存器，算术强度直接顶到理论上限。OpenBLAS这类库的核心，就是手写汇编或intrinsic精心排布的这种micro-kernel，外面再一层层套上缓存分块。

矩阵分块：大问题拆成小问题的艺术

这就是系列2提过一句的分块（blocking），这次把道理说透。分块矩阵乘法在数学上和普通矩阵乘法长得一模一样，只是"元素"换成了"子矩阵"：

$$ C_{ij} = \sum_{k=1}^{s} A_{ik} B_{kj} $$

妙处在于每个子矩阵都足够小，能整个蹲进缓存。$A_{ik}$ 从内存搬进缓存一次，要被行上所有的 $C_{ij}$ 反复使用；不分块时这种复用隔着几层循环，早被冲走了。块的大小按L1/L2/L3的容量像俄罗斯套娃一样逐层设计，就是BLAS库里实际的干活方式；还有更潇洒的cache-oblivious递归分块，不做任何参数调优，一路对半切到缓存装得下为止，任意层级的缓存都自然受益。

分块解决的是"搬数据"的问题，属于把 $O(n^3)$ 的常数做小。还有一路更狠的剑走偏锋：Strassen算法，直接把指数往下砍。1969年Strassen发现，两个2×2分块的矩阵相乘，常规要8次子矩阵乘法，通过7次巧妙的加减组合，7次子矩阵乘法就够：

$$ T(n) = 7\,T(n/2) + O(n^2) \; \Rightarrow \; T(n) = O(n^{\log_2 7}) \approx O(n^{2.807}) $$

每递归一层省1/8的乘法，n=2048时理论运算量能省两成多。听起来美好，实际工程里主流BLAS默认并不用它：加减的额外内存搬运、数值稳定性变差、还有和分块优化的纠缠，导致它只在矩阵特别大时才划算，属于知道有这回事、面试能加分、生产谨慎用的角色。

组合拳成绩单

把这一路优化叠起来，还是系列2那个N=2048的实验，账目大概是：

ijk朴素版(指令串行):        约 9 秒   —— 教科书基线
ikj重排(缓存友好):          约 1.5 秒 —— 访存局部性
+SIMD与寄存器分块(单线程):   约 0.4 秒 —— 指令并行, 吃掉理论峰值
OpenBLAS(分块+SIMD+多线程):  约 0.05 秒 —— 组合拳全开

每一步都不是玄学，分别对应：少搬数据、搬一次多用几回、一条指令算一批、多核一起上。所谓高性能计算，很大程度上就是把这几条朴素的道理一层层叠到极致。

BTW

验证自己代码有没有被向量化，gcc加 -fopt-info-vec 参数，编译时会打出哪些循环被向量化了，哪条循环没被向化和为什么，一目了然。另外现代编译器在 -O3 -march=native 下自动向化的能力已经不错，简单循环根本不用手写intrinsic，真正需要动手的场合，是像micro-kernel那种要把寄存器分配和指令排布都攥在自己手里的硬仗。

标签: none

添加新评论