c - 如何提高以下循环的性能

Question

我在 C 中有一个简单的循环，我将 magnitude 和 angle 转换为 real 和 虚数部分。我有两个版本的循环。 版本 1 是一个简单的 for 循环，我使用以下代码执行转换

for(k = 0; k < n; k++){
    xReal[k] = Mag[k] * cos(Angle[k]);
    xImag[k] = Mag[k] * sin(Angle[k]);
}

版本 2，其中 Intrinsics 用于矢量化循环。

__m256d cosVec, sinVec;
__m256d resultReal, resultImag;
__m256d angVec, voltVec;
for(k = 0; k < SysData->totNumOfBus; k+=4){

    voltVec = _mm256_loadu_pd(volt + k);
    angVec = _mm256_loadu_pd(theta + k);

    sinVec = _mm256_sincos_pd(&cosVec, angVec);

    resultImag = _mm256_mul_pd(voltVec, sinVec);
    resultReal = _mm256_mul_pd(voltVec, cosVec);

    _mm256_store_pd(xReal+k, resultReal);
    _mm256_store_pd(xImag+k, resultImag);

}

在 Core i7 2600k @3.4GHz 处理器上，这些循环给出以下结果:

Version 1: n = 18562320, Time: 0.2sec
Version 2: n = 18562320, Time: 0.16sec

使用这些值进行的简单计算表明，在版本 1 中，每次迭代几乎需要 36 个周期才能完成，而它需要 117 个周期版本 2 完成。考虑到 sine 和 cosine 函数的计算自然是昂贵的，这些数字似乎并不可怕。然而，这个循环是我的函数的一个严重瓶颈，因为分析显示几乎 1/3 的时间花在了循环内。所以，我想知道是否有任何方法可以加快这个循环(例如，以不同方式计算 sine 和 cosine 函数)。如果能帮助我解决这个问题并让我知道是否有改进此循环性能的空间，我们将不胜感激。

预先感谢您的帮助

PS:我正在使用icc 来编译代码。另外，我应该提到数据没有对齐(也不能对齐)。但是，对齐数据只会带来较小的性能提升(不到 1%)。

Answer 1

我建议制作基于泰勒级数的 sin/cos 函数和 _mm256_stream_pd() 来存储数据。这是基本示例代码。

    __m256d sin_req[10];
    __m256d cos_req[10];
    __m256d one_pd =  _mm256_set1_pd(1.0);

    for(int i=0; i<10; ++i)
    {
        sin_req[i] = i%2 == 0 ? _mm256_set1_pd(-1.0/Factorial((i+1)*2+1) ) : _mm256_set1_pd(+1.0/Factorial((i+1)*2+1) );
        cos_req[i] = i%2 == 0 ? _mm256_set1_pd(-1.0/Factorial((i+1)*2+0) ) : _mm256_set1_pd(+1.0/Factorial((i+1)*2+0) );
    }

    for(int i=0; i<count; i+=4)
    {
            __m256d voltVec = _mm256_load_pd(volt + i);
            __m256d angVec = _mm256_load_pd(theta + i);

            // sin/cos by taylor series
            __m256d angleSq = angVec * angVec;
            __m256d sinVec = angVec;
            __m256d cosVec = one_pd;
            __m256d sin_serise = sinVec;
            __m256d cos_serise = one_pd;
            for(int j=0; j<10; ++j)
            {
                sin_serise = sin_serise * angleSq; // [1]
                cos_serise = cos_serise * angleSq;
                sinVec = sinVec + sin_serise * sin_req[j];
                cosVec = cosVec + cos_serise * cos_req[j];
            }

            __m256d resultReal = voltVec * sinVec;
            __m256d resultImag = voltVec * cosVec;

            _mm256_store_pd(xReal + i, resultReal);
            _mm256_store_pd(xImag + i, resultImag );
    }

我可以获得 57~58 个 CPU 周期来计算 4 个分量。

我搜索了谷歌并运行了一些测试来确定我的 sin/cos 的准确性。一些文章说 10 次迭代是 double 精确的，而 -M_PI/2 < angle < +M_PI/2。我的测试结果表明它在 -M_PI < angle < +M_PI 范围内比 math.h 的 sin/cos 更准确。如果需要，您可以增加迭代以提高大角度值的准确性。

不过，我会更深入地优化这段代码。此代码有延迟问题计算泰尔级数。 AVX 的乘法延迟是 5 个 CPU 周期，这意味着我们不能以快于 5 个周期的速度运行一次迭代，因为 [1] 使用的是前一次迭代结果的结果。

我们可以像这样简单地展开它。

    for(int i=0; i<count; i+=8)
    {
        __m256d voltVec0 = _mm256_load_pd(volt + i + 0);
        __m256d voltVec1 = _mm256_load_pd(volt + i + 4);
        __m256d angVec0  = _mm256_load_pd(theta + i + 0);
        __m256d angVec1  = _mm256_load_pd(theta + i + 4);
        __m256d sinVec0;
        __m256d sinVec1;
        __m256d cosVec0;
        __m256d cosVec1;

        __m256d angleSq0 = angVec0 * angVec0;
        __m256d angleSq1 = angVec1 * angVec1;
        sinVec0 = angVec0;
        sinVec1 = angVec1;
        cosVec0 = one_pd;
        cosVec1 = one_pd;
        __m256d sin_serise0 = sinVec0;
        __m256d sin_serise1 = sinVec1;
        __m256d cos_serise0 = one_pd;
        __m256d cos_serise1 = one_pd;

        for(int j=0; j<10; ++j)
        {
            sin_serise0 = sin_serise0 * angleSq0;
            cos_serise0 = cos_serise0 * angleSq0;
            sin_serise1 = sin_serise1 * angleSq1;
            cos_serise1 = cos_serise1 * angleSq1;
            sinVec0 = sinVec0 + sin_serise0 * sin_req[j];
            cosVec0 = cosVec0 + cos_serise0 * cos_req[j];
            sinVec1 = sinVec1 + sin_serise1 * sin_req[j];
            cosVec1 = cosVec1 + cos_serise1 * cos_req[j];
        }

        __m256d realResult0 = voltVec0 * sinVec0;
        __m256d imagResult0 = voltVec0 * cosVec0;
        __m256d realResult1 = voltVec1 * sinVec1;
        __m256d imagResult1 = voltVec1 * cosVec1;

        _mm256_store_pd(xReal + i + 0, realResult0);
        _mm256_store_pd(xImag + i + 0, imagResult0);
        _mm256_store_pd(xReal + i + 4, realResult1);
        _mm256_store_pd(xImag + i + 4, imagResult1);
    }

这个结果51~51.5 cycles for 4 component calculation。 (8个分量102~103个循环)

它消除了泰勒计算循环中的乘法延迟，并使用了 85% 的 AVX 乘法单元。展开将解决很多延迟问题，同时它不会将寄存器交换到内存。编译时生成 asm 文件，看看你的编译器如何处理你的代码。我尝试展开更多，但结果很糟糕，因为它无法容纳 16 个 AVX 寄存器。

现在我们进行内存优化。将 _mm256_store_ps() 替换为 _mm256_stream_ps()。

    _mm256_stream_pd(xReal + i + 0, realResult0);
    _mm256_stream_pd(xImag + i + 0, imagResult0);
    _mm256_stream_pd(xReal + i + 4, realResult1);
    _mm256_stream_pd(xImag + i + 4, imagResult1);

替换内存写入代码结果 48 个周期用于 4 分量计算。

如果您不打算回读它，_mm256_stream_pd() 总是更快。它跳过缓存系统并将数据直接发送到内存 Controller 并且不会污染您的缓存。通过使用 _mm256_stream_pd()，您将获得更多的数据总线/缓存空间来读取数据。

让我们试试预取。

    for(int i=0; i<count; i+=8)
    {
    _mm_prefetch((const CHAR *)(volt + i + 5 * 8), _MM_HINT_T0);
    _mm_prefetch((const CHAR *)(theta + i + 5 * 8), _MM_HINT_T0);

            // calculations here.
    }

现在每次计算我得到 45.6~45.8 个 CPU 周期。 AVX 乘法单元忙于 94%。

Prefech 提示缓存以加快读取速度。我建议根据物理内存的 RAS-CAS 延迟在 400~500 个 CPU 周期之前预加载。在最坏的情况下，物理内存延迟最多可能需要 300 个周期。可能因硬件配置而异，即使您使用昂贵的低 RAS-CAS 延迟内存，也不会小于 200 个周期。

0.064 秒(计数 = 18562320)

sin/cos 优化结束。 :-)