AVX2：将 4 个复数值与 4 个双精度值相乘和求和的最佳方法是什么？

xnec2c项目中电磁模拟的一大热点采用这种形式，并且在整个计算过程中以各种方式重复相同的形式：

*dst += (*a1) * (*a2) + (*b1) * (*b2) + (*c1) * (*c2);  // + (*d1) * (*d2)

（

这个Github链接是matrix.c的真实示例）：

通常 a1、b1、c1 是复数。有时 a2、b2、c2 是复数，有时它们是实数双精度数。（我评论了 d1,d2，因为虽然它们在 EM 模拟中不存在，但我猜我们在 AVX 计算中免费获得它们，因此下面的代码编写了 4 个双精度数）。

以下是 Peter Cordes 回答的

修改版本：

*dst

static inline void avx_fma_4zd4d_sum(double complex * restrict dst,
  const double complex * restrict A, 
  double *b0, double *b1, double *b2, double *b3)
{
    // 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 ]

    // Q: b0-b3 are spread out in memory. Is this the best way to load them?
    __m256d B0 = _mm256_set_pd(*b1, *b1, *b0, *b0); // reverse from realA order because set_pd is.
    __m256d B2 = _mm256_set_pd(*b3, *b3, *b2, *b2); 

    // desired: real=rArB, imag=rBiA
    A0 = _mm256_mul_pd(A0, B0);
    A2 = _mm256_mul_pd(A2, B2);

    // sum A0-A3
    A0 = _mm256_add_pd(A0, A2);

    // Q: Now A0 has 2x complex doubles in a single vector.  How do I sum them?

    // Q: Once I do, how do I add the result to dst?

    //_mm256_storeu_pd((double *) dst, A0);
    //_mm256_storeu_pd((double *) (dst + 2), A2);
}

你可以在上面的评论中看到我的问题。我正在看

this answer，但它对所有双精度数的向量行求和，我的是复杂双精度数，所以它看起来像 [A0r, A0i, A1r, A1i] 并且总和需要交错。如果我跳出内在函数，结果将是这样，但我想了解此操作的内在细节：

dst += (A0[0]+A0[2]) + (A0[1]+A0[3]) * I;

请注意，

*dst

 不在内存中靠近

*A

，因此需要单独加载。

我没有 FMA 硬件，但我希望它能够编译，以便 gcc/clang 将减少为 FMA（如果可用）。

相关问题：

a2, b2, c2, d2

 很复杂怎么办？ （在向量函数中命名为b0-b3）

dst
```
 是否更容易？
```

a1, b1, c1

 放在同一个向量中，但是

a2, b2, c2

 都超出了内存。有没有比使用

_mm256_set_pd

 更好的方法来加载它们？

感谢您的帮助！

4
投票

这是您需要的一些部件。未经测试。

// 4 FP64 complex numbers in 2 AVX SIMD vectors
struct Complex4 { __m256d vec1, vec2; };

// Change the macro if you have FMA and optimizing for Intel
#define USE_FMA 0

Complex4 add( Complex4 a, Complex4 b )
{
    Complex4 res;
    res.vec1 = _mm256_add_pd( a.vec1, b.vec1 );
    res.vec2 = _mm256_add_pd( a.vec2, b.vec2 );
    return res;
}

// Multiply vectors component-wise.
// This is not a complex multiplication, just 2 vmulpd instructions
Complex4 mul_cwise( Complex4 a, Complex4 b )
{
    __m256d r0, r1;
    Complex4 res;
    // Compute the product
    res.vec1 = _mm256_mul_pd( a.vec1, b.vec1 );
    res.vec2 = _mm256_mul_pd( a.vec2, b.vec2 );
    return res;
}

// Multiply + accumulate vectors vectors component-wise
// This is not a complex multiplication.
Complex4 fma_cwise( Complex4 a, Complex4 b, Complex4 acc )
{
#if USE_FMA
    Complex4 res;
    res.vec1 = _mm256_fmadd_pd( a.vec1, b.vec1, acc.vec1 );
    res.vec2 = _mm256_fmadd_pd( a.vec2, b.vec2, acc.vec2 );
    return res;
#else
    return add( mul_cwise( a, b ), acc );
#endif
}

// Load 4 scalars from memory, and duplicate them into 2 AVX vectors
// This is for the complex * real use case
Complex4 load_duplicating( double* rsi )
{
    Complex4 res;
#ifdef __AVX2__
    __m256d tmp = _mm256_loadu_pd( rsi );
    res.vec1 = _mm256_permute4x64_pd( tmp, _MM_SHUFFLE( 1, 1, 0, 0 ) );
    res.vec2 = _mm256_permute4x64_pd( tmp, _MM_SHUFFLE( 3, 3, 2, 2 ) );
#else
    res.vec1 = _mm256_setr_m128d( _mm_loaddup_pd( rsi ), _mm_loaddup_pd( rsi + 1 ) );
    res.vec2 = _mm256_setr_m128d( _mm_loaddup_pd( rsi + 2 ), _mm_loaddup_pd( rsi + 3 ) );
#endif
    return res;
}

// Load 4 complex numbers
Complex4 load_c4( double* rsi )
{
    Complex4 res;
    res.vec1 = _mm256_loadu_pd( rsi );
    res.vec2 = _mm256_loadu_pd( rsi + 4 );
    return res;
}

// Multiply 2 complex numbers by another 2 complex numbers
__m256d mul_CxC_2( __m256d a, __m256d b )
{
    // The formula is [ a.x*b.x - a.y*b.y,  a.x*b.y + a.y*b.x ]
    // Computing it as a.xy * b.xx -+ a.yx * b.yy

    __m256d ayx = _mm256_permute_pd( a, _MM_SHUFFLE2( 0, 1 ) ); // a.yx
    __m256d byy = _mm256_permute_pd( b, _MM_SHUFFLE2( 1, 1 ) ); // b.yy
    __m256d p1 = _mm256_mul_pd( ayx, byy ); // a.yx * b.yy

    __m256d bxx = _mm256_permute_pd( b, _MM_SHUFFLE2( 0, 0 ) ); // b.xx

#if USE_FMA
    return _mm256_fmaddsub_pd( a, bxx, p1 );
#else
    return _mm256_addsub_pd( _mm256_mul_pd( a, bxx ), p1 );
#endif
}

还有更多注释。

考虑使用 C++ 来编写此类代码。其一，运算符和重载函数会有所帮助。此外，您的整个热点可以使用

Eigen 库用几行代码编写，其性能通常可与手动编写的内在函数相媲美。

使用编译器标志或 GCC 特定的函数属性来确保所有这些函数始终是内联的。

将 #2 的总和添加到 *dst

为这样的代码存储中间值从来都不是一个好主意。如果您有 6 个输入要相乘/相加，请将它们全部计算出来，并且只将结果存储一次。一般来说，内存访问效率很低，尤其是先写后读。

这样 gcc/clang 将减少为 FMA（如果可用）。

将事物简化为 FMA 对于 Intel 来说通常是个好主意，但对于 AMD 来说不一定。与 100% FMA 相比，AMD 芯片的 50% 加法/50% 乘法指令混合吞吐量高出两倍。与 Intel 不同，AMD 浮点加法和乘法不会竞争执行单元。我所说的“吞吐量”是指每秒指令数；营销人员决定 1 次 FMA 运算 = 2 次浮点运算，因此 FLOP 数是相同的。

问题描述投票：0回答：1

1个回答

最新问题

AVX2：将 4 个复数值与 4 个双精度值相乘和求和的最佳方法是什么？

问题描述 投票：0回答：1

1个回答

最新问题

问题描述投票：0回答：1