优化快速乘法但缓慢添加:FMA和doubledouble
Z b*_*son 9 floating-point x86 assembly mandelbrot fma
当我第一次使用Haswell处理器时,我尝试使用FMA来确定Mandelbrot集.主要算法是这样的:
intn = 0;
for(int32_t i=0; i<maxiter; i++) {
floatn x2 = square(x), y2 = square(y); //square(x) = x*x
floatn r2 = x2 + y2;
booln mask = r2<cut; //booln is in the float domain non integer domain
if(!horizontal_or(mask)) break; //_mm256_testz_pd(mask)
n -= mask
floatn t = x*y; mul2(t); //mul2(t): t*=2
x = x2 - y2 + cx;
y = t + cy;
}
这确定n像素是否在Mandelbrot集中.因此对于双浮点,它运行超过4个像素(floatn = __m256d,intn = __m256i).这需要4个SIMD浮点乘法和4个SIMD浮点加法.
然后我修改了这个就像这样使用FMA
intn n = 0;
for(int32_t i=0; i<maxiter; i++) {
floatn r2 = mul_add(x,x,y*y);
booln mask = r2<cut;
if(!horizontal_or(mask)) break;
add_mask(n,mask);
floatn t = x*y;
x = mul_sub(x,x, mul_sub(y,y,cx));
y = mul_add(2.0f,t,cy);
}
其中mul_add调用_mm256_fmad_pd和mul_sub调用_mm256_fmsub_pd.该方法使用4个FMA SIMD操作和两个SIMD乘法,这是没有FMA的两个算术运算.此外,FMA和乘法可以使用两个端口,只添加一个.
为了减少我的测试偏差,我放大了一个完全在Mandelbrot集中的区域,所以所有的值都是maxiter.在这种情况下,使用FMA的方法快约27%.这肯定是一个改进,但从SSE到AVX的性能翻了一番,所以我希望FMA可能有另外两个因素.
但后来我发现这个在问候FMA答案在那里说
融合乘加指令的重要方面是中间结果的(虚拟)无限精度.这有助于提高性能,但不是因为两个操作在单个指令中编码 - 它有助于提高性能,因为中间结果的几乎无限精度有时很重要,并且通过普通乘法和加法来恢复非常昂贵精确度正是程序员追求的目标.
然后给出一个double*double到double-double乘法的例子
high = a * b; /* double-precision approximation of the real product */
low = fma(a, b, -high); /* remainder of the real product */
由此,我得出结论,我正在非优化地实施FMA,因此我决定实施SIMD双倍.我基于GPU计算的扩展精度浮点数实现了双倍.这篇论文用于双浮动,所以我修改它为双倍.此外,不是在SIMD寄存器中打包一个双倍值,而是将4个双倍值打包到一个AVX高位寄存器和一个AVX低位寄存器中.
对于Mandelbrot集合我真正需要的是双倍乘法和加法.在那篇论文中,这些是df64_add和df64_mult功能.下图显示了我df64_mult的软件FMA(左)和硬件FMA(右)功能的组件.这清楚地表明硬件FMA是双倍乘法的重大改进.
那么硬件FMA如何在双倍Mandelbrot集合计算中执行？答案是它只比软件FMA快15%左右.这比我希望的要少得多.两双曼德尔布罗计算需要4两双添加和四个双双乘法(x*x,y*y,x*y,和2*(x*y)).然而,2*(x*y)乘法对于双倍是微不足道的,因此可以在成本中忽略这种乘法.因此,我认为使用硬件FMA的改进如此之小的原因是计算主要是慢速双倍加法(见下面的装配).
过去,乘法比加法慢(并且程序员使用了几个技巧来避免乘法)但是对于Haswell来说,它似乎是另一种方式.不仅是因为FMA,还因为乘法可以使用两个端口但只添加一个.
所以我的问题(最后)是:
- 与乘法相比,当加法缓慢时,如何优化？
- 是否有一种代数方法来改变我的算法以使用更多乘法和更少的加法？我知道有相反的方法,例如
(x+y)*(x+y) - (x*x+y*y) = 2*x*y,使用另外两个加法来减少一次乘法. - 有没有办法简单地使用df64_add函数(例如使用FMA)？
如果有人想知道双重方法比双重方法慢十倍.这并不是那么糟糕,我认为好像有一个硬件四精度类型,它可能至少是double的两倍慢,所以我的软件方法比我预期的硬件慢五倍(如果它存在的话).
df64_add 部件
vmovapd 8(%rsp), %ymm0
movq %rdi, %rax
vmovapd 72(%rsp), %ymm1
vmovapd 40(%rsp), %ymm3
vaddpd %ymm1, %ymm0, %ymm4
vmovapd 104(%rsp), %ymm5
vsubpd %ymm0, %ymm4, %ymm2
vsubpd %ymm2, %ymm1, %ymm1
vsubpd %ymm2, %ymm4, %ymm2
vsubpd %ymm2, %ymm0, %ymm0
vaddpd %ymm1, %ymm0, %ymm2
vaddpd %ymm5, %ymm3, %ymm1
vsubpd %ymm3, %ymm1, %ymm6
vsubpd %ymm6, %ymm5, %ymm5
vsubpd %ymm6, %ymm1, %ymm6
vaddpd %ymm1, %ymm2, %ymm1
vsubpd %ymm6, %ymm3, %ymm3
vaddpd %ymm1, %ymm4, %ymm2
vaddpd %ymm5, %ymm3, %ymm3
vsubpd %ymm4, %ymm2, %ymm4
vsubpd %ymm4, %ymm1, %ymm1
vaddpd %ymm3, %ymm1, %ymm0
vaddpd %ymm0, %ymm2, %ymm1
vsubpd %ymm2, %ymm1, %ymm2
vmovapd %ymm1, (%rdi)
vsubpd %ymm2, %ymm0, %ymm0
vmovapd %ymm0, 32(%rdi)
vzeroupper
ret
为了回答我的第三个问题,我找到了一个更快的双倍加法解决方案.我在图形硬件上实现float-float操作符的文章中找到了另一种定义.
Theorem 5 (Add22 theorem) Let be ah+al and bh+bl the float-float arguments of the following
algorithm:
Add22 (ah ,al ,bh ,bl)
1 r = ah ? bh
2 if | ah | ? | bh | then
3 s = ((( ah ? r ) ? bh ) ? b l ) ? a l
4 e l s e
5 s = ((( bh ? r ) ? ah ) ? a l ) ? b l
6 ( rh , r l ) = add12 ( r , s )
7 return (rh , r l)
这是我实现这个(伪代码)的方式:
static inline doubledoublen add22(doubledoublen const &a, doubledouble const &b) {
doublen aa,ab,ah,bh,al,bl;
booln mask;
aa = abs(a.hi); //_mm256_and_pd
ab = abs(b.hi);
mask = aa >= ab; //_mm256_cmple_pd
// z = select(cut,x,y) is a SIMD version of z = cut ? x : y;
ah = select(mask,a.hi,b.hi); //_mm256_blendv_pd
bh = select(mask,b.hi,a.hi);
al = select(mask,a.lo,b.lo);
bl = select(mask,b.lo,a.lo);
doublen r, s;
r = ah + bh;
s = (((ah - r) + bh) + bl ) + al;
return two_sum(r,s);
}
Add22的这个定义使用11个加法而不是20个,但它需要一些额外的代码来确定是否|ah| >= |bh|. 以下是有关如何实现SIMD minmag和maxmag函数的讨论.幸运的是,大多数附加代码不使用端口1.现在只有12条指令转到端口1而不是20.
以下是新的Add22的IACA吞吐量分析表
Throughput Analysis Report
--------------------------
Block Throughput: 12.05 Cycles Throughput Bottleneck: Port1
Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
---------------------------------------------------------------------------------------
| Cycles | 0.0 0.0 | 12.0 | 2.5 2.5 | 2.5 2.5 | 2.0 | 10.0 | 0.0 | 2.0 |
---------------------------------------------------------------------------------------
| Num Of | Ports pressure in cycles | |
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | |
---------------------------------------------------------------------------------
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm3, ymmword ptr [rip]
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm0, ymmword ptr [rdx]
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm4, ymmword ptr [rsi]
| 1 | | | | | | 1.0 | | | | vandpd ymm2, ymm4, ymm3
| 1 | | | | | | 1.0 | | | | vandpd ymm3, ymm0, ymm3
| 1 | | 1.0 | | | | | | | CP | vcmppd ymm2, ymm3, ymm2, 0x2
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm3, ymmword ptr [rsi+0x20]
| 2 | | | | | | 2.0 | | | | vblendvpd ymm1, ymm0, ymm4, ymm2
| 2 | | | | | | 2.0 | | | | vblendvpd ymm4, ymm4, ymm0, ymm2
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm0, ymmword ptr [rdx+0x20]
| 2 | | | | | | 2.0 | | | | vblendvpd ymm5, ymm0, ymm3, ymm2
| 2 | | | | | | 2.0 | | | | vblendvpd ymm0, ymm3, ymm0, ymm2
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm3, ymm1, ymm4
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm2, ymm1, ymm3
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm2, ymm4
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm1, ymm0
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm0, ymm1, ymm5
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm2, ymm3, ymm0
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm1, ymm2, ymm3
| 2^ | | | | | 1.0 | | | 1.0 | | vmovapd ymmword ptr [rdi], ymm2
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm0, ymm0, ymm1
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm1, ymm2, ymm1
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm3, ymm3, ymm1
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm0, ymm3, ymm0
| 2^ | | | | | 1.0 | | | 1.0 | | vmovapd ymmword ptr [rdi+0x20], ymm0
这是旧的吞吐量分析
Throughput Analysis Report
--------------------------
Block Throughput: 20.00 Cycles Throughput Bottleneck: Port1
Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
---------------------------------------------------------------------------------------
| Cycles | 0.0 0.0 | 20.0 | 2.0 2.0 | 2.0 2.0 | 2.0 | 0.0 | 0.0 | 2.0 |
---------------------------------------------------------------------------------------
| Num Of | Ports pressure in cycles | |
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | |
---------------------------------------------------------------------------------
| 1 | | | 1.0 1.0 | | | | | | | vmovapd ymm0, ymmword ptr [rsi]
| 1 | | | | 1.0 1.0 | | | | | | vmovapd ymm1, ymmword ptr [rdx]
| 1 | | | 1.0 1.0 | | | | | | | vmovapd ymm3, ymmword ptr [rsi+0x20]
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm4, ymm0, ymm1
| 1 | | | | 1.0 1.0 | | | | | | vmovapd ymm5, ymmword ptr [rdx+0x20]
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm2, ymm4, ymm0
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm1, ymm1, ymm2
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm2, ymm4, ymm2
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm0, ymm0, ymm2
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm2, ymm0, ymm1
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm3, ymm5
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm6, ymm1, ymm3
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm5, ymm5, ymm6
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm6, ymm1, ymm6
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm2, ymm1
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm3, ymm3, ymm6
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm2, ymm4, ymm1
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm3, ymm3, ymm5
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm4, ymm2, ymm4
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm1, ymm1, ymm4
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm0, ymm1, ymm3
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm2, ymm0
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm2, ymm1, ymm2
| 2^ | | | | | 1.0 | | | 1.0 | | vmovapd ymmword ptr [rdi], ymm1
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm0, ymm0, ymm2
| 2^ | | | | | 1.0 | | | 1.0 | | vmovapd ymmword ptr [rdi+0x20], ymm0
如果除了FMA之外还有三个操作数单舍入模式指令,则更好的解决方案.在我看来应该有单一的舍入模式指令
a + b + c
a * b + c //FMA - this is the only one in x86 so far
a * b * c
| 归档时间： |
|
| 查看次数： |
944 次 |
| 最近记录： |