英特尔AVX:256位版本的点积,用于双精度浮点变量

Question

英特尔高级矢量扩展指令集(AVX)在256位版本(YMM寄存器)中不提供双精度浮点变量的点积."为什么？" 问题已在另一个论坛(此处)和Stack Overflow(此处)进行了简要处理.但我面临的问题是如何以有效的方式用其他AVX指令替换这条缺失的指令？

对于单精度浮点变量,存在256位版本的点积(此处参考):

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

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

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

更具体地说,我想从__m128(四个浮点数)转换为__m256d(4个双精度数)的代码使用以下指令:

   __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,val4和val0.

也许这可以给出建议的提示？

Answer 1

我会使用4*双乘法,然后a hadd(不幸的是,在上半部分和下半部分只添加2*2个浮点数),提取上半部分(一个shuffle应该同样工作,可能更快)并将其添加到下半部分.

结果是低64位dotproduct.

__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 );

Answer 2

我会延长drhirsch的答案,同时执行两个点产品,节省一些工作:

__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.

Answer 3

对于单个点积，它只是垂直乘法和水平和（请参阅在 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 的单向量答案在英特尔上 1 shuffle uop 更好，并且在 AMD Jaguar / Bulldozer-family / Ryzen 上获得更大的胜利，因为它立即缩小到 128b 而不是做一个一堆 256b 的东西。AMD 将 256b ops 分成两个 128b uops。

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

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

haddpd内部 shufflexy和zw两种不同的方式放在一起，并将其提供给一个垂直addpd，无论如何我们都会手工做。如果我们保持xy并zw分开，我们需要 2 次 shuffle + 2 次添加以获得点积（在单独的寄存器中）。因此，hadd作为第一步，将它们混洗在一起，我们节省了混洗的总数，仅节省了添加和总 uop 计数。

/*  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
*/

但是对于vextractf128非常便宜的AMD 来说，256b 的hadd成本hadd是 128b 的 2 倍，将每个 256b 产品分别缩小到 128b，然后再与 128b hadd 结合是有意义的。

实际上，根据Agner Fog 的表格，haddpd xmm,xmm在 Ryzen 上是 4 uops。（而 256b ymm 版本是 8 uops）。因此，如果数据正确，在 Ryzen 上手动使用 2x vshufpd+实际上更好vaddpd。可能不是：他的 Piledriver 数据有 3 uop haddpd xmm,xmm，并且只有 4 uop 带有内存操作数。对我来说，他们不能hadd仅实现3 个（或 ymm 的 6 个）uop是没有意义的。

为了dot将结果打包为一个4秒，所__m256d问的确切问题，我认为@hirschhornsalz 的答案对于英特尔 CPU 来说看起来非常好。我没有仔细研究它，但成对组合hadd很好。 vperm2f128在 Intel 上效率很高（但在 AMD 上很糟糕：Ryzen 上 8 uops，每 3c 吞吐量一个）。