快速矢量化rsqrt和SSE/AVX的倒数取决于精度

Question

假设有必要计算打包浮点数据的倒数或倒数平方根.两者都可以轻松完成:

__m128 recip_float4_ieee(__m128 x) { return _mm_div_ps(_mm_set1_ps(1.0f), x); }
__m128 rsqrt_float4_ieee(__m128 x) { return _mm_div_ps(_mm_set1_ps(1.0f), _mm_sqrt_ps(x)); }

这种方法效果很好但很慢:根据指南,它们在Sandy Bridge上进行了14次和28次循环(吞吐量).对应的AVX版本在Haswell上几乎占用相同的时间.

另一方面,可以使用以下版本:

__m128 recip_float4_half(__m128 x) { return _mm_rcp_ps(x); }
__m128 rsqrt_float4_half(__m128 x) { return _mm_rsqrt_ps(x); }

它们只需要一到两个时间周期(吞吐量),从而大大提升了性能.但是,它们非常接近:它们产生的结果相对误差小于1.5*2 ^ -12.鉴于单精度浮点数的机器epsilon是2 ^？24,我们可以说这种近似具有大约一半的精度.

似乎可以添加Newton-Raphson迭代以产生具有单精度的结果(可能不像IEEE标准所要求的那样精确),参见GCC,ICC,LLVM上的讨论.理论上,相同的方法可用于双精度值,产生半精度或单精度或双精度.

我有兴趣为float和double数据类型以及所有(half,single,double)精度实现此方法的实现.处理特殊情况(除以零,sqrt(-1),inf/nan等)不是必需的.此外,我不清楚这些例程中的哪一个比普通的IEEE编译解决方案更快,哪个更慢.

以下是对答案的一些小限制,请:

1. 在代码示例中使用内在函数.程序集依赖于编译器,因此不太有用.
2. 对函数使用类似的命名约定.
3. 实现例程,将单个SSE/AVX寄存器包含密集打包的float/double值作为输入.如果有相当大的性能提升,你也可以发布几个寄存器作为输入的例程(两个reg可能是可行的).
4. 如果两个SSE/AVX版本绝对等于将_mm更改为_mm256,则不要发布它们,反之亦然.

欢迎任何性能评估,测量和讨论.

摘要

以下是具有一次NR迭代的单精度浮点数的版本:

__m128 recip_float4_single(__m128 x) {
  __m128 res = _mm_rcp_ps(x);
  __m128 muls = _mm_mul_ps(x, _mm_mul_ps(res, res));
  return res =  _mm_sub_ps(_mm_add_ps(res, res), muls);
}
__m128 rsqrt_float4_single(__m128 x) {
  __m128 three = _mm_set1_ps(3.0f), half = _mm_set1_ps(0.5f);
  __m128 res = _mm_rsqrt_ps(x); 
  __m128 muls = _mm_mul_ps(_mm_mul_ps(x, res), res); 
  return res = _mm_mul_ps(_mm_mul_ps(half, res), _mm_sub_ps(three, muls)); 
}

Peter Cordes给出的答案解释了如何创建其他版本,并包含全面的理论性能分析.

您可以在此处找到所有已实施的基准测试解决方案: recip_rsqrt_benchmark.

Ivy Bridge获得的吞吐量结果如下所示.只有单注册SSE实现已经过基准测试.花费的时间以每次呼叫的周期给出.第一个数字用于半精度(无NR),第二个用于单精度(1个NR迭代),第三个用于2个NR迭代.

1. RECIP上浮动取1,4次循环相对于7个周期.
2. rsqrt上浮动取1,6个循环对14个周期.
3. RECIP上双取3,6,9个周期相对于14个循环.
4. rsqrt on double需要3个,8个,13个周期而不是28个周期.

警告:我必须创造性地将原始结果舍入...

Answer 1

在实践中有很多算法的例子.例如:

带有SSE2的Newton Raphson - 有人可以解释我这3行有一个答案,解释了英特尔的一个例子所使用的迭代.

对于井上分析,让我们说Haswell(它可以在两个执行端口上进行FP,与以前的设计不同),我将在这里重现代码(每行一个操作).有关指令吞吐量和延迟的表格,以及有关如何理解更多背景的文档,请参见http://agner.org/optimize/.

// execute (aka dispatch) on cycle 1, results ready on cycle 6
nr = _mm_rsqrt_ps( x );

// both of these execute on cycle 6, results ready on cycle 11
xnr = _mm_mul_ps( x, nr );         // dep on nr
half_nr = _mm_mul_ps( half, nr );  // dep on nr

// can execute on cycle 11, result ready on cycle 16
muls = _mm_mul_ps( xnr , nr );     // dep on xnr

// can execute on cycle 16, result ready on cycle 19
three_minus_muls = _mm_sub_ps( three, muls );  // dep on muls

// can execute on cycle 19, result ready on cycle 24
result = _mm_mul_ps( half_nr, three_minus_muls ); // dep on three_minus_muls
// result is an approximation of 1/sqrt(x), with ~22 to 23 bits of precision in the mantissa.

如果它不是依赖链的一部分,那么这里有很多其他计算重叠的空间.但是,如果你的代码的每一次迭代的数据依赖于从以前的数据,你可能会用的11个周期的延迟更好sqrtps.或者即使每个循环迭代足够长,乱序执行也无法通过重叠独立迭代来隐藏它们.

为了得到sqrt(x)而不是1/sqrt(x)乘以x(并且如果x可以为零则_mm_andn_ps修复,例如通过使用CMPPS的结果使用掩蔽()来抵消0.0).最佳的方法是更换half_nr用half_xnr = _mm_mul_ps( half, xnr );.这不会延长dep链,因为half_xnr可以在第11周期开始,但直到结束(第19周期)才需要.与可用的FMA相同:总延迟没有增加.

如果精度11位是足够(无牛顿迭代),英特尔公司的优化手册(第11.12.3)建议使用rcpps(rsqrt(X)),这是比由原始x相乘,至少与AVX更糟.它也许可以节省128位SSE一个MOVDQA指令,但256B rcpps慢于256B MUL或FMA.(并且它允许您在最后一步使用FMA而不是mul免费添加sqrt结果.)

此循环的SSE版本(不考虑任何移动指令)为6 uops.这意味着它应该在Haswell上具有每3个周期一个的吞吐量(假设两个执行端口可以处理FP mul,而rsqrt在FP add/sub的相对端口上).在SnB/IvB(可能还有Nehalem)上,它应该具有每5个周期一个的吞吐量,因为mulps和rsqrtps竞争端口0.subps在port1上,这不是瓶颈.

对于Haswell,我们可以使用FMA将减法与mul结合起来.但是,分频器/ sqrt单元不是256b宽,因此与其他所有器件不同,ymm regs上的divps/sqrtps/rsqrtps/rcpps需要额外的uops和额外的周期.

// vrsqrtps ymm has higher latency
// execute on cycle 1, results ready on cycle 8
nr = _mm256_rsqrt_ps( x );

// both of can execute on cycle 8, results ready on cycle 13
xnr = _mm256_mul_ps( x, nr );         // dep on nr
half_nr = _mm256_mul_ps( half, nr );  // dep on nr

// can execute on cycle 13, result ready on cycle 18
three_minus_muls = _mm256_fnmadd_ps( xnr, nr, three );  // -(xnr*nr) + 3

// can execute on cycle 18, result ready on cycle 23
result = _mm256_mul_ps( half_nr, three_minus_muls ); // dep on three_minus_muls

我们用FMA节省了3个周期.这可以通过使用2-cyle-slow 256b rsqrt来抵消,净延迟减少1个周期(相当于两倍宽).SnB/IvB AVX将是128b的24c延迟,256b的26c延迟.

256b FMA版本总共使用7次.(VRSQRTPS3个uop,2个用于p0,1个用于p1/5.)256b mulps和fma都是单uop指令,并且都可以在端口0或端口1上运行.(p0仅在Haswell之前).因此,吞吐量应该是:每3c一个,如果OOO引擎将uop调度到最佳执行端口.(即来自rsqrt的shuffle uop总是进入p5,永远不会进入p1,它将占用mul/fma带宽.)至于与其他计算重叠(不仅仅是自身的独立执行),它再次非常轻量级.只有7个具有23个循环dep链的uops为其他东西留下了大量空间,而那些uops位于重新排序缓冲区中.

如果这是巨型dep链中的一个步骤而没有其他任何事情发生(甚至不是独立的下一次迭代),则256b vsqrtps是19个周期延迟,吞吐量为每14个周期一个.(Haswell的).如果你仍然需要倒数,那么256b vdivps也有18-21c延迟,每14c吞吐量一个.因此对于普通sqrt,指令具有较低的延迟.对于recip sqrt,近似的迭代是较少的延迟.(如果它与自身重叠,那么FAR会增加吞吐量.如果与其他没有分割单元的东西重叠,sqrtps则不是问题.)

sqrtps可以是吞吐量获胜与rsqrt+牛顿迭代,如果它是循环体的一部分,其中有足够的其他非分割和非sqrt工作,则除法单元未饱和.

如果你需要sqrt(x),尤其如此1/sqrt(x).例如,在带有AVX2的Haswell上,复制+ arcsinh循环在一个适合L3缓存的浮点数组中,实现为fastlog(v + sqrt(v*v + 1)),与真实VSQRTPS或VRSQRTPS + Newton-Raphson迭代以大约相同的吞吐量运行.(即使对log()进行非常快速的近似,所以总循环体大约是9 FMA/add/mul/convert操作,2布尔值加上VSQRTPS ymm.使用just的速度很快fastlog(v2_plus_1 * rsqrt(v2_plus_1) + v2_plus_1),所以它没有瓶颈关于内存带宽,但它可能是延迟的瓶颈(因此无序执行无法利用独立迭代的所有并行性).

其他精度

对于半精度,没有指令在半浮点数上进行计算.您打算在使用转换说明加载/存储时动态转换.

对于双精度,没有rsqrtpd.据推测,我们的想法是,如果你不需要完全精确,你应该首先使用浮动.所以你可以转换为float和back,然后执行完全相同的算法,但pd不是ps内在函数.或者,您可以将数据保持浮动一段时间.例如,将两个ymm的双打寄存器转换成一个ymm的单个寄存器.

不幸的是,没有一条指令需要两个ymm的双精度寄存器并输出一个单独的ymm.你必须去ymm-> xmm两次,然后_mm256_insertf128_ps一个xmm到另一个的高128.但是,您可以将256b ymm向量提供给相同的算法.

如果你要在之后转换回加倍,那么分别对两个单独的128b寄存器进行rsqrt + Newton-Raphson迭代是有意义的.额外的2个微指令插入/提取物,以及额外的2个微指令的256B rsqrt,开始增加,更不用提的3个周期的延迟vinsertf128/ vextractf128.计算将重叠,两个结果将分开几个周期.(或者相隔1个周期,如果你有一个特殊版本的代码,可以同时对2个输入进行交错操作).

请记住,单精度的指数范围小于double,因此转换可以溢出到+ Inf或下溢到0.0.如果你不确定,一定只是使用正常_mm_sqrt_pd.

有了AVX512F,就有了_mm512_rsqrt14_pd( __m512d a).使用AVX512ER(KNL但不是SKX或Cannonlake),有_mm512_rsqrt28_pd. _ps版本当然也存在.在更多情况下,14位尾数精度可能足以在没有牛顿迭代的情况下使用.

牛顿迭代仍然不会像常规sqrt 那样为您提供正确舍入的单精度浮点数.