密集矩阵向量乘法最有效的实现是什么？
What's the most efficient implementation of dense-matrix-vector multiplication?
如果M是一个稠密的m x n矩阵，v是一个n分量向量，那么乘积u = Mv是u[i] = sum(M[i,j] * v[j], 1 <= j <= n)给出的m分量向量.这种乘法的一个简单实现是
allocate m-component vector u of zeroes
for i = 1:m
for j = 1:n
u[i] += M[i,j] * v[j]
end
end
一次构建向量 u 一个元素。另一种实现是交换循环：
allocate m-component vector u of zeroes
for j = 1:n
for i = 1:m
u[i] += M[i,j] * v[j]
end
end
其中整个向量是一起构建的。
这些实现中的哪一个（如果有）通常用于 C 和 Fortran 等专为高效数值计算而设计的语言中？我的猜测是，像 C 这样在内部以行优先顺序存储矩阵的语言使用前一种实现，而像 Fortran 这样使用列优先顺序的语言使用后者，以便内循环访问矩阵 M 的连续内存位置。是这个正确吗？
前者的实现似乎更高效，因为写入的内存位置只改变m次，而在后一种实现中，写入的内存位置改变m*n次，即使只m 曾经写入过唯一位置。 （当然，按照相同的逻辑，后一种实现对于行向量矩阵乘法会更有效，但这不太常见。）另一方面，我相信 Fortran 通常在密集矩阵向量上更快乘法而不是 C，所以也许我要么 (a) 猜错了它们的实现，要么 (b) 错误判断了两种实现的相对效率。
可能使用已建立的 BLAS 实现是最常见的。除此之外，还有一些简单实现的问题可能值得关注。例如，在 C（或 C++ 中），指针别名通常会阻止大量优化，因此 for example
void matvec(double *M, size_t n, size_t m, double *v, double * u)
{
for (size_t i = 0; i < m; i++) {
for (size_t j = 0; j < n; j++) {
u[i] += M[i * n + j] * v[j];
}
}
}
被Clang 5转成这个（内循环摘录）
.LBB0_4: # Parent Loop BB0_3 Depth=1
vmovsd xmm1, qword ptr [rcx + 8*rax] # xmm1 = mem[0],zero
vfmadd132sd xmm1, xmm0, qword ptr [r13 + 8*rax - 24]
vmovsd qword ptr [r8 + 8*rbx], xmm1
vmovsd xmm0, qword ptr [rcx + 8*rax + 8] # xmm0 = mem[0],zero
vfmadd132sd xmm0, xmm1, qword ptr [r13 + 8*rax - 16]
vmovsd qword ptr [r8 + 8*rbx], xmm0
vmovsd xmm1, qword ptr [rcx + 8*rax + 16] # xmm1 = mem[0],zero
vfmadd132sd xmm1, xmm0, qword ptr [r13 + 8*rax - 8]
vmovsd qword ptr [r8 + 8*rbx], xmm1
vmovsd xmm0, qword ptr [rcx + 8*rax + 24] # xmm0 = mem[0],zero
vfmadd132sd xmm0, xmm1, qword ptr [r13 + 8*rax]
vmovsd qword ptr [r8 + 8*rbx], xmm0
add rax, 4
cmp r11, rax
jne .LBB0_4
看着真的很痛，执行起来会更痛。编译器 "had to" 这样做是因为 u 可能与 M and/or v 互为别名，因此存储到 u 会受到极大的怀疑（"had to" 用引号引起来，因为编译器可以插入一个别名测试，并为好的情况提供快速路径）。在 Fortran 中，默认情况下过程参数不能别名，所以这个问题不会存在。这是一个典型的原因，为什么在没有特殊技巧的情况下随机输入的代码在 Fortran 中比在 C 中更快——我的其余答案不会是关于这个的，我只是想让 C 代码更快一点（最后我回到专栏主要 M)。在 C 中，别名问题可以 be fixed 的一种方式是 restrict，但它唯一要做的就是它不是侵入性的（使用显式累加器而不是求和到 u[i]也可以解决问题，但不依赖于魔术关键字）
void matvec(double *M, size_t n, size_t m, double *v, double * restrict u)
{
for (size_t i = 0; i < m; i++) {
for (size_t j = 0; j < n; j++) {
u[i] += M[i * n + j] * v[j];
}
}
}
现在发生这种情况：
.LBB0_8: # Parent Loop BB0_3 Depth=1
vmovupd ymm5, ymmword ptr [rcx + 8*rbx]
vmovupd ymm6, ymmword ptr [rcx + 8*rbx + 32]
vmovupd ymm7, ymmword ptr [rcx + 8*rbx + 64]
vmovupd ymm8, ymmword ptr [rcx + 8*rbx + 96]
vfmadd132pd ymm5, ymm1, ymmword ptr [rax + 8*rbx - 224]
vfmadd132pd ymm6, ymm2, ymmword ptr [rax + 8*rbx - 192]
vfmadd132pd ymm7, ymm3, ymmword ptr [rax + 8*rbx - 160]
vfmadd132pd ymm8, ymm4, ymmword ptr [rax + 8*rbx - 128]
vmovupd ymm1, ymmword ptr [rcx + 8*rbx + 128]
vmovupd ymm2, ymmword ptr [rcx + 8*rbx + 160]
vmovupd ymm3, ymmword ptr [rcx + 8*rbx + 192]
vmovupd ymm4, ymmword ptr [rcx + 8*rbx + 224]
vfmadd132pd ymm1, ymm5, ymmword ptr [rax + 8*rbx - 96]
vfmadd132pd ymm2, ymm6, ymmword ptr [rax + 8*rbx - 64]
vfmadd132pd ymm3, ymm7, ymmword ptr [rax + 8*rbx - 32]
vfmadd132pd ymm4, ymm8, ymmword ptr [rax + 8*rbx]
add rbx, 32
add rbp, 2
jne .LBB0_8
它不再是标量了，这很好。但并不理想。虽然这里有 8 个 FMA，但它们排列成四对相关的 FMA。纵观整个循环，实际上只有4条独立的FMA依赖链。 FMA 通常具有较长的延迟和不错的吞吐量，例如在 Skylake 上它的延迟为 4，吞吐量为 2/周期，因此需要 8 个独立的 FMA 链来利用所有计算吞吐量。 Haswell 更糟，FMA 的延迟为 5，而吞吐量已经是 2/周期，因此它需要 10 个独立的链。另一个问题是很难实际喂养所有这些 FMA，上面的结构甚至没有真正尝试：每个 FMA 使用 2 个负载，而负载实际上具有与 FMA 相同的吞吐量，因此它们的比率应该是 1:1.
可以通过展开 outer 循环来提高 load:FMA 比率，这让我们可以将 v 的负载重新用于几行M（这甚至不是缓存考虑因素，但它也有帮助）。展开外环也有助于实现拥有更多独立 FMA 链的目标。编译器通常不喜欢展开除内部循环之外的任何内容，因此这需要一些 manual work。 "Tail" 迭代省略（或：假设 m 是 4 的倍数）。
void matvec(double *M, size_t n, size_t m, double *v, double * restrict u)
{
size_t i;
for (i = 0; i + 3 < m; i += 4) {
for (size_t j = 0; j < n; j++) {
size_t it = i;
u[it] += M[it * n + j] * v[j];
it++;
u[it] += M[it * n + j] * v[j];
it++;
u[it] += M[it * n + j] * v[j];
it++;
u[it] += M[it * n + j] * v[j];
}
}
}
不幸的是，Clang 仍然决定错误地展开内部循环，"wrong" 是那个天真的串行展开。只要仍然只有4条独立链，意义不大：
.LBB0_8: # Parent Loop BB0_3 Depth=1
vmovupd ymm5, ymmword ptr [rcx + 8*rdx]
vmovupd ymm6, ymmword ptr [rcx + 8*rdx + 32]
vfmadd231pd ymm4, ymm5, ymmword ptr [r12 + 8*rdx - 32]
vfmadd231pd ymm3, ymm5, ymmword ptr [r13 + 8*rdx - 32]
vfmadd231pd ymm2, ymm5, ymmword ptr [rax + 8*rdx - 32]
vfmadd231pd ymm1, ymm5, ymmword ptr [rbx + 8*rdx - 32]
vfmadd231pd ymm4, ymm6, ymmword ptr [r12 + 8*rdx]
vfmadd231pd ymm3, ymm6, ymmword ptr [r13 + 8*rdx]
vfmadd231pd ymm2, ymm6, ymmword ptr [rax + 8*rdx]
vfmadd231pd ymm1, ymm6, ymmword ptr [rbx + 8*rdx]
add rdx, 8
add rdi, 2
jne .LBB0_8
如果我们停止偷懒并最终做出一些 explicit accumulators:
，这个问题就会消失
void matvec(double *M, size_t n, size_t m, double *v, double *u)
{
size_t i;
for (i = 0; i + 3 < m; i += 4) {
double t0 = 0, t1 = 0, t2 = 0, t3 = 0;
for (size_t j = 0; j < n; j++) {
size_t it = i;
t0 += M[it * n + j] * v[j];
it++;
t1 += M[it * n + j] * v[j];
it++;
t2 += M[it * n + j] * v[j];
it++;
t3 += M[it * n + j] * v[j];
}
u[i] += t0;
u[i + 1] += t1;
u[i + 2] += t2;
u[i + 3] += t3;
}
}
现在 Clang 这样做：
.LBB0_6: # Parent Loop BB0_3 Depth=1
vmovupd ymm8, ymmword ptr [r10 - 32]
vmovupd ymm9, ymmword ptr [r10]
vfmadd231pd ymm6, ymm8, ymmword ptr [rdi]
vfmadd231pd ymm7, ymm9, ymmword ptr [rdi + 32]
lea rax, [rdi + r14]
vfmadd231pd ymm4, ymm8, ymmword ptr [rdi + 8*rsi]
vfmadd231pd ymm5, ymm9, ymmword ptr [rdi + 8*rsi + 32]
vfmadd231pd ymm1, ymm8, ymmword ptr [rax + 8*rsi]
vfmadd231pd ymm3, ymm9, ymmword ptr [rax + 8*rsi + 32]
lea rax, [rax + r14]
vfmadd231pd ymm0, ymm8, ymmword ptr [rax + 8*rsi]
vfmadd231pd ymm2, ymm9, ymmword ptr [rax + 8*rsi + 32]
add rdi, 64
add r10, 64
add rbp, -8
jne .LBB0_6
这是体面的。 load:FMA 比率为 10:8 并且 Haswell 的累加器太少，因此仍有可能进行一些改进。其他一些有趣的展开组合是（外部 x 内部）4x3（12 个累加器、3 个临时累加器、5/4 load:FMA）、5x2（10、2、6/5）、7x2（14、2、8/7） , 15x1（15、1、16/15）。这使得它看起来通过展开外循环更好，但是有太多不同的流（即使不是 "streaming" 在 "streaming loads" 的意义上）不利于自动预取，并且在实际流式传输时它可能超过填充缓冲区的数量是不好的（实际细节很少）。手动预取也是一种选择。获得一个真正好的 MVM 程序需要做更多的工作，尝试很多这样的事情。
将商店保存到 u 以供内部循环外意味着不再需要 restrict。我认为，最令人印象深刻的是，不需要 SIMD 内在函数就可以做到这一点——如果没有可怕的潜在别名，Clang 就可以很好地做到这一点。 GCC 和 ICC 确实尝试过，但展开得不够多，但更多的手动展开可能会成功..
循环平铺也是一个选项，但这是 MVM。平铺对于 MMM 来说是极其必要的，但是 MMM 具有几乎无限量的数据重用，这是 MVM 所没有的。只有向量被重用，矩阵只流过一次。流式传输巨大矩阵的内存吞吐量可能比向量不适合缓存的问题更大。
对于列优先 M，它是不同的，没有明显的循环携带依赖性。通过记忆存在依赖性，但它有很多时间。 load:FMA 比率仍然需要降低，所以它仍然需要展开一些外循环，但总体来说似乎更容易处理。它可以重新排列以主要使用加法，但无论如何 FMA 具有高吞吐量（在 HSW 上，高于加法！）。它不需要水平和，这很烦人，但无论如何它们都发生在内部循环之外。在 return 中，内循环中有商店。如果不尝试，我不认为这些方法之间存在很大的内在差异，这两种方法似乎都应该可以调整到计算吞吐量的 80% 到 90%（对于可缓存大小）。 "annoying extra load" 本质上可以防止任意接近 100%。
如果M是一个稠密的m x n矩阵，v是一个n分量向量，那么乘积u = Mv是u[i] = sum(M[i,j] * v[j], 1 <= j <= n)给出的m分量向量.这种乘法的一个简单实现是
allocate m-component vector u of zeroes
for i = 1:m
for j = 1:n
u[i] += M[i,j] * v[j]
end
end
一次构建向量 u 一个元素。另一种实现是交换循环：
allocate m-component vector u of zeroes
for j = 1:n
for i = 1:m
u[i] += M[i,j] * v[j]
end
end
其中整个向量是一起构建的。
这些实现中的哪一个（如果有）通常用于 C 和 Fortran 等专为高效数值计算而设计的语言中？我的猜测是，像 C 这样在内部以行优先顺序存储矩阵的语言使用前一种实现，而像 Fortran 这样使用列优先顺序的语言使用后者，以便内循环访问矩阵 M 的连续内存位置。是这个正确吗？
前者的实现似乎更高效，因为写入的内存位置只改变m次，而在后一种实现中，写入的内存位置改变m*n次，即使只m 曾经写入过唯一位置。 （当然，按照相同的逻辑，后一种实现对于行向量矩阵乘法会更有效，但这不太常见。）另一方面，我相信 Fortran 通常在密集矩阵向量上更快乘法而不是 C，所以也许我要么 (a) 猜错了它们的实现，要么 (b) 错误判断了两种实现的相对效率。
可能使用已建立的 BLAS 实现是最常见的。除此之外，还有一些简单实现的问题可能值得关注。例如，在 C（或 C++ 中），指针别名通常会阻止大量优化，因此 for example
void matvec(double *M, size_t n, size_t m, double *v, double * u)
{
for (size_t i = 0; i < m; i++) {
for (size_t j = 0; j < n; j++) {
u[i] += M[i * n + j] * v[j];
}
}
}
被Clang 5转成这个（内循环摘录）
.LBB0_4: # Parent Loop BB0_3 Depth=1
vmovsd xmm1, qword ptr [rcx + 8*rax] # xmm1 = mem[0],zero
vfmadd132sd xmm1, xmm0, qword ptr [r13 + 8*rax - 24]
vmovsd qword ptr [r8 + 8*rbx], xmm1
vmovsd xmm0, qword ptr [rcx + 8*rax + 8] # xmm0 = mem[0],zero
vfmadd132sd xmm0, xmm1, qword ptr [r13 + 8*rax - 16]
vmovsd qword ptr [r8 + 8*rbx], xmm0
vmovsd xmm1, qword ptr [rcx + 8*rax + 16] # xmm1 = mem[0],zero
vfmadd132sd xmm1, xmm0, qword ptr [r13 + 8*rax - 8]
vmovsd qword ptr [r8 + 8*rbx], xmm1
vmovsd xmm0, qword ptr [rcx + 8*rax + 24] # xmm0 = mem[0],zero
vfmadd132sd xmm0, xmm1, qword ptr [r13 + 8*rax]
vmovsd qword ptr [r8 + 8*rbx], xmm0
add rax, 4
cmp r11, rax
jne .LBB0_4
看着真的很痛，执行起来会更痛。编译器 "had to" 这样做是因为 u 可能与 M and/or v 互为别名，因此存储到 u 会受到极大的怀疑（"had to" 用引号引起来，因为编译器可以插入一个别名测试，并为好的情况提供快速路径）。在 Fortran 中，默认情况下过程参数不能别名，所以这个问题不会存在。这是一个典型的原因，为什么在没有特殊技巧的情况下随机输入的代码在 Fortran 中比在 C 中更快——我的其余答案不会是关于这个的，我只是想让 C 代码更快一点（最后我回到专栏主要 M)。在 C 中，别名问题可以 be fixed 的一种方式是 restrict，但它唯一要做的就是它不是侵入性的（使用显式累加器而不是求和到 u[i]也可以解决问题，但不依赖于魔术关键字）
void matvec(double *M, size_t n, size_t m, double *v, double * restrict u)
{
for (size_t i = 0; i < m; i++) {
for (size_t j = 0; j < n; j++) {
u[i] += M[i * n + j] * v[j];
}
}
}
现在发生这种情况：
.LBB0_8: # Parent Loop BB0_3 Depth=1
vmovupd ymm5, ymmword ptr [rcx + 8*rbx]
vmovupd ymm6, ymmword ptr [rcx + 8*rbx + 32]
vmovupd ymm7, ymmword ptr [rcx + 8*rbx + 64]
vmovupd ymm8, ymmword ptr [rcx + 8*rbx + 96]
vfmadd132pd ymm5, ymm1, ymmword ptr [rax + 8*rbx - 224]
vfmadd132pd ymm6, ymm2, ymmword ptr [rax + 8*rbx - 192]
vfmadd132pd ymm7, ymm3, ymmword ptr [rax + 8*rbx - 160]
vfmadd132pd ymm8, ymm4, ymmword ptr [rax + 8*rbx - 128]
vmovupd ymm1, ymmword ptr [rcx + 8*rbx + 128]
vmovupd ymm2, ymmword ptr [rcx + 8*rbx + 160]
vmovupd ymm3, ymmword ptr [rcx + 8*rbx + 192]
vmovupd ymm4, ymmword ptr [rcx + 8*rbx + 224]
vfmadd132pd ymm1, ymm5, ymmword ptr [rax + 8*rbx - 96]
vfmadd132pd ymm2, ymm6, ymmword ptr [rax + 8*rbx - 64]
vfmadd132pd ymm3, ymm7, ymmword ptr [rax + 8*rbx - 32]
vfmadd132pd ymm4, ymm8, ymmword ptr [rax + 8*rbx]
add rbx, 32
add rbp, 2
jne .LBB0_8
它不再是标量了，这很好。但并不理想。虽然这里有 8 个 FMA，但它们排列成四对相关的 FMA。纵观整个循环，实际上只有4条独立的FMA依赖链。 FMA 通常具有较长的延迟和不错的吞吐量，例如在 Skylake 上它的延迟为 4，吞吐量为 2/周期，因此需要 8 个独立的 FMA 链来利用所有计算吞吐量。 Haswell 更糟，FMA 的延迟为 5，而吞吐量已经是 2/周期，因此它需要 10 个独立的链。另一个问题是很难实际喂养所有这些 FMA，上面的结构甚至没有真正尝试：每个 FMA 使用 2 个负载，而负载实际上具有与 FMA 相同的吞吐量，因此它们的比率应该是 1:1.
可以通过展开 outer 循环来提高 load:FMA 比率，这让我们可以将 v 的负载重新用于几行M（这甚至不是缓存考虑因素，但它也有帮助）。展开外环也有助于实现拥有更多独立 FMA 链的目标。编译器通常不喜欢展开除内部循环之外的任何内容，因此这需要一些 manual work。 "Tail" 迭代省略（或：假设 m 是 4 的倍数）。
void matvec(double *M, size_t n, size_t m, double *v, double * restrict u)
{
size_t i;
for (i = 0; i + 3 < m; i += 4) {
for (size_t j = 0; j < n; j++) {
size_t it = i;
u[it] += M[it * n + j] * v[j];
it++;
u[it] += M[it * n + j] * v[j];
it++;
u[it] += M[it * n + j] * v[j];
it++;
u[it] += M[it * n + j] * v[j];
}
}
}
不幸的是，Clang 仍然决定错误地展开内部循环，"wrong" 是那个天真的串行展开。只要仍然只有4条独立链，意义不大：
.LBB0_8: # Parent Loop BB0_3 Depth=1
vmovupd ymm5, ymmword ptr [rcx + 8*rdx]
vmovupd ymm6, ymmword ptr [rcx + 8*rdx + 32]
vfmadd231pd ymm4, ymm5, ymmword ptr [r12 + 8*rdx - 32]
vfmadd231pd ymm3, ymm5, ymmword ptr [r13 + 8*rdx - 32]
vfmadd231pd ymm2, ymm5, ymmword ptr [rax + 8*rdx - 32]
vfmadd231pd ymm1, ymm5, ymmword ptr [rbx + 8*rdx - 32]
vfmadd231pd ymm4, ymm6, ymmword ptr [r12 + 8*rdx]
vfmadd231pd ymm3, ymm6, ymmword ptr [r13 + 8*rdx]
vfmadd231pd ymm2, ymm6, ymmword ptr [rax + 8*rdx]
vfmadd231pd ymm1, ymm6, ymmword ptr [rbx + 8*rdx]
add rdx, 8
add rdi, 2
jne .LBB0_8
如果我们停止偷懒并最终做出一些 explicit accumulators:
，这个问题就会消失void matvec(double *M, size_t n, size_t m, double *v, double *u)
{
size_t i;
for (i = 0; i + 3 < m; i += 4) {
double t0 = 0, t1 = 0, t2 = 0, t3 = 0;
for (size_t j = 0; j < n; j++) {
size_t it = i;
t0 += M[it * n + j] * v[j];
it++;
t1 += M[it * n + j] * v[j];
it++;
t2 += M[it * n + j] * v[j];
it++;
t3 += M[it * n + j] * v[j];
}
u[i] += t0;
u[i + 1] += t1;
u[i + 2] += t2;
u[i + 3] += t3;
}
}
现在 Clang 这样做：
.LBB0_6: # Parent Loop BB0_3 Depth=1
vmovupd ymm8, ymmword ptr [r10 - 32]
vmovupd ymm9, ymmword ptr [r10]
vfmadd231pd ymm6, ymm8, ymmword ptr [rdi]
vfmadd231pd ymm7, ymm9, ymmword ptr [rdi + 32]
lea rax, [rdi + r14]
vfmadd231pd ymm4, ymm8, ymmword ptr [rdi + 8*rsi]
vfmadd231pd ymm5, ymm9, ymmword ptr [rdi + 8*rsi + 32]
vfmadd231pd ymm1, ymm8, ymmword ptr [rax + 8*rsi]
vfmadd231pd ymm3, ymm9, ymmword ptr [rax + 8*rsi + 32]
lea rax, [rax + r14]
vfmadd231pd ymm0, ymm8, ymmword ptr [rax + 8*rsi]
vfmadd231pd ymm2, ymm9, ymmword ptr [rax + 8*rsi + 32]
add rdi, 64
add r10, 64
add rbp, -8
jne .LBB0_6
这是体面的。 load:FMA 比率为 10:8 并且 Haswell 的累加器太少，因此仍有可能进行一些改进。其他一些有趣的展开组合是（外部 x 内部）4x3（12 个累加器、3 个临时累加器、5/4 load:FMA）、5x2（10、2、6/5）、7x2（14、2、8/7） , 15x1（15、1、16/15）。这使得它看起来通过展开外循环更好，但是有太多不同的流（即使不是 "streaming" 在 "streaming loads" 的意义上）不利于自动预取，并且在实际流式传输时它可能超过填充缓冲区的数量是不好的（实际细节很少）。手动预取也是一种选择。获得一个真正好的 MVM 程序需要做更多的工作，尝试很多这样的事情。
将商店保存到 u 以供内部循环外意味着不再需要 restrict。我认为，最令人印象深刻的是，不需要 SIMD 内在函数就可以做到这一点——如果没有可怕的潜在别名，Clang 就可以很好地做到这一点。 GCC 和 ICC 确实尝试过，但展开得不够多，但更多的手动展开可能会成功..
循环平铺也是一个选项，但这是 MVM。平铺对于 MMM 来说是极其必要的，但是 MMM 具有几乎无限量的数据重用，这是 MVM 所没有的。只有向量被重用，矩阵只流过一次。流式传输巨大矩阵的内存吞吐量可能比向量不适合缓存的问题更大。
对于列优先 M，它是不同的，没有明显的循环携带依赖性。通过记忆存在依赖性，但它有很多时间。 load:FMA 比率仍然需要降低，所以它仍然需要展开一些外循环，但总体来说似乎更容易处理。它可以重新排列以主要使用加法，但无论如何 FMA 具有高吞吐量（在 HSW 上，高于加法！）。它不需要水平和，这很烦人，但无论如何它们都发生在内部循环之外。在 return 中，内循环中有商店。如果不尝试，我不认为这些方法之间存在很大的内在差异，这两种方法似乎都应该可以调整到计算吞吐量的 80% 到 90%（对于可缓存大小）。 "annoying extra load" 本质上可以防止任意接近 100%。