如何进一步优化矩阵乘法的性能？
How to further optimize performance of Matrix Multiplication?
我正在尝试在单核上优化我的矩阵乘法代码 运行。我怎样才能进一步提高循环展开方面的性能，FMA/SSE？我也很想知道如果在内循环中使用四个而不是两个总和，为什么性能不会提高。
问题大小是 1000x1000 矩阵乘法。 gcc 9 和 icc 19.0.5 都可用。 Intel Xeon @ 3.10GHz，32K L1d 缓存，Skylake 架构。使用 gcc -O3 -mavx.
编译
void mmult(double* A, double* B, double* C)
{
const int block_size = 64 / sizeof(double);
__m256d sum[2], broadcast;
for (int i0 = 0; i0 < SIZE_M; i0 += block_size) {
for (int k0 = 0; k0 < SIZE_N; k0 += block_size) {
for (int j0 = 0; j0 < SIZE_K; j0 += block_size) {
int imax = i0 + block_size > SIZE_M ? SIZE_M : i0 + block_size;
int kmax = k0 + block_size > SIZE_N ? SIZE_N : k0 + block_size;
int jmax = j0 + block_size > SIZE_K ? SIZE_K : j0 + block_size;
for (int i1 = i0; i1 < imax; i1++) {
for (int k1 = k0; k1 < kmax; k1++) {
broadcast = _mm256_broadcast_sd(A+i1*SIZE_N+k1);
for (int j1 = j0; j1 < jmax; j1+=8) {
sum[0] = _mm256_load_pd(C+i1*SIZE_K+j1+0);
sum[0] = _mm256_add_pd(sum[0], _mm256_mul_pd(broadcast, _mm256_load_pd(B+k1*SIZE_K+j1+0)));
_mm256_store_pd(C+i1*SIZE_K+j1+0, sum[0]);
sum[1] = _mm256_load_pd(C+i1*SIZE_K+j1+4);
sum[1] = _mm256_add_pd(sum[1], _mm256_mul_pd(broadcast, _mm256_load_pd(B+k1*SIZE_K+j1+4)));
_mm256_store_pd(C+i1*SIZE_K+j1+4, sum[1]);
// doesn't improve performance.. why?
// sum[2] = _mm256_load_pd(C+i1*SIZE_K+j1+8);
// sum[2] = _mm256_add_pd(sum[2], _mm256_mul_pd(broadcast, _mm256_load_pd(B+k1*SIZE_K+j1+8)));
// _mm256_store_pd(C+i1*SIZE_K+j1+8, sum[2]);
// sum[3] = _mm256_load_pd(C+i1*SIZE_K+j1+12);
// sum[3] = _mm256_add_pd(sum[3], _mm256_mul_pd(broadcast, _mm256_load_pd(B+k1*SIZE_K+j1+12)));
// _mm256_store_pd(C+i1*SIZE_K+j1+4, sum[3]);
}
}
}
}
}
}
}
此代码每个 FMA 有 2 个负载（如果发生 FMA 收缩），但 Skylake 理论上只支持 最多 每个 FMA 一个负载（如果你想最大化 2 /clock FMA 吞吐量），即使这样在实践中通常也太多了。 （峰值是每个时钟 2 个负载 + 1 个存储，但它通常不能完全维持）。请参阅 Intel 的优化指南和 https://agner.org/optimize/
循环开销不是最大的问题，正文本身迫使代码以一半的速度运行。
如果 k-loop 是内部循环，则可以 链接 大量积累，而不必 load/store 来回 C。这有一个缺点：对于像这样的循环携带的依赖链，需要明确地编码确保有足够的独立工作要做。
为了负载少但有足够的独立工作，内部循环的主体可以计算来自A的小列向量和来自B的小行向量之间的乘积，对于使用 4 个标量广播加载列和从 B 加载 2 个法向向量的示例，导致 8 个独立 FMA 仅加载 6 次（甚至可能更低的比率），这足以让 Skylake 保持满意的独立 FMA 而不是太多许多负载。即使是 3x4 足迹也是可能的，它也有足够的独立 FMA 来让 Haswell 满意（它至少需要 10 个）。
我正好有一些示例代码，它是针对单精度和 C++ 的，但你会明白要点的：
sumA_1 = _mm256_load_ps(&result[i * N + j]);
sumB_1 = _mm256_load_ps(&result[i * N + j + 8]);
sumA_2 = _mm256_load_ps(&result[(i + 1) * N + j]);
sumB_2 = _mm256_load_ps(&result[(i + 1) * N + j + 8]);
sumA_3 = _mm256_load_ps(&result[(i + 2) * N + j]);
sumB_3 = _mm256_load_ps(&result[(i + 2) * N + j + 8]);
sumA_4 = _mm256_load_ps(&result[(i + 3) * N + j]);
sumB_4 = _mm256_load_ps(&result[(i + 3) * N + j + 8]);
for (size_t k = kk; k < kk + akb; k++) {
auto bc_mat1_1 = _mm256_set1_ps(*mat1ptr);
auto vecA_mat2 = _mm256_load_ps(mat2 + m2idx);
auto vecB_mat2 = _mm256_load_ps(mat2 + m2idx + 8);
sumA_1 = _mm256_fmadd_ps(bc_mat1_1, vecA_mat2, sumA_1);
sumB_1 = _mm256_fmadd_ps(bc_mat1_1, vecB_mat2, sumB_1);
auto bc_mat1_2 = _mm256_set1_ps(mat1ptr[N]);
sumA_2 = _mm256_fmadd_ps(bc_mat1_2, vecA_mat2, sumA_2);
sumB_2 = _mm256_fmadd_ps(bc_mat1_2, vecB_mat2, sumB_2);
auto bc_mat1_3 = _mm256_set1_ps(mat1ptr[N * 2]);
sumA_3 = _mm256_fmadd_ps(bc_mat1_3, vecA_mat2, sumA_3);
sumB_3 = _mm256_fmadd_ps(bc_mat1_3, vecB_mat2, sumB_3);
auto bc_mat1_4 = _mm256_set1_ps(mat1ptr[N * 3]);
sumA_4 = _mm256_fmadd_ps(bc_mat1_4, vecA_mat2, sumA_4);
sumB_4 = _mm256_fmadd_ps(bc_mat1_4, vecB_mat2, sumB_4);
m2idx += 16;
mat1ptr++;
}
_mm256_store_ps(&result[i * N + j], sumA_1);
_mm256_store_ps(&result[i * N + j + 8], sumB_1);
_mm256_store_ps(&result[(i + 1) * N + j], sumA_2);
_mm256_store_ps(&result[(i + 1) * N + j + 8], sumB_2);
_mm256_store_ps(&result[(i + 2) * N + j], sumA_3);
_mm256_store_ps(&result[(i + 2) * N + j + 8], sumB_3);
_mm256_store_ps(&result[(i + 3) * N + j], sumA_4);
_mm256_store_ps(&result[(i + 3) * N + j + 8], sumB_4);
这意味着 j 循环和 i 循环被展开，但 k 循环没有展开，即使它现在是内部循环。稍微展开 k 循环确实对我的实验有所帮助。
请参阅@harold 的回答以获得实际改进。这主要是为了转载我在评论里写的东西
four instead of two sums in the inner loop. (Why doesn't unrolling help?)
sum[i] 没有循环携带的依赖。下一次迭代分配 sum[0] = _mm256_load_pd(C+i1*SIZE_K+j1+0);，它不依赖于先前的值。
因此，将同一架构寄存器重命名为不同的 物理 寄存器足以避免可能使流水线停止的写后写危险。无需使用多个 tmp 变量使源复杂化。请参阅 （在那个问题中，2 个数组的一个点积， 是 一个循环通过累加器进行依赖。在那里，使用多个累加器 是 对隐藏 FP FMA 延迟很有价值，因此我们的瓶颈是 FMA 吞吐量，而不是延迟。）
没有 寄存器重命名的管道 （大多数有序 CPU）将从 "software pipelining" 中受益，以静态调度乱序 exec 可以在fly：加载到不同的寄存器中，因此每个负载和消耗它的 FMA 之间存在距离（充满独立工作）。然后在那和商店之间。
但是所有现代 x86 CPU 都是 OoO；甚至 Knight's Landing 对于 SIMD 向量也有一些有限的 OoO exec。 （Silvermont 不支持 AVX，但按顺序 运行 SIMD 指令，仅对整数执行 OoO）。
在没有任何多累加器情况来隐藏延迟的情况下，展开的好处（明确地在源代码中或使用 -fprofile-use 启用的 -funroll-loop，或默认情况下在 clang 中）是：
- 减少前端带宽以使循环开销进入后端。每个循环开销更有用的工作微指令。因此，如果您的 "useful work" 在前端接近瓶颈，它会有所帮助。
- 减少 运行 循环开销的后端执行单元需求。通常在 Haswell 及更高版本或 Zen 上不是问题；当指令组合包含一些整数内容和一些纯加载指令时，后端基本上可以跟上前端。
- 每次完成的工作总 uops 更少意味着 OoO exec 可以 "see" 更远的内存 loads/stores。
- 有时对短 运行ning 循环的分支预测更好：迭代次数越少意味着分支预测学习的模式越短。因此，对于短行程计数，当执行脱离循环时，有更好的机会正确预测最后一次迭代未采用的次数。
- 有时在更复杂的情况下保存
mov reg,reg，因为编译器更容易在不同的 reg 中生成新结果。同一个变量可以在两个 regs 之间交替，而不是需要 moved 回到同一个 regs 为下一次迭代做好准备。特别是如果你有一个以依赖的方式使用 a[i] 和 a[i+1] 的循环，或者类似 Fibonacci.
循环中有 2 个加载 + 1 个存储，这可能是瓶颈，而不是 FMA 或前端带宽。按 2 展开可能有助于避免前端瓶颈，但更重要的是只有来自另一个超线程的争用才重要。
评论中提出了一个有趣的问题：展开不需要很多寄存器才有用吗？
哈罗德评论道：
16 is not a huge number of registers, but it's enough to have 12
accumulators and 3 pieces of row vector from B and the broadcasted
scalar from A, so it works out to just about enough. The loop from OP
above barely uses any registers anyway. The 8 registers in 32bit are
indeed too few.
当然，由于问题中的代码在循环迭代中的寄存器中没有 "accumulators"，仅添加到内存中，编译器可以优化所有 sum[0..n] 以重用相同的寄存器汇编；存储后是"dead"。所以实际套印压力很低。
是的，x86-64 有点寄存器不足，这就是为什么 AVX512 将矢量寄存器 (zmm0..31) 的数量和宽度加倍的原因。是的，许多 RISC 有 32 int / 32 fp regs，包括 AArch64，ARM 有 16 个。
x86-64有16个标量整型寄存器（包括栈指针，不包括程序计数器），所以普通函数可以用15个。还有16个向量寄存器，xmm0..15。 （对于 AVX，它们的宽度是 ymm0..15 的两倍）。
（其中一些是在我注意到 sum[0..n] 毫无意义，不是循环携带之前写的。）
在这种情况下，将寄存器重命名到一个大的物理寄存器文件就足够了。在其他情况下，拥有更多架构寄存器会有所帮助，尤其是对于更高的 FP 延迟，因此 AVX512 具有 32 个 zmm 寄存器。但是对于整数 16 已经足够了。 RISC CPU 通常设计为有序，无需重新命名，需要 SW 管道。
对于 OoO exec，就减少 spill/reloads 而言，从 8 到 16 个架构 GP 整数 reg 的跳跃比从 16 到 32 的跳跃更重要。 （我看过一篇论文，用不同数量的架构寄存器测量了 SPECint 的总动态指令数。我没有再次查找，但 8->16 可能总节省了 10%，而 16->32 只是几个 %。类似的东西）。
但是这个具体问题不需要很多 FP 寄存器，sum[0..3] 只需要 4 个向量（如果它们是循环携带的），也许 1 个是临时的； x86 可以使用内存源 mul/add/FMA。寄存器重命名消除了任何 WAW 危害，因此我们可以重复使用相同的临时寄存器而不需要软件流水线。 （并且 OoO exec 还隐藏了负载和 ALU 延迟。）
当存在循环携带的依赖项时，您需要多个累加器。这段代码是添加到内存中，而不是添加到一些向量累加器中，所以任何依赖都是通过 store/reload。但这只有大约 7 个周期的延迟，所以任何理智的缓存阻塞因素都会隐藏它。
我正在尝试在单核上优化我的矩阵乘法代码 运行。我怎样才能进一步提高循环展开方面的性能，FMA/SSE？我也很想知道如果在内循环中使用四个而不是两个总和，为什么性能不会提高。
问题大小是 1000x1000 矩阵乘法。 gcc 9 和 icc 19.0.5 都可用。 Intel Xeon @ 3.10GHz，32K L1d 缓存，Skylake 架构。使用 gcc -O3 -mavx.
void mmult(double* A, double* B, double* C)
{
const int block_size = 64 / sizeof(double);
__m256d sum[2], broadcast;
for (int i0 = 0; i0 < SIZE_M; i0 += block_size) {
for (int k0 = 0; k0 < SIZE_N; k0 += block_size) {
for (int j0 = 0; j0 < SIZE_K; j0 += block_size) {
int imax = i0 + block_size > SIZE_M ? SIZE_M : i0 + block_size;
int kmax = k0 + block_size > SIZE_N ? SIZE_N : k0 + block_size;
int jmax = j0 + block_size > SIZE_K ? SIZE_K : j0 + block_size;
for (int i1 = i0; i1 < imax; i1++) {
for (int k1 = k0; k1 < kmax; k1++) {
broadcast = _mm256_broadcast_sd(A+i1*SIZE_N+k1);
for (int j1 = j0; j1 < jmax; j1+=8) {
sum[0] = _mm256_load_pd(C+i1*SIZE_K+j1+0);
sum[0] = _mm256_add_pd(sum[0], _mm256_mul_pd(broadcast, _mm256_load_pd(B+k1*SIZE_K+j1+0)));
_mm256_store_pd(C+i1*SIZE_K+j1+0, sum[0]);
sum[1] = _mm256_load_pd(C+i1*SIZE_K+j1+4);
sum[1] = _mm256_add_pd(sum[1], _mm256_mul_pd(broadcast, _mm256_load_pd(B+k1*SIZE_K+j1+4)));
_mm256_store_pd(C+i1*SIZE_K+j1+4, sum[1]);
// doesn't improve performance.. why?
// sum[2] = _mm256_load_pd(C+i1*SIZE_K+j1+8);
// sum[2] = _mm256_add_pd(sum[2], _mm256_mul_pd(broadcast, _mm256_load_pd(B+k1*SIZE_K+j1+8)));
// _mm256_store_pd(C+i1*SIZE_K+j1+8, sum[2]);
// sum[3] = _mm256_load_pd(C+i1*SIZE_K+j1+12);
// sum[3] = _mm256_add_pd(sum[3], _mm256_mul_pd(broadcast, _mm256_load_pd(B+k1*SIZE_K+j1+12)));
// _mm256_store_pd(C+i1*SIZE_K+j1+4, sum[3]);
}
}
}
}
}
}
}
此代码每个 FMA 有 2 个负载（如果发生 FMA 收缩），但 Skylake 理论上只支持 最多 每个 FMA 一个负载（如果你想最大化 2 /clock FMA 吞吐量），即使这样在实践中通常也太多了。 （峰值是每个时钟 2 个负载 + 1 个存储，但它通常不能完全维持）。请参阅 Intel 的优化指南和 https://agner.org/optimize/
循环开销不是最大的问题，正文本身迫使代码以一半的速度运行。
如果 k-loop 是内部循环，则可以 链接 大量积累，而不必 load/store 来回 C。这有一个缺点：对于像这样的循环携带的依赖链，需要明确地编码确保有足够的独立工作要做。
为了负载少但有足够的独立工作，内部循环的主体可以计算来自A的小列向量和来自B的小行向量之间的乘积，对于使用 4 个标量广播加载列和从 B 加载 2 个法向向量的示例，导致 8 个独立 FMA 仅加载 6 次（甚至可能更低的比率），这足以让 Skylake 保持满意的独立 FMA 而不是太多许多负载。即使是 3x4 足迹也是可能的，它也有足够的独立 FMA 来让 Haswell 满意（它至少需要 10 个）。
我正好有一些示例代码，它是针对单精度和 C++ 的，但你会明白要点的：
sumA_1 = _mm256_load_ps(&result[i * N + j]);
sumB_1 = _mm256_load_ps(&result[i * N + j + 8]);
sumA_2 = _mm256_load_ps(&result[(i + 1) * N + j]);
sumB_2 = _mm256_load_ps(&result[(i + 1) * N + j + 8]);
sumA_3 = _mm256_load_ps(&result[(i + 2) * N + j]);
sumB_3 = _mm256_load_ps(&result[(i + 2) * N + j + 8]);
sumA_4 = _mm256_load_ps(&result[(i + 3) * N + j]);
sumB_4 = _mm256_load_ps(&result[(i + 3) * N + j + 8]);
for (size_t k = kk; k < kk + akb; k++) {
auto bc_mat1_1 = _mm256_set1_ps(*mat1ptr);
auto vecA_mat2 = _mm256_load_ps(mat2 + m2idx);
auto vecB_mat2 = _mm256_load_ps(mat2 + m2idx + 8);
sumA_1 = _mm256_fmadd_ps(bc_mat1_1, vecA_mat2, sumA_1);
sumB_1 = _mm256_fmadd_ps(bc_mat1_1, vecB_mat2, sumB_1);
auto bc_mat1_2 = _mm256_set1_ps(mat1ptr[N]);
sumA_2 = _mm256_fmadd_ps(bc_mat1_2, vecA_mat2, sumA_2);
sumB_2 = _mm256_fmadd_ps(bc_mat1_2, vecB_mat2, sumB_2);
auto bc_mat1_3 = _mm256_set1_ps(mat1ptr[N * 2]);
sumA_3 = _mm256_fmadd_ps(bc_mat1_3, vecA_mat2, sumA_3);
sumB_3 = _mm256_fmadd_ps(bc_mat1_3, vecB_mat2, sumB_3);
auto bc_mat1_4 = _mm256_set1_ps(mat1ptr[N * 3]);
sumA_4 = _mm256_fmadd_ps(bc_mat1_4, vecA_mat2, sumA_4);
sumB_4 = _mm256_fmadd_ps(bc_mat1_4, vecB_mat2, sumB_4);
m2idx += 16;
mat1ptr++;
}
_mm256_store_ps(&result[i * N + j], sumA_1);
_mm256_store_ps(&result[i * N + j + 8], sumB_1);
_mm256_store_ps(&result[(i + 1) * N + j], sumA_2);
_mm256_store_ps(&result[(i + 1) * N + j + 8], sumB_2);
_mm256_store_ps(&result[(i + 2) * N + j], sumA_3);
_mm256_store_ps(&result[(i + 2) * N + j + 8], sumB_3);
_mm256_store_ps(&result[(i + 3) * N + j], sumA_4);
_mm256_store_ps(&result[(i + 3) * N + j + 8], sumB_4);
这意味着 j 循环和 i 循环被展开，但 k 循环没有展开，即使它现在是内部循环。稍微展开 k 循环确实对我的实验有所帮助。
请参阅@harold 的回答以获得实际改进。这主要是为了转载我在评论里写的东西
four instead of two sums in the inner loop. (Why doesn't unrolling help?)
sum[i] 没有循环携带的依赖。下一次迭代分配 sum[0] = _mm256_load_pd(C+i1*SIZE_K+j1+0);，它不依赖于先前的值。
因此，将同一架构寄存器重命名为不同的 物理 寄存器足以避免可能使流水线停止的写后写危险。无需使用多个 tmp 变量使源复杂化。请参阅
没有 寄存器重命名的管道 （大多数有序 CPU）将从 "software pipelining" 中受益，以静态调度乱序 exec 可以在fly：加载到不同的寄存器中，因此每个负载和消耗它的 FMA 之间存在距离（充满独立工作）。然后在那和商店之间。
但是所有现代 x86 CPU 都是 OoO；甚至 Knight's Landing 对于 SIMD 向量也有一些有限的 OoO exec。 （Silvermont 不支持 AVX，但按顺序 运行 SIMD 指令，仅对整数执行 OoO）。
在没有任何多累加器情况来隐藏延迟的情况下，展开的好处（明确地在源代码中或使用 -fprofile-use 启用的 -funroll-loop，或默认情况下在 clang 中）是：
- 减少前端带宽以使循环开销进入后端。每个循环开销更有用的工作微指令。因此，如果您的 "useful work" 在前端接近瓶颈，它会有所帮助。
- 减少 运行 循环开销的后端执行单元需求。通常在 Haswell 及更高版本或 Zen 上不是问题；当指令组合包含一些整数内容和一些纯加载指令时，后端基本上可以跟上前端。
- 每次完成的工作总 uops 更少意味着 OoO exec 可以 "see" 更远的内存 loads/stores。
- 有时对短 运行ning 循环的分支预测更好：迭代次数越少意味着分支预测学习的模式越短。因此，对于短行程计数，当执行脱离循环时，有更好的机会正确预测最后一次迭代未采用的次数。
- 有时在更复杂的情况下保存
mov reg,reg，因为编译器更容易在不同的 reg 中生成新结果。同一个变量可以在两个 regs 之间交替，而不是需要moved 回到同一个 regs 为下一次迭代做好准备。特别是如果你有一个以依赖的方式使用a[i]和a[i+1]的循环，或者类似 Fibonacci.
循环中有 2 个加载 + 1 个存储，这可能是瓶颈，而不是 FMA 或前端带宽。按 2 展开可能有助于避免前端瓶颈，但更重要的是只有来自另一个超线程的争用才重要。
评论中提出了一个有趣的问题：展开不需要很多寄存器才有用吗？
哈罗德评论道：
16 is not a huge number of registers, but it's enough to have 12 accumulators and 3 pieces of row vector from B and the broadcasted scalar from A, so it works out to just about enough. The loop from OP above barely uses any registers anyway. The 8 registers in 32bit are indeed too few.
当然，由于问题中的代码在循环迭代中的寄存器中没有 "accumulators"，仅添加到内存中，编译器可以优化所有 sum[0..n] 以重用相同的寄存器汇编；存储后是"dead"。所以实际套印压力很低。
是的，x86-64 有点寄存器不足，这就是为什么 AVX512 将矢量寄存器 (zmm0..31) 的数量和宽度加倍的原因。是的，许多 RISC 有 32 int / 32 fp regs，包括 AArch64，ARM 有 16 个。
x86-64有16个标量整型寄存器（包括栈指针，不包括程序计数器），所以普通函数可以用15个。还有16个向量寄存器，xmm0..15。 （对于 AVX，它们的宽度是 ymm0..15 的两倍）。
（其中一些是在我注意到 sum[0..n] 毫无意义，不是循环携带之前写的。）
在这种情况下，将寄存器重命名到一个大的物理寄存器文件就足够了。在其他情况下，拥有更多架构寄存器会有所帮助，尤其是对于更高的 FP 延迟，因此 AVX512 具有 32 个 zmm 寄存器。但是对于整数 16 已经足够了。 RISC CPU 通常设计为有序，无需重新命名，需要 SW 管道。
对于 OoO exec，就减少 spill/reloads 而言，从 8 到 16 个架构 GP 整数 reg 的跳跃比从 16 到 32 的跳跃更重要。 （我看过一篇论文，用不同数量的架构寄存器测量了 SPECint 的总动态指令数。我没有再次查找，但 8->16 可能总节省了 10%，而 16->32 只是几个 %。类似的东西）。
但是这个具体问题不需要很多 FP 寄存器，sum[0..3] 只需要 4 个向量（如果它们是循环携带的），也许 1 个是临时的； x86 可以使用内存源 mul/add/FMA。寄存器重命名消除了任何 WAW 危害，因此我们可以重复使用相同的临时寄存器而不需要软件流水线。 （并且 OoO exec 还隐藏了负载和 ALU 延迟。）
当存在循环携带的依赖项时，您需要多个累加器。这段代码是添加到内存中，而不是添加到一些向量累加器中，所以任何依赖都是通过 store/reload。但这只有大约 7 个周期的延迟，所以任何理智的缓存阻塞因素都会隐藏它。