如何在M1芯片上优化C语言1024x1024矩阵乘法以匹配NumPy性能

我正在优化C语言实现的1024x1024矩阵乘法，目标是匹配NumPy在M1芯片上单线程约90 GFlop/s的性能。目前我的C实现仅能达到约55 GFlop/s，已采用ARM NEON intrinsics进行向量化，并使用多种分块大小提升缓存利用率。

以下是我的C代码：

#define N 1024
#define BLOCK_X 16 
#define BLOCK_K 8 
#define BLOCK_Y 8 
void matmul(float* const A, float* const B, float* const C){
    for(int by = 0; by < N; by+=BLOCK_Y) {
        for(int bx = 0; bx < N; bx+=BLOCK_X) {
            for(int bk = 0; bk < N; bk+=BLOCK_K) {
                for (int y=by; y<by+BLOCK_Y; y++){
                    for (int x=bx; x<bx+BLOCK_X; x+=4){
                        float32x4_t sum_vec = vld1q_f32(&C[y * N + x]);
                        __builtin_prefetch(&C[y * N + x + 4], 1, 0);
                        for (int k=bk; k<bk+BLOCK_K; k++){
                            float32x4_t a_vec = vdupq_n_f32(A[y * N + k]);
                            float32x4_t b_vec = vld1q_f32(&B[k*N+x]);
                            sum_vec = vmlaq_f32(sum_vec, a_vec, b_vec);
                        }
                        vst1q_f32(&C[y * N + x], sum_vec);
                    }
                }
            }
        }
    }
}

我原本预期性能更接近NumPy，但目前差距较大，请问当前实现中哪些瓶颈或问题限制了性能？

• B矩阵的非连续内存访问：当前对B矩阵的访问是B[k*N+x]，属于按列访问（stride=N），M1的缓存系统对连续内存访问的带宽利用率远高于非连续访问，这种模式会导致大量缓存miss，严重拖慢性能。NumPy通常会在乘法前转置B矩阵，将列访问转为连续的行访问，最大化缓存命中率。
• 分块策略与缓存匹配度不足：M1核心的L1数据缓存为32KB，当前分块大小（BLOCK_Y=8、BLOCK_X=16、BLOCK_K=8）的内存总占用远低于L1缓存容量，没有充分利用缓存带宽。同时循环顺序（by -> bx -> bk）没有最大化数据复用，应该调整为bk -> by -> bx或类似顺序，让A和B的分块数据在缓存中停留更久，减少重复加载。
• NEON指令并行度未充分发挥：内循环中每次仅处理A的单个元素（通过vdupq_n_f32广播），与B的4个元素做乘加。M1的NEON支持同时处理更多向量操作，比如将BLOCK_K调整为16，一次加载4个A元素组成向量，配合B的4x4连续块进行批量乘加，能大幅提升指令吞吐量。
• 预取策略错误：当前预取的是C矩阵的后续地址，但性能瓶颈主要在A、B的读访问，而非C的写回。应该针对后续要访问的A、B分块进行预取，比如在处理当前bk块时，预取下一个bk块的A和B数据，提前将数据加载到缓存中。
• 循环展开不足：内循环k仅循环8次，未手动展开会保留循环控制开销，且无法充分利用M1的指令流水线并行能力。手动展开该循环可减少分支指令，让编译器生成更紧凑的并行代码。
• C矩阵的冗余内存操作：每次循环都从C加载初始值sum_vec，若C初始化为全零，可直接将sum_vec初始化为零向量（vdupq_n_f32(0.0f)），省去加载步骤，减少一次内存读操作。

内容的提问来源于stack exchange，提问作者Steven Daniel Anderson