c++ - 将FMA指令用于FFT算法

我有一些C++代码，随着时间的流逝，它已成为一些有用的FFT库，并且使用SSE和AVX指令使其运行得相当快。当然，它们全都基于radix-2算法，但仍然有效。我最新的尝试是使蝶形计算与FMA指令一起使用。基本的基数2蝶形由4个乘法和6个加法或减法组成。一种简单的方法是用2条FMA指令替换2个加法和减法以及2个乘法，从而产生数学上相同的蝶形图，但是显然有更好的方法:

https://books.google.com/books?id=2HG0DwAAQBAJ&pg=PA56&lpg=PA56&dq=radix+2+fft+fma&source=bl&ots=R5XDWyYBVv&sig=ACfU3U0S2n1hcgiP63LTKMxI5Oc85eEZaQ&hl=en&sa=X&ved=2ahUKEwiz_I3PsrToAhVoHzQIHYmVDGIQ6AEwDXoECAoQAQ#v=onepage&q=radix%202%20fft%20fma&f=false

ci1 = ci1 / cr1
u0 = zinr(0)
v0 = zini(0)
r = zinr(1)
s = sini(1)
u1 = r - s * ci1
v1 = r * ci1 + s
zoutr(0) = u0 + u1 * cr1
zouti(0) = v0 + v1 * cr1
zoutr(1) = u0 - u1 * cr1
zouti(1) = v0 - v1 * cr1

如果将旋转因子的虚部除以实部，则作者用6个FMA替换了所有10个添加项，子项和杂项。文本的部分内容为“请注意，cr1!= 0”。简而言之，这实际上是我的问题。对于所有旋转因子，数学似乎都像广告中所述的那样工作，除了当实际旋转为零时，在这种情况下，我们最终除以零。在这里效率至关重要的地方，当cr1 == 0时将代码分支到另一只蝶形不是一个好选择，尤其是当我们使用SIMD一次处理多个旋转和蝶形时，其中cr1 ==的一个元素0。我的直觉告诉我应该是这种情况，即当cr1 == 0时，cr1和ci1应该完全是其他一些值，并且FMA代码仍将得出正确的答案，但是我似乎无法弄清楚。如果我能弄清楚的话，修改FMA蝴蝶的预先计算的旋转因子将是相对简单的事情，当然，我们也可以避免在蝴蝶开始时进行除法运算。

最佳答案

这本书似乎暗示cr1 != 0总是正确的。但不幸的是，情况并非总是如此(旋转角度为PI/2时)。

我认为您无法通过调整旋转因子来解决此问题。我看到的唯一选择是使用一些很小的数字而不是零。它可以工作，但是很丑陋，并且在某些情况下可能会导致错误。

可能的解决方案:

将循环拆分为两个，并处理这种中心情况(发生零除)，特别是

而不是除以cr1，而除以ci1，然后相应地修改论坛。这种情况下的除数仍然为零，但是会在循环的第一次迭代中发生。因此，您必须专门处理第一次迭代(而不是中心)(因此只需要一个循环)。

使用不同的FMA公式:

注意，

zoutr(1) = u0 - u1
         = u0 - u1 - (u0 + u1) + (u0 + u1)
         = u0 - u1 - zoutr(0) + u0 + u1
         = 2*u0 - zoutr(0)

因此，可以在1 FMA中完成此操作。

并且如果您将u1替换为zoutr(0)的表达式:

zoutr(0) = u0 + u1
         = u0 + r*cr1 - s*ci1

这可以通过2个FMA来完成。

可以按照与zouti相同的方式来计算zoutr。因此，您需要使用6个FMA操作，这与本书所进行的操作相同。

(请注意，这并不意味着此变体将自动运行得更快，因为它具有不同的数据依赖链)

关于c++ - 将FMA指令用于FFT算法，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/60862508/