Mar*_*ari 6 c++ optimization avx matrix-multiplication cpu-cache
我的任务是实现一个优化的矩阵乘法微内核，该微内核C = A*B从以下代码片段开始用 C++ 进行计算。我遇到了一些违反直觉的行为，我需要一些帮助来更好地理解正在发生的事情。
void mat_mul(double* A, double* B, double* C) {\n for (int n = 0; n < N; ++n) {\n for (int k = 0; k < K; ++k) {\n for (int m = 0; m < M; ++m) {\n C[m + n*M] += A[m + k*M] * B[k + n*K];\n }\n }\n }\n}\n挑战条件如下：
\n通过查看循环的顺序，当我们线性迭代缓冲区时，内存访问顺序似乎已经最大限度地减少了缓存未命中。
\n我的第一个想法是尝试对代码进行矢量化。这就是我想出的：
\nvoid mat_mul(double* A, double* B, double* C) {\n for (int n = 0; n < N; ++n) {\n for (int k = 0; k < K; ++k) {\n\n __m256d B_v = _mm256_broadcast_sd(B + k + n*K);\n \n for (int m = 0; m < M; m+=4) {\n\n __m256d A_v = _mm256_load_pd(A + m + k*M);\n __m256d C_v = _mm256_load_pd(C + m + n*M);\n __m256d rax = _mm256_fmadd_pd(A_v, B_v, C_v);\n _mm256_store_pd(C + m + n*M, rax);\n\n }\n }\n }\n}\n\n当 M=N=32 时，最高可达 23 GFLOP 左右。当减少 N 和 M 时，性能会下降。
\n想了一会儿，我意识到EPYC 7742的L1d缓存是32KiB，相当于4096个双精度。理想情况下，我希望将所有三个矩阵加载到 L1 缓存中。
\n这会产生以下优化问题：
\n最大化 N>0,M>0，使得 N*M + 128*N + 128 * M \xe2\x89\xa4 4096。
\n我注意到 M = 4, N = 8 是一个使用二次幂值的足够好的解决方案。
\n假设 M=4，我可以通过保留单个向量作为累加变量来消除 m 循环。
\n这个想法产生了以下代码：
\nvoid mat_mul(double* A, double* B, double* C) {\n \n __m256d A_v, B_v;\n register __m256d C_v;\n double* aPtr; \n double* bPtr;\n double* cPtr;\n\n for (int n = 0; n < N; ++n) {\n cPtr = C + n*M;\n C_v = _mm256_load_pd(cPtr);\n\n for (int k = 0; k < K; ++k) {\n aPtr = A + k*M;\n bPtr = B + k + n*K;\n\n B_v = _mm256_broadcast_sd(bPtr);\n A_v = _mm256_load_pd(aPtr);\n\n C_v = _mm256_fmadd_pd(A_v, B_v, C_v);\n }\n\n _mm256_store_pd(cPtr, C_v);\n }\n}\n我认为这会快得多，但是我获得的性能约为 4 GFLOP，这与以次优 M=4、N=8 运行第一个代码相同。
\n鉴于在第一种情况下，矩阵太大而无法完全放入 L1d 缓存，这似乎表明代码的第二个版本正在读取数据并将数据写入 L2 缓存，即使矩阵应该适合 L1d 缓存。
\n我的怀疑正确吗？如果是这样，这种行为是由我的思维过程中的某些错误引起的，还是我在这背后遗漏了某些原因？
\n请给出一些解释，而不仅仅是发布代码的更好优化版本，因为我真的很想更好地理解概念上发生的事情。理想情况下，如果你们能给我一些关于我可以自我审视的事情的提示，那就太好了。
\n谢谢 ：）
\n我尝试遵循@doug 和@ElderBug 评论中建议的一些技巧。
\n@doug 特别建议尝试转置 B，但是考虑到矩阵采用列主格式，我找不到实现他们的想法的方法。相反，我所做的是转置 A 并将乘积累加到 tmp 变量中。
\n这是我想出的：
\nvoid mat_mul(double* A, double* B, double* C) {\n \n __m256d B_v, A_v, rax;\n\n // Transpose A\n double AT[K*M];\n\n for(int k=0; k < K; ++k){\n for(int m=0; m<M; ++m){\n AT[m*K + k] = A[m + k*M];\n }\n }\n\n // Compute product\n for (int n = 0; n < N; ++n) {\n for (int m = 0; m < M; ++m) {\n \n double tmp = 0.;\n\n for (int k = 0; k < K; k+=4) {\n B_v = _mm256_load_pd(B + n*K + k);\n A_v = _mm256_load_pd(AT + m*K + k);\n rax = _mm256_mul_pd(A_v, B_v);\n\n double* pv = (double*)&rax;\n tmp += pv[0] + pv[1] + pv[2] + pv[3];\n }\n\n C[n*M + m] = tmp;\n }\n }\n}\n当 M=4、N=8 时，其运行速度仍约为 4 GFLOP。房间里的大象似乎正在减少 rax 向量。我无法找到一种不同的方法来更有效地做到这一点。
\n如果我消除tmp += pv[0] + pv[1] + pv[2] + pv[3];性能加倍，但我在 M=N=32 的情况下达到 14 GFLOP 的峰值，所以这仍然比我的简单矢量化方法更糟糕。
如果有人有任何进一步的建议/意见，他们将不胜感激。
\n我忘了提及我正在使用icc (ICC) 2021.5.0 20211109以下优化标志来编译代码：\n-O3 -funroll-loops -std=c++11 -mavx2 -march=native -fma -ftz -fomit-frame-pointer
目标是以尽可能最好的方式实现这个串行微内核，以便我可以重新使用它来进行分块并行矩阵乘法。根据我的计算，理论峰值应该是 46 GFLOPS，所以我得到了大约 50%。OpenBLAS 大约有 40 GFLOP，所以 32-35 左右就不错了。
\n如果有人能够深入了解为什么我尝试的某些选项不起作用，以便我可以考虑如何解决它们，那就太好了。
\n感谢 @ElderBug\ 的评论，我对如何管理第一个矢量化实现中的存储操作进行了一些额外的思考，并且通过进行一些巧妙的展开，我能够达到 40GFLOps，这与 OpenBLAS 相当！
\n阅读本文也有帮助：https://github.com/flame/blis/blob/master/docs/KernelsHowTo.md#gemm-microkernel
\n现代 CPU 是巨大的、令人难以置信的超标量和无序野兽。您的 Zen 2 CPU 可以运行 100 多条指令。这意味着当您的第一个 FMA 完成时，CPU 可能已经提前 10 个循环发出负载。这足以隐藏 L1 甚至 L2 延迟。当然，这种方法的实现是有条件的，例如，必须不存在可能阻止提前计算负载的依赖关系。分支预测还必须具有高预测率。矩阵乘法是一个“不错”的算法，一切都很好。
当然，这并不意味着应该忽略缓存。局部性是最重要的，以确保您使用在缓存中获取的所有数据。重用是很好的，可以减少带宽的使用。可预测的加载模式是可取的，以便硬件预取器可以在您需要它之前填充缓存。
在您的情况下，缓存预取是完美的，因为您的内部访问是连续的。这可能就是为什么你的矩阵是否适合 L1 并不重要：大多数访问都会命中，因为它们是预测的。其他访问可能会丢失，但由于它们的数量少得多，所以影响较小（并且您不能排除它们仍然是预取的）。
TL;DR：虽然您可能仍然会发现缓存方面的一些小改进，但无论矩阵有多大，它都不太可能成为问题。
现在怎么办 ？这意味着提高性能的主要方法是尽可能保持管道由独立的 FMA 供电。这意味着展开循环和巧妙地使用寄存器。这很无聊，但事实就是如此。您的第一个版本展开为仅依赖于一次广播的 8 个并行 FMA。我的第一个建议是只采用您的第一个版本，并将 C mat 访问移出循环，存储在 a 中__m256d[8]，以便删除仍在该循环中的无用代码。
需要注意的一件重要事情是，如果您认为 L1 不相关，则意味着您正在从 L2 中获取 FMA。这不一定是所解释的问题，但 L2 的带宽仍然小于 L1。L2带宽够用吗？
Zen 及更高版本的 CPU 每周期 L2 带宽为 32B。每个 FMA 大约有一个 256 位负载，因此这将 FMA 的带宽消耗限制为每个 FMA 32B。Zen FMA 吞吐量为每个周期 1，因此这非常适合 L2。然而，对于 Zen 2（您的 CPU），FMA 吞吐量加倍到每个周期 2，这意味着 L2 无法以最大带宽为 FMA 提供数据。
所以最后，矩阵大小对于获得后半部分的性能实际上很重要。更大的矩阵将趋向于最大 FLOPS 的一半。在你的情况下，它们可能仍然足够小，即使它们并不全部适合，你仍然可以获得 L1 命中。
要解决这个问题，您可以更改算法以保证一定程度的 L1 重用。在迭代 K 时，可能会交错 2 或 4 N 列，从而将 M 减少那么多，以便它们仍然适合寄存器。当然，让矩阵足够小以适合缓存也是一种解决方案。