英特尔 AVX：双精度浮点变量的256位点积版本

Question

英特尔高级矢量扩展（AVX）在256位版本（YMM寄存器）中不提供双精度浮点变量的点积。关于“为什么”问题已经在另一个论坛（here）和Stack Overflow（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（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之间的点积结果。也许这可以为建议提供一些提示？

- gleeen.gould

---

谢谢您的建议，但我应该在我的应用程序中保持双精度。 - gleeen.gould

此外，转换+浮点数点积所需的时间比双倍点积更长。 - Gunther Piez

Answer 1

我会使用4倍精度浮点数乘法，然后使用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 );

- Gunther Piez

---

谢谢！你能解释一下最后一行吗？我不太确定是否理解得很好。难道不是 _mm256_add_pd 吗？ - gleeen.gould

我仍然很难理解它。我会使用以下代替最后两行：__m256d swapped1 = _mm256_permute2f128_pd( temp01, temp23, 0x30 ); __m256d swapped2 = _mm256_permute2f128_pd( temp01, temp23, 0x21 ); __m256d dotproduct = _mm256_add_pd( swapped1, swapped2 );代码未经测试（我的机器上没有 AVX :-)） - gleeen.gould

@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个周期。(希望我得到了正确的常数) - Norbert P.

@gleeen.gould：AVX2已经发布了吗？我以为它会在2013年的Haswell上推出。我是在谈论Sandy Bridge上当前一代AVX：请参见Agner Fog的指令表，第129页。 - Norbert P.

注意：即使使用 AVX2 和 vpermpd，使用 vextractf128 和 addpd 的解决方案的组合延迟仍然比连续应用 vpermpd 和 vhaddpd 生成水平和要低。 - Pixelchemist

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)在dotproduct的low double中，而dot(z,w)则在high double中。

- Norbert P.

---

我正在尝试使用Visual Studio 2022，但遇到了一个错误：error C2440: 'type cast': cannot convert from '__m256d' to '__m128d'。修复它的最佳方法是什么？ - Jepessen

Answer 3

对于单个点积，它只是一个垂直乘法和水平求和（参见在x86上执行水平浮点向量求和的最快方法）。hadd需要2次洗牌和1次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个洗牌uop，在AMD Jaguar / Bulldozer-family / Ryzen上获得了更大的胜利，因为它立即收缩到128b，而不是做一堆256b的事情。 AMD将256b操作拆分为两个128b uop。

---

在并行进行2或4个点积时，使用hadd可能是值得的，其中您正在将其与2个不同的输入向量一起使用。 Norbert的两对向量点乘在想要打包结果时看起来最优。即使使用AVX2 vpermpd作为跨越车道的洗牌，我也看不到任何更好的方法。 当然，如果您真的希望有一个更大的dot（8个或更多个double），请使用垂直add（使用多个累加器来隐藏vaddps延迟），并在最后进行水平求和。 如果可用，也可以使用fma。

---

haddpd内部以两种不同的方式将xy和zw混合在一起，并将其馈送到垂直的addpd，这正是我们手动完成的操作。如果我们保持xy和zw分开，则需要为每个点积进行2次洗牌+ 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
*/

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

实际上，根据Agner Fog的表格，haddpd xmm，xmm在Ryzen上有4个uop。 （256b ymm版本为8个uop）。因此，在Ryzen上手动使用2x vshufpd + vaddpd比使用hadd更好，如果数据正确的话。但这可能不是真的：他的Piledriver数据具有3 uop的haddpd xmm，xmm，并且使用内存操作数只有4个uops。我觉得他们无法将hadd实现为仅3（或ymm为6）个uop。

---

要将4个dot结果打包到一个__m256d中，正是所要求的精确问题，我认为@ hirschhornsalz的答案看起来非常适用于英特尔CPU。 我还没有仔细研究过，但使用hadd成对组合是很好的。 vperm2f128在英特尔上效率很高（但在AMD上非常糟糕：每3c吞吐量为1）。

- Peter Cordes