4

我目前正在尝试最有效地对一个复数数组（内存对齐方式与 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 并允许它使用重载运算符进行 scalar * complex 进行乘法运算，但这最终产生了相同的结果。

那么，有任何想法的人吗？非常感谢您花时间考虑这一点！目标平台是 Windows。我使用的是 Visual Studio 2008。产品也不能包含 GPL 代码！非常感谢。

4

4 回答 4

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 个标量值，然后将其解包为两个向量：

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);
    }
}
于 2011-07-28T19:42:52.737 回答
1

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

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

我会考虑使用不同的展开大小，并使用和的对齐方式scalar（values内存访问模式可能会对缓存效果产生很大影响）。

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

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;

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

于 2011-07-28T19:57:34.927 回答
1

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

于 2011-07-28T18:57:12.897 回答
0

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

于 2011-07-28T18:53:47.297 回答