CPU 矩阵乘法优化

1. 结论

2. 实现原理

2.1 矩阵乘法基础

矩阵乘法 C = A × B，其中 A[M×K], B[K×N], C[M×N]：

2.2 为什么循环顺序影响性能？

CPU 缓存是按行（连续内存）加载的。不同循环顺序导致不同的内存访问模式：

2.3 分块优化原理

将大矩阵分成小块，使每个块能完全放入 L1/L2 缓存：

2.4 OpenMP 并行化

使用 #pragma omp parallel for 将工作分配给多个线程：

2.5 OpenBLAS 为什么这么快？

3. 编译与运行

3.1 编译命令

3.2 运行方式

4. 核心源代码

4.1 朴素版本 (ijk)

void matmul_naive(const float* A, const float* B, float* C, int M, int K, int N) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];  // B 按列访问，缓存不友好
            }
            C[i * N + j] = sum;
        }
    }
}

4.2 循环重排版本 (ikj)

void matmul_ikj(const float* A, const float* B, float* C, int M, int K, int N) {
    memset(C, 0, M * N * sizeof(float));
    for (int i = 0; i < M; i++) {
        for (int k = 0; k < K; k++) {
            float a_ik = A[i * K + k];  // 提取到寄存器
            for (int j = 0; j < N; j++) {
                C[i * N + j] += a_ik * B[k * N + j];  // B 按行访问，缓存友好
            }
        }
    }
}

4.3 分块版本

#define BLOCK_SIZE 64  // 块大小，通常选择 L1 缓存大小的 1/3

void matmul_blocked(const float* A, const float* B, float* C, int M, int K, int N) {
    memset(C, 0, M * N * sizeof(float));

    for (int i0 = 0; i0 < M; i0 += BLOCK_SIZE) {
        for (int k0 = 0; k0 < K; k0 += BLOCK_SIZE) {
            for (int j0 = 0; j0 < N; j0 += BLOCK_SIZE) {
                // 块内乘法
                int i_end = (i0 + BLOCK_SIZE < M) ? i0 + BLOCK_SIZE : M;
                int k_end = (k0 + BLOCK_SIZE < K) ? k0 + BLOCK_SIZE : K;
                int j_end = (j0 + BLOCK_SIZE < N) ? j0 + BLOCK_SIZE : N;

                for (int i = i0; i < i_end; i++) {
                    for (int k = k0; k < k_end; k++) {
                        float a_ik = A[i * K + k];
                        for (int j = j0; j < j_end; j++) {
                            C[i * N + j] += a_ik * B[k * N + j];
                        }
                    }
                }
            }
        }
    }
}

4.4 OpenMP 并行 + 分块

#include <omp.h>

void matmul_omp(const float* A, const float* B, float* C, int M, int K, int N) {
    memset(C, 0, M * N * sizeof(float));

    #pragma omp parallel for schedule(dynamic)  // 动态调度，负载均衡
    for (int i0 = 0; i0 < M; i0 += BLOCK_SIZE) {
        for (int k0 = 0; k0 < K; k0 += BLOCK_SIZE) {
            for (int j0 = 0; j0 < N; j0 += BLOCK_SIZE) {
                // ... 块内计算同上
            }
        }
    }
}

4.5 OpenBLAS 版本

#include <cblas.h>

void matmul_openblas(const float* A, const float* B, float* C, int M, int K, int N) {
    // C = alpha * A * B + beta * C
    cblas_sgemm(CblasRowMajor,   // 行优先存储
                CblasNoTrans,    // A 不转置
                CblasNoTrans,    // B 不转置
                M, N, K,         // 矩阵维度
                1.0f,            // alpha = 1
                A, K,            // A 矩阵
                B, N,            // B 矩阵
                0.0f,            // beta = 0
                C, N);           // C 矩阵
}

5. 运行结果与分析

5.1 测试环境

5.2 性能测试结果

5.3 结果分析

单线程优化效果：

• 循环重排和分块性能相近 (~31 GFLOPs)，都实现了良好的缓存利用
• 对于 1024×1024 矩阵，分块优势不明显，更大矩阵差距会更大

OpenMP 简单并行为什么更慢？

OpenMP + 分块为什么有效？

• 每个线程处理独立的块，减少缓存竞争
• 分块保证了缓存局部性，多线程不会互相干扰
• 20核实现 2.4x 加速（未达到理想的 20x，因为内存带宽限制）

OpenBLAS 为什么快这么多？

• AVX-512: 一条指令处理 16 个 float（我们的代码每次只处理 1 个）
• FMA 指令: 融合乘加，一条指令完成 a*b+c
• 寄存器分块: 充分利用 32 个 AVX 寄存器
• 精确的缓存调优: 块大小针对 L1/L2/L3 精确优化

6. 更多加速技术

6.1 GPU 加速示例 (CUDA cuBLAS)

#include <cublas_v2.h>

// GPU 矩阵乘法
void matmul_gpu(const float* A, const float* B, float* C, int M, int K, int N) {
    cublasHandle_t handle;
    cublasCreate(&handle);

    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A, M * K * sizeof(float));
    cudaMalloc(&d_B, K * N * sizeof(float));
    cudaMalloc(&d_C, M * N * sizeof(float));

    cudaMemcpy(d_A, A, M * K * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, B, K * N * sizeof(float), cudaMemcpyHostToDevice);

    float alpha = 1.0f, beta = 0.0f;
    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N);

    cudaMemcpy(C, d_C, M * N * sizeof(float), cudaMemcpyDeviceToHost);

    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    cublasDestroy(handle);
}

6.2 性能层级总结

7. 各硬件平台编程生态

不同加速硬件有各自独立的编程接口，CUDA 仅适用于 NVIDIA GPU，不能用于 TPU 或 NPU：

测试环境: Ubuntu 22 / GCC / 20核 CPU

生成时间: 2026-01-29