如何让矩阵乘法快10倍？一个真实高性能计算优化案例

你有没有遇到过这样的场景：训练一个深度学习模型，光是前向传播就卡了几十秒；做一次图像卷积，等结果等到泡了三杯咖啡；跑个科学模拟，一晚上都算不完？

背后元凶之一，可能就是那个看似简单的操作——矩阵乘法。

别小看它。在 $ N=2048 $ 的规模下，两个方阵相乘需要超过80亿次浮点运算。如果用最朴素的三重循环串行实现，哪怕你的CPU主频高达3GHz，也得跑好几秒钟。而这还只是单次调用。

怎么破局？答案是：并行 + 缓存优化 + 向量化三管齐下。

今天我们就来拆解一个真实的高性能计算优化案例，带你一步步把一个“教科书式”的慢速矩阵乘法，变成接近硬件极限的高速引擎。这不是理论推演，而是你在OpenBLAS、Intel MKL这些工业级库中每天都在用的技术组合。

从零开始：先写个“正确但很慢”的版本

我们先从最基础的串行实现出发：

void matmul_basic(double A[N][N], double B[N][N], double C[N][N]) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double sum = 0.0; for (int k = 0; k < N; k++) { sum += A[i][k] * B[k][j]; } C[i][j] = sum; } } }

逻辑清晰，代码简洁，考试满分。但在性能上……几乎是“灾难”。

为什么这么慢？三个关键问题：
1.只用了单核，现代CPU动辄8核16线程，白白浪费；
2.内存访问不友好，对B[k][j]是列访问，缓存命中率极低；
3.没有利用SIMD指令，每个周期只能算一次乘加。

这就像开着拖拉机跑F1赛道——车没问题，路也没问题，但你没踩油门。

第一步加速：多核并行 —— 让所有核心动起来

既然有多核，那就别闲着。我们可以把外层循环i拆开，每个线程处理一部分行。

借助 OpenMP，一行预处理指令就能完成：

#include <omp.h> void matmul_parallel(double A[N][N], double B[N][N], double C[N][N]) { #pragma omp parallel for collapse(2) for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double sum = 0.0; for (int k = 0; k < N; k++) { sum += A[i][k] * B[k][j]; } C[i][j] = sum; } } }

加上collapse(2)是为了让双重循环被整体调度，避免负载不均。

编译时打开优化：

gcc -O3 -fopenmp -march=native matrix_mult.c -o matmul

效果如何？在一台16核服务器上，$ N=1024 $ 时，速度直接提升6~8倍。听起来不错？其实还有很大空间。

因为你会发现，随着核心数增加，加速比不再线性上升——瓶颈已经转移到内存子系统了。

第二步突破：分块优化（Tiling）—— 把数据“搬进”缓存

现在的问题是：虽然我们并行了，但每个线程还是频繁地从主存读取数据，而现代CPU的缓存带宽比主存高一个数量级。

举个例子：当你访问B[k][j]时，如果j固定、k变化，相当于按列访问二维数组。而C语言中数组是行优先存储的，这意味着每次访问都不是连续内存，导致缓存行利用率极低。

解决方案是什么？把大矩阵切成小块，一块一块地算。

这就是所谓的分块矩阵乘法（Blocked Matrix Multiplication 或 Tiling）。

分块的核心思想

我们将矩阵划分为若干 $ B_s \times B_s $ 的小块（tile），使得每个块能完整放入L1缓存。然后按块遍历：

for ii ← 0 to N step Bs for jj ← 0 to N step Bs for kk ← 0 to N step Bs // 计算 C[ii:ii+Bs, jj:jj+Bs] += A[ii:ii+Bs, kk:kk+Bs] × B[kk:kk+Bs, jj:jj+Bs]

这样，在内层计算中，A、B、C的子块都能被重复使用，大大提升数据局部性。

实际代码实现

#define BLOCK_SIZE 64 void matmul_tiled(double A[N][N], double B[N][N], double C[N][N]) { for (int ii = 0; ii < N; ii += BLOCK_SIZE) for (int jj = 0; jj < N; jj += BLOCK_SIZE) for (int kk = 0; kk < N; kk += BLOCK_SIZE) // 内部小块乘加 for (int i = ii; i < ii + BLOCK_SIZE && i < N; i++) for (int j = jj; j < jj + BLOCK_SIZE && j < N; j++) { double temp = 0.0; for (int k = kk; k < kk + BLOCK_SIZE && k < N; k++) { temp += A[i][k] * B[k][j]; } C[i][j] += temp; } }

注意：这里初始C[i][j]应为0，或改用累加模式。

结合 OpenMP 并行最外层两个块循环：

#pragma omp parallel for collapse(2) for (int ii = 0; ii < N; ii += BLOCK_SIZE) for (int jj = 0; jj < N; jj += BLOCK_SIZE) ...

性能提升有多大？

实测表明，在 $ N=1024 $ 场景下，仅靠分块即可再提速2~4倍，尤其在非NUMA均衡架构上更为显著。缓存命中率从不足40%提升至85%以上。

第三步压榨：SIMD向量化 —— 单指令多数据流

到现在为止，我们已经解决了“并行”和“缓存”两大难题。接下来要挑战的是指令级并行。

现代x86 CPU支持 AVX/AVX2 指令集，可以一次性处理4个双精度浮点数（256位）。如果你能让编译器生成这些指令，就能实现“一拍四算”。

可惜，不是所有循环都能自动向量化。比如原始的内层点积：

for (k = 0; k < N; k++) sum += A[i][k] * B[k][j];

由于B[k][j]是跨步访问，编译器通常不敢向量化。

但我们可以在分块的基础上，对内部小块启用显式向量操作。

使用内在函数（Intrinsics）手动向量化

#include <immintrin.h> void block_multiply_vectorized(double *a_block, double *b_block, double *c_block, int bs) { for (int i = 0; i < bs; i++) { for (int j = 0; j < bs; j += 4) { __m256d c_vec = _mm256_loadu_pd(&c_block[i*bs + j]); __m256d b_col0 = _mm256_loadu_pd(&b_block[0*bs + j]); __m256d b_col1 = _mm256_loadu_pd(&b_block[1*bs + j]); __m256d b_col2 = _mm256_loadu_pd(&b_block[2*bs + j]); __m256d b_col3 = _mm256_loadu_pd(&b_block[3*bs + j]); for (int k = 0; k < bs; k++) { __m256d a_val = _mm256_set1_pd(a_block[i*bs + k]); __m256d b_vals = _mm256_set_pd( b_block[k*bs + j+3], b_block[k*bs + j+2], b_block[k*bs + j+1], b_block[k*bs + j+0] ); c_vec = _mm256_add_pd(c_vec, _mm256_mul_pd(a_val, b_vals)); } _mm256_storeu_pd(&c_block[i*bs + j], c_vec); } } }

当然，上面只是示意。实际更高效的做法是采用GEMM 分块算法 + Register Blocking + SIMD + 多线程的完整链条。

不过好消息是：你不需要自己写这么多底层代码。

工业级方案参考：为什么 BLAS 库这么快？

像 Intel MKL、OpenBLAS、BLIS 这些库之所以能做到极致性能，正是融合了上述所有技术：

它们甚至会根据 CPU 型号自动选择最优块大小和线程策略。

所以当你调用cblas_dgemm时，背后是一整套精密协作的高性能引擎在工作。

性能对比：优化前后差距有多大？

我们在一台 Intel Xeon Gold 6230（20核40线程）上测试 $ N=1024 $ 的双精度矩阵乘法：

看到没？最终性能相差三十多倍。这还不包括更高级的流水线重叠、预取优化等技巧。

调优实战建议：五个必须知道的坑

1. 块大小不是越大越好
   - 太大会溢出L1缓存，太小则开销占比高。
   - 推荐范围：32~64，可通过实验绘制性能曲线确定最优值。

2. 内存对齐很重要
   c double *A = (double*)aligned_alloc(32, sizeof(double)*N*N);
   使用32字节对齐有助于AVX加载，避免性能降级。

3. 别盲目开启超线程
   - 矩阵乘法是计算密集型任务，通常设为物理核心数即可。
   - 可通过OMP_NUM_THREADS=20控制。

4. 编译选项决定下限
   必须启用：
   bash -O3 -march=native -ffast-math -funroll-loops

5. NUMA系统要小心
   在双路服务器上，若内存绑定不当，远程访问延迟可达本地2倍。
   运行时使用：
   bash numactl --interleave=all ./matmul

结语：掌握这套方法，你能优化的不只是矩阵乘法

今天我们走完了从“教科书代码”到“接近极限性能”的全过程。总结一下关键技术栈：

并行化 × 数据局部性 × 向量化 = 高性能计算三大支柱

这套方法不仅适用于矩阵乘法，还能迁移到：
- 卷积神经网络中的 im2col + GEMM
- FFT 中的蝶形运算并行
- 稀疏矩阵与向量乘法（SpMV）
- 动态规划类算法（如序列比对）

更重要的是，它教会我们一种思维方式：不要只关注算法复杂度，更要关心数据如何流动、指令如何执行、缓存如何工作。

下次当你觉得“程序太慢”的时候，不妨问自己三个问题：
1. 我的代码用满所有核心了吗？
2. 数据是不是一直在“长途跋涉”访问内存？
3. CPU的SIMD单元是不是在“摸鱼”？

只要答好这三个问题，你就离写出真正高效的代码不远了。

如果你正在实现自己的数值计算库，或者想深入理解BLAS背后的原理，欢迎留言交流。也可以分享你在项目中做过哪些令人印象深刻的性能优化！