c - 如何用256位AVX vector 对两个复数 double 平方？

Question

Matt Scarpino提供了一个很好的解释(尽管他承认他不确定这是否是最佳算法，但我要感谢他)，以了解如何将两个复杂的double与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

因此，我确实有两个问题:是否有更好的(更简单和/或更快速)将两个复数 double 函数与AVX内在函数平方的方法？如果没有，我该如何修改Matt的代码呢？

Answer 1

有关将不同的复数相乘而不是平方的一般情况，请参见我的其他答案。

TL:DR :只需在我的其他答案中使用代码，两个输入都相同。编译器在冗余方面做得很好。

平方稍微简化了数学运算:rAiB和rBiA不需要相同的叉积，而无需两个不同的叉积。但是它仍然需要加倍，因此基本上我们最终得到2 mul +1 FMA +1添加，而不是2 mul + 2 FMA。

使用SIMD不友好的交错存储格式，由于只有一个输入可以随机播放，因此大大提高了解交错方法。 Matt的方法完全没有好处，因为它使用相同的 vector 乘积来计算两个叉积。

使用我的其他答案中的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

可能有不同的指令排序会更好，以减少资源冲突。例如因为先解压缩了实 vector ，所以它加倍了，因此VADDPD可以在imag * imag VMULPD之前更快地开始一个周期。但是，C语言中的重新排序行通常不会直接转换为asm重新排序，因为现代编译器是复杂的野兽。 (IIRC，gcc并未特别尝试为x86安排指令，因为乱序执行通常会隐藏这些效果。)

无论如何，每4个复杂正方形:

2个加载项(从4个减少)+ 2个商店，显而易见的原因

4洗牌(从6减少)，再次明显

2 VMULPD(相同)

1 FMA + 1 VADDPD(低于2FMA。VADDPD的延迟比Haswell/Broadwell的FMA更低，与Skylake相同)。

马特(Matt)的版本仍将是6个改组，其他都相同。