使用 SIMD 提高数组的浮点点积的性能

Question

我有这个函数来计算一个双精度数组：

void avx2_mul_64_block(double& sum, double* lhs_arr, double* rhs_arr) noexcept
{
    __m256i accumulator = _mm256_setzero_pd();

    for (std::size_t block = 0; block < 64; block += 4)
    {
        __m256i lhs = _mm256_set_pd(
            lhs_arr[block    ],
            lhs_arr[block + 1],
            lhs_arr[block + 2],
            lhs_arr[block + 3]);

        __m256i rhs = _mm256_set_pd(
            rhs_arr[block    ],
            rhs_arr[block + 1],
            rhs_arr[block + 2],
            rhs_arr[block + 3]);

        accumulator = _mm256_add_pd(accumulator, _mm256_mul_pd(lhs, rhs));
    }
    double* res = reinterpret_cast<double*>(&accumulator);
    sum += res[0] + res[1] + res[2] + res[3];
}

，这段代码的性能不是我想要的。我相信让它成功 - 避免为其所有元素创建双数组，但我不知道如何做到这一点。

顺便说一句，与_mm256_setzero_si256相比，_mm256_setzero_pd 将整个功能减慢了一半。

我的标志：-O3 -ftree-vectorize -march=native

附：这不是真正的问题，只是设计问题。

Answer 1

一些建议已经在 cmets 中提到，但我会尝试提出一些建议。

Haswell 和更新的 CPU 上的大多数 SIMD 浮点指令的倒数吞吐量都小于其延迟，这意味着如果并行执行多条指令，则可以提高性能。例如，根据 Agner Fog 的instruction tables，vaddpd 在 Haswell 上的延迟为 3，倒数吞吐量为 1 个时钟周期，这意味着 CPU 可以潜在地并行执行 3 条指令。更多的vmulpd 指令可以在其 5 和 0.5 个时钟周期内并行执行，以实现延迟和互惠吞吐量。

您的代码可能没有利用这种指令级并行性 (ILP)，因为循环体依赖于在上一次循环迭代中更新的 accumulator 值。这是因为不允许编译器执行许多优化，例如对 FP 数的数学运算重新排序，因为这可能会导致数学上不同的结果。因此，您的算法会受到延迟限制。

您可以通过使用特定于编译器的选项来缓解这种情况，例如用于 gcc 和兼容编译器的 -ffast-math，但考虑到这种并行性，只重写您的算法更便于移植。

我还将在下面的代码中加入其他建议，例如：

• 修复了错误的向量类型，应该是__m256d。
• 使用专用指令从内存中加载整个向量，而不是依赖编译器优化_mm256_set_pd。
• 使用 FMA 内在函数而不是依赖编译器优化 _mm256_add_pd+_mm256_mul_pd 对。 FMA 指令减少了计算的延迟，使加法有效地免费。 FMA 还产生更精确的结果，因为乘法和加法之间没有舍入。请注意，FMA 需要 AVX2，它不适用于仅支持 AVX 的 CPU。
• 使用适当的内在函数从向量中提取最终总和（由于 double 无论如何都存储在向量寄存器中，这可能会在最终汇编程序中被优化掉）。

void avx2_mul_64_block(double& sum, double* lhs_arr, double* rhs_arr) noexcept
{
    __m256d accumulator1 = _mm256_setzero_pd();
    __m256d accumulator2 = _mm256_setzero_pd();

    for (std::size_t block = 0; block < 64; block += 4 * 2)
    {
        __m256d lhs1 = _mm256_loadu_pd(lhs_arr + block);
        __m256d lhs2 = _mm256_loadu_pd(lhs_arr + block + 4);
        __m256d rhs1 = _mm256_loadu_pd(rhs_arr + block);
        __m256d rhs2 = _mm256_loadu_pd(rhs_arr + block + 4);

        accumulator1 = _mm256_fmadd_pd(lhs1, rhs1, accumulator1);
        accumulator2 = _mm256_fmadd_pd(lhs2, rhs2, accumulator2);
    }

    accumulator1 = _mm256_add_pd(accumulator1, accumulator2);

    __m128d accumulator = _mm_add_pd(_mm256_castpd256_pd128(accumulator1),
        _mm256_extractf128_pd(accumulator1, 1));
    accumulator = _mm_add_pd(accumulator,
        _mm_unpackhi_pd(accumulator, accumulator));

    sum += _mm_cvtsd_f64(accumulator);
}

在上面的代码中，我使用了两个单独的累加器，因此 CPU 现在能够并行执行两条累加链。进一步提高并行性可能是有益的（参见上面提到的性能数据），但如果块长度不能被累加器的数量乘以向量中的元素数量整除，则问题可能会更大。您可能需要设置尾部处理，这可能会产生一些轻微的性能开销。

请注意，如前所述，由于算术运算和 FMA 的不同顺序以及因此数学误差的不同累积，此算法可能会产生不严格等于原始结果的结果。但是，这通常不是问题，尤其是double 的高精度。

---

上面代码中没有用到的一些建议：

• 未使用水平加法 (_mm256_hadd_pd) 来累积最终总和，因为在当前的 Intel（最高 Coffee Lake）和 AMD（最高 Zen 2）处理器上，vhadd 指令的延迟比vunpckhpd 稍长+vaddpd 对，即使后者有依赖链。这可能会在未来的处理器中发生变化，并且使用水平加法可能会变得有益。不过，在当前 CPU 中，水平相加可能有助于节省一些代码大小。
• 代码使用来自内存的未对齐负载，而不是_mm256_load_pd，因为支持 AVX2 的现代 CPU 上的未对齐负载不会有性能损失，只要负载实际上在运行时对齐或至少不对齐跨缓存线边界。跨越缓存线，尤其是页面边界时会有开销，但通常这仍然不会降低现代 CPU 上的性能（在 Peter Cordes 的this 帖子中有一些性能数据，如以及那里链接的材料）。

Answer 2

正如 Marc Glisse 在上面的 cmets 中提到的，您可能希望在编译器标志中设置 -ffast-math。这是编译器很容易优化的函数之一，最好直接用 C++ 编写代码。

void mul_64_block(double& sum, double* lhs_arr, double* rhs_arr) {
    double res = 0;
    for(int i = 0; i < 64; ++i) {
        res += lhs_arr[i] * rhs_arr[i];
    }
    sum += res;
}

此 C++ 代码产生与您的 simd 代码相同的输出。

https://godbolt.org/z/cddPMd