如何用 256 位 AVX 向量对两个复数进行平方？
How to square two complex doubles with 256-bit AVX vectors?
Matt Scarpino 对如何使用 Intel 的 AVX 内在函数将两个复数双精度相乘给出了很好的解释（尽管他承认他不确定这是最佳算法，但我向他表示感谢）。这是我验证过的他的方法：
__m256d vec1 = _mm256_setr_pd(4.0, 5.0, 13.0, 6.0);
__m256d vec2 = _mm256_setr_pd(9.0, 3.0, 6.0, 7.0);
__m256d neg = _mm256_setr_pd(1.0, -1.0, 1.0, -1.0);
/* Step 1: Multiply vec1 and vec2 */
__m256d vec3 = _mm256_mul_pd(vec1, vec2);
/* Step 2: Switch the real and imaginary elements of vec2 */
vec2 = _mm256_permute_pd(vec2, 0x5);
/* Step 3: Negate the imaginary elements of vec2 */
vec2 = _mm256_mul_pd(vec2, neg);
/* Step 4: Multiply vec1 and the modified vec2 */
__m256d vec4 = _mm256_mul_pd(vec1, vec2);
/* Horizontally subtract the elements in vec3 and vec4 */
vec1 = _mm256_hsub_pd(vec3, vec4);
/* Display the elements of the result vector */
double* res = (double*)&vec1;
printf("%lf %lf %lf %lf\n", res[0], res[1], res[2], res[3]);
我的问题是我想对两个复数双打求平方。我试着像这样使用 Matt 的技术：
struct cmplx a;
struct cmplx b;
a.r = 2.5341;
a.i = 1.843;
b.r = 1.3941;
b.i = 0.93;
__m256d zzs = squareZ(a, b);
double* res = (double*) &zzs;
printf("\nA: %f + %f, B: %f + %f\n", res[0], res[1], res[2], res[3]);
使用Haskell的复数运算，我验证了结果是正确的除了，如你所见，B的实部:
A: 3.025014 + 9.340693, B: 0.000000 + 2.593026
所以我真的有两个问题：是否有更好（更简单 and/or 更快）的方法来用 AVX 内在函数对两个复数双精度数求平方？如果没有，我该如何修改 Matt 的代码来做到这一点？
这个答案涵盖了将两个复数数组相乘的一般情况
理想情况下，将数据存储在单独的实部和虚部数组中，这样您就可以只加载实部和虚部的连续向量。这使得可以自由地进行交叉乘法（只需使用不同的寄存器/变量），而不必在向量中随机移动。
您可以相当便宜地在交错 double complex 样式和 SIMD 友好的分离向量样式之间即时转换，受 AVX 通道内洗牌的变幻莫测的影响。例如如果您不关心临时向量中数据的实际顺序，使用 unpacklo / unpackhi 洗牌可以非常便宜地去交错或在通道内重新交错。
执行这种洗牌实际上非常便宜，以至于在运行中为单个复数乘法执行它比 Matt 的代码（甚至是经过调整的版本）稍微领先，尤其是在支持 FMA 的 CPU 上。这需要以 4 个复数双打（2 个结果向量）为一组生成结果。
如果您一次只需要生成一个结果向量，我还想出了一种替代 Matt 算法的方法，它可以使用 FMA（实际上是 FMADDSUB）并避免单独的符号更改 insn.
gcc 自动将简单的复数乘法标量循环矢量化为非常好的代码，只要您使用 -ffast-math。它像我建议的那样去交错。
#include <complex.h>
// even with -ffast-math -ffp-contract=fast, clang doesn't manage to use vfmaddsubpd, instead using vmulpd and vaddsubpd :(
// gcc does use FMA though.
// auto-vectorizes with a lot of extra shuffles
void cmul(double complex *restrict dst,
const double complex *restrict A, const double complex *restrict B)
{ // clang and gcc change strategy slightly for i<1 or i<2, vs. i<4
for (int i=0; i<4 ; i++) {
dst[i] = A[i] * B[i];
}
}
参见Godbolt compiler explorer. I'm not sure how good clang's asm is; it uses a lot of 64b->128b VMODDDUP broadcast-loads. This form is handled purely in the load ports on Intel CPUs (see Agner Fog's insn tables)上的asm，不过还是有很多操作。如前所述，gcc 在乘法 / FMA 之前使用 4 个 VPERMPD 洗牌在通道内重新排序，然后在将结果与 VSHUFPD 组合之前使用另一个 4 VPERMPD 对结果重新排序。这是 4 次复数乘法的 8 次额外洗牌。
将 gcc 的版本转换回内在函数并删除冗余的 shuffle 得到最佳代码。 （gcc 显然希望其临时对象按 A B C D 顺序排列，而不是 VUNPCKLPD (_mm256_unpacklo_pd) 的车道内行为导致的 A C B D 顺序）。
我将代码 on Godbolt 与 Matt 代码的调整版本一起放置。因此，您可以使用不同的编译器选项以及不同的编译器版本。
// multiplies 4 complex doubles each from A and B, storing the result in dst[0..3]
void cmul_manualvec(double complex *restrict dst,
const double complex *restrict A, const double complex *restrict B)
{
// low element first, little-endian style
__m256d A0 = _mm256_loadu_pd((double*)A); // [A0r A0i A1r A1i ] // [a b c d ]
__m256d A2 = _mm256_loadu_pd((double*)(A+2)); // [e f g h ]
__m256d realA = _mm256_unpacklo_pd(A0, A2); // [A0r A2r A1r A3r ] // [a e c g ]
__m256d imagA = _mm256_unpackhi_pd(A0, A2); // [A0i A2i A1i A3i ] // [b f d h ]
// the in-lane behaviour of this interleaving is matched by the same in-lane behaviour when we recombine.
__m256d B0 = _mm256_loadu_pd((double*)B); // [m n o p]
__m256d B2 = _mm256_loadu_pd((double*)(B+2)); // [q r s t]
__m256d realB = _mm256_unpacklo_pd(B0, B2); // [m q o s]
__m256d imagB = _mm256_unpackhi_pd(B0, B2); // [n r p t]
// desired: real=rArB - iAiB, imag=rAiB + rBiA
__m256d realprod = _mm256_mul_pd(realA, realB);
__m256d imagprod = _mm256_mul_pd(imagA, imagB);
__m256d rAiB = _mm256_mul_pd(realA, imagB);
__m256d rBiA = _mm256_mul_pd(realB, imagA);
// gcc and clang will contract these nto FMA. (clang needs -ffp-contract=fast)
// Doing it manually would remove the option to compile for non-FMA targets
__m256d real = _mm256_sub_pd(realprod, imagprod); // [D0r D2r | D1r D3r]
__m256d imag = _mm256_add_pd(rAiB, rBiA); // [D0i D2i | D1i D3i]
// interleave the separate real and imaginary vectors back into packed format
__m256d dst0 = _mm256_shuffle_pd(real, imag, 0b0000); // [D0r D0i | D1r D1i]
__m256d dst2 = _mm256_shuffle_pd(real, imag, 0b1111); // [D2r D2i | D3r D3i]
_mm256_storeu_pd((double*)dst, dst0);
_mm256_storeu_pd((double*)(dst+2), dst2);
}
Godbolt asm output: gcc6.2 -O3 -ffast-math -ffp-contract=fast -march=haswell
vmovupd ymm0, YMMWORD PTR [rsi+32]
vmovupd ymm3, YMMWORD PTR [rsi]
vmovupd ymm1, YMMWORD PTR [rdx]
vunpcklpd ymm5, ymm3, ymm0
vunpckhpd ymm3, ymm3, ymm0
vmovupd ymm0, YMMWORD PTR [rdx+32]
vunpcklpd ymm4, ymm1, ymm0
vunpckhpd ymm1, ymm1, ymm0
vmulpd ymm2, ymm1, ymm3
vmulpd ymm0, ymm4, ymm3
vfmsub231pd ymm2, ymm4, ymm5 # separate mul/sub contracted into FMA
vfmadd231pd ymm0, ymm1, ymm5
vunpcklpd ymm1, ymm2, ymm0
vunpckhpd ymm0, ymm2, ymm0
vmovupd YMMWORD PTR [rdi], ymm1
vmovupd YMMWORD PTR [rdi+32], ymm0
vzeroupper
ret
对于（2 对输入向量的）4 个复数乘法，我的代码使用：
- 4 次加载（每次 32B）
- 2 家商店（每家 32B）
- 6 次通道内洗牌（每个输入向量一个，每个输出一个）
- 2 VMULPD
- 2 VFMA...某事
- （如果我们可以在分离的实数和图像向量中使用结果，则只有 4 次打乱，或者如果输入也已经采用这种格式，则为 0 次打乱）
- Intel Skylake 上的延迟（不计算 loads/stores）：14 个周期 = 4c 进行 4 次随机播放，直到第二个 VMULPD 可以启动 + 4 个周期（第二个 VMULPD）+ 4c（第二个 vfmadd231pd）+ 1c（首先随机播放result vector ready 1c earlier) + 1c (shuffle second result vector)
因此对于吞吐量，这完全是随机端口上的瓶颈。 （每个时钟吞吐量 1 次洗牌，而在 Intel Haswell 及更高版本上每个时钟总共 2 MUL/FMA/ADD）。这就是打包存储很糟糕的原因：洗牌的吞吐量有限，花更多的指令洗牌而不是做数学是不好的。
Matt Scarpino's 代码经过我的细微调整（重复执行 4 次复数乘法）。 （请参阅下文，了解我如何更有效地一次生成一个向量的重写）。
- 同6个loads/stores
- 6 次通道内洗牌（HSUBPD 是 2 次洗牌，是当前 Intel 和 AMD CPU 的减法）
- 4 倍
- 2 次减法（不能与 muls 组合成 FMA）
- 一个额外的指令（+一个常量）来翻转虚数元素的符号。 Matt 选择乘以 1.0 或 -1.0，但有效的选择是对符号位进行异或（即 XORPD 与
-0.0）。
- 第一个结果向量在 Intel Skylake 上的延迟：11 个周期。 1c（同一周期中的 vpermilpd 和 vxorpd）+ 4c（第二个 vmulpd）+ 6c（vhsubpd）。第一个 vmulpd 与其他操作重叠，与 shuffle 和 vxorpd 在同一周期开始。第二个结果向量的计算应该很好地交错。
Matt 代码的主要优点是它一次只处理一个矢量宽度的复数乘法，而不是要求您有 4 个输入数据矢量。它的延迟稍低。但请注意，我的版本不需要 2 对输入向量来自连续的内存，或者完全相互关联。它们在处理时混合在一起，但结果是 2 个独立的 32B 向量。
我调整后的 Matt 代码版本在没有 FMA 的 CPU 上几乎和 4-at-a-time 版本一样好（只是花费额外的 VXORPD），但当它阻止我们利用FMA。此外，它永远不会以非压缩形式提供结果，因此您不能将分离形式用作另一个乘法的输入并跳过改组。
一次一个矢量结果，使用 FMA：
如果您有时计算平方而不是乘以两个不同的复数，请不要使用此方法。这就像 Matt 的算法一样，公共子表达式消除不会简化。
我没有为此输入 C 内在函数，只是算出了数据移动。由于所有洗牌都是在车道内进行的，因此我只会展示低车道。使用相关指令的 256b 版本在两个通道中进行相同的洗牌。他们保持分开。
// MULTIPLY: for each AVX lane: am-bn, an+bm
r i r i
a b c d // input vectors: a + b*i, etc.
m n o p
算法:
用 movshdup(a b) + mulpd
创建 bm bn
用 shufpd 在上一个结果上创建 bn bm。 （或创建 n m 并在 mul 之前进行随机播放）
用 movsldup(a b)
创建 a a
使用 fmaddsubpd 产生最终结果：[a|a]*[m|n] -/+ [bn|bm].
是的，SSE/AVX 有 ADDSUBPD to do alternating subtract/add in even/odd elements (in that order, presumably because of this use-case). FMA includes FMADDSUB132PD 减法和加法（相反，FMSUBADD 加法和减法）。
每 4 个结果：6x 随机播放、2x mul、2xfmaddsub。因此，除非我弄错了，否则 它与去交错方法 一样有效（当不对相同的数字求平方时）。 Skylake 延迟 = 10c = 1+4+1 以创建 bn bm（与 1 个周期重叠以创建 a a），+ 4 (FMA)。所以它比 Matt 的延迟低一个周期。
在 Bulldozer 系列上，将两个输入混洗到第一个 mul 将是一个胜利，因此 mul->fmaddsub 关键路径保持在 FMA 域内（延迟低 1 个周期）。以另一种方式进行操作有助于阻止愚蠢的编译器通过过早执行 movsldup(a b) 并延迟 mulpd 来造成资源冲突。 （不过，在一个循环中，许多迭代将在随机端口上运行和瓶颈。）
这仍然比 Matt 的平方更好（仍然保存 XOR，并且可以使用 FMA），但我们不保存任何混洗：
// SQUARING: for each AVX lane: aa-bb, 2*ab
// ab bb // movshdup + mul
// bb ab // ^ -> shufpd
// a a // movsldup
// aa-bb ab+ab // fmaddsubpd : [a|a]*[a|b] -/+ [bb|ab]
// per 4 results: 6x shuffle, 2x mul, 2xfmaddsub
我也尝试过一些可能性，例如 (a+b) * (a+b) = aa+2ab+bb 或 (r-i)*(r+i) = rr - ii，但没有得到任何结果。步骤之间的舍入意味着 FP 数学不能完全抵消，所以做这样的事情甚至不会产生完全相同的结果。
请参阅我对不同复数相乘的一般情况的其他回答，而不是平方。
TL:DR：只需使用我的其他答案中的代码，两个输入都相同。编译器在冗余方面做得很好。
平方稍微简化了数学运算：不需要两个不同的叉积，rAiB 和 rBiA 是相同的。但它仍然需要加倍，所以基本上我们最终得到 2 mul + 1 FMA + 1 add，而不是 2 mul + 2 FMA。
对于 SIMD 不友好的交错存储格式，它极大地促进了去交错方法，因为只有一个输入要洗牌。 Matt 的方法根本没有任何好处，因为它计算两个具有相同向量乘法的叉积。
使用我的其他答案中的 cmul_manualvec()：
// squares 4 complex doubles from A[0..3], storing the result in dst[0..3]
void csquare_manual(double complex *restrict dst,
const double complex *restrict A) {
cmul_manualvec(dst, A, A);
}
gcc 和 clang 足够聪明，可以优化使用相同输入两次的冗余，因此无需制作带有内在函数的自定义版本。 clang 在标量自动矢量化版本上做得不好，所以不要使用它。我看不出这个 asm 输出有什么好处 (from Godbolt):
clang3.9 -O3 -ffast-math -ffp-contract=fast -march=haswell
vmovupd ymm0, ymmword ptr [rsi]
vmovupd ymm1, ymmword ptr [rsi + 32]
vunpcklpd ymm2, ymm0, ymm1
vunpckhpd ymm0, ymm0, ymm1 # doing this shuffle first would let the first multiply start a cycle earlier. Silly compiler.
vmulpd ymm1, ymm0, ymm0 # imag*imag
vfmsub231pd ymm1, ymm2, ymm2 # real*real - imag*imag
vaddpd ymm0, ymm0, ymm0 # imag+imag = 2*imag
vmulpd ymm0, ymm2, ymm0 # 2*imag * real
vunpcklpd ymm2, ymm1, ymm0
vunpckhpd ymm0, ymm1, ymm0
vmovupd ymmword ptr [rdi], ymm2
vmovupd ymmword ptr [rdi + 32], ymm0
vzeroupper
ret
可能不同的指令顺序会更好，以减少资源冲突。例如将实际向量加倍，因为它先解包，所以 VADDPD 可以在 imag*imag VMULPD 之前更快地开始一个循环。但是在 C 源代码中重新排序行通常不会直接转换为 asm 重新排序，因为现代编译器是复杂的野兽。 （IIRC，gcc 并没有特别尝试为 x86 安排指令，因为乱序执行主要隐藏了这些影响。）
总之，每 4 个复方块：
- 2 次装载（从 4 次减少）+ 2 家商店，原因显而易见
- 4 次洗牌（从 6 次减少），同样明显
- 2 VMULPD（相同）
- 1 FMA + 1 VADDPD（低于 2 FMA。VADDPD 在 Haswell/Broadwell 上的延迟低于 FMA，在 Skylake 上也是如此）。
马特的版本仍然是 6 次洗牌，其他一切都一样。
Matt Scarpino 对如何使用 Intel 的 AVX 内在函数将两个复数双精度相乘给出了很好的解释（尽管他承认他不确定这是最佳算法，但我向他表示感谢）。这是我验证过的他的方法：
__m256d vec1 = _mm256_setr_pd(4.0, 5.0, 13.0, 6.0);
__m256d vec2 = _mm256_setr_pd(9.0, 3.0, 6.0, 7.0);
__m256d neg = _mm256_setr_pd(1.0, -1.0, 1.0, -1.0);
/* Step 1: Multiply vec1 and vec2 */
__m256d vec3 = _mm256_mul_pd(vec1, vec2);
/* Step 2: Switch the real and imaginary elements of vec2 */
vec2 = _mm256_permute_pd(vec2, 0x5);
/* Step 3: Negate the imaginary elements of vec2 */
vec2 = _mm256_mul_pd(vec2, neg);
/* Step 4: Multiply vec1 and the modified vec2 */
__m256d vec4 = _mm256_mul_pd(vec1, vec2);
/* Horizontally subtract the elements in vec3 and vec4 */
vec1 = _mm256_hsub_pd(vec3, vec4);
/* Display the elements of the result vector */
double* res = (double*)&vec1;
printf("%lf %lf %lf %lf\n", res[0], res[1], res[2], res[3]);
我的问题是我想对两个复数双打求平方。我试着像这样使用 Matt 的技术：
struct cmplx a;
struct cmplx b;
a.r = 2.5341;
a.i = 1.843;
b.r = 1.3941;
b.i = 0.93;
__m256d zzs = squareZ(a, b);
double* res = (double*) &zzs;
printf("\nA: %f + %f, B: %f + %f\n", res[0], res[1], res[2], res[3]);
使用Haskell的复数运算，我验证了结果是正确的除了，如你所见，B的实部:
A: 3.025014 + 9.340693, B: 0.000000 + 2.593026
所以我真的有两个问题：是否有更好（更简单 and/or 更快）的方法来用 AVX 内在函数对两个复数双精度数求平方？如果没有，我该如何修改 Matt 的代码来做到这一点？
这个答案涵盖了将两个复数数组相乘的一般情况
理想情况下，将数据存储在单独的实部和虚部数组中，这样您就可以只加载实部和虚部的连续向量。这使得可以自由地进行交叉乘法（只需使用不同的寄存器/变量），而不必在向量中随机移动。
您可以相当便宜地在交错 double complex 样式和 SIMD 友好的分离向量样式之间即时转换，受 AVX 通道内洗牌的变幻莫测的影响。例如如果您不关心临时向量中数据的实际顺序，使用 unpacklo / unpackhi 洗牌可以非常便宜地去交错或在通道内重新交错。
执行这种洗牌实际上非常便宜，以至于在运行中为单个复数乘法执行它比 Matt 的代码（甚至是经过调整的版本）稍微领先，尤其是在支持 FMA 的 CPU 上。这需要以 4 个复数双打（2 个结果向量）为一组生成结果。
如果您一次只需要生成一个结果向量，我还想出了一种替代 Matt 算法的方法，它可以使用 FMA（实际上是 FMADDSUB）并避免单独的符号更改 insn.
gcc 自动将简单的复数乘法标量循环矢量化为非常好的代码，只要您使用 -ffast-math。它像我建议的那样去交错。
#include <complex.h>
// even with -ffast-math -ffp-contract=fast, clang doesn't manage to use vfmaddsubpd, instead using vmulpd and vaddsubpd :(
// gcc does use FMA though.
// auto-vectorizes with a lot of extra shuffles
void cmul(double complex *restrict dst,
const double complex *restrict A, const double complex *restrict B)
{ // clang and gcc change strategy slightly for i<1 or i<2, vs. i<4
for (int i=0; i<4 ; i++) {
dst[i] = A[i] * B[i];
}
}
参见Godbolt compiler explorer. I'm not sure how good clang's asm is; it uses a lot of 64b->128b VMODDDUP broadcast-loads. This form is handled purely in the load ports on Intel CPUs (see Agner Fog's insn tables)上的asm，不过还是有很多操作。如前所述，gcc 在乘法 / FMA 之前使用 4 个 VPERMPD 洗牌在通道内重新排序，然后在将结果与 VSHUFPD 组合之前使用另一个 4 VPERMPD 对结果重新排序。这是 4 次复数乘法的 8 次额外洗牌。
将 gcc 的版本转换回内在函数并删除冗余的 shuffle 得到最佳代码。 （gcc 显然希望其临时对象按 A B C D 顺序排列，而不是 VUNPCKLPD (_mm256_unpacklo_pd) 的车道内行为导致的 A C B D 顺序）。
我将代码 on Godbolt 与 Matt 代码的调整版本一起放置。因此，您可以使用不同的编译器选项以及不同的编译器版本。
// multiplies 4 complex doubles each from A and B, storing the result in dst[0..3]
void cmul_manualvec(double complex *restrict dst,
const double complex *restrict A, const double complex *restrict B)
{
// low element first, little-endian style
__m256d A0 = _mm256_loadu_pd((double*)A); // [A0r A0i A1r A1i ] // [a b c d ]
__m256d A2 = _mm256_loadu_pd((double*)(A+2)); // [e f g h ]
__m256d realA = _mm256_unpacklo_pd(A0, A2); // [A0r A2r A1r A3r ] // [a e c g ]
__m256d imagA = _mm256_unpackhi_pd(A0, A2); // [A0i A2i A1i A3i ] // [b f d h ]
// the in-lane behaviour of this interleaving is matched by the same in-lane behaviour when we recombine.
__m256d B0 = _mm256_loadu_pd((double*)B); // [m n o p]
__m256d B2 = _mm256_loadu_pd((double*)(B+2)); // [q r s t]
__m256d realB = _mm256_unpacklo_pd(B0, B2); // [m q o s]
__m256d imagB = _mm256_unpackhi_pd(B0, B2); // [n r p t]
// desired: real=rArB - iAiB, imag=rAiB + rBiA
__m256d realprod = _mm256_mul_pd(realA, realB);
__m256d imagprod = _mm256_mul_pd(imagA, imagB);
__m256d rAiB = _mm256_mul_pd(realA, imagB);
__m256d rBiA = _mm256_mul_pd(realB, imagA);
// gcc and clang will contract these nto FMA. (clang needs -ffp-contract=fast)
// Doing it manually would remove the option to compile for non-FMA targets
__m256d real = _mm256_sub_pd(realprod, imagprod); // [D0r D2r | D1r D3r]
__m256d imag = _mm256_add_pd(rAiB, rBiA); // [D0i D2i | D1i D3i]
// interleave the separate real and imaginary vectors back into packed format
__m256d dst0 = _mm256_shuffle_pd(real, imag, 0b0000); // [D0r D0i | D1r D1i]
__m256d dst2 = _mm256_shuffle_pd(real, imag, 0b1111); // [D2r D2i | D3r D3i]
_mm256_storeu_pd((double*)dst, dst0);
_mm256_storeu_pd((double*)(dst+2), dst2);
}
Godbolt asm output: gcc6.2 -O3 -ffast-math -ffp-contract=fast -march=haswell
vmovupd ymm0, YMMWORD PTR [rsi+32]
vmovupd ymm3, YMMWORD PTR [rsi]
vmovupd ymm1, YMMWORD PTR [rdx]
vunpcklpd ymm5, ymm3, ymm0
vunpckhpd ymm3, ymm3, ymm0
vmovupd ymm0, YMMWORD PTR [rdx+32]
vunpcklpd ymm4, ymm1, ymm0
vunpckhpd ymm1, ymm1, ymm0
vmulpd ymm2, ymm1, ymm3
vmulpd ymm0, ymm4, ymm3
vfmsub231pd ymm2, ymm4, ymm5 # separate mul/sub contracted into FMA
vfmadd231pd ymm0, ymm1, ymm5
vunpcklpd ymm1, ymm2, ymm0
vunpckhpd ymm0, ymm2, ymm0
vmovupd YMMWORD PTR [rdi], ymm1
vmovupd YMMWORD PTR [rdi+32], ymm0
vzeroupper
ret
对于（2 对输入向量的）4 个复数乘法，我的代码使用：
- 4 次加载（每次 32B）
- 2 家商店（每家 32B）
- 6 次通道内洗牌（每个输入向量一个，每个输出一个）
- 2 VMULPD
- 2 VFMA...某事
- （如果我们可以在分离的实数和图像向量中使用结果，则只有 4 次打乱，或者如果输入也已经采用这种格式，则为 0 次打乱）
- Intel Skylake 上的延迟（不计算 loads/stores）：14 个周期 = 4c 进行 4 次随机播放，直到第二个 VMULPD 可以启动 + 4 个周期（第二个 VMULPD）+ 4c（第二个 vfmadd231pd）+ 1c（首先随机播放result vector ready 1c earlier) + 1c (shuffle second result vector)
因此对于吞吐量，这完全是随机端口上的瓶颈。 （每个时钟吞吐量 1 次洗牌，而在 Intel Haswell 及更高版本上每个时钟总共 2 MUL/FMA/ADD）。这就是打包存储很糟糕的原因：洗牌的吞吐量有限，花更多的指令洗牌而不是做数学是不好的。
Matt Scarpino's 代码经过我的细微调整（重复执行 4 次复数乘法）。 （请参阅下文，了解我如何更有效地一次生成一个向量的重写）。
- 同6个loads/stores
- 6 次通道内洗牌（HSUBPD 是 2 次洗牌，是当前 Intel 和 AMD CPU 的减法）
- 4 倍
- 2 次减法（不能与 muls 组合成 FMA）
- 一个额外的指令（+一个常量）来翻转虚数元素的符号。 Matt 选择乘以 1.0 或 -1.0，但有效的选择是对符号位进行异或（即 XORPD 与
-0.0）。 - 第一个结果向量在 Intel Skylake 上的延迟：11 个周期。 1c（同一周期中的 vpermilpd 和 vxorpd）+ 4c（第二个 vmulpd）+ 6c（vhsubpd）。第一个 vmulpd 与其他操作重叠，与 shuffle 和 vxorpd 在同一周期开始。第二个结果向量的计算应该很好地交错。
Matt 代码的主要优点是它一次只处理一个矢量宽度的复数乘法，而不是要求您有 4 个输入数据矢量。它的延迟稍低。但请注意，我的版本不需要 2 对输入向量来自连续的内存，或者完全相互关联。它们在处理时混合在一起，但结果是 2 个独立的 32B 向量。
我调整后的 Matt 代码版本在没有 FMA 的 CPU 上几乎和 4-at-a-time 版本一样好（只是花费额外的 VXORPD），但当它阻止我们利用FMA。此外，它永远不会以非压缩形式提供结果，因此您不能将分离形式用作另一个乘法的输入并跳过改组。
一次一个矢量结果，使用 FMA：
如果您有时计算平方而不是乘以两个不同的复数，请不要使用此方法。这就像 Matt 的算法一样，公共子表达式消除不会简化。
我没有为此输入 C 内在函数，只是算出了数据移动。由于所有洗牌都是在车道内进行的，因此我只会展示低车道。使用相关指令的 256b 版本在两个通道中进行相同的洗牌。他们保持分开。
// MULTIPLY: for each AVX lane: am-bn, an+bm
r i r i
a b c d // input vectors: a + b*i, etc.
m n o p
算法:
用 movshdup(a b) + mulpd
创建bm bn用 shufpd 在上一个结果上创建
bn bm。 （或创建n m并在 mul 之前进行随机播放）用 movsldup(a b)
创建a a使用 fmaddsubpd 产生最终结果：
[a|a]*[m|n] -/+ [bn|bm].是的，SSE/AVX 有 ADDSUBPD to do alternating subtract/add in even/odd elements (in that order, presumably because of this use-case). FMA includes FMADDSUB132PD 减法和加法（相反，FMSUBADD 加法和减法）。
每 4 个结果：6x 随机播放、2x mul、2xfmaddsub。因此，除非我弄错了，否则 它与去交错方法 一样有效（当不对相同的数字求平方时）。 Skylake 延迟 = 10c = 1+4+1 以创建 bn bm（与 1 个周期重叠以创建 a a），+ 4 (FMA)。所以它比 Matt 的延迟低一个周期。
在 Bulldozer 系列上，将两个输入混洗到第一个 mul 将是一个胜利，因此 mul->fmaddsub 关键路径保持在 FMA 域内（延迟低 1 个周期）。以另一种方式进行操作有助于阻止愚蠢的编译器通过过早执行 movsldup(a b) 并延迟 mulpd 来造成资源冲突。 （不过，在一个循环中，许多迭代将在随机端口上运行和瓶颈。）
这仍然比 Matt 的平方更好（仍然保存 XOR，并且可以使用 FMA），但我们不保存任何混洗：
// SQUARING: for each AVX lane: aa-bb, 2*ab
// ab bb // movshdup + mul
// bb ab // ^ -> shufpd
// a a // movsldup
// aa-bb ab+ab // fmaddsubpd : [a|a]*[a|b] -/+ [bb|ab]
// per 4 results: 6x shuffle, 2x mul, 2xfmaddsub
我也尝试过一些可能性，例如 (a+b) * (a+b) = aa+2ab+bb 或 (r-i)*(r+i) = rr - ii，但没有得到任何结果。步骤之间的舍入意味着 FP 数学不能完全抵消，所以做这样的事情甚至不会产生完全相同的结果。
请参阅我对不同复数相乘的一般情况的其他回答，而不是平方。
TL:DR：只需使用我的其他答案中的代码，两个输入都相同。编译器在冗余方面做得很好。
平方稍微简化了数学运算：不需要两个不同的叉积，rAiB 和 rBiA 是相同的。但它仍然需要加倍，所以基本上我们最终得到 2 mul + 1 FMA + 1 add，而不是 2 mul + 2 FMA。
对于 SIMD 不友好的交错存储格式，它极大地促进了去交错方法，因为只有一个输入要洗牌。 Matt 的方法根本没有任何好处，因为它计算两个具有相同向量乘法的叉积。
使用我的其他答案中的 cmul_manualvec()：
// squares 4 complex doubles from A[0..3], storing the result in dst[0..3]
void csquare_manual(double complex *restrict dst,
const double complex *restrict A) {
cmul_manualvec(dst, A, A);
}
gcc 和 clang 足够聪明，可以优化使用相同输入两次的冗余，因此无需制作带有内在函数的自定义版本。 clang 在标量自动矢量化版本上做得不好，所以不要使用它。我看不出这个 asm 输出有什么好处 (from Godbolt):
clang3.9 -O3 -ffast-math -ffp-contract=fast -march=haswell
vmovupd ymm0, ymmword ptr [rsi]
vmovupd ymm1, ymmword ptr [rsi + 32]
vunpcklpd ymm2, ymm0, ymm1
vunpckhpd ymm0, ymm0, ymm1 # doing this shuffle first would let the first multiply start a cycle earlier. Silly compiler.
vmulpd ymm1, ymm0, ymm0 # imag*imag
vfmsub231pd ymm1, ymm2, ymm2 # real*real - imag*imag
vaddpd ymm0, ymm0, ymm0 # imag+imag = 2*imag
vmulpd ymm0, ymm2, ymm0 # 2*imag * real
vunpcklpd ymm2, ymm1, ymm0
vunpckhpd ymm0, ymm1, ymm0
vmovupd ymmword ptr [rdi], ymm2
vmovupd ymmword ptr [rdi + 32], ymm0
vzeroupper
ret
可能不同的指令顺序会更好，以减少资源冲突。例如将实际向量加倍，因为它先解包，所以 VADDPD 可以在 imag*imag VMULPD 之前更快地开始一个循环。但是在 C 源代码中重新排序行通常不会直接转换为 asm 重新排序，因为现代编译器是复杂的野兽。 （IIRC，gcc 并没有特别尝试为 x86 安排指令，因为乱序执行主要隐藏了这些影响。）
总之，每 4 个复方块：
- 2 次装载（从 4 次减少）+ 2 家商店，原因显而易见
- 4 次洗牌（从 6 次减少），同样明显
- 2 VMULPD（相同）
- 1 FMA + 1 VADDPD（低于 2 FMA。VADDPD 在 Haswell/Broadwell 上的延迟低于 FMA，在 Skylake 上也是如此）。
马特的版本仍然是 6 次洗牌，其他一切都一样。