将 FMA 指令用于 FFT 算法

Question

我有一些 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=假

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

如果旋转因子的虚部除以实部，作者将所有 10 个加法、减法和乘法替换为 6 个 FMA。部分文字为“注意 cr1 != 0”。简而言之，这基本上是我的问题。数学似乎对所有旋转因子都有效，除非真正的旋转因子为零，在这种情况下，我们最终除以零。在这里效率绝对至关重要，当 cr1 == 0 时将代码分支到不同的蝴蝶不是一个好的选择，尤其是当我们使用 SIMD 一次处理多个旋转和蝴蝶时，其中可能只有 cr1 == 的一个元素0. 我的直觉告诉我应该是这样，当 cr1 == 0 时，cr1 和 ci1 应该完全是其他一些值，而 FMA 代码仍然会产生正确的答案，但我似乎无法弄清楚这一点。如果我能弄清楚，修改 FMA 蝴蝶的预先计算的旋转因子将是一件相对简单的事情，我们当然也可以避免蝴蝶开始时的除法运算。

Answer 1

这本书似乎表明这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运算，这和书本上的运算量是一样的。

（请注意，这并不意味着该变体会自动运行得更快，因为它具有不同的数据依赖链）