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

问题描述

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

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_pd 相比，_mm256_setzero_si256 将整个函数减慢了一半。

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

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

解决方法

评论中已经提到了一些建议，但除此之外我会尝试提出一些建议。

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 帖子中有一些性能数据，以及那里链接的材料）。

,

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

void mul_64_block(double& sum,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