AVX2:计算512个浮点阵的点积

Question

在此之前，我将说我是SIMD Intrinsics的一个完全的初学者。

实际上，我有一个支持AVX2 instrinsic的CPU(Intel(R)Core(TM)i5-7500T CPU@2.70GHz)。我想知道计算两个std::vector 大小512的点积的最快方法。

我在网上做了一些挖掘，发现了这个和这个，这个堆栈溢出问题建议使用以下函数__m256_mm256_dp_ps(__m256 m1,__m256 m2,const int mask)；，但是，这些都建议了不同的方法来执行点积，我不确定什么是正确的（也是最快的）方法。

编辑1：我对-mavx2gcc标志也有点困惑。如果我使用这些AVX2函数，我需要在编译时添加标志吗？另外，如果我编写一个简单的点产品实现，gcc是否能够为我进行这些优化（例如，如果我使用-ofastgcc标志）？

编辑2如果有人有时间和精力，我将非常感谢如果你能写一个完整的实现。我相信其他初学者也会重视这些信息。

Answer 1

_mm256_dp_ps只对2到4个元素的点积有用；对于更长的向量，在循环中使用垂直SIMD，并在结束时减少到标量。在循环中使用_mm256_dp_ps和_mm256_add_ps会慢得多。

与MSVC和ICC不同，GCC和clang要求您启用（使用命令行选项）使用intrinsics的ISA扩展。

下面的代码可能接近CPU的理论性能限制。未经测试。

#include <immintrin.h>
#include <vector>
#include <assert.h>

// CPUs support RAM access like this: "ymmword ptr [rax+64]"
// Using templates with offset int argument to make easier for compiler to emit good code.

// Multiply 8 floats by another 8 floats.
template<int offsetRegs>
inline __m256 mul8( const float* p1, const float* p2 )
{
    constexpr int lanes = offsetRegs * 8;
    const __m256 a = _mm256_loadu_ps( p1 + lanes );
    const __m256 b = _mm256_loadu_ps( p2 + lanes );
    return _mm256_mul_ps( a, b );
}

// Returns acc + ( p1 * p2 ), for 8-wide float lanes.
template<int offsetRegs>
inline __m256 fma8( __m256 acc, const float* p1, const float* p2 )
{
    constexpr int lanes = offsetRegs * 8;
    const __m256 a = _mm256_loadu_ps( p1 + lanes );
    const __m256 b = _mm256_loadu_ps( p2 + lanes );
    return _mm256_fmadd_ps( a, b, acc );
}

// Compute dot product of float vectors, using 8-wide FMA instructions.
float dotProductFma( const std::vector<float>& a, const std::vector<float>& b )
{
    assert( a.size() == b.size() );
    assert( 0 == ( a.size() % 32 ) );
    if( a.empty() )
        return 0.0f;

    const float* p1 = a.data();
    const float* const p1End = p1 + a.size();
    const float* p2 = b.data();

    // Process initial 32 values. Nothing to add yet, just multiplying.
    __m256 dot0 = mul8<0>( p1, p2 );
    __m256 dot1 = mul8<1>( p1, p2 );
    __m256 dot2 = mul8<2>( p1, p2 );
    __m256 dot3 = mul8<3>( p1, p2 );
    p1 += 8 * 4;
    p2 += 8 * 4;

    // Process the rest of the data.
    // The code uses FMA instructions to multiply + accumulate, consuming 32 values per loop iteration.
    // Unrolling manually for 2 reasons:
    // 1. To reduce data dependencies. With a single register, every loop iteration would depend on the previous result.
    // 2. Unrolled code checks for exit condition 4x less often, therefore more CPU cycles spent computing useful stuff.
    while( p1 < p1End )
    {
        dot0 = fma8<0>( dot0, p1, p2 );
        dot1 = fma8<1>( dot1, p1, p2 );
        dot2 = fma8<2>( dot2, p1, p2 );
        dot3 = fma8<3>( dot3, p1, p2 );
        p1 += 8 * 4;
        p2 += 8 * 4;
    }

    // Add 32 values into 8
    const __m256 dot01 = _mm256_add_ps( dot0, dot1 );
    const __m256 dot23 = _mm256_add_ps( dot2, dot3 );
    const __m256 dot0123 = _mm256_add_ps( dot01, dot23 );
    // Add 8 values into 4
    const __m128 r4 = _mm_add_ps( _mm256_castps256_ps128( dot0123 ), _mm256_extractf128_ps( dot0123, 1 ) );
    // Add 4 values into 2
    const __m128 r2 = _mm_add_ps( r4, _mm_movehl_ps( r4, r4 ) );
    // Add 2 lower values into the final result
    const __m128 r1 = _mm_add_ss( r2, _mm_movehdup_ps( r2 ) );
    // Return the lowest lane of the result vector.
    // The intrinsic below compiles into noop, modern compilers return floats in the lowest lane of xmm0 register.
    return _mm_cvtss_f32( r1 );
}

确保两个输入向量对齐，例如，使用自定义分配器，在msvc上调用_aligned_malloc/_aligned_free，或者在gcc&clang上调用aligned_alloc/free。然后将_mm256_loadu_ps替换为_mm256_load_ps。

要自动向量化一个简单的标量点积，还需要OpenMP SIMD或-ffast-math（由-ofast)表示），以便编译器将FP math视为关联，即使它不是关联的（因为舍入）。但是GCC在自动向量化时不会使用多个累加器，即使它确实展开，所以你会在FMA延迟上遇到瓶颈，而不是负载吞吐量。

（每个FMA加载2次意味着此代码的吞吐量瓶颈是向量加载，而不是实际的FMA操作。）