Использование инструкций FMA для алгоритма FFT

У меня есть немного кода на C++, который со временем стал довольно полезной библиотекой БПФ, и он был сделан для достаточно быстрой работы с использованием инструкций SSE и AVX. Конечно, все это основано только на алгоритме счисления по основанию 2, но он все еще работает. Последнее, что мне не терпится поцарапать, — заставить вычисления бабочки работать с инструкциями FMA. Базовая бабочка по основанию-2 состоит из 4 умножений и 6 сложений или вычитаний. Простой подход будет включать замену 2 операций сложения и вычитания и 2 операций умножения двумя инструкциями FMA, в результате чего получится математически идентичная бабочка, но, по-видимому, есть и лучшие способы сделать это:

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

Автор заменяет все 10 добавлений, сабвуферов и мультов на 6 FMA при условии, что мнимая часть фактора поворота делится на действительную часть. Часть текста гласит: «Обратите внимание, что cr1 != 0». Это, по сути, моя проблема в двух словах. Кажется, что математика работает точно так же, как рекламируется, для всех коэффициентов оборота, за исключением случаев, когда реальный оборот равен нулю, и в этом случае мы заканчиваем делением на ноль. Там, где эффективность здесь абсолютно критична, код ветвления, когда cr1 == 0, к другой бабочке не является хорошим вариантом, особенно когда мы используем SIMD для одновременной обработки нескольких твиддлов и бабочек, где, возможно, только один элемент cr1 == 0. Моя интуиция говорит мне, что это должно иметь место, так это то, что когда cr1 == 0, cr1 и ci1 должны быть совершенно другими значениями, и код FMA все равно приведет к правильному ответу, но я не могу понять это . Если бы я мог понять это, было бы относительно просто изменить предварительно вычисленные коэффициенты поворота для бабочек FMA, и мы также, конечно, могли бы избежать операции деления в начале бабочки.


person Kumputer    schedule 26.03.2020    source источник
comment
Я не могу ничего прочитать по указанной вами ссылке. Частично связано: если вас интересует эффективность, вы тестировали с помощью radix-4?   -  person Damien    schedule 01.04.2020


Ответы (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, что равно количеству операций, которое есть в книге.

(Обратите внимание, это не означает, что этот вариант будет работать быстрее автоматически, так как он имеет другую цепочку зависимостей данных)

person geza    schedule 28.03.2020