使用 SSE 缩放复杂向量

Posted 2023-02-16

【中文标题】使用 SSE 缩放复杂向量

【英文标题】：Scaling of a complex vector using SSE

【发布时间】：2016-05-04 15:53:38

【问题描述】：

我想将 SSE 指令应用于包含复数的向量。如果没有 SSE 说明，我可以使用以下代码来完成。但是，当我应用 SSE 指令时，我不知道如何将计算出的实部和虚部返回到数组中。我该如何解决这个问题？

complex double * complexScaling(complex double *input, double c, int length) 

  for(int i=0; i<length; i++) 
    input[i] = c*input[i];
  
  return input;


complex double * complexScalingSSE(complex double *input, double c, int length) 

  __m128d multiplier,real,imag;
  multiplier = _mm_set1_pd(c);
  for(int i=0; i<length; i+=2) 
    real = _mm_set_pd(creal(input[i]),creal(input[i+1]));
    real = _mm_mul_pd(real, multiplier);
    imag = _mm_set_pd(cimag(input[i]),cimag(input[i+1]));
    imag = _mm_mul_pd(imag, multiplier);
  
  return input;

【问题讨论】：

您的 SSE 级别是多少？ CPUID 指示了几个指令集：1, 2, 3, 4, 4.1, 4.2, AVX, AVX2, FMA3, FMA4, BMI, BMI2, ... 我问这个，因为用 FMAx 计算复数是一个蛋糕...

使用_mm_store*，_mm_storeu_pd 或 _mm_store_pd

必须兼容SSE2

任何看起来像set(var, var) 的东西都会产生非常烦人的代码。唯一安全的集合是带有常量参数。

【参考方案1】：

你可能想要这样的东西：

complex double * complexScalingSSE(complex double *input, double c, int length)

    const __m128d vc = _mm_set1_pd(c);

    for (int i = 0; i < length; i++)
    
        __m128d v = _mm_loadu_pd((double *)&input[i]);  // load one complex double
        v = _mm_mul_pd(v, vc);                          // scale it
        _mm_storeu_pd((double *)&input[i], v);          // store it
    
    return input;

您可能想尝试将循环展开 2 倍（但请注意，对于大型向量，此例程可能会受到 I/O 限制）：

complex double * complexScalingSSE(complex double *input, double c, int length)

    const __m128d vc = _mm_set1_pd(c);
    int i;

    for (i = 0; i <= length - 2; i += 2)
    
        __m128d v1 = _mm_loadu_pd((double *)&input[i]); // load two complex doubles
        __m128d v2 = _mm_loadu_pd((double *)&input[i + 1]);
        v1 = _mm_mul_pd(v1, vc);                        // scale them
        v2 = _mm_mul_pd(v2, vc);
        _mm_storeu_pd((double *)&input[i], v1);         // store them
        _mm_storeu_pd((double *)&input[i + 1], v2); 
    
    if (i < length)                                     // handle odd element at end
    
        __m128d v = _mm_loadu_pd((double *)&input[i]);  // load one complex double
        v = _mm_mul_pd(v, vc);                          // scale it
        _mm_storeu_pd((double *)&input[i], v);          // store it
    
    return input;

另请注意，体面的编译器应该自动矢量化您的原始例程，在这种情况下，您不会看到手动编码的 SSE 有任何好处。

【讨论】：

是的，gcc5 auto-vectorizes this just fine with SSE2，即使没有-ffast-math。不过，即使 gcc 4.9.2 似乎也无法管理它。 Clang 像往常一样将其展开合理的数量。 （仅 clang 3.6 和更新的自动矢量化。）

gcc 4.9.2 似乎需要 -ffast-math 才能对其进行矢量化（但尚不清楚为什么会这样）。

也许它只是害怕，因为有浮动，浮动优化很可怕

我将如何展开 2 倍？

@Nils：每次循环迭代只需执行两次加载、两次乘法和两次存储（使用两个单独的临时变量，例如 v1、v2），并将 i 增加 2 而不是 1。