SIMD и векторные расширения: параллелизм внутри ядра

Возьмите цикл, который складывает два массива по миллиону чисел одинарной точности. Работы в нём — миллион сложений. Но процессор, исполняя наивную скалярную версию, выполнит миллион загрузок пары значений, миллион сложений, миллион записей, миллион увеличений индекса, миллион сравнений и миллион переходов. Полезная арифметика — примерно шестая часть команд, и за каждую из этих команд, включая полезную, ядро платит полную цену: чтение из кэша команд, декодирование, переименование регистров, место в буфере переупорядочивания, планирование, фиксация.

Внеочередное ядро умеет исполнять эти команды параллельно — это тема главы 07. Но оно не умеет их не выполнять. Каждая команда всё равно проедет через фронтенд и всю бухгалтерию, потратив на это энергию и место в очередях.

SIMD — другой ответ на ту же задачу. Не «исполним много команд одновременно», а «сделаем одну команду, которая обрабатывает много данных». Базовая механика команды, регистра и АЛУ здесь предполагается известной (Как работает процессор). Формально это single instruction, multiple data, но полезнее держать в голове экономическую формулировку: вся фиксированная стоимость команды делится на число элементов, которые она обработала.

Три вида параллелизма, которые нельзя путать

Трек про железо всё время говорит о параллелизме, но это три разных ресурса: их добывают разные механизмы, и упираются они в разные пределы.

Кто находит параллелизм. ILP извлекается железом динамически: окно внеочередного исполнения смотрит на несколько десятков или сотен команд и сам решает, что можно запустить раньше. DLP объявляется статически: в момент компиляции кто-то решил, что эти восемь операций независимы, и записал их одной командой. Железу больше нечего доказывать — гарантия приехала в кодировке команды.

Во что упирается. ILP упирается в длину цепочки зависимостей и в размер окна; удвоить ширину внеочередного ядра стоит примерно квадратичного роста сложности планировщика и сети обходов. DLP упирается в то, лежат ли данные так, чтобы их можно было взять пачкой, — а дальше в пропускную способность памяти. TLP упирается в синхронизацию.

Как сочетается. Эти три вида параллелизма перемножаются, а не заменяют друг друга. Восемь полос × четыре независимые цепочки накопления × шестнадцать ядер — это три разных множителя, и потерять любой из них означает потерять свою долю. Обзорно эта таксономия разбиралась в Параллелизме и конкурентности; здесь нас интересует только средний множитель.

Почему вектор дешевле по энергии, а не только по времени

Самый частый способ объяснить SIMD — «восемь сложений вместо одного». Это верно, но объясняет только половину выигрыша, причём менее важную.

Разложим стоимость одной скалярной команды на высокопроизводительном внеочередном ядре 2020-х годов:

Полезная работа:      одна операция АЛУ над 32 битами
Фиксированные траты:  чтение строки кэша команд и её декодирование
                      переименование архитектурного регистра в физический
                      выделение записи в буфере переупорядочивания
                      постановка в очередь планировщика и пробуждение по готовности
                      чтение операндов из физического регистрового файла
                      запись результата и фиксация в порядке программы

Ключевой факт, который стоит запомнить: на универсальном ядре энергия этой бухгалтерии сопоставима с энергией самой арифметики, а часто заметно больше её. Это не оценка «на глаз»: измерения энергетического бюджета универсальных процессоров показывают, что на простых целочисленных операциях доля собственно вычисления в общем расходе — единицы процентов, всё остальное съедают выборка, декодирование и работа с регистрами (Hameed et al., «Understanding Sources of Inefficiency in General-Purpose Chips», ISCA 2010). Порядок величины, разумеется, зависит от класса ядра: у простого микроконтроллера пропорции другие, там фронтенд гораздо дешевле.

Теперь то же самое, но команда работает над 512-битным регистром из 16 значений float32. Фиксированные траты не изменились — запись в буфере переупорядочивания одна, переименование одно, планирование одно. Изменился знаменатель: они делятся на 16. А для int8 в том же регистре — на 64. К этому добавляется вторая, чисто временная экономия: накладные расходы самого цикла. Увеличение индекса, сравнение с границей и условный переход выполняются один раз на вектор, а не один раз на элемент. Для короткого тела цикла это заметная доля всех команд. Отсюда практический вывод, который переворачивает интуицию: векторизация часто выгодна даже там, где она не даёт ускорения в разы. Мобильное или встраиваемое ядро, работающее в жёстком бюджете мощности, выигрывает от SIMD прежде всего джоулями на обработанный элемент, и лишь во вторую очередь секундами (Мощность и пределы).

Устройство векторного блока

Векторный регистр — это просто широкий регистр, который интерпретируется как массив элементов. Сами биты не знают, сколько их в элементе: разбиение задаёт команда.

Ширина регистра — сколько бит помещается. 128 у SSE и NEON, 256 у AVX, 512 у AVX-512; у SVE и RISC-V V это свойство конкретной реализации, а не контракта.

Тип элемента — как эти биты нарезаны. Одна и та же 256-битная строка данных — это 4 значения float64, 8 float32, 16 int16 или 32 int8. Команда vaddps складывает как float32, vpaddb — как int8 с разрядом переноса, обрезанным на границе каждого байта. Никакого переноса между полосами: полосы изолированы, и это главное упрощение, которое делает векторный блок дешёвым.

Число полос исполнения — сколько АЛУ физически стоит рядом. Вот здесь начинается тонкость, о которой обычно молчат: ширина в контракте не обязана совпадать с шириной в кремнии. Реализация имеет полное право взять 256-битную команду и прогнать её за два такта через 128-битный тракт. Код при этом работает без изменений, но пропускная способность вдвое ниже ожидаемой. Так делали при первом появлении новых расширений в младших ядрах: сначала совместимость, потом полная ширина. Это ровно та линия «архитектура против микроархитектуры», которую мы проводили в главе 03.

Что значит «одна операция за такт при ширине 256 бит»: исполнительный порт способен принимать новую 256-битную операцию каждый такт. Латентность при этом остаётся многотактовой — для сложения или умножения с плавающей точкой на высокопроизводительных ядрах последнего десятилетия это порядок 3–5 тактов. Различие латентности и пропускной способности здесь работает ровно так же, как в главе 06: чтобы загрузить порт полностью, нужно несколько независимых цепочек накопления. Отдельно стоит FMA — слитное умножение со сложением. Одна команда считает a*b + c с единственным округлением в конце. Это удваивает пиковую арифметику на такт и одновременно слегка меняет результат по сравнению с раздельными операциями — деталь, к которой мы вернёмся в разделе про редукции.

Две модели: фиксированная ширина и переменная длина

Фиксированная ширина: SSE, AVX, AVX-512, NEON

Между двумя моделями проходит главный архитектурный водораздел, и разбирать его надо как инженерный компромисс, а не как соревнование. Ширина зашита в кодировку команды и в тип данных языка. __m256 — это ровно 256 бит, всегда. Компилятор точно знает, сколько элементов обрабатывается за итерацию, и может на этом строить раскрутку, планирование и константные перестановки.

Цена системная и состоит из трёх частей. Хвост цикла существует всегда: массив из 1003 элементов при ширине 8 оставляет три элемента, которые надо обработать отдельно. Каждая новая ширина — новый набор команд: MMX, SSE, SSE2, AVX, AVX2, AVX-512 — не расширения одного механизма, а последовательные наборы с собственными мнемониками и кодировкой, под каждый из которых программу приходится как минимум пересобирать. И каждое поколение навсегда занимает место в кодовом пространстве команд (глава 03).

Переменная длина: RISC-V V и ARM SVE

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

В RISC-V V это команда vsetvli: она получает желаемое число элементов, тип элемента и коэффициент группировки регистров, а возвращает фактическую длину vl — минимум из запрошенного и того, что реализация может за раз.

  # saxpy на RVV 1.0: y[i] = a * x[i] + y[i]
  # a0 = сколько элементов осталось, a1 = указатель x, a2 = указатель y, fa0 = a
loop:
    vsetvli   t0, a0, e32, m1, ta, ma   # прошу a0 элементов по 32 бита;
                                        # железо вернуло t0 — сколько взяло
    vle32.v   v0, (a1)                  # загрузка ровно t0 элементов x
    vle32.v   v1, (a2)                  # загрузка ровно t0 элементов y
    vfmacc.vf v1, fa0, v0               # v1 += a * v0 по всем активным полосам
    vse32.v   v1, (a2)                  # запись ровно t0 элементов
    slli      t1, t0, 2                 # t0 элементов = t0 * 4 байта
    add       a1, a1, t1
    add       a2, a2, t1
    sub       a0, a0, t0                # осталось обработать
    bnez      a0, loop                  # хвоста нет: последний проход просто короче

Обратите внимание, чего в этом коде нет: числа 4, 8 или 16. Нет проверки «хватает ли на полный вектор», нет второго скалярного цикла, нет константы ширины. Тот же двоичный код корректно работает и на реализации со 128-битными регистрами, и на реализации с 1024-битными, причём на второй он выполнит меньше итераций. Это и называется vector-length agnostic — код, безразличный к длине вектора.

Последняя итерация в этой схеме автоматически оказывается короткой: vsetvli вернёт остаток, и все команды тела отработают ровно по нему. Хвост цикла как отдельная конструкция исчезает, превращаясь в обычную итерацию.

ARM SVE решает ту же задачу иначе: длина фиксирована для конкретной реализации, но неизвестна на этапе компиляции. Цикл строится вокруг предиката, который команда whilelo формирует из текущего индекса и границы: активны те полосы, чей индекс меньше границы. Приращение индекса делается командами вида incw/incb, которые прибавляют текущую ширину вектора в элементах или байтах. Результат тот же: одна кодовая ветка, никаких констант ширины.

Инженерная цена модели переменной длины тоже вполне конкретна. Компилятор теряет константу: раскрутка цикла на известное число итераций, статическое планирование под конкретную ширину и константные перестановки либо усложняются, либо становятся невозможными. Векторный тип перестаёт быть обычным типом языка: в SVE типы вроде svfloat32_t объявлены «безразмерными» — их нельзя положить в структуру, нельзя завести их массив, нельзя вычислить sizeof, а компилятор обязан уметь выделять под них стек динамически. Появляется новое состояние, которое надо сохранять: регистр длины вектора и настройки типа элемента становятся частью контекста, и переключение контекста дорожает. Наконец, алгоритмы с перестановками между полосами усложняются: пока код не знает ширину, он не может сказать «переставь элементы 0..7 вот в таком порядке».

Идея не новая: это возврат к векторным машинам 1970-х, где регистр длины вектора был штатным механизмом. Её похоронили в 1990-х вместе с суперкомпьютерами на векторных процессорах и достали обратно, когда стало ясно, во что обходится фрагментация наборов фиксированной ширины: в универсальных ядрах она шла непрерывно с середины 1990-х — MMX, SSE, SSE2, AVX, AVX2, AVX-512, — и попытка собрать её обратно в уровни началась только в 2020-х.

Маскирование: как векторизуется if

Условие внутри цикла — главная головная боль векторизации. Полосы изолированы и исполняют одну и ту же команду; «пойти в разные ветки» они физически не могут. Решение — предикация: вычислить обе ветки для всех полос, а затем выбрать результат по маске. Маска — это вектор битов, по одному на полосу, полученный из сравнения.

// Скалярный вариант: одно ветвление на элемент
for (size_t i = 0; i < n; ++i)
    out[i] = (a[i] > 0.0f) ? a[i] * k : a[i] * m;

// Векторный вариант, по шагам: маска = сравнить(вектор a, ноль) — один бит на полосу;
// t1 = вектор a * k и t2 = вектор a * m — оба посчитаны для ВСЕХ полос без исключения;
// результат = выбрать(маска, t1, t2) — смешивание по маске.

Отсюда — правило стоимости, которое надо держать в голове всегда: цена ветвления в векторном коде равна сумме обеих веток, а не максимуму из них. Скалярный код с хорошо предсказуемым условием выполняет одну ветку и почти ничего не платит (глава 08); векторный выполняет обе всегда. Практическое следствие: если одна ветка дешёвая, а другая очень дорогая и берётся редко, векторизация может проиграть скалярному коду. При этом сами маски прошли эволюцию от «вектора байт с нулями и единицами» до первоклассных регистров.

Регистры масок дают три вещи, которых не было раньше. Первая — экономия векторных регистров. Вторая — выбор семантики неактивной полосы: обнулить или сохранить старое значение; второе позволяет накапливать результат по частям без лишних смешиваний. Третья, самая важная и самая недооценённая — подавление отказов. Маскированная загрузка не вызывает исключения для неактивных полос, даже если их адрес не отображён. Именно это делает законным то, что иначе было бы неопределённым поведением: можно взять полный вектор, выходящий за конец массива, замаскировать лишние полосы — и не получить отказ страницы. Без этого механизма компилятор обязан заканчивать цикл скалярно, и это одна из главных причин, по которой автовекторизатор так часто отказывается работать с короткими массивами (Неопределённое поведение).

Родственный механизм — загрузки «только первый отказ» (vle8ff в RVV, first-faulting loads в SVE). Они позволяют векторизовать поиск по строке неизвестной длины: загрузить сколько получится, а если часть вектора попала на неотображённую страницу, вернуть только успешную часть и сообщить об этом. Так векторизуют strlen и разбор текста.

Память: где на самом деле теряется выигрыш

Выравнивание

Векторный блок бесполезен, если данные не приходят достаточно быстро; практически всегда именно память, а не арифметика, определяет итоговое ускорение. Исторически SIMD-загрузка требовала, чтобы адрес был кратен ширине регистра, а иначе — исключение. На современных высокопроизводительных x86 и ARM-ядрах это в прошлом: невыровненная загрузка стоит столько же, сколько выровненная, пока она не пересекает границу строки кэша. Пересечение превращает одно обращение в два и требует склейки — на ядрах последнего десятилетия это порядок нескольких дополнительных тактов. Пересечение границы страницы дороже ещё заметно.

Отсюда практическое правило: выравнивать данные стоит не на ширину вектора, а на строку кэша — порядка 64 байт на большинстве архитектур, 128 на некоторых. Тогда 512-битная загрузка по выровненному адресу никогда не расщепляется. Как устроены строки и почему их длина именно такая — глава 09.

Gather и scatter

Команда gather берёт вектор индексов и собирает элементы из произвольных мест. Выглядит как магия, а на деле это одна команда, которая внутри разворачивается в столько обращений, сколько разных строк кэша она задела. Кэш первого уровня имеет один-два порта загрузки; восемь элементов из восьми разных строк — это восемь последовательных обращений плюс сборка результата.

Порядок величины для настольных и серверных x86 конца 2010-х: собрать вектор из 8 элементов по 32 бита с произвольными индексами стоит примерно столько же, сколько 8 обычных загрузок, плюс накладные расходы самой команды. Выигрыш gather — в компактности кода и в том, что цикл остаётся векторным целиком; выигрыша по времени, как правило, нет. Если индексы случайно оказались в одной строке кэша, некоторые реализации это используют, а некоторые нет — рассчитывать на это нельзя. scatter, обратная операция разброса, ещё хуже: к тем же проблемам добавляется разбор конфликтов, когда две полосы пишут по одному адресу.

Раскладка данных решает больше, чем интринсики

Это главный практический тезис главы. Массив структур почти всегда убивает векторизацию, а структура массивов почти всегда её включает — и никакие интринсики этого не компенсируют.

// Массив структур (AoS): координаты вперемешку. Чтобы собрать 8 значений x, нужен
// доступ с шагом 12 байт — gather либо три загрузки и цепочка перестановок.
// Плюс в кэш притащены y и z, которые в этом цикле не нужны.
struct ParticleAoS { float x, y, z; } pa[N];

// Структура массивов (SoA): 8 значений x — одна загрузка 32 байт, ни одного лишнего байта.
struct ParticlesSoA { float *x, *y, *z; };   // по N элементов каждый

// Гибрид (AoSoA): блоки по ширине вектора. Локальность AoS для одной частицы
// плюс непрерывность SoA внутри блока.
#define W 8
struct ParticleBlock { float x[W], y[W], z[W]; } blocks[(N + W - 1) / W];

Разница здесь не только в числе команд. При AoS обработка 8 частиц по координате x затрагивает 96 байт памяти вместо 32 — то есть в три раза больше строк кэша и в три раза больше трафика к памяти, из которого две трети выбрасывается. При векторе шириной 512 бит пропорция сохраняется, а абсолютные потери растут. Подробный разбор локальности — Кэш и локальность. Правило, которое стоит выучить: сначала раскладка, потом векторизация. Если данные лежат правильно, автовекторизатор чаще всего справится сам; если неправильно — не спасут и интринсики.

Редукции: почему сумма по вектору ломается о плавающую точку

Полосы изолированы, а сумма всех элементов вектора — операция поперёк полос. Такие операции называются горизонтальными, и они устроены принципиально дороже.

Сложение 8 значений в 256-битном регистре — лог-дерево из перестановок и сложений:
  шаг 1: верхние 128 бит + нижние  → 4 значения
  шаг 2: верхние  64 бита + нижние → 2 значения
  шаг 3: верхние  32 бита + нижние → 1 значение
Каждый шаг зависит от предыдущего: log2(8) = 3 пары команд подряд, без перекрытия.

Отсюда единственно правильный шаблон: не сворачивать вектор на каждой итерации. Держим вектор частичных сумм всё время цикла и сворачиваем один раз в конце.

// Тело цикла: O(n / W) векторных FMA, независимых внутри итерации.
// Эпилог: O(log2 W) операций ровно один раз. Память: O(1) сверх входных массивов.
float dot(const float *restrict a, const float *restrict b, size_t n) {
    float acc[8] = {0};                      // концептуально — вектор частичных сумм
    size_t i = 0;
    for (; i + 8 <= n; i += 8)
        for (int j = 0; j < 8; ++j)
            acc[j] += a[i + j] * b[i + j];   // 8 независимых цепочек накопления
    float s = 0;
    for (int j = 0; j < 8; ++j) s += acc[j]; // свёртка один раз, в самом конце
    for (; i < n; ++i) s += a[i] * b[i];     // хвост
    return s;
}

Сложность: O(n) по времени с константой, уменьшенной примерно в W раз на счётной части, и O(1) дополнительной памяти. Класс сложности векторизация не меняет никогда — она меняет только константу. Это стоит помнить, когда обсуждение производительности сползает в интринсики: алгоритм с лучшей асимптотикой побьёт векторизованный переборный вариант на достаточно больших данных при любой ширине регистра.

А теперь самое важное. Сложение чисел с плавающей точкой не ассоциативно: (a+b)+c и a+(b+c) дают разные результаты, потому что каждое сложение округляет. Восемь частичных сумм складывают элементы в другом порядке, чем последовательный цикл, — значит, дают другой ответ в младших битах.

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

gcc -O3 -ffast-math prog.c          # пушка по площадям: переассоциация, игнор NaN
                                    # и знака нуля, денормалы выключены — прощай IEEE 754
gcc -O3 -fassociative-math -fno-signed-zeros -fno-trapping-math prog.c   # точечнее
gcc -O3 -fopenmp-simd prog.c        # разрешение выдаётся в исходнике, на конкретный цикл,
                                    # директивой "#pragma omp simd reduction(+:s)"

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

Хвост цикла

Массив редко делится на ширину вектора нацело. Есть четыре стратегии.

Скалярный хвост — не мелочь. При ширине 16 и массиве из 20 элементов четыре пятых работы делает одна векторная итерация, а остаток обрабатывается скалярно, то есть в шестнадцать раз медленнее на элемент. Именно поэтому компиляторы вставляют перед векторным циклом проверку вида «если элементов меньше порога — сразу в скалярную версию»: на коротких массивах пролог, эпилог и проверки перекрывают весь выигрыш.

Автовекторизация: что мешает компилятору

Первый инструмент векторизации — не интринсики, а отчёт компилятора; второй — профилировщик, который показывает, стоит ли вообще трогать этот цикл (Профилирование CPU). Оптимизатор уже пытался векторизовать ваш цикл и, если не смог, готов сказать почему (Оптимизации компилятора).

gcc -O3 -march=x86-64-v3 -fopt-info-vec -fopt-info-vec-missed -c hot.c   # что и почему НЕ
clang -O3 -march=x86-64-v3 -Rpass=loop-vectorize \
      -Rpass-missed=loop-vectorize -Rpass-analysis=loop-vectorize -c hot.c
clang -O3 -fsave-optimization-record -c hot.c   # машиночитаемый отчёт hot.opt.yaml

Типичные диагнозы и что с ними делать.

Алиасинг указателей. Компилятор обязан считать, что out и in могут указывать на пересекающиеся области, — тогда порядок обращений менять нельзя. Обычно он выкручивается версионированием: генерирует и векторный, и скалярный вариант плюс проверку перекрытия во время работы. Это стоит кода и одной ветки на входе в цикл. restrict снимает вопрос:

void scale(float *out, const float *in, float k, size_t n);  // перекрытие возможно →
                                                             // проверка во время работы
// Обещание «эти области не пересекаются»: векторизуется без проверок.
// Нарушение обещания — неопределённое поведение, а не просто ошибка.
void scale(float *restrict out, const float *restrict in, float k, size_t n) {
    for (size_t i = 0; i < n; ++i) out[i] = in[i] * k;
}

Зависимость, переносимая по итерациям. a[i] = a[i-1] + b[i] — настоящая зависимость с расстоянием 1, векторизации не поддаётся вообще. А вот a[i] = a[i-8] + b[i] — расстояние 8, и цикл прекрасно векторизуется при ширине до 8 элементов включительно. Компилятор считает расстояния сам; ваша задача — не создавать зависимостей там, где их могло не быть. Отдельный случай той же болезни — редукция по плавающей точке: она разобрана выше и лечится директивой #pragma omp simd reduction.

Вызов функции внутри цикла. Если функция не встроилась, векторизация невозможна: непонятно, что она делает с памятью. Исключения — математические функции, для которых существуют векторные реализации (libmvec в glibc, подключается вместе с -ffast-math или через OpenMP SIMD), и собственные функции, помеченные #pragma omp declare simd.

Непредсказуемое число итераций само по себе не мешает: компилятор сгенерирует проверку и хвост. Мешает выход из цикла посередине — break, return, исключение; такие циклы векторизуются только специальными средствами вроде загрузок с подавлением отказов. Переменный или неизвестный шаг — a[idx[i]] превращается в gather, и компилятор часто решает, что оно того не стоит; обычно он прав.

Знаковость индекса. Тонкий момент на стыке с языком: переполнение знакового int — это неопределённое поведение, поэтому компилятор вправе считать, что i не переполнится, и доказать число итераций. С unsigned переполнение определено как заворот, доказательство ломается, и цикл может не векторизоваться. Один из редких случаев, когда UB работает на пользу (Неопределённое поведение).

Практика: интринсики и что из них получается

Когда отчёт исчерпан, а цикл действительно того стоит, пишут интринсики. Вот честный saxpy под AVX2 с FMA.

#include <immintrin.h>
#include <stddef.h>

// y[i] = a * x[i] + y[i], восемь элементов float32 за итерацию.
// Время O(n) с константой примерно в 8 раз меньше на счётной части, память O(1).
void saxpy_avx2(float a, const float *restrict x, float *restrict y, size_t n) {
    const __m256 va = _mm256_set1_ps(a);         // скаляр разослан по всем восьми полосам
    size_t i = 0;
    for (; i + 8 <= n; i += 8) {
        __m256 vx = _mm256_loadu_ps(x + i);      // невыровненная загрузка 32 байт
        __m256 vy = _mm256_loadu_ps(y + i);
        vy = _mm256_fmadd_ps(va, vx, vy);        // одно округление на пару операций
        _mm256_storeu_ps(y + i, vy);
    }
    for (; i < n; ++i) y[i] = a * x[i] + y[i];   // скалярный хвост: до 7 элементов
}

Собирается как gcc -O3 -mavx2 -mfma. Ядро цикла превращается примерно в это:

.L4:
        vmovups ymm1, YMMWORD PTR [rsi+rax]     ; 8 x float32 из x
        vmovups ymm2, YMMWORD PTR [rdx+rax]     ; 8 x float32 из y
        vfmadd132ps ymm1, ymm2, ymm0            ; ymm1 = ymm1 * ymm0 + ymm2
        vmovups YMMWORD PTR [rdx+rax], ymm1     ; запись результата
        add     rax, 32                         ; шаг цикла — 32 байта, а не 4
        cmp     rax, rcx
        jne     .L4

Семь команд на восемь элементов против примерно шести команд на один элемент в скалярной версии. Как читать такой листинг — Основы ассемблера; удобнее всего смотреть на Compiler Explorer, переключая целевую архитектуру. А теперь о том, что здесь принципиально плохо и о чём стоит подумать до того, как писать такой код: восьмёрка зашита в трёх местах, тип регистра зашит в сигнатуре, а на машине с 512-битными регистрами этот код навсегда останется вдвое медленнее возможного. Компилятор, встретив интринсики, уже не станет ничего перепланировать — вы взяли ответственность на себя. Промежуточный вариант — библиотеки-абстракции над разными наборами: Google Highway, Vector Class Library, ISPC, в перспективе std::simd в стандартной библиотеке C++.

Тот же механизм этажом ниже: NumPy

«Векторизация» в высокоуровневых библиотеках — не метафора, а ровно тот же механизм, просто вызванный из другого места.

import numpy as np

n = 1 << 24                                # 16 777 216 элементов по 4 байта = 64 МиБ
x = np.random.rand(n).astype(np.float32)
y = np.random.rand(n).astype(np.float32)
a = np.float32(2.5)

z = a * x + y              # три прохода по памяти и два временных массива по 64 МиБ:
                           # t1 = a * x, затем z = t1 + y

z = np.empty_like(x)       # тот же результат за два прохода и ноль аллокаций в цикле
np.multiply(x, a, out=z)   # проход 1: читаем x, пишем z
np.add(z, y, out=z)        # проход 2: читаем z и y, пишем z

Каждая операция NumPy — это цикл на C внутри универсальной функции, и в современных сборках этот цикл собран с SIMD-ядрами под конкретный набор расширений, выбираемыми во время работы (документация NumPy по SIMD). То есть между вашим a * x + y и командой vfmadd132ps ровно один слой. Практический вывод для Python отсюда неочевидный и очень полезный: на больших массивах выигрыш даёт не столько сама векторизация — она уже есть, — сколько сокращение числа проходов по памяти и числа временных массивов. Три прохода по 64 МиБ гарантированно вылетают из кэша последнего уровня, и время определяется трафиком к оперативной памяти, а не арифметикой. Для сложных выражений это лечат слиянием операций — numexpr, numba, JIT в других библиотеках; подробнее — Производительность Python. Та же логика управляет производительностью вывода нейросетей, где почти вся работа сводится к плотным матричным операциям (Инференс и деплой).

Инженерные пределы

Потолок пропускной способности памяти

Самая частая причина разочарования в SIMD: векторизованный цикл требует данных в W раз быстрее, а память отдаёт их с той же скоростью, что и раньше. Полезная величина — арифметическая интенсивность: сколько операций приходится на байт, прочитанный из памяти. У y[i] += a * x[i] она чудовищно низкая: два байта чтения и четыре байта записи-чтения на две операции. У плотного умножения матриц с правильной блокировкой — высокая, потому что каждый загруженный блок переиспользуется многократно.

Соотношение простое: высокопроизводительное ядро 2020-х выполняет за такт порядка десятков операций с плавающей точкой одинарной точности, а устойчивая пропускная способность памяти в пересчёте на ядро — порядок единиц байт за такт (Память и производительность). Значит, чтобы упереться в счёт, а не в байты, нужны десятки операций на каждый прочитанный байт. Это очень много, и большинство реальных циклов до этого не дотягивают. Модель, формализующая рассуждение, называется roofline (Williams, Waterman, Patterson, CACM 2009).

Снижение частоты при широких векторах

Физика простая и уже разобранная в главе 13: динамическая мощность пропорциональна числу переключающихся вентилей, ёмкости, квадрату напряжения и частоте. 512-битная FMA переключает на порядок больше транзисторов за такт, чем скалярное сложение, и делает это на большой площади кристалла. Чип обязан остаться в тепловом и токовом бюджете — значит, при интенсивной широкой векторной нагрузке он снижает частоту. На серверных x86 конца 2010-х это было формализовано в виде «лицензий» на частоту: устойчивая тяжёлая работа с 512-битными операциями опускала частоту всех ядер на величину порядка десятков процентов. Числа зависят от поколения, модели и числа активных ядер, и приводить их как константу нельзя — важен сам механизм.

Практическое следствие контринтуитивно: микробенчмарк векторного ядра может показать отличный результат, а приложение целиком — просесть, потому что снизилась частота всех остальных потоков. Отсюда правило: узкие векторные вставки в общем коде часто не окупаются, а измерять нужно систему целиком, а не цикл (Измерения).

Фрагментация наборов и диспетчеризация во время работы

Расширений много, они не образуют цепочку вложенности, и целевая машина заранее неизвестна. Индустрия отвечает уровнями и профилями: x86-64-v2/v3/v4 у x86, профили RVA22/RVA23 у RISC-V — фиксированные наборы, под которые можно собирать дистрибутив (глава 03, глава 05).

Когда одного уровня мало, применяют диспетчеризацию во время работы — вручную или силами компилятора.

// 1. Ручная проверка: один раз на старте, дальше вызовы через указатель.
typedef void (*saxpy_fn)(float, const float *, float *, size_t);

static saxpy_fn pick_saxpy(void) {
    if (__builtin_cpu_supports("avx512f")) return saxpy_avx512;
    if (__builtin_cpu_supports("avx2"))    return saxpy_avx2;
    return saxpy_scalar;
}

// 2. Function multiversioning: компилятор сам собирает варианты и сам
//    разрешает вызов через ifunc при загрузке программы.
__attribute__((target_clones("default", "avx2", "avx512f")))
void scale(float *restrict out, const float *restrict in, float k, size_t n) {
    for (size_t i = 0; i < n; ++i) out[i] = in[i] * k;
}

Механизм разрешения на Linux — ifunc: динамический компоновщик один раз при загрузке вызывает функцию-резолвер и подставляет нужный адрес в таблицу переходов. Так устроены memcpy и strlen в glibc: одна и та же библиотека содержит несколько реализаций под разные наборы. Цена — косвенный вызов и невозможность встроить функцию; поэтому диспетчеризовать надо крупные ядра, а не мелкие операции внутри горячего цикла. Отсюда правило, которое стоит принять как безусловное: никогда не раздавайте пользователям бинарники, собранные с -march=native. На машине сборки они работают, на машине пользователя падают с недопустимой командой, и падают в произвольном месте без внятного сообщения. Правильно: выбрать явный базовый уровень (-march=x86-64-v2, например), отдельно указать -mtune под типичную машину и добавить диспетчеризацию там, где выигрыш это оправдывает (Выбор железа).

Наконец, о границах. Векторный блок — это параллелизм по данным внутри одного ядра и внутри одного потока команд. Как только требуется на порядки больше полос, специализированное матричное железо или отдельная иерархия памяти, начинается другая архитектура: GPU, NPU, матричные блоки и систолические массивы. Это глава 12. Как только речь заходит о нескольких потоках команд, согласованности кэшей и барьерах памяти — это глава 11.

Типичные заблуждения

«SIMD ускоряет всё в W раз». Только счётную часть и только если цикл в неё упирался. Если он упирался в память, ускорение будет околонулевым, и это самый частый исход.

«Интринсики всегда быстрее автовекторизации». Часто медленнее. Они замораживают ширину, блокируют дальнейшие преобразования компилятора и устаревают вместе с расширением, под которое написаны. Начинать надо с отчёта векторизатора.

«Расширение шире — значит быстрее». Не обязательно: играют роль снижение частоты, число исполнительных портов, ширина физического тракта и то, упирается ли код в память вообще. Единственный способ узнать — измерить на целевом железе.

«Векторизация не меняет результат». Для целых — не меняет. Для плавающей точки меняет: другой порядок суммирования, а при FMA ещё и другое число округлений. Если результат сравнивается побитово, это ломается при смене машины.

«Выравнивание больше не важно». Важно, но на уровне строки кэша, а не ширины регистра: расщеплённые обращения и пересечения страниц никуда не делись. Из той же оперы — «gather дешевле цикла загрузок»: по коду он компактнее, по времени обычно нет, потому что обращений в кэш ровно столько, сколько задето разных строк.

«Компилятор не векторизовал — значит, надо писать интринсики». Сначала прочитайте, почему. В подавляющем большинстве случаев это алиасинг, зависимость по итерациям или запрещённая переассоциация — и всё это чинится в исходнике, без единого интринсика.

«SIMD заменяет многопоточность». Это ортогональные множители: векторный код без потоков использует одно ядро, потоки без векторизации — долю каждого ядра.

«Маска — это способ пропустить работу». Это способ пропустить запись и отказ, но не вычисление: неактивные полосы всё равно занимают исполнительный блок и потребляют энергию.

Мини-итог

• Параллелизм по данным объявляется статически и стоит дёшево: фиксированная стоимость команды — фронтенд, переименование, планирование, фиксация — делится на число элементов. Именно поэтому SIMD выигрывает джоули, а не только секунды.
• Полосы изолированы; всё, что требует движения данных поперёк полос — редукции, перестановки, gather, — принципиально дороже поэлементных операций.
• Ширина в контракте не равна ширине в кремнии: одна команда может исполняться за такт или за несколько, и код об этом не узнает.
• Модель переменной длины убирает хвост цикла и перекомпиляцию под каждую ширину ценой потери константы на этапе компиляции и усложнения типов языка. Модель фиксированной ширины даёт компилятору больше знаний ценой фрагментации наборов.
• Маски превращают ветвление в вычисление обеих веток; цена условия становится суммой, а не максимумом. Подавление отказов на неактивных полосах — то, что делает законным маскированный хвост.
• Раскладка данных важнее интринсиков: AoS против SoA меняет и число команд, и объём трафика, а gather не спасает.
• Редукция по плавающей точке требует явного разрешения на переассоциацию, потому что меняет результат. Лучше выдавать это разрешение на конкретный цикл, а не всей программе.
• Потолок почти всегда — пропускная способность памяти: считайте операции на байт прежде, чем считать полосы. Широкие векторы вдобавок снижают частоту, а наборы расширений фрагментированы — отсюда уровни, профили и диспетчеризация во время работы вместо -march=native.

Источники

• Hennessy J., Patterson D. Computer Architecture: A Quantitative Approach — глава 4 целиком посвящена параллелизму по данным: векторные архитектуры, SIMD-расширения и GPU разобраны в одной системе понятий.
• Hameed R. et al. «Understanding Sources of Inefficiency in General-Purpose Chips», ISCA 2010 — https://dl.acm.org/doi/10.1145/1815961.1815968: откуда берётся энергия в универсальном ядре и сколько её достаётся собственно арифметике.
• Williams S., Waterman A., Patterson D. «Roofline: An Insightful Visual Performance Model», CACM 2009 — https://dl.acm.org/doi/10.1145/1498765.1498785: как понять, упирается ли код в счёт или в память, до того, как что-то оптимизировать.
• Intel Intrinsics Guide — справочник по интринсикам x86 с семантикой и требуемым набором; Arm Neon и SVE intrinsics — то же для AArch64, с фильтром по расширению.
• The RISC-V Instruction Set Manual — том с расширением V: vsetvl, группировка регистров, маскирование, загрузки с подавлением отказов.
• Auto-vectorization in GCC и Auto-Vectorization in LLVM — что векторизаторы умеют, какие формы циклов распознают и какие диагностики выдают.
• Fog A. Optimizing software in C++ и таблицы латентностей — https://www.agner.org/optimize/; измеренные пропускные способности векторных команд — https://uops.info/.
• NumPy SIMD optimizations — как устроен слой универсальных функций и выбор ядер во время работы.
• Google Highway — переносимая абстракция над SIMD, включая наборы переменной длины; хороший пример того, как выглядит код без зашитой ширины.

Что дальше

Векторный блок выжимает параллелизм из данных внутри одного потока команд, и его потолок мы теперь понимаем: изолированные полосы, изолированный регистровый файл, одна последовательность команд и общая для всех иерархия памяти. Следующий множитель производительности лежит в другой плоскости — несколько независимых потоков команд, работающих над общей памятью. И тут же возникают вопросы, которых у SIMD не было вовсе: как разные ядра видят записи друг друга, в каком порядке эти записи становятся видимыми, чем за это платит железо и какие гарантии оно обязано дать программисту.

Многоядерность: согласованность, барьеры, модели памяти