为何GCC不对8×250×4矩阵乘法做向量化及寄存器累加优化？

在计算8×250矩阵与250×4矩阵相乘时，GCC生成的代码性能表现极差，以下是具体分析：

原始C实现 #

对应的矩阵乘法C代码如下：

#define MR 8
#define NR 4
#define KC 250

void matmul(double *a, double *b, double *c)
{
    for (int i = 0; i < MR; i++) {
        for (int j = 0; j < NR; j++) {
            for (int k = 0; k < KC; k++) {
                c[i * NR + j] += a[i * KC + k] * b[k * NR + j];
            }
        }
    }
}

GCC生成的热点循环汇编 #

GCC编译后，热点循环的汇编代码如下：

.L3:
        vmovsd  xmm0, QWORD PTR [rdx]
        add     rax, 32
        add     rdx, 8
        vmulsd  xmm0, xmm0, QWORD PTR [rax-32]
        vaddsd  xmm1, xmm1, xmm0
        vmovsd  QWORD PTR [rcx], xmm1
        cmp     rax, rsi
        jne     .L3

从这段汇编可以看出两个关键性能问题：

• 未针对j循环实现向量化：理论上j循环（NR=4）可以通过SIMD指令一次性处理4个double元素，但当前代码仅用标量指令vmulsd/vaddsd逐元素计算。
• 累加操作未在寄存器中完成：每次循环都将累加结果写回内存（vmovsd QWORD PTR [rcx], xmm1），频繁的内存读写大幅降低了性能。

理想的内循环汇编实现 #

更高效的手写内循环汇编如下，充分利用了AVX2指令集的广播和FMA（融合乘加）能力：

.L12:
        vmovupd ymm0, YMMWORD PTR [rsi]
        vbroadcastsd    ymm10, QWORD PTR [rdi]
        add     rdi, 8
        add     rsi, 32
        vbroadcastsd    ymm9, QWORD PTR [rdi+1992]
        vfmadd231pd     ymm8, ymm0, ymm10
        vbroadcastsd    ymm10, QWORD PTR [rdi+3992]
        vfmadd231pd     ymm7, ymm0, ymm9
        vbroadcastsd    ymm9, QWORD PTR [rdi+5992]
        vfmadd231pd     ymm6, ymm0, ymm10
        vbroadcastsd    ymm10, QWORD PTR [rdi+7992]
        vfmadd231pd     ymm5, ymm0, ymm9
        vbroadcastsd    ymm9, QWORD PTR [rdi+9992]
        vfmadd231pd     ymm4, ymm0, ymm10
        vbroadcastsd    ymm10, QWORD PTR [rdi+11992]
        vfmadd231pd     ymm3, ymm0, ymm9
        vbroadcastsd    ymm9, QWORD PTR [rdi+13992]
        vfmadd231pd     ymm2, ymm0, ymm10
        vfmadd231pd     ymm1, ymm0, ymm9
        cmp     rax, rdi
        jne     .L12

性能对比 #

在Zen2处理器上，手写汇编的性能比GCC生成的代码快5倍；Clang虽然对代码做了向量化处理，但生成的代码性能反而比GCC慢2倍。

GCC无法完成优化的原因 #

1. 内存访问模式识别不足：原代码中b[k * NR + j]的访问模式是跨 stride 的（NR=4，每个double占8字节，stride为32字节），GCC未能识别这种模式可以通过vbroadcastsd广播单个a元素，再结合FMA指令一次性完成4个j元素的乘加操作，尤其是KC=250这种非2的幂次，进一步限制了编译器的循环展开与向量化启发式触发。
2. 别名分析与寄存器提升限制：编译器无法确定c指针指向的内存与a/b是否存在别名，因此不敢将c[i * NR + j]的累加操作完全放在寄存器中完成，只能每次写回内存，导致额外的内存开销。
3. 循环顺序未自动调整：原代码的循环顺序为i→j→k，而更适合向量化的顺序是i→k→j（让j循环的连续访问适配SIMD宽度），GCC的循环变换优化未自动调整这一顺序以适配向量化需求。

内容的提问来源于stack exchange，提问作者asdfldsfdfjjf