密集矩阵向量乘法的最有效实现是什么？答案

【问题标题】：What's the most efficient implementation of dense-matrix-vector multiplication?密集矩阵向量乘法的最有效实现是什么？
【发布时间】：2017-11-10 20:02:14
【问题描述】：

如果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) 误判了两种实现的相对效率。

【问题讨论】：

在线性代数运算中，缓冲区加载和卸载是瓶颈。算法已经饱和en.wikipedia.org/wiki/…，所以你不会在任何方法上获得太大的性能差异。这不是您可以轻松实现并期望性能的东西。数十年的优化已经到位。为此，您可以安全地选择专为此类操作优化的 MKL 或 OpenBLAS 实现。

标签： performance matrix language-agnostic linear-algebra matrix-multiplication

【解决方案1】：

可能使用已建立的 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

看着真的很痛，执行起来会更痛。编译器“必须”这样做是因为u 可能与M 和/或v 有别名，因此对u 的存储会受到极大怀疑（“必须”在引号中，因为编译器可以插入测试混叠并为好案例提供快速路径）。在 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/cycle，因此它需要 10 个独立的链。另一个问题是，实际上很难为所有这些 FMA 提供数据，上面的结构甚至没有真正尝试过：每个 FMA 使用 2 个负载，而负载实际上与 FMA 具有相同的吞吐量，因此它们的比率应该是 1:1。

可以通过展开外部循环来提高负载：FMA 比率，这让我们可以将来自v 的负载重新用于M 的几行（这甚至不是缓存考虑，但它也有帮助）。展开外循环也有助于实现拥有更多独立 FMA 链的目标。编译器通常不喜欢展开除内部循环之外的任何内容，因此这需要一些 manual work。省略“尾”迭代（或：假设 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 仍然决定错误地展开内部循环，“错误”是那种幼稚的串行展开。只要仍然只有 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

这是体面的。负载：FMA 比为 10:8，Haswell 的蓄电池太少，因此仍有一些改进的可能。其他一些有趣的展开组合是（外部 x 内部）4x3（12 个累加器，3 个临时，5/4 负载：FMA）、5x2（10、2、6/5）、7x2（14、2、8/7）、15x1 (15, 1, 16/15)。这使得通过展开外部循环看起来更好，但是有太多不同的流（即使不是“流加载”意义上的“流”）对于自动预取是不利的，并且在实际流时它可能是不好的超过填充缓冲区的数量（实际细节很少）。手动预取也是一种选择。要获得一个真正好的 MVM 过程需要更多的工作，尝试很多这些事情。

将商店保存到u 以供在内循环之外意味着不再需要restrict。我认为，最令人印象深刻的是，不需要 SIMD 内在函数就可以走到这一步——如果没有可怕的潜在混叠，Clang 就很好。 GCC 和 ICC 确实尝试过，但展开不够充分，但更多的手动展开可能会奏效..

循环平铺也是一种选择，但这是 MVM。平铺对于 MMM 来说是非常必要的，但是 MMM 具有几乎无限量的数据重用，而 MVM 没有。仅重用向量，矩阵仅流过一次。流式传输巨大矩阵的内存吞吐量可能比不适合缓存的向量更大。

使用列优先的 M，情况就不同了，没有明显的循环携带依赖性。通过内存存在依赖关系，但它有很多时间。负载：FMA 比率仍然必须降低，因此仍然需要展开外部循环，但总体而言似乎更容易处理。它可以重新排列以主要使用加法，但无论如何 FMA 具有高吞吐量（在 HSW 上，高于加法！）。它不需要水平总和，这很烦人，但无论如何它们都发生在内循环之外。作为回报，在内循环中有商店。如果不尝试，我预计这些方法之间不会有很大的内在差异，似乎两种方法都应该可以调整到计算吞吐量的 80% 到 90% 之间（对于可缓存大小）。 “烦人的额外负载”本质上可以防止任意接近 100%。

【讨论】：

流式加载 movntdqa 在回写内存上没有任何特殊作用，仅在 USWC 上（例如从视频内存复制回来）。但是 NT 预取在 WB 内存上仍然是特殊的。
离题：stackoverflow.com/questions/47461547/… 已被 OP 编辑以删除失败的尝试。您删除的答案现在将回答问题。因为我无法对已删除的答案发表评论，所以在此处向您发送最近的其他答案。