如何优化矩阵乘法的简单for循环以实现最快执行速度

你的核心需求是对内存中连续存储的1960个8×8单精度复数矩阵，计算所有元素的模平方累加和，以下优化按收益从高到低落地，最终速度可以做到接近内存带宽上限：

• 第一优先级：修正内存访问顺序，消除跨步访问，最大化缓存利用率
  原代码最大的性能问题是内存访问顺序完全错误：固定外层i的前提下，循环顺序是j从0到7、k从0到7，对应索引i*64 +k*8 +j的访问步长是8个复数元素（64字节），属于典型的大跨步访问，CPU硬件预取完全失效，缓存行利用率仅1/8，这种写法的性能通常会比顺序访问低5~10倍。
  由于整个Q数组总大小约1MB，完全可以放进L1缓存，只需要交换j和k的循环顺序，就能实现完全顺序的内存访问；更进一步可以直接用指针步进遍历，完全消除循环内的乘法索引计算，三层嵌套循环可以直接简化为单循环遍历——因为模平方累加满足交换律，遍历顺序不影响计算结果。

• 第二优先级：循环展开+多路累加，打破指令依赖链
  原代码所有累加操作都依赖同一个sq变量，每一次加法都需要等待上一次计算完成，CPU的乱序执行流水线完全无法填满。你可以设置4~8个独立的累加变量，每次并行处理多个元素，最后再把所有累加变量求和，让CPU可以同时调度多条乘加指令，IPC（每周期指令数）可以拉到CPU理论上限。
  由于内层8×8的循环尺寸是编译期固定常量，编译器也会自动配合做展开优化。

• 第三优先级：对齐内存+开启编译器向量化优化
  编译时开启-O3 -march=native（GCC/Clang）或者/O2 /arch:AVX2（MSVC），让编译器自动生成SSE/AVX/AVX-512向量指令，单条指令即可完成4~16组单精度浮点数的乘加运算。
  内存分配优先使用对齐接口分配32/64字节对齐的内存，避免非对齐向量访问的额外惩罚。

优化后参考代码 #

#include <stdlib.h>

typedef struct
{
  float real;
  float img;
} my_complex;

// 分配32字节对齐内存，适配AVX向量指令要求
my_complex* Q = (my_complex*)aligned_alloc(32, 8 * 8 * 1960 * sizeof(my_complex));

float calc_sum_sq(const my_complex* Q, size_t total_cnt)
{
    // 4路独立累加，打破指令依赖链
    float sq0 = 0.0f, sq1 = 0.0f, sq2 = 0.0f, sq3 = 0.0f;
    const my_complex* p = Q;
    size_t i = 0;

    // 每次处理4个复数，循环步长为4
    for (; i + 3 < total_cnt; i += 4) {
        sq0 += p[0].real * p[0].real + p[0].img * p[0].img;
        sq1 += p[1].real * p[1].real + p[1].img * p[1].img;
        sq2 += p[2].real * p[2].real + p[2].img * p[2].img;
        sq3 += p[3].real * p[3].real + p[3].img * p[3].img;
        p += 4;
    }

    // 处理剩余不足4个的尾元素
    for (; i < total_cnt; i++) {
        sq0 += p->real * p->real + p->img * p->img;
        p++;
    }

    return sq0 + sq1 + sq2 + sq3;
}

// 调用：总元素数固定为1960*8*8=125440
float sq = calc_sum_sq(Q, 125440);

注：如果你不需要手动写向量intrinsics指令，只需要把循环改成顺序单循环、开启O3和对应架构优化选项，编译器自动生成的向量代码性能和手动优化版本几乎没有差距。

额外可选优化 #

如果对精度要求不极端，不要用双精度做累加，单精度累加在这个数据规模下误差可以忽略，且向量指令处理单精度的吞吐量是双精度的2倍。
不要在循环内做任何分支判断，保证循环是单一的直出控制流，避免分支预测失败的开销。

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