5

我有一些 C++ 代码，随着时间的推移，它已经成为一个有用的 FFT 库，并且使用 SSE 和 AVX 指令使其运行得相当快。当然，这一切都只是基于 radix-2 算法，但它仍然成立。我最近想抓挠的是让蝴蝶计算与 FMA 指令一起工作。基本的 radix-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

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

4

1 回答 1

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运算，和书上的运算量是一样的。

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

于 2020-03-28T19:09:30.017 回答