当我第一次使用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 . 在这种情况下 the method using FMA is about 27% faster. 这肯定是一个改进，但从SSE到AVX的性能翻了一倍，所以我希望FMA可能有另外两个因素 .
但后来我发现了有关FMA的答案.2904328_答案
融合乘加指令的重要方面是中间结果的（虚拟）无限精度 . 这有助于提高性能，但不是因为两个操作在一条指令中编码 - 它有助于提高性能，因为中间结果的几乎无限精度有时很重要，并且通过普通乘法和加法来恢复非常昂贵精确度正是程序员追求的目标 .
然后给出一个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双倍 . 我根据论文Extended-Precision Floating-Point Numbers for GPU Computation实施了双倍 . 这篇论文用于双浮动，所以我修改它为双倍 . 此外，不是在SIMD寄存器中打包一个双倍值，而是将4个双倍值打包到一个AVX高位寄存器和一个AVX低位寄存器中 .
对于Mandelbrot集合我真正需要的是双倍乘法和加法 . 在那篇论文中，这些是 df64_add 和 df64_mult 函数 . 下图显示software FMA（左）和硬件FMA（右）的 df64_mult 函数的程序集 . 这清楚地表明硬件FMA是双倍乘法的重大改进 .
那么硬件FMA如何在双倍Mandelbrot集合计算中执行？ The answer is that's only about 15% faster than with software FMA. That's much less than I hoped for. 双倍Mandelbrot计算需要4次双重加法和4次双重乘法（ x*x ， y*y ， x*y 和 2*(x*y) ） . 但是，2*(x*y) multiplication is trivial for double-double所以这个乘法可以在成本中忽略 . 因此，我认为使用硬件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
3 回答
为了回答我的第三个问题，我找到了一个更快的双倍加法解决方案 . 我在论文Implementation of float-float operators on graphics hardware中找到了另一种定义 .
这是我实现这个（伪代码）的方式：
Add22的这个定义使用11个加法而不是20个，但它需要一些额外的代码来确定
|ah| >= |bh|. Here is a discussion on how to implement SIMD minmag and maxmag functions . 幸运的是，大多数附加代码不使用端口1.现在只有12条指令端口1而不是20 .以下是新Add22的吞吐量分析表单IACA
这是旧的吞吐量分析
如果除了FMA之外还有三个操作数单舍入模式指令，则更好的解决方案 . 在我看来应该有单一的舍入模式指令
为了加快算法速度，我使用了基于2 fma，1 mul和2 add的简化版本 . 我以这种方式处理8次迭代 . 然后计算转义半径并在必要时回滚最后8次迭代 .
使用x86内部函数编写的以下关键循环X = X ^ 2 C很好地由编译器展开，并且在展开之后您将发现2 FMA操作彼此之间并非严重依赖 .
然后我计算逃逸半径（| X | <threshold），它只花费其他fma和另一个乘法，每8次迭代 .
你提到“加速很慢”，这不完全正确，但你说得对，乘法吞吐量随着时间的推移在最近的架构上变得越来越高 .
乘法延迟和依赖关系是关键 . FMA的吞吐量为1个周期，延迟为5个周期 . 独立FMA指令的执行可以重叠 .
基于乘法结果的加法得到完整的延迟命中 .
因此，您必须通过执行“代码拼接”并在同一循环中计算2个点来打破这些直接依赖关系，并在检查IACA之前只交错代码将会发生什么 . 下面的代码有两组变量（后缀为0和1代表X0 = X0 ^ 2 C0，X1 = X1 ^ 2 C1）并开始填充FMA孔
总结一下，
你可以在你的关键循环中 halve the number of instructions
你可以 add more independent instructions 并获得高吞吐量与乘法的低延迟以及融合乘法和加法的优势 .
你提到以下代码：
如果仔细检查，这些主要是依赖操作，并且不符合关于延迟/吞吐效率的基本规则 . 大多数指令取决于前一个指令的结果，或之前的2个指令 . 该序列包含30个循环的关键路径（关于“3个循环延迟”/“1个循环吞吐量”约9或10个指令） .
您的IACA在关键路径中报告“CP”=>指令，评估的成本为20个周期的吞吐量 . 您应该获得延迟报告，因为如果您对执行速度感兴趣，那么它就是重要的 .
为了消除这个关键路径的成本，如果编译器不能这样做，你必须交错大约20个类似的指令（例如，因为你的双重代码是在没有-flto优化的情况下编译的单独库中，并且在函数入口和出口处到处都是vzeroupper ，vectorizer只适用于内联代码） .
一种可能性是并行运行2次计算（参见前一篇文章中的代码拼接以改进流水线操作）
如果我认为你的双重代码看起来像这个“标准”实现
然后你必须考虑以下代码，其中指令以非常精细的方式交错 .
然后，这将创建类似下面的代码（关于延迟报告中的相同关键路径，以及大约35条指令） . 有关运行时的详细信息，乱序执行应该在没有停顿的情况下飞越 .
摘要：
inline 您的双重双重源代码：由于ABI约束，编译器和向量化程序无法跨函数调用进行优化，并且由于担心别名而无法跨内存访问进行优化 .
stitch 代码可以 balancer 吞吐量和延迟并最大化CPU端口使用（并且还可以最大化每个周期的指令），只要编译器不会将过多的寄存器溢出到内存中 .
您可以使用perf实用程序（包linux-tools-generic和linux-cloud-tools-generic）跟踪优化影响，以获取执行的指令数和每个周期的指令数 .