快速增殖但加法慢的优化: FMA和双倍

Question

当我第一次得到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运算和2次SIMD乘法运算，比不使用FMA的算法运算少2次。另外，FMA和乘法可以使用两个端口，而加法只能使用一个端口。

为了减少测试的偏倚性，我放大到一个完全位于Mandelbrot集中的区域，所以所有的值都是maxiter。在这种情况下，使用FMA的方法大约快了27%，当然是一个改进，但是从SSE到AVX，我的性能翻了一番，所以我希望使用FMA的另一个因素是两个。

但是我找到了关于FMA的这答案，上面写着

融合乘法加指令的重要方面是中间结果的(实质上)无限精度。这有助于提高性能，但这并不是因为两个操作是在一条指令中编码的--这有助于性能，因为中间结果的几乎无限精度有时很重要，而且在程序员真正追求这种精度的情况下，用普通乘法和加法恢复这种精度非常昂贵。

然后给出了一个双倍到双倍乘法的例子。

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%左右。这比我所期望的要少得多。的双双Mandelbrot计算需要4个双倍加法和4个双倍乘法(x*x、y*y、x*y和2*(x*y))。然而，乘法对于双倍来说是微不足道的。使这种乘法可以忽略在成本上。因此，我认为使用硬件FMA的改进是如此小的原因是，计算是以缓慢的双倍加法为主(见下面的组装)。

过去，乘法比加法慢(编程人员用了几个技巧来避免乘法)，但对于Haswell，情况似乎正好相反。这不仅是由于FMA，也是因为乘法可以使用两个端口，但只有一个。

所以我的问题(最后)是：

1. 当加法比乘法慢时，如何优化？
2. 有没有一种代数的方法来改变我的算法来使用更多的乘法和更少的加法？我知道有一些方法可以做相反的事情，例如(x+y)*(x+y) - (x*x+y*y) = 2*x*y，它使用两个加法来减少一个乘法。
3. 是否有一种方法可以简化df64_add函数(例如使用FMA)？

如果有人想知道双倍的方法比双倍慢十倍。这并不是很糟糕，我认为如果有硬件四精度类型，它可能至少是双倍速度的两倍，所以我的软件方法比我期望的硬件慢五倍。

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

Answer 1

为了回答我的第三个问题，我找到了一个双倍加法的更快的解决方案。我在论文浮点操作符在图形硬件上的实现中找到了另一个定义。

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个。

下面是新的IACA的吞吐量分析表单Add22

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

Answer 2

为了加快算法的速度，我使用了基于2 fma、1 mul和2 add的简化版本。我以这种方式处理8次迭代。然后计算转义半径，必要时回滚最后8次迭代。

下面的关键循环X= X^2 +C是用x86内部函数编写的，编译器很好地展开了它，并且在展开之后，您会发现两个FMA操作彼此并不严重依赖。

//  IACA_START;
for (j = 0; j < 8; j++) {
    Xrm = _mm256_mul_ps(Xre, Xim);
    Xtt = _mm256_fmsub_ps(Xim, Xim, Cre);
    Xrm = _mm256_add_ps(Xrm, Xrm);
    Xim = _mm256_add_ps(Cim, Xrm);
    Xre = _mm256_fmsub_ps(Xre, Xre, Xtt);
}       // for
//  IACA_END;

然后我计算出逃逸半径(x，x，<阈值)，它需要一个其他的fma和另一个乘法，每8次迭代一次。

cmp = _mm256_mul_ps(Xre, Xre);
cmp = _mm256_fmadd_ps(Xim, Xim, cmp);
cmp = _mm256_cmp_ps(cmp, vec_threshold, _CMP_LE_OS);
if (_mm256_testc_si256((__m256i) cmp, vec_one)) {
    i += 8;
    continue;
}

您提到“加法是缓慢的”，这不是完全正确的，但您是对的，乘法吞吐量随着时间的推移在最近的架构上越来越高。

乘法延迟和依赖关系是关键。FMA的吞吐量为1个周期，延迟为5个周期。独立FMA指令的执行可能重叠。

基于乘法结果的加法得到了完全的延迟命中。

因此，您必须通过执行“代码拼接”并在同一个循环中计算两个点来打破这些直接的依赖关系，并且只需在与IACA检查将要发生的事情之前将代码插入。下面的代码有2组变量( X0=X0^2+C0和X1=X1^2+C1的后缀为0和1)，并开始填充FMA漏洞

for (j = 0; j < 8; j++) {
    Xrm0 = _mm256_mul_ps(Xre0, Xim0);
    Xrm1 = _mm256_mul_ps(Xre1, Xim1);
    Xtt0 = _mm256_fmsub_ps(Xim0, Xim0, Cre);
    Xtt1 = _mm256_fmsub_ps(Xim1, Xim1, Cre);
    Xrm0 = _mm256_add_ps(Xrm0, Xrm0);
    Xrm1 = _mm256_add_ps(Xrm1, Xrm1);
    Xim0 = _mm256_add_ps(Cim0, Xrm0);
    Xim1 = _mm256_add_ps(Cim1, Xrm1);
    Xre0 = _mm256_fmsub_ps(Xre0, Xre0, Xtt0);
    Xre1 = _mm256_fmsub_ps(Xre1, Xre1, Xtt1);
}       // for

总之，

• 您可以将关键循环中的指令数()减半。
• 您可以添加更多独立的指令，，并获得高吞吐量和低延迟的乘法和融合乘法和添加的优势。

Answer 3

您提到以下代码：

vsubpd  %ymm0, %ymm4, %ymm2
vsubpd  %ymm2, %ymm1, %ymm1  <-- immediate dependency ymm2
vsubpd  %ymm2, %ymm4, %ymm2
vsubpd  %ymm2, %ymm0, %ymm0  <-- immediate dependency ymm2
vaddpd  %ymm1, %ymm0, %ymm2  <-- immediate dependency ymm0
vaddpd  %ymm5, %ymm3, %ymm1
vsubpd  %ymm3, %ymm1, %ymm6  <-- immediate dependency ymm1
vsubpd  %ymm6, %ymm5, %ymm5  <-- immediate dependency ymm6
vsubpd  %ymm6, %ymm1, %ymm6  <-- dependency ymm1, ymm6
vaddpd  %ymm1, %ymm2, %ymm1
vsubpd  %ymm6, %ymm3, %ymm3  <-- dependency ymm6
vaddpd  %ymm1, %ymm4, %ymm2 
vaddpd  %ymm5, %ymm3, %ymm3  <-- dependency ymm3
vsubpd  %ymm4, %ymm2, %ymm4 
vsubpd  %ymm4, %ymm1, %ymm1  <-- immediate dependency ymm4
vaddpd  %ymm3, %ymm1, %ymm0  <-- immediate dependency ymm1, ymm3
vaddpd  %ymm0, %ymm2, %ymm1  <-- immediate dependency ymm0
vsubpd  %ymm2, %ymm1, %ymm2  <-- immediate dependency ymm1

如果仔细检查，这些操作大多是依赖的操作，并且没有满足关于延迟/吞吐量效率的基本规则。大多数指令取决于上一条指令的结果，或之前的2条指令。该序列包含30周期的关键路径(关于"3周期延迟“/”1周期吞吐量“的大约9或10条指令)。

您的IACA在关键路径中报告"CP“=>指令，估计的成本是20个周期的吞吐量。您应该得到延迟报告，因为如果您对执行速度感兴趣，那么它才是重要的。

为了消除这条关键路径的成本，如果编译器不能这样做，您必须再插入大约20条类似的指令(例如，因为您的双双代码在一个单独的库中编译，没有-flto优化，函数输入和退出时到处都有vzeroupper，向量程序只能很好地处理内联代码)。

一种可能是并行运行2次计算(参见前一篇文章中关于代码拼接以改进流水线)。

如果我假设您的双双代码看起来像这样的“标准”实现

// (r,e) = x + y
#define two_sum(x, y, r, e) 
    do { double t; r = x + y; t = r - x; e = (x - (r - t)) + (y - t); } while (0)
#define two_difference(x, y, r, e) \
    do { double t; r = x - y; t = r - x; e = (x - (r - t)) - (y + t); } while (0)
.....

然后，您必须考虑以下代码，其中指令以相当精细的方式交织在一起。

// (r1, e1) = x1 + y1, (r2, e2) x2 + y2
#define two_sum(x1, y1, x2, y2, r1, e1, r2, e2) 
    do { double t1, t2 \
    r1 = x1 + y1; r2 = x2 + y2; \
    t1 = r1 - x1; t2 = r2 - x2; \
    e1 = (x1 - (r1 - t1)) + (y1 - t1); e2 = (x2 - (r2 - t2)) + (y2 - t2);  \
} while (0)
....

然后创建如下代码(延迟报告中的相同关键路径，以及大约35条指令)。有关运行时、无序执行的详细信息，请在不延迟的情况下执行。

vsubsd  %xmm2, %xmm0, %xmm8
vsubsd  %xmm3, %xmm1, %xmm1
vaddsd  %xmm4, %xmm4, %xmm4
vaddsd  %xmm5, %xmm5, %xmm5
vsubsd  %xmm0, %xmm8, %xmm9
vsubsd  %xmm9, %xmm8, %xmm10
vaddsd  %xmm2, %xmm9, %xmm2
vsubsd  %xmm10, %xmm0, %xmm0
vsubsd  %xmm2, %xmm0, %xmm11
vaddsd  %xmm14, %xmm4, %xmm2
vaddsd  %xmm11, %xmm1, %xmm12
vsubsd  %xmm4, %xmm2, %xmm0
vaddsd  %xmm12, %xmm8, %xmm13
vsubsd  %xmm0, %xmm2, %xmm11
vsubsd  %xmm0, %xmm14, %xmm1
vaddsd  %xmm6, %xmm13, %xmm3
vsubsd  %xmm8, %xmm13, %xmm8
vsubsd  %xmm11, %xmm4, %xmm4
vsubsd  %xmm13, %xmm3, %xmm15
vsubsd  %xmm8, %xmm12, %xmm12
vaddsd  %xmm1, %xmm4, %xmm14
vsubsd  %xmm15, %xmm3, %xmm9
vsubsd  %xmm15, %xmm6, %xmm6
vaddsd  %xmm7, %xmm12, %xmm7
vsubsd  %xmm9, %xmm13, %xmm10
vaddsd  16(%rsp), %xmm5, %xmm9
vaddsd  %xmm6, %xmm10, %xmm15
vaddsd  %xmm14, %xmm9, %xmm10
vaddsd  %xmm15, %xmm7, %xmm13
vaddsd  %xmm10, %xmm2, %xmm15
vaddsd  %xmm13, %xmm3, %xmm6
vsubsd  %xmm2, %xmm15, %xmm2
vsubsd  %xmm3, %xmm6, %xmm3
vsubsd  %xmm2, %xmm10, %xmm11
vsubsd  %xmm3, %xmm13, %xmm0

摘要：

• 内联您的双双源代码:编译器和向量器不能优化跨函数调用，因为ABI约束，和跨内存访问担心混叠。
• 编写代码，以平衡吞吐量和延迟，并最大限度地利用CPU端口(并且每个周期的指令也最大化)，只要编译器不将太多寄存器泄漏到内存中。

您可以使用perf实用程序(包linux generic)跟踪优化影响，以获得每个周期执行的指令数量和指令数量。