FMA в Raspberrypi

Разбираясь с одним интересным багом в коде DSP фильтра, я с удивлением обнаружил целую область неизведанного в области чисел с плавающий точкой. Есть замечательная статья о том, что каждый программист должен знать об операциях с плавающей точкой. В ней очень подробно описывается проблема и её решение в стандарте IEEE 754. Однако, там нет явного упоминания такой вещи, как FMA (Fused multiply-add). Мне показалось интересным разобраться, как эта технология работает.

FMA

FMA - это оптимизация на уровне процессора следующей операции:

result += A * B

Выбор такой операции не случаен. Это ни что иное, как операция свёртки, которая очень часто используется в анализе сигналов (DSP), различных алгоритмах матанализа и пр. В общем виде её можно записать так:

result[i] = A[i]*B[i] + A[i+1]*B[i+1] + ...

Собственно оптимизация заключается в двух вещах:

  • умножение и сложение делается одной ассемблерной инструкцией вместо двух
  • округление происходит только после сложения.

На последнем я хотел бы немного остановится. Дело в том, что регистры процессора имеют конечную длину. Записать туда бесконечное число не получится. Поэтому результат умножения должен быть округлён, сложен и округлён ещё раз:

round(round(A[i] * B[i]) + result[i])

Однако, из-за того, что FMA - это одна ассемблерная инструкция, то можно сделать небольшую оптимизацию и произвести умножение на гораздо большем регистре, сложить и результат округлить:

round(A[i]*B[i] + result[i])

Это позволит существенно повысить точность, так как округление происходит только один раз. Понятное дело, что полностью от ошибки округления избавиться нельзя, но вот значительно сократить ошибку вполне возможно.

BCM2837

Для проверки FMA я совершенно неслучайно выбрал Raspberrypi 3. На нём установлен SoC Broadcom BCM2837. Это чип, в котором одновременно есть и CPU и GPU. А поскольку они оба могут производить операции с плавающей точкой, то, конечно же, интересно было бы узнать как работает FMA на каждом из них.

Хотя на самом деле нет. Я просто пытаюсь найти багу в своём GPU коде.

rpi-fma

Я написал небольшой проект - rpi-fma, в котором явным образом с помощью ассемблерных инструкций считаю простое выражение:

1.0000001F * 1.0000001F - 1.0000002F

Я взял его не случайно, а нашёл в одной статье от Nvidia, в которой описывалась работа их GPU процессора. Она достаточно обширна, но написана была весной 2011. В частности, там утверждается, что:

At the time this paper is written (Spring 2011) there are no commercially available x86 CPUs which offer hardware FMA

Если использовать одни и те же числа, то можно сравнить их с результатами инженеров Nvidia.

Я реализовал несколько разных вычислений одного и того же выражения:

  • С использованием двух инструкций: VMUL и VADD
  • С использованием FMA инструкции: VFMA
  • С использованием MLA инструкции: VMLA
  • OpenCL код для запуска на GPU, который в итоге компилировался в две инструкции: FMUL и FADD

MLA - это одна инструкция для умножения и сложения с результатом. Её основное отличие от FMA заключается в том, что не делается “Fuse” - оптимизация округления.

Результаты в таблице ниже:

Тест Результат
VMUL + VADD 0
VFMA 1.42108547e-14
VMLA 0
FMUL + FADD 0

Выводы

Выводов неожиданно несколько:

  • FMA существует и на процессорах.
  • FMA работает и даёт результаты отличные от наивной реализации через 2 инструкции.
  • В GPU Raspberrypi нет поддержки FMA. Я дополнительно залез в спецификацию процессора и не обнаружил там ничего похожего.