程序性能优化之矩阵乘法(1)

一：实现简单矩阵乘法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

// C=A*B, A=M*N, B=N*P  ==> C=M*P
void matmul_naive(int M, int N, int P,
		              double *A, double *B, double *C) {
// 初始化C为0
  for(int i = 0; i < M*P; i++)
    C[i] = 0;
  for(int i = 0; i < M; i++)
    for (int j = 0; j < P; j++)
      for (int k = 0; k < N; k++)
        C[i * P + j] = A[i * N + k] * B[k * P + j];
}

int main() {
  int M = 512, N = 512, P = 512;
  double *A = malloc(M * N * sizeof(double));
  double *B = malloc(N * P * sizeof(double));
  double *C = malloc(M * P * sizeof(double));

  // 随机初始化
  for (int i = 0; i < M * N; i++)
    A[i] = (double)rand() / RAND_MAX;
  for (int i = 0; i < N * P; i++)
    B[i] = (double)rand() / RAND_MAX;
  clock_t start = clock();
  matmul_naive(M, N, P, A, B, C);
  printf("Naive: %.3f s\n", (double)(clock() - start) / CLOCKS_PER_SEC);
  free(A);
  free(B);
  free(C);
  return 0;
}

问题:行优先（Row-Major）vs 列优先（Column-Major）

1
2
3
C[i][j] = C[i * P + j]   // 行优先：第i行第j列
             ↑      ↑
         行偏移   列偏移
方式 公式 含义 使用场景
行优先 i * 列数 + j 一行一行存 C/C++/Python/NumPy默认
列优先 j * 行数 + i 一列一列存 Fortran/MATLAB默认

直观图解

1
2
3
4
5
6
7
8
9
10
11
12
13
内存地址增长方向 → → → → → → → → → → → →

行优先存储（C语言）：
[0][0] [0][1] [0][2] [0][3] | [1][0] [1][1] [1][2] [1][3] | [2][0] [2][1] [2][2] 
  0      1      2      3        4      5      6      7        8      9     10     
[2][3]
  11

列优先存储（Fortran）：
[0][0] [1][0] [2][0] | [0][1] [1][1] [2][1] | [0][2] [1][2] [2][2] | [0][3] [1][3] 
  0      1      2        3      4      5        6      7      8        9     10   
[2][3]
  11

为什么c语言行优先

代码直观性

1
2
3
4
5
6
7
8
// 二维数组在C中的声明
double C[3][4];  // 3行4列

// 访问方式完全匹配直觉
C[i][j]  // 第i行，第j列

// 编译器实际转换的公式（行优先）
*((double*)C + i * 4 + j)

关键观察：行的元素在内存中连续

1
2
3
4
5
6
7
// 遍历一行（快！缓存友好）
for (int j = 0; j < 4; j++)
    C[i][j] = ...;  // 访问 0,1,2,3 连续地址

// 遍历一列（慢！缓存不友好）  
for (int i = 0; i < 3; i++)
    C[i][j] = ...;  // 访问 0,4,8 跳跃地址

二： 理解性能瓶颈—缓存

内存层次(速度差异巨大)

层级 访问时间 容量
CPU寄存器 ~0.3 ns 几十个
L1缓存 ~1 ns 32-64 KB
L2缓存 ~4 ns 256-512 KB
L3缓存 ~10 ns 几-几十 MB
内存(DRAM) ~100 ns 几-几十 GB

缓存行（Cache Line）

CPU一次从内存取64字节（8个double）到缓存。

关键问题：如果访问的内存不连续，每次都要重新从内存加载！

三：优化1—循环交换

问题分析：

1
2
3
4
// 简单版本： 内层循环k, 访问B[k*P+j]
// k增加1,B的访问跳跃P个元素(不连续的！)
C[i*P+j] += A[i*N+k] * B[k*p+j];
//          ^连续         ^跳跃

优化：让j在最内层（连续访问）

1
2
3
4
5
6
7
8
void matmul_opt1(int M, int N, int P, double *A, double *B, double *C) {
  for (int i = 0; i < M * P; i++)
    C[i] = 0;
  for (int i = 0; i < M; i++)
    for (int k = 0; k < N; k++)
      for (int j = 0; j < P; j++)
        C[i * P + j] += A[i * N + k] * B[k * P + j];
}

四： 优化2—分块(Tiling/Blocking)

核心思想

缓存装不下整个矩阵，但装得下一个小块。

把大矩阵乘法拆成小矩阵乘法，让小矩阵能留在缓存中复用。

1
2
3
4
5
6
7
8
9
10
┌─────────────┐    ┌─────────────┐
│ A11 │ A12   │    │ B11 │ B12   │
├─────┼───────┤ ×  ├─────┼───────┤
│ A21 │ A22   │    │ B21 │ B22   │
└─────────────┘    └─────────────┘

= 计算 2×2 = 4 个小矩阵乘法
  C11 = A11×B11 + A12×B21
  C12 = A11×B12 + A12×B22
  ...
1

五： 优化3—SIMD向量化

六：优化4—分块 + SIMD + 展开

七： 参考OpenBLAS工业级实现