AVX 内在澄清，4x4 矩阵乘法奇数

AVX Intrinsic Clarification, 4x4 Matrix Multiplication Oddities

本文关键字：4x4 AVX 更新时间：2023-10-16

在纸上，我画出了这个算法的长形式，在纸上它应该可以正常工作。我是否遇到了寄存器转换(256/128/256)的微妙之处，还是我实际上在某处搞砸了算法结构？

为了方便起见，我将原版代码和 AVX 代码放在 Godbolt 查看器上，以便您可以随意查看生成的程序集。

标准代码 https://godbolt.org/g/v47RKH

我的 AVX 尝试 1： https://godbolt.org/g/oH1DpO

我的AVX尝试2： https://godbolt.org/g/QFtdKr(减少 5 个循环，减少铸造需求，更易于阅读)

奇怪的是，SSE 代码使用的是标量操作，这让我难以置信，因为这绝对可以通过水平广播、muls 和添加来加速。我想做的是把这个概念提升一个层次。

RHS永远不需要改变，但本质上如果LHS是{a， b， ...， p}， LHS 是 {1， 2， ...， 16}，那么我们只需要 2 个寄存器来保存 RHS 的 2 个半部分，然后需要 2 个寄存器来保存给定的 LHS 行，形式为 {a， a， a， a， b， b， b， b} 和 {c， c， c， c， d， d， d， d}。这是通过 2 次广播和 256/128/256 演员表实现的。

我们得到的中间结果

{a*1， a*2， a*3， a*4， b*5， b*6， b*

7， b*8} =>行[0]

和

{C*9， C*10， C*11， C*12， D*13， D*14， D*

15， D*16} =>行[1]

一旦 w.r.t LHS 展开，我们就会生成

{e*1， ...f*8}， {g*9， ...h*16} => 行[2]，行[3]

接下来将 r0，r1 和 r2，r3 加在一起(保留 r0 和 r2 作为当前中间体)

最后，将行 [0] 的高半部分提取到 resHalf 的下半部分，将行 [2] 的低半部分插入 resHalf 的高半部分，将行的高半部分 [2] 插入行的高半部分 [0]，然后将行 [0] 添加到 resHalf。

按所有权利，这应该给我们留下 resHalf[0] 等于迭代结束时的以下内容 i = 0

{a*1 + b*2 + c*3 + d*

4， a*5 + b*6 + c*7 + d*8，

a*9 + b*10 + c*11 + d*12， a*13 + b*14 + c*15 + d*16，

e*1 + ... + h*4， e*5 + ... + h*8，

e*9 + ... + h*12， e*13 + ... + h*16}

但是，我的算法产生如下：

2x {a*1 + c*3， a*5 + c*7， a*9 + c*

11， a*13 + c*15}，

2x {e*1 + g*3， e*5 + g*7， e*9 + g*11， e*13 + g*15}

更可怕的是，如果我在三元条件中交换 rhsHolders[0/1]，它根本不会改变结果。就好像编译器忽略了其中一个交换和添加。Clang 4 和 GCC 7 都这样做了，所以我在哪里搞砸了？

编辑：输出应该是 4 行 {10， 26， 42， 58}，但我得到 {4， 12， 20， 28}

奇怪的是，SSE 代码使用的是标量操作，这让我难以置信，因为这绝对可以通过水平广播、muls 和 add 来加速。

是指编译器生成的汇编代码吗？ clang4.0 和 gcc7.1 输出中MatMul()的所有 AVX 指令都在 ymm 向量上运行。除了 clang 愚蠢的广播负载：它执行标量加载，然后执行单独的 AVX2 广播指令，这是特别糟糕的，因为英特尔 CPU 将广播负载作为单 uop ALU 指令处理。加载端口本身可以进行广播。但是，如果源是寄存器，则需要一个用于随机端口的 ALU uop。

vmovss  xmm5, dword ptr [rdi + 24] # xmm5 = mem[0],zero,zero,zero
vbroadcastss    xmm5, xmm5

与 GCC 使用的 AVX1vbroadcastss xmm5, [rdi + 24]相比，Clang 的实际输出(上图)真的很愚蠢。

在main()中，clang 确实发出标量运算。

由于你的输入矩阵都是编译时常量，唯一的谜团是为什么它没有优化到cout << "a long string with the numbers already formattedn";，或者至少优化去掉所有的数学，只把double的结果准备好打印。 (是的，它们正在打印循环中从float转换为vcvtss2sddouble。

它通过一些内在的随机和数学进行优化，在编译时进行优化。我猜 clang 在洗牌的某个地方迷路了，仍然发出一些数学运算。它们是标量的事实可能表明它在编译时没有做太多工作，但它没有重新排序以对其进行矢量化。

请注意，某些常量不会出现在源中，并且它们在内存中不按升序排列。

...
.LCPI1_5:
.long   1092616192              # float 10
.LCPI1_6:
.long   1101004800              # float 20
.LCPI1_7:
.long   1098907648              # float 16
...

chlang如何将浮点值放在位模式的整数表示之后的注释中，这真的很好。

还是我实际上在某处搞砸了算法结构？

好吧，这部分实现看起来完全是假的。从rows[j]初始化lowerHalf，然后在下一条语句中覆盖该值。

__m128 lowerHalf = _mm256_castps256_ps128(rows[j]);
lowerHalf = _mm_broadcast_ss(&lhs[offset+2*j]);

然后你做一个 256b 乘法，rows[j]未定义的上 128b 车道。

rows[j] = _mm256_castps128_ps256(lowerHalf);
rows[j] = _mm256_mul_ps(rows[j], (chooser) ? rhsHolders[0] : rhsHolders[1]);

在 gcc 和 clang 的 asm 中，上通道全部为零(因为它们显然选择了使用标量 -> xmm 广播最后写入的 ymm 寄存器，该寄存器隐式地零扩展到最大矢量宽度)。请注意，_mm256_castps128_ps256不保证零扩展。除非__m128本身是从 256b 或更宽的向量中提取/强制转换的结果，否则很有可能，但它是不确定的。请参阅如何清除__m256值的上 128 位？适用于在矢量中需要归零的上车道的情况。

无论如何，这意味着你会从 128b 向量乘法 (vmulps xmm, xmm, xmm) 得到相同的结果：在这些指令之后，上面的 4 个元素都将为零(或 NaN)

vbroadcastss    xmm0, DWORD PTR [rdi+40]
vmulps  ymm0, ymm2, ymm0

这种 asm 输出(来自 gcc7.1)极不可能成为正确 matmul 实现的一部分。

我没有仔细研究你到底想在源代码中做什么，但我认为它不完全是这个。

更可怕的是，如果我在三元条件中交换 rhsHolders[0/1]，它根本不会改变结果。就好像编译器忽略了其中一个交换和添加。

当更改源代码中的某些内容不会在 asm 输出中产生您期望的更改时，这表明您可能弄错了源代码，并且某些内容正在优化。有时我复制/粘贴一个内部函数，忘记在新行中更改输入变量，所以我的函数忽略了它的一些计算结果并使用另一个结果两次。

它几乎是从我昨天在 SO 上的答案中复制和粘贴的:)

试试这个

void MatMul(const float* __restrict lhs , const float* __restrict rhs , float* __restrict out ) 
{
lhs = reinterpret_cast<float*>(__builtin_assume_aligned (lhs, 32));
rhs = reinterpret_cast<float*>(__builtin_assume_aligned (rhs, 32));
out = reinterpret_cast<float*>(__builtin_assume_aligned (out, 32));
for(int i = 0; i < 4; i++){
for(int j = 0; j < 4; j++){
for (int k = 0; k < 4; k++){
out[i*4 + j] += lhs[i*4 + k]*rhs[k*4 + i];
}
}     
}     
}

使用以下之一进行编译(衡量哪一个对您最快)

-O3 -mavx
-O3 -mavx2
-O3 -mavx2 -mfma
-O3 -mavx2 -mfma -ffast-math

这在 GCC 下有效(我的意思是矢量化)，cLANG 由于某种原因未能这样做。GCC 也将展开循环。