Генерация фракталов с помощью SSE/SSE2 - Метод MOV-ADD

ОГЛАВЛЕНИЕ

Метод MOV-ADD

Профилирование показало, что функция PaintMandelbrotFPU занимает около 95%-99% полного времени выполнения. Рисование занимает только 5% с GDI или 1% при использовании DirectDraw. Оптимизировать данную функцию гораздо важнее, чем любую иную часть программы.

Один из лучших способов оптимизации - использование SSE (потоковые SIMD-расширения). Эти команды позволяют вычислять значения 4 пикселей одновременно. Следовательно, хотя бы в теории, скорость возрастет в 4 раза (реальные цифры - около 3.2-3.3). Дополнительную информацию о SSE можно найти в Руководстве разработчика ПО под архитектуру Intel, том 1, глава 10, и в Руководстве программиста под архитектуру AMD64, том 1, глава 4. Дальнейшие объяснения предполагают, что вы знаете основы программирования на ассемблере под SSE и x86.

Расположить команды SSE для этой функции можно разными способами, и некоторые из них быстрее остальных. Ниже приведена очень быстрая последовательность для процессора Pentium IV, называемая методом MOV-ADD:

; xmm0 = zx  xmm1 = zy  xmm2 = tmp  xmm3 = tmp  
;            xmm4 = px  xmm5 = py  xmm6 = результат  xmm7 = 4.0
; eax = tmp    ecx = счётчик цикла
   MOVAPS xmm0,xmm4
   XORPS  xmm6,xmm6
   MOVAPS xmm1,xmm5
   mov ecx, ITER
ILOOP:
   ; xmm0 = zx             xmm1 = zy
   MOVAPS xmm2, xmm0
   MULPS  xmm0, xmm0
   MOVAPS xmm3, xmm1
   ADDPS  xmm1, xmm1
   ; xmm0 = zx^2           xmm1 = 2 * zy     xmm2 = zx           xmm3 = zy
   MULPS  xmm1, xmm2
   MOVAPS xmm2, xmm0
   MULPS  xmm3, xmm3
   ; xmm0 = zx^2           xmm1 = 2*zy*zx    xmm2 = zx^2         xmm3 = zy^2
   ADDPS  xmm1, xmm5
   SUBPS  xmm0, xmm3
   ADDPS  xmm2, xmm3
   ; xmm0 = zx^2 - zy^2    xmm1=2*zy*zx+py   xmm2 = zx^2 + zy^2  xmm3 = zy^2
   CMPLEPS xmm2, xmm7
   ADDPS   xmm0, xmm4
   MOVMSKPS eax, xmm2
   test eax, eax
   jz EXIT
   ANDPS   xmm2, xmm7      ; xmm6 += (xmm2 < 4.0) ? 4.0 : 0.0;
   ADDPS   xmm6, xmm2
   sub ecx, 1
   jnz ILOOP
EXIT:

Здесь показан только внутренний цикл. К началу цикла xmm0 и xmm1 содержат px и py, начальные значения для zx и zy. Регистры xmm4 и xmm5 хранят px и py (две постоянные цикла). Регистр xmm7 содержит константу 4.0 во всех двойных словах.

Программа начинает работу с сохранения xmm0 и xmm1 в буферные регистры xmm2 и xmm3. Затем zx возводится в квадрат, и zy удваивается путем сложения xmm1 с самим собой. Остальные шаги расчета показаны в комментариях выше.

Регистр xmm6 хранит число полных итераций. Заметьте, что число в xmm6 должно считаться независимо для каждого из 4 пикселей. Команда CMPLEPS помещает битовую маску в регистр xmm2, и последовательность MovMskPS-test-jz проверяет, вышли ли все 4 пикселя за пределы круга радиусом 2.0. Если это происходит, то цикл прерывается. В иных случаях программа должна увеличить число полных итераций для пикселей, все еще пребывающих внутри круга. Затем производится операция "логическое и"над маской и 4.0, чтобы получить 4.0 в двойных словах, соответствующих таким пикселям, и 0.0 в других двойных словах. Это значение прибавляется к xmm6. Короче говоря, когда цикл заканчивается, xmm6 будет содержать 4.0 * число полных итераций для каждого из 4 пикселей. Значение позже будет использоваться как сдвиг в палитре.

Теперь вычисляется время выполнения для Pentium IV и находятся узкие места в этом коде. Руководство от Anger Fog описывает, как провести такой анализ; если вы впервые знакомитесь с оптимизацией для Pentium IV, прочитайте руководство. Также смотрите официальные документы от Intel: Справочник по оптимизации для архитектуры Intel и Микроархитектура процессора Pentium 4.

Во-первых, обращение к памяти не является узким местом, потому что вычисление занимает гораздо больше времени, чем запись в память. Для обычного разрешения 1024x768 программа записывает 1024*768*4 байта = 3 Мб. Типичная производительность современной памяти DDR составляет 1-2 Гб/с. Грубая оценка показывает, что обращение к памяти займет 1.5-3 миллисекунды, тогда как реальное время выполнения этого кода - 70-80 мс. Анализ кода дает такой же результат: происходит всего 2 обращения к памяти для 1 пикселя (смотрите код в следующем разделе), но внутренний цикл выполняется ITER раз (обычно около 64 раз), и он включает в себя команды SSE с долгими задержками.

Поэтому разные приемы оптимизации производительности кэша бесполезны для этого кода. То же самое верно для расшифровки команд: циклы столь малы, что они полностью вмещаются в трассировочный кэш и расшифровываются лишь один раз. Реальное узкое место – долгие задержки и ограниченная производительность команд SSE.

Для простоты предположим, что конвейер вначале пустой, и весь код находится в трассировочном кэше. Имена регистров xmm1, xmm2, и т.д. являются сокращениями для x1, x2 чтобы они вошли на рисунок.

Первые три команды (MOVAPS x2, x0, MULPS x0, x0, и MOVAPS x3, x1) отправляются из трассировочного кэша в очередь операций. Кластер выполнения FP/SSE имеет два порта: один - для операций add, mul, div, названный блоком fp в руководстве Agner Fog, а второй порт - для пересылок между регистрами и для сохранений (блок выполнения mov). Эти два блока могут начать выполнять команды MOVAPS и MULPS одновременно.

Во втором тактовом цикле следующие три команды уходят в очередь. Из них MOVAPS x3, x1 и ADDPS x1, x1 могут выполняться параллельно. Взаимная производительность для команды MOVAPS равна 1, следовательно, в этом цикле может запуститься вторая команда.
Задержка команды MOVAPS составляет 6 тактовых циклов, и MULPS x1, x2 зависит от результатов MOVAPS x2, x0 (зависимость показана стрелками на схеме). Следовательно, MULPS может запуститься в шестом цикле (к этому времени очередь должна быть полна операций, ожидающих выполнения).

Следующая команда MOVAPS x2, x0 зависит от MULPS x0, x0. Команда MULPS также имеет дополнительную задержку 1, а значит, ей надо ждать 1 лишний тактовый цикл, если следующая операция уходит в другой блок выполнения. MULPS и MOVAPS родом из разных блоков выполнения fp и mov, и между ними пришлось вставить задержку в один цикл.
Теперь команда MULPS x3, x3. Подблок умножения (умножитель FP в понятиях Intel) может производить одно умножение каждые два цикла, поэтому второй команде MULPS придется ждать один цикл. Она запустится в восьмом цикле.

Следующие пять команд уходят в подблок суммирования FP (сумматор FP). Все они имеют задержку в 4 тактовых цикла. ADDPS x1, x5 может запуститься на 12-м цикле, когда будут готовы результаты умножения MUL x1, x2. То же самое касается SUBPS x0, x3; она запустится в 14-м цикле. Команда ADDPS x2, x3 может запуститься раньше, чем на 16-м цикле, но производительность сумматора FP ограничивает число одновременно выполняющихся команд до 2. Последовательность SUBPS x0, x3; ADDPS x2, x3; ADDPS x0, x4; CMPPS x2, x7 идеально соответствует этому требованию сумматора FP. Две соседних команды взяты из разных цепей зависимости, следовательно, они полностью независимы и могут выполняться параллельно.
Задача – завершить вычисление x0 = x0 ^ 2 - x1 ^ 2 + x4, x1 = 2 * x0 * x1 + x5 как можно быстрее. На графике команды, вычисляющие новые значения для x0 и x1, отмечены разными оттенками серого, а неважные команды белые. Надо включить в график серые команды раньше остальных, потому что следующая итерация цикла не может начаться без новых значений для x0 и x1. Команды ANDPS и MOVMSKPS из предыдущей итерации могут выполняться параллельно со следующей итерацией, при условии, что переход будет предсказан правильно. Стало быть, их можно включить в график позже серых команд.

Проведенные измерения времени подтвердили эти соображения. Сейчас есть долгая задержка между MOVAPS x2, x0 и ADDPS x2, x3 (смотрите схему). Но если поместить ADDPS x2, x3 раньше SUBPS x0, x3, то полное время выполнения увеличится. Значит, задержка не важна, при условии что x0 и x1 вычисляются рано.

Быстродействие кода зависит от задержек команд и производительности, потому что есть три независимых цепи, способных выполняться одновременно. Критический путь отмечен темно-серым цветом. MOVAPS добавляет лишнюю задержку (может быть, если бы архитектура x86 имела трехадресные команды, как RISC(вычисления с сокращённым набором команд), выполнялся бы быстрее, как знать). Подблоки не используются полностью. На первой стадии цикла (циклы 0-12) большая нагрузка ложится на блок MOV и на умножитель FP; позже, в циклах 14-22, перегружается сумматор FP, а умножитель FP вообще не используется. Можно взять команды из двух итераций и перемешать их, чтобы получить более равномерную нагрузку, но для этого надо больше регистров SSE. В AMD64 (или Intel EM64T) имеется 8 дополнительных регистров, следовательно, код можно дополнительно оптимизировать в 64-битном режиме.
После рисования схемы команды были упорядочены по их времени запуска. Например, команды MOVAPS x2, x0 и MULPS x0, x0 запускаются раньше всех, поэтому поставлены в коде первыми. Затем запускаются MOVAPS x3, x1 и ADDPS x1, x1, поэтому они являются следующими командами в исходном коде. Упорядочивание команд сокращает число операций, ждущих в очереди. Теперь процессору не надо переупорядочивать команды при выполнении: они уже упорядочены согласно их задержкам, производительности и зависимостям.

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

Были случайно выбраны два увеличенных изображения наборов, затем было вычислено распределение числа итераций для них и для одного большого изображения:

 

Голубые полосы обозначают большое изображение набора Мандельброта; красные и желтые полосы показывают два увеличенных изображения. ITER в программе равно 64, поэтому диаграмма показывает частоты точек от 0 до 5 итераций, от 6 до 11 итераций, и так далее до 60 итераций и более. Минимальное и максимальное числа итераций составляют крупнейшие части изображения. Для голубой полосы минимальное значение равно 0-5 итераций (59% всех пикселей); для красной полосы - 12-17 (43%); и для желтой полосы - 30-35 (8%) и 36-41 (15%). Минимальное число итераций соответствует огромной фоновой области, закрашенной одним цветом, а максимальное число итераций означает черные точки, принадлежащие набору. Таким образом, 26% пикселей принадлежат набору Мандельброта на большом изображении, а 21% и 62% точек принадлежат набору на двух увеличенных изображениях.

Из этого анализа ясно, что если цикл прерывается, то велика вероятность, что он будет прерван рано, на минимальном числе итераций. Но если цикл не прерывается рано, обычно он вообще не прерывается. Эта статистика интересна, но ей не нашлось применения в коде. Может быть, читатель сможет сделать это.

Ниже сделано несколько выводов о методе MOV-ADD:

•    Общая стратегия – как можно раньше вычислить значения, требуемые в следующей итерации (в данном случае, x0 и x1). Эти вычисления составляют критический путь программы, и полное время выполнения зависит от их времени завершения.
•    Планирование одновременного выполнения независимых команд улучшает параллелизм на уровне команд (ILP). Например, первые четыре команды выполняют сложение, умножение и две пересылки параллельно.
•    Команды упорядочены по их времени запуска, что уменьшает длину очереди операций.