c++ - 有效地将大复数 vector 乘以标量 C++

Question

我目前正在尝试通过一组标量值最有效地执行复数数组的就地乘法(内存对齐方式与 std::complex 相同，但目前使用我们自己的 ADT)与复数数组大小相同。

该算法已经并行化，即调用对象将工作拆分为多个线程。此计算是在数以百万计的数组上完成的 - 因此，可能需要一些时间才能完成。 CUDA 不是该产品的解决方案，尽管我希望它是。我确实可以使用 boost，因此有可能使用 BLAS/uBLAS。

不过，我认为 SIMD 可能会产生更好的结果，但我对如何使用复数执行此操作还不够熟悉。我现在的代码如下(请记住，它被分成与目标机器上的内核数量相对应的线程)。目标机器也是未知的。因此，通用方法可能是最好的。

void cmult_scalar_inplace(fcomplex *values, const int start, const int end, const float *scalar)
{
    for (register int idx = start; idx < end; ++idx)
    {
        values[idx].real *= scalar[idx];
        values[idx].imag *= scalar[idx];
    }
}

fcomplex定义如下:

struct fcomplex
{
    float real;
    float imag;
};

我已经尝试手动展开循环，因为我的 finally 循环计数将始终是 2 的幂，但编译器已经为我这样做了(我已经展开到 32)。我已经尝试了对标量的 const float 引用——我想我会保存一次访问——结果证明它等于编译器已经在做的事情。我已经尝试过 STL 和变换，游戏结果很接近，但更糟。我还尝试转换为 std::complex 并允许它使用重载运算符进行标量 * 复数的乘法运算，但这最终产生了相同的结果。

那么，有人有什么想法吗？非常感谢您花时间考虑这个问题!目标平台是 Windows。我正在使用 Visual Studio 2008。产品也不能包含 GPL 代码!非常感谢。

Answer 1

您可以使用 SSE 很容易地做到这一点，例如

void cmult_scalar_inplace(fcomplex *values, const int start, const int end, const float *scalar)
{
    for (int idx = start; idx < end; idx += 2)
    {
        __m128 vc = _mm_load_ps((float *)&values[idx]);
        __m128 vk = _mm_set_ps(scalar[idx + 1], scalar[idx + 1], scalar[idx], scalar[idx]);
        vc = _mm_mul_ps(vc, vk);
        _mm_store_ps((float *)&values[idx], vc);
    }
}

请注意，values 和 scalar 需要 16 字节对齐。

或者您可以只使用英特尔 ICC 编译器，让它为您完成繁重的工作。

---

更新

这是一个改进版本，它将循环展开 2 倍，并使用单个加载指令获取 4 个标量值，然后将其解压缩为两个 vector :

void cmult_scalar_inplace(fcomplex *values, const int start, const int end, const float *scalar)
{
    for (int idx = start; idx < end; idx += 4)
    {
        __m128 vc0 = _mm_load_ps((float *)&values[idx]);
        __m128 vc1 = _mm_load_ps((float *)&values[idx + 2]);
        __m128 vk = _mm_load_ps(&scalar[idx]);
        __m128 vk0 = _mm_shuffle_ps(vk, vk, 0x50);
        __m128 vk1 = _mm_shuffle_ps(vk, vk, 0xfa);
        vc0 = _mm_mul_ps(vc0, vk0);
        vc1 = _mm_mul_ps(vc1, vk1);
        _mm_store_ps((float *)&values[idx], vc0);
        _mm_store_ps((float *)&values[idx + 2], vc1);
    }
}