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

Question

如果M是一个密集的m×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）错误地判断两个实现的相对效率。

Answer 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]; 
     } 
    } 
} 

是由锵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]也可以做到这一点，但不依赖于magic关键字）

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/cycle，因此需要8个独立的FMA链来利用所有的计算吞吐量。 Haswell甚至更糟，FMA的延迟为5，并且吞吐量已达每周2次，因此需要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]; 
     } 
    } 
} 

不幸的是，铿锵仍然决定展开内循环错误，与“错误”是天真的串行展开。

.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; 
    } 
} 

现在锵做到这一点：

只要仍然只有4个独立的连锁店有没有什么点

.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就非常好。海湾合作委员会和国际商会做的尝试，但不展开足够的，但更多的手动展开可能会做的伎俩..

循环平铺也是一个选项，但这是MVM。平铺对于MMM来说是非常必要的，但是MMM具有几乎无限量的数据重用，而MVM不具备这一点。只有矢量被重用，矩阵只是流过一次。对于流式处理大型矩阵而言，内存吞吐量可能会比不适合缓存的向量更大。

对于列主M，它是不同的，没有显着的循环携带依赖性。通过记忆存在依赖性，但它有很多时间。负载：FMA比率仍然需要降低，所以它仍然需要一些外循环的展开，但总体来看似乎更容易处理。它可以重新安排使用大多数添加，但FMA无论如何都具有高吞吐量（在HSW上，高于添加！）。它不需要水平的总和，这是令人讨厌的，但它们发生在内部循环之外。作为回报，内部循环中有商店。没有尝试它，我不认为这些方法之间存在巨大的内在差异，似乎两种方法都应该调整到计算吞吐量的80％到90％之间（对于可缓存大小）。 “恼人的额外负荷”固有地防止任意方式接近100％。