（如何）我可以使用 openmp 向量化 `std::complex<double>` 吗？

Question

我想使用矢量化优化我的应用程序。更具体地说，我想对 std::complex<double> 类型的数学运算进行矢量化。然而，这似乎相当困难。考虑以下示例：

#define TEST_LEN 100
#include <algorithm>
#include <complex>
typedef std::complex<double> cmplx;
using namespace std::complex_literals;

#pragma omp declare simd
cmplx add(cmplx a, cmplx b)
{
     return a + b;
}

#pragma omp declare simd
cmplx mult(cmplx a, cmplx b)
{
     return a * b;
}

void k(cmplx *x, cmplx *&y, int i0, int N)
{
    #pragma omp for simd
    for (int i = i0; i < N; i++)
        y[i] = add(mult(-(1i + 1.0), x[i]), 1i);
}

int main(int argc, char **argv)
{
    cmplx *x = new cmplx[TEST_LEN];
    cmplx *y = new cmplx[TEST_LEN];

    for (int i = 0; i < TEST_LEN; i++)
        x[i] = 0;

    for (int i = 0; i < TEST_LEN; i++)
    {
        int N = std::min(4, TEST_LEN - i);
        k(x, y, i, N);
    }

    delete[] x;
    delete[] y;

    return 1;
}

我正在使用 g++ 编译器。对于这段代码，编译器给出以下警告：

warning: unsupported return type 'cmplx' {aka 'std::complex'} for simd

包含 mult 和 add 函数的行。似乎无法像这样矢量化 std::complex<double> 类型。

有没有其他方法可以实现这一点？

Answer 1

不容易。当您在接下来的 N 个步骤中具有相同行为方式的值时，SIMD 工作得很好。因此，例如考虑一个二维向量数组：

X Y X Y X Y X Y

如果我们在这里做一个向量加法运算，

X Y X Y X Y X Y
+ + + + + + + +
X Y X Y X Y X Y

编译器会很好地向量化该操作。然而，如果我们想要对 X 和 Y 值做一些不同的事情，内存布局就会成为 SIMD 的问题：

X Y X Y X Y X Y
+ / + / + / + /
X Y X Y X Y X Y

如果你考虑乘法的例子：

(a + bi) (c + di) = (ac - bd)  (ad + bc)i

操作突然在 SIMD 通道之间跳跃，这几乎会扼杀任何合适的矢量化。

快速看一下这个神螺栓：https://godbolt.org/z/rnVVgl 加法归结为一些 vaddps 指令（一次处理 8 个浮点数）。 Multiply 最终使用 vfmadd231ss 和 vmulss（它们一次都在 1 个浮点数上工作）。

唯一简单自动向量化复杂代码的方法是将实部和虚部分成 2 个数组：

struct ComplexArray {
  float* real;
  float* imaginary;
};

在这个 godbolt 中，您可以看到编译器现在正在使用 vfmadd213ps 指令（所以再次回到一次处理 8 个浮点数）。

https://godbolt.org/z/Ostaax

（如何）我可以使用 openmp 向量化 `std::complex<double>` 吗？

(How) Can I vectorize `std::complex<double>` using openmp?

c++

vectorization

openmp