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

Question

如果M是稠密的mxn矩阵，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

其一次一个元素地建立向量CCD_ 5。另一个实现是交换循环：

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的唯一位置被写入，或者(b)错误判断两种实现方式的相对效率。

Answer 1

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

看着真的很痛苦，执行起来会更痛苦。编译器"不得不"这样做，因为u可能与M和/或v别名，因此对u中的存储非常怀疑("不得不"用引号括起来，因为编译器可以插入别名测试，并为良好的情况提供快速路径)。在Fortran中，过程参数默认情况下不能别名，所以这个问题不会存在。这就是为什么在Fortran中随机键入而不使用特殊技巧的代码比在C中更快的典型原因——我的其余答案将不涉及此，我只想让C代码更快一点(最后我回到M专栏)。在C中，可以解决混叠问题的一种方法是使用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。

可以通过展开外部循环来提高负载：FMA比率，这使我们可以对M的几行重用v中的负载(这甚至不是缓存的考虑因素，但也有助于此)。展开外环也有助于实现拥有更多独立的FMA链的目标。编译器通常不喜欢展开除内部循环之外的任何内容，因此这需要一些手动操作。省略了"尾"迭代(或者：假设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

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

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%。