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

可能使用已建立的 BLAS 实现是最常见的。除此之外，简单实现也存在一些问题，这些问题可能很有趣。例如，在 C（或 C++）中，指针别名通常会阻止大量优化，因此例如

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

这看起来真的很痛苦，执行起来会更痛苦。编译器“必须”这样做，因为可能用 and/或 进行别名，因此 store into 会受到极大的怀疑（“不得不”在引号中，因为编译器可以插入别名测试，并为好情况提供快速路径）。在 Fortran 中，过程参数默认不能别名，因此不存在此问题。这就是为什么在没有特殊技巧的情况下随机输入的代码在 Fortran 中比在 C 中更快的典型原因 - 我的其余答案不会是关于这个的，我只是要让 C 代码更快一点（最后我回到列专业）。在 C 语言中，可以解决混叠问题的一种方法是 ，但它唯一能做的就是它不是侵入性的（使用显式累加器而不是求和也可以解决问题，但不依赖于魔术关键字）uMvuMrestrictu[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 链来利用所有这些计算吞吐量。更糟糕的是，FMA 的延迟为 5，而吞吐量已经达到 2/周期，因此它需要 10 条独立的链。另一个问题是，很难真正为所有这些 FMA 供电，上面的结构甚至没有真正尝试过：它每个 FMA 使用 2 个负载，而负载实际上具有与 FMA 相同的吞吐量，因此它们的比例应该是 1：1。

提高负载：FMA 比率可以通过展开外部循环来完成，这让我们可以重用几行的负载（这甚至不是缓存方面的考虑，但它也有助于这一点）。展开外循环也有助于实现拥有更多独立的 FMA 链的目标。编译器通常不喜欢展开任何内容，只喜欢内部循环，因此这需要一些手动工作。省略“尾部”迭代（或：假设是 4 的倍数）。vMm

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

如果我们停止懒惰并最终制作一些显式累加器，这个问题就会消失：

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程序，需要更多的工作，尝试很多这些东西。

将存储保存到内部循环之外意味着不再需要这样做。我认为，最令人印象深刻的是，不需要 SIMD 内部函数就可以走到这一步——如果没有可怕的潜在混叠，Clang 在这方面做得非常好。GCC 和 ICC 确实尝试过，但展开不够，但更多的手动展开可能会奏效。urestrict

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

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