Intel AVX：用于双精度浮点变量的 256 位点积版本答案

【问题标题】：Intel AVX: 256-bits version of dot product for double precision floating point variablesIntel AVX：用于双精度浮点变量的 256 位点积版本
【发布时间】：2012-05-14 07:29:33
【问题描述】：

英特尔高级矢量扩展 (AVX) 在 256 位版本（YMM 寄存器）中不为双精度浮点变量提供点积。 “为什么？”这个问题已经在另一个论坛 (here) 和 StackOverflow (here) 上得到了非常简短的处理。但我面临的问题是如何以有效的方式用其他 AVX 指令替换这条缺失的指令？

256位版本的点积存在单精度浮点变量（reference here）：

 __m256 _mm256_dp_ps(__m256 m1, __m256 m2, const int mask);

我们的想法是为这个缺失的指令找到一个有效的等价物：

 __m256d _mm256_dp_pd(__m256d m1, __m256d m2, const int mask);

更具体地说，我想从__m128（四个浮点数）转换为__m256d（四个双精度数）的代码使用以下说明：

   __m128 val0 = ...; // Four float values
   __m128 val1 = ...; //
   __m128 val2 = ...; //
   __m128 val3 = ...; //
   __m128 val4 = ...; //

   __m128 res = _mm_or_ps( _mm_dp_ps(val1,  val0,   0xF1),
                _mm_or_ps( _mm_dp_ps(val2,  val0,   0xF2),
                _mm_or_ps( _mm_dp_ps(val3,  val0,   0xF4),
                           _mm_dp_ps(val4,  val0,   0xF8) )));

此代码的结果是一个由四个浮点数组成的_m128 向量，其中包含val1 和val0、val2 和val0、val3 和val0、@ 之间的点积结果987654336@和val0。

也许这可以为建议提供提示？

【问题讨论】：

感谢您的想法，但我应该在我的应用程序中保持双精度。
此外，转换+浮点乘积比双点乘积需要更多时间。

标签： c++ performance simd avx

【解决方案1】：

我会使用 4*double 乘法，然后是 hadd（不幸的是，它只在上半部分和下半部分添加了 2*2 浮点数），提取上半部分（随机播放应该同样有效，也许更快）并添加它到下半部分。

结果在dotproduct的低64位。

__m256d xy = _mm256_mul_pd( x, y );
__m256d temp = _mm256_hadd_pd( xy, xy );
__m128d hi128 = _mm256_extractf128_pd( temp, 1 );
__m128d dotproduct = _mm_add_pd( (__m128d)temp, hi128 );

编辑：
根据 Norbert P 的想法，我扩展了这个版本，一次做 4 个点积。

__m256d xy0 = _mm256_mul_pd( x[0], y[0] );
__m256d xy1 = _mm256_mul_pd( x[1], y[1] );
__m256d xy2 = _mm256_mul_pd( x[2], y[2] );
__m256d xy3 = _mm256_mul_pd( x[3], y[3] );

// low to high: xy00+xy01 xy10+xy11 xy02+xy03 xy12+xy13
__m256d temp01 = _mm256_hadd_pd( xy0, xy1 );   

// low to high: xy20+xy21 xy30+xy31 xy22+xy23 xy32+xy33
__m256d temp23 = _mm256_hadd_pd( xy2, xy3 );

// low to high: xy02+xy03 xy12+xy13 xy20+xy21 xy30+xy31
__m256d swapped = _mm256_permute2f128_pd( temp01, temp23, 0x21 );

// low to high: xy00+xy01 xy10+xy11 xy22+xy23 xy32+xy33
__m256d blended = _mm256_blend_pd(temp01, temp23, 0b1100);

__m256d dotproduct = _mm256_add_pd( swapped, blended );

【讨论】：

感谢您的建议，这很好。为了更具体，我已经编辑了我的问题。
谢谢！你能解释一下最后一行吗？我不敢很好地理解它。不是_mm256_add_pd吗？
@drhirsch：好主意。但是 gleeen.gould 是对的，你需要一个额外的洗牌。我推荐：__m256d swapped = _mm256_permute2f128_pd( temp01, temp23, 0x21 ); __m256d mixed = _mm256_blend_pd(temp01, temp23, 12); __m256d dotproduct = _mm256_add_pd( swapped, mixed );。唯一的原因是 VPERM2F128 需要 2 个周期，而 VBLENDPD 需要 1 个周期。（希望我得到了正确的常量）
@gleeen.gould：AVX2 已经出来了吗？我以为是coming to Haswell in 2013。我在谈论 Sandy Bridge 上的当前一代 AVX：请参阅 Agner Fog's instruction tables，第 1 页。 129.跨度>
注意：即使使用AVX2 和vpermpd，使用vextractf128 和addpd 的解决方案也比随后应用vpermpd 和vhaddpd 来生成水平总和。

【解决方案2】：

我会扩展 drhirsch's answer 以同时执行两个点积，节省一些工作：

__m256d xy = _mm256_mul_pd( x, y );
__m256d zw = _mm256_mul_pd( z, w );
__m256d temp = _mm256_hadd_pd( xy, zw );
__m128d hi128 = _mm256_extractf128_pd( temp, 1 );
__m128d dotproduct = _mm_add_pd( (__m128d)temp, hi128 );

那么dot(x,y) 是低双，dot(z,w) 是dotproduct 的高双。

【讨论】：

【解决方案3】：

对于单个点积，它只是一个垂直乘法和水平和（参见Fastest way to do horizontal float vector sum on x86）。 hadd 花费 2 次随机播放 + 一个 add。当使用两个输入 = 相同的向量时，吞吐量几乎总是次优。

// both elements = dot(x,y)
__m128d dot1(__m256d x, __m256d y) {
    __m256d xy = _mm256_mul_pd(x, y);

    __m128d xylow  = _mm256_castps256_pd128(xy);   // (__m128d)cast isn't portable
    __m128d xyhigh = _mm256_extractf128_pd(xy, 1);
    __m128d sum1 =   _mm_add_pd(xylow, xyhigh);

    __m128d swapped = _mm_shuffle_pd(sum1, sum1, 0b01);   // or unpackhi
    __m128d dotproduct = _mm_add_pd(sum1, swapped);
    return dotproduct;
}

如果您只需要一个点积，这比 @hirschhornsalz 的单向量答案要好，在 Intel 上 1 shuffle uop，并且在 AMD Jaguar / Bulldozer-family / Ryzen 上取得更大的胜利，因为它立即缩小到 128b做一堆256b的东西。 AMD 将 256b 的 ops 拆分为两个 128b uop。

hadd 在并行处理 2 个或 4 个点积的情况下是值得使用的，其中您将它与 2 个不同的输入向量一起使用。如果您想要打包结果，Norbert 的两对向量的dot 看起来是最佳的。即使将 AVX2 vpermpd 用作车道交叉洗牌，我也看不出有任何方法可以做得更好。

当然，如果您真的想要更大的dot（8 个或更多doubles），请使用垂直add（使用多个累加器来隐藏vaddps 延迟）并在以下位置进行水平求和结束。如果可用，您也可以使用fma。

haddpd 在内部以两种不同的方式将xy 和zw 混洗在一起，并将其提供给垂直的addpd，无论如何，这就是我们要手动执行的操作。如果我们将 xy 和 zw 分开，我们需要 2 次洗牌 + 2 次加法来获得一个点积（在单独的寄存器中）。因此，第一步将它们与hadd 一起洗牌，我们节省了洗牌的总数，仅减少了添加数和微指令总数。

/*  Norbert's version, for an Intel CPU:
    __m256d temp = _mm256_hadd_pd( xy, zw );   // 2 shuffle + 1 add
    __m128d hi128 = _mm256_extractf128_pd( temp, 1 ); // 1 shuffle (lane crossing, higher latency)
    __m128d dotproduct = _mm_add_pd( (__m128d)temp, hi128 ); // 1 add
     // 3 shuffle + 2 add
*/

但对于 AMD，vextractf128 非常便宜，而 256b hadd 的成本是 128b hadd 的 2 倍，因此将每个 256b 产品分别缩小到 128b 然后与 128b hadd 组合是有意义的.

实际上，根据Agner Fog's tables，haddpd xmm,xmm 在 Ryzen 上是 4 uop。（而 256b ymm 版本是 8 uops）。因此，如果数据正确，实际上最好在 Ryzen 上手动使用 2x vshufpd + vaddpd。可能不是：他的 Piledriver 数据有 3 uop haddpd xmm,xmm，而且只有 4 uop 和内存操作数。对我来说，他们无法将hadd 实现为只有 3 个（或 6 个 ymm）微指令。

为了完成 4 个 dots 并将结果打包到一个 __m256d 中，所问的确切问题，我认为 @hirschhornsalz 的答案对于 Intel CPU 来说看起来非常好。我还没有仔细研究过它，但是与hadd 成对组合是很好的。 vperm2f128 在 Intel 上效率很高（但在 AMD 上相当糟糕：Ryzen 上 8 微指令，每 3c 吞吐量一个）。

【讨论】：