有效地将大复数向量乘以标量 C++

Question

我目前正在尝试最有效地对复数数组（内存对齐方式与 std::complex 相同，但当前使用我们自己的 ADT）与相同的标量值数组进行就地乘法size 作为复数数组。

该算法已经并行化，即调用对象将工作分成线程。此计算是在数百个数百万的数组上完成的 - 因此，可能需要一些时间才能完成。 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和transform，这两种游戏的结果很接近，但仍然更糟。我还尝试转换为 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);
    }
}

请注意，值和标量需要 16 字节对齐。

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

---

更新

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

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);
    }
}

Answer 2

最好的选择是使用优化的 BLAS，它将利用目标平台上可用的任何内容。

Answer 3

我看到的一个问题是，在函数中，编译器很难理解标量指针实际上并未指向复数数组的中间（标量理论上可能指向复数或实数）综合体的一部分）。
这实际上强制了评估的顺序。

我看到的另一个问题是，这里的计算非常简单，其他因素会影响原始速度，因此，如果您真的关心性能，我认为唯一的解决方案是实现几种变体并在用户计算机上运行时测试它们以发现什么是最快的。

我会考虑使用不同的展开大小，并调整标量和值的对齐（内存访问模式会对缓存效果产生很大影响）。

对于不需要的序列化问题，一个选项是查看生成的代码是什么，

float r0 = values[i].real, i0 = values[i].imag, s0 = scalar[i];
float r1 = values[i+1].real, i1 = values[i+1].imag, s1 = scalar[i+1];
float r2 = values[i+2].real, i2 = values[i+2].imag, s2 = scalar[i+2];
values[i].real = r0*s0; values[i].imag = i0*s0;
values[i+1].real = r1*s1; values[i+1].imag = i1*s1;
values[i+2].real = r2*s2; values[i+2].imag = i2*s2;

因为理论上优化器有更多的自由度。

Answer 4

您可以访问英特尔的集成性能基元吗？
集成性能原语 它们有许多函数可以处理诸如这具有相当不错的性能。您可能会在解决特定问题上取得一些成功，但如果您的编译器已经在优化代码方面做得很好，我不会感到惊讶。