Есть ли шанс ускорить повторяющийся код с SIMD?

Question

Рассмотрим следующий код, где a - это массив параметров float, а s - изначально неинициализированный массив результатов float:

s[n - 1] = mu * a[n - 1];
for (int j = n - 2; j >= 0; j--)
    s[j] = mu * (a[j] + s[j + 1]);
return s;

Есть ли шанс улучшитьпроизводительность такого повторяющегося кода с SIMD (AVX2)?

РЕДАКТИРОВАТЬ: я выяснил позже, что эта формула / алгоритм называется «сумма со скидкой», но не смог найти его параллельную версию в Интернете.

Answer 1

Связано: Можно ли использовать SIMD для последовательной зависимости в вычислениях, например, для экспоненциального фильтра скользящего среднего? - Если есть формула закрытой формы для n шагов вперед, вы можете использовать ее длябоковые последовательные зависимости.Но я не думаю, что это так.

Это похоже на префиксную сумму типа последовательной зависимости, поверх вертикального добавления с a[j].Есть способы ускорить это, получив ускорение примерно как
O( SIMD_width / log(SIMD_width) ).

• Сумма префикса SIMD на процессоре Intel
• параллельная префиксная (накопительная) сумма с SSE

Answer 2

Иногда помогает писать выражения в виде матрично-векторных произведений.Предполагая, что вы уже знаете sₖ₊₈, вы можете вычислить sₖ до sₖ₊₇ с aₖ до aₖ₊₇, используя

 [ µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷ µ⁸]   [aₖ₊₀     ]
 [ 0  µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷]   [aₖ₊₁     ]
 [ 0  0  µ  µ² µ³ µ⁴ µ⁵ µ⁶]   [aₖ₊₂     ]
 [ 0  0  0  µ  µ² µ³ µ⁴ µ⁵]   [aₖ₊₃     ]
 [ 0  0  0  0  µ  µ² µ³ µ⁴] * [aₖ₊₄     ]
 [ 0  0  0  0  0  µ  µ² µ³]   [aₖ₊₅     ]
 [ 0  0  0  0  0  0  µ  µ²]   [aₖ₊₆     ]
 [ 0  0  0  0  0  0  0  µ ]   [aₖ₊₇+sₖ₊₈]

Поскольку sₖ₊₈, вероятно, будет иметь некоторую задержку при вычислении,имеет смысл вывести его из продукта.Это можно рассчитать с помощью одной широковещательной рассылки и одного сложного-кратного сложения:

 [ µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷ µ⁸]   [aₖ₊₀]   [ µ⁸]
 [ 0  µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷]   [aₖ₊₁]   [ µ⁷]
 [ 0  0  µ  µ² µ³ µ⁴ µ⁵ µ⁶]   [aₖ₊₂]   [ µ⁶]
 [ 0  0  0  µ  µ² µ³ µ⁴ µ⁵]   [aₖ₊₃]   [ µ⁵]
 [ 0  0  0  0  µ  µ² µ³ µ⁴] * [aₖ₊₄] + [ µ⁴] * sₖ₊₈
 [ 0  0  0  0  0  µ  µ² µ³]   [aₖ₊₅]   [ µ³]
 [ 0  0  0  0  0  0  µ  µ²]   [aₖ₊₆]   [ µ²]
 [ 0  0  0  0  0  0  0  µ ]   [aₖ₊₇]   [ µ ] 

И первую матрицу можно разложить на три матрицы, которые можно рассчитать, используя одну случайную последовательность и одну FMA каждую:

 [ 1  0  0  0  µ⁴ 0  0  0 ]   [ 1  0  µ² 0  0  0  0  0 ]   [ µ  µ² 0  0  0  0  0  0 ]   [aₖ₊₀]   [ µ⁸]
 [ 0  1  0  0  µ³ 0  0  0 ]   [ 0  1  µ  0  0  0  0  0 ]   [ 0  µ  0  0  0  0  0  0 ]   [aₖ₊₁]   [ µ⁷]
 [ 0  0  1  0  µ² 0  0  0 ]   [ 0  0  1  0  0  0  0  0 ]   [ 0  0  µ  µ² 0  0  0  0 ]   [aₖ₊₂]   [ µ⁶]
 [ 0  0  0  1  µ  0  0  0 ]   [ 0  0  0  1  0  0  0  0 ]   [ 0  0  0  µ  0  0  0  0 ]   [aₖ₊₃]   [ µ⁵]
 [ 0  0  0  0  1  0  0  0 ] * [ 0  0  0  0  1  0  µ² 0 ] * [ 0  0  0  0  µ  µ² 0  0 ] * [aₖ₊₄] + [ µ⁴] * sₖ₊₈
 [ 0  0  0  0  0  1  0  0 ]   [ 0  0  0  0  0  1  µ  0 ]   [ 0  0  0  0  0  µ  0  0 ]   [aₖ₊₅]   [ µ³]
 [ 0  0  0  0  0  0  1  0 ]   [ 0  0  0  0  0  0  1  0 ]   [ 0  0  0  0  0  0  µ  µ²]   [aₖ₊₆]   [ µ²]
 [ 0  0  0  0  0  0  0  1 ]   [ 0  0  0  0  0  0  0  1 ]   [ 0  0  0  0  0  0  0  µ ]   [aₖ₊₇]   [ µ ]

Самый правый матричный вектор-продукт на самом деле на одно умножение больше.

В целом, для 8 элементов вам нужно 4FMA, одно умножение и 4 тасования / широковещания ($$$$ означает что угодно (конечно) может быть здесь - в качестве альтернативы, если они гарантированно равны 0, векторы µ могут быть частично разделены. Все векторы обозначаются как наименее значимые-сначала, все умножения поэлементно):

bₖₖ₊₇  = [aₖ₊₀, aₖ₊₁, aₖ₊₂, aₖ₊₃, aₖ₊₄, aₖ₊₅, aₖ₊₆, aₖ₊₇] * [µ  µ  µ  µ  µ  µ  µ  µ ]     vmulps
bₖₖ₊₇ += [aₖ₊₁, $$$$, aₖ₊₃, $$$$, aₖ₊₅, $$$$, aₖ₊₆, $$$$] * [µ² 0  µ² 0  µ² 0  µ² 0 ]     vshufps (or vpsrlq) + vfmadd 

cₖₖ₊₇  =  bₖₖ₊₇ 
cₖₖ₊₇ += [bₖ₊₂, bₖ₊₂, $$$$, $$$$, bₖ₊₆, bₖ₊₆, $$$$, $$$$] * [µ² µ  0  0  µ² µ  0  0 ]     vshufps + vfmadd

dₖₖ₊₇  =  cₖₖ₊₇ 
dₖₖ₊₇ += [cₖ₊₄, cₖ₊₄, cₖ₊₄, cₖ₊₄, $$$$, $$$$, $$$$, $$$$] * [µ⁴ µ³ µ² µ  0  0  0  0 ]     vpermps + vfmadd

sₖₖ₊₇ = dₖₖ₊₇ 
       + [sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈] * [µ⁸ µ⁷ µ⁶ µ⁵ µ⁴ µ³ µ² µ ]     vbroadcastss + vfmadd

Если я проанализировал это правильно, вычисление кратного dₖ может чередоваться, что сведет на нет задержки.И единственный горячий путь будет заключительным vbroadcastss + vfmadd для вычисления sₖₖ₊₇ из dₖₖ₊₇ и sₖ₊₈. может иметь смысл вычислять блоки 16 sₖ и fmadd µⁱ * sₖ₊₁₆.Кроме того, FMA для расчета dₖ используют только половину элементов.При некотором изощренном стремлении можно вычислить два блока с одинаковым количеством FMA (я полагаю, это не стоит усилий, но вы можете попробовать это).

Для сравнения: для чисто скалярной реализации требуется 8сложения и 8 умножений для 8 элементов, и каждая операция зависит от предыдущего результата.

---

Примечание: вы можете сохранить одно умножение, если вместо формулы рассчитали:

sₖ = aₖ₊₁ + µ*sₖ₊₁

Кроме того, в скалярной версии вы бы использовали Fused-Multiple-Adds вместо того, чтобы сначала добавлять и умножать потом.Результат будет отличаться только в µ.

раза

Answer 3

Если вы имеете в виду векторизацию, то нет, не напрямую.Поскольку для вычисления следующей итерации используется значение предыдущей итерации, она не является тривиальной векторизацией.

Кроме того, использование s без выравнивания может вызвать проблемы.Цикл с конца может быть тоже.

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