Pull to refresh

Comments 72

Не помешала бы сводная таблица в конце статьи.
А дон знает толк в ветряных мельницах! И результат весьма впечатляет, почти как у пресловутого дона — о нем мир помнит, а мельницы сгинули.

Можно поинтересоваться судьбой софтового рендера?
Конечно. Пока посмотреть можно тут: www.youtube.com/watch?v=gBp1PoWzeO8
Сначала писал статью. Потом понял, что одна не получится. Потом понял, что сразу всё писать сложно. Эта статья как раз вышла из неё. Решил посмотреть как пойдёт и будет ли интерес.
Если кратко, то: Quake2 1600x1200, MT (multythreaded), SSE3, i7-3770, 30 fps.
Переделываю с нуля под AVX. Цикл статьей планируется по SSE варианту.
Минутку. А где unit test из 100500 повторений вызова функции, замером времени, и расчётом реального количества инструкций на функцию?
Вы настолько верите IACA? Тогда вас может ждать много удивительных открытий при погружении в реальный мир.
1. Да, верю. По опыту. Результаты как минимум коррелируют. Ну будут у меня цифры немного другие, сути это не меняет. А IACA разбирается в архитектуре Intel лучше меня уж точно.
2. Про реальный мир вы зря: я давно уже использую SIMD и тестирую и так, и так. Тут решил обойтись.
3. Не могу протестировать всё. У меня максимум AVX (i7-3770).
4. Это будет синтетика. Я уже отмечал, что тысячами матрицы я не умножаю.
5. Однако вы правы: пункт 4 имеет место быть, но IACA не может учитывать реальную работу с данными, памятью. Поэтому стоит написать. Интересно будет не только вам.
3. Не могу протестировать всё. У меня максимум AVX (i7-3770).
Однако как раз AVX512 это особенно важно. Дело в том, что увидев AVX512-инструкцию процессор «пугается» и снижает частоту. Сразу, превентивно. Потому скорость в тактах становится куда менее важной чем для предыдущих SIMD-команд.
> На каждый элемент мы делаем 4 умножения и 3 сложения. В сумме это 256 умножений и 192 сложения.

В матрице 16 элементов, в сумме это 64 умножения и 48 сложений.
Спасибо!!! Как проглядел непонятно. Исправлю.
Кстати я понял почему так вышло. Я взял те самые 4*4*4 из текста, предназначенные для одного умножения во внутреннем цикле, и перенёс их на mulps, в котором 4 умножения.
Хорошее велосипедостроительство. Но реализация libxsmm мне нравиться больше:
Small matrix multiplication kernels (SMMs) are generated for Intel SSE, Intel AVX, Intel AVX2, IMCI (KNCni) for Intel Xeon Phi coprocessors (KNC), and Intel AVX‑512 as found in the Intel Xeon Phi processor family (KNL, KNM) and Intel Xeon processors (SKX). Highly optimized code for small convolutions is targeting Intel AVX2 and Intel AVX‑512, whereas other targets can automatically leverage specialized SMMs to perform convolutions.
Универсальные решения как правило проигрывают специализированным а libxsmm даже не самая быстрая реализация из существующих. Если нужно постоянно работать с матрицами 4x4 а не зоопарком где фигурируют перемножения матриц 5x12 на 12x7 то использование кастомного кода выглядит разумнее. Но было бы конечно любопытно взглянуть на код который libxsmm сгенерирует для 4x4 матрицы из примера автора
А можно утверждение про «не самая быстрая» чем то подкрепить?
Незачем тащить с собой монструозную некроссплатформенную библиотеку с нетривиальной установкой и подключением для решения мелкой проблемы, если нет в этом особой необходимости. К тому же, насколько я понял, libxsmm принуждает к использованию собственных типов для хранения матриц, что может потребовать либо переписывание кода, либо ненужное копирование данных.

Интересно, удалось ли обогнать оптимизатор GCC с флагами "-m64 -O3 -march=native -funroll-loops -ffast-math"? Перемножение двух целочисленных uint32_t матриц 4096x4096 у меня заняло ~6.8 сек на одном ядре Intel Xeon E3-1220 v5, задача здесь. Алгоритм — слегка оптимизированный O(n^3), отданный на растерзание компилятору.

GCC плохо/либо же ужасно умеет использовать такие инструкции. Хотя по времени выигрыш может оказаться и не очень большим.

Статью стоит рассматривать как мини введение в векторные инструкции, а не как конкретную реализацию умножения матриц либо еще какой конкретной задачи.
Согласен. Если бы компиляторы так умели, я бы точно не стал заморачиваться.
А у Вас нет случаем бенчмарков для сравнения Вашего кода и того, что нагенерировал GCC/Clang/ICC для этого случая?

Потому что мне кажется, что тот же ICC очень неплохо занимается векторизацией. Было бы интересно посмотреть на результаты бенчмарка.
Бенчмарки делаю, но пока не очень получается. Пока погоду показывают. Ни линукса, ни ICC (он вроде платный) у меня под руками нет. Потом когда будут бенчмарки, затестю на том, что подвернётся.
А пока разве только разные сайты использовать, где можно онлайн скомпилить код любым компилятором, и посмотреть результат. Но кому интересно, уже наверняка это сделали. Код то я выложил.
Какие именно проблемы с бенчмарками? Может что подсказать стоит?
ICC компилятор имеет уже давно бесплатную(возобновляемую каждые 90 дней) даже коммерческую лицензию. Просто заходите и качаете. К тому же Intel System Studio работает и на Windows, так что можете провести измерения msvc vs icc.

Как уже отметили, иногда сложно предсказать, какие способы будут работать быстрее, глядя просто на то, что там нагенерировал компилятор — а вот по бенчмаркам уже проще смотреть.
i7-3770:
x86: тормозят все реализации раза в 2 от ожидаемых
x64: все SSE быстрые, все AVX медленнее в 20-30 от ожидаемых
На 8700K сказали не проявляется.
AVX_v1 быстрее без стриминга. AVX_v2 быстрее со стримингом.
Справедливости ради стоит отметить, что в GCC 8.2, который вышел прям совсем недавно, улучшили генерацию кода с AVX-512
4096x4096 это по сути 4k экран. Тут уже SIMD + MT (multithreading). А у меня речь про чистый SIMD пока. А сама по себе задача интересная.
AVX+FMA вариант (финальный)… В данном случае нам требуется вычислить a*b + c*d + e*f + g*h.

Выполняется одной командой _mm_dp_ps из SSE4 — она гораздо быстрее
В моём рабочем коде классов классе пишу вроде
__forceinline __m128 __vectorcall operator + (__m128 const) const noexcept;

Однако здесь всё итак работает хорошо. Всегда проверяю сгенерированный ассемблерный код. Насчёт *assume_aligned я посмотрю что это. Скорее всего то, чего мне сильно и давно не хватает. Но если оно linux/intel_compiler/… specific, то мне, к сожалению, не подойдёт.
Посмотрел. Да, это что мне нужно. Увы, ICC specific.

__assume_aligned — это Intel, __builtin_assume_aligned — это GCC.

а почему бы не использовать DPPS _mm_dp_ps (SSE4.1/AVX)? готовое скалярное произведение пары векторов
Очень разумный вопрос. Статья по факту представляет чистый положительный опыт. Если собрать весь отрицательный, то она будет минимум в 3 раза больше. Конечно dpps я пробовал когда то. Не оправдала от слова совсем.

Когда SSE4.1 только рекламировали, я уже надеялся, что такая классная инструкция dpps всё упростит. И ждал. Но увы…

1. Она дорогая: 2 такта (и судя по всему не параллелится как mulps), то есть даже без учёта транспонирования, их надо 16 по 2 такта, получается 32!!! (IACA согласна)
1.1 dpps — 4 умножения, 3 сложения
1.2 mulps — 4 умножения, параллелится, это 8 умножений за такт, 16 за два такта, или, грубо, 8 умножений + 3 сложения за 2.5 такта
2. Одну матрицу надо предварительно трансформировать: dpps(строка M, столбец N)
3.1 дополнительная инструкция movss для извлечения результата и работа с памятью на уровне 1-го float
ИЛИ
3.2 дополнительный код для компоновки результата из 4 float в одном xmm

Не стоит думать об инструкции, как вещи в себе. Инструкции нужны данные, они должны быть в нужном представлении как на входе, так и на выходе.

По совокупности факторов её даже Intel не использует в своих примерах матричного умножения.
… Очень разумный вопрос. Статья по факту представляет чистый положительный опыт. Если собрать весь отрицательный, то она будет минимум в 3 раза больше…
А может будет это вторая статья?

По моему опыту, dpps имеет приличный выигрыш, но случаи бывают всякие. И подготовка данных, полностью согласен, стоит времени ЦПУ
Как вариант приложить в коде. Но на статью, пожалуй, это не тянет. Простой набросок без транспонирования дал уже 32 такта.
Я уверен, что всё равно найдутся несогласные. Но для меня это пройденный опыт, и я не хочу на нём останавливаться. Если вы считаете, что можете лучше с dpps, киньте сюда свою реализацию вычисления умножения с ней, которая будет лучше хотя бы 20 тактов (с учётом чтения матриц m/n и записью в результирующую матрицу r), и я публично признаю свою ошибку.
Спасибо за информацию. Погуглил, пишут что снижение частоты зависит от степени нагрева процессора. При равной тепловой нагрузке AVX-код будет выполняться на частоте до ~20% ниже (т.е. 8 реальных тактов автора превратятся в 10 «эквивалентных»), это все равно намного быстрее чем 16 в SSE-варианте.
Подстава заключается в том, что частота падает не только во время выполнения этих инструкций, но и некоторое время после. В реальности отдельная функция может начать работать быстрее, но программа в целом замедлится. Этот эффект наблюдается с новым OpenSSL — по всем бенчмаркам самой библиотеки стало быстрее, однако веб-сервер начинает обслуживать https запросы медленнее.
Это вы ещё AVX512 не видели. Там вообще нет смысла в использовании, если приложение ускоряется меньше чем на 20% по тактам, т.к. частота на эти же 20% может сбрасываться. Относительно версии с AVX/AVX2.
Дополнительно ещё стоит учесть работу с памятью.
Например, разницы при сложении двух очень длинных векторов с использованием SSE и AVX при максимальном распараллеливании не будет, т.к. всё упрётся в пропускную способность памяти.
А AMD не снижают. Стало быть корректно?
Зачем быть заложником одной компании?
У современных процессоров AMD в два раза меньше вычислительных блоков AVX, чем у процессоров Intel. Так что AMD пока в проигрыше.
Смотря каких.
www.agner.org/optimize/blog/read.php?i=838

AMD has four 128-bit units for floating point and vector operations. Two of these can do addition and two can do multiplication. Intel has two 256-bit units, both of which can do addition as well as multiplication. This means that floating point code with scalars or vectors of up to 128 bits will execute on the AMD processor at a maximum rate of four instructions per clock (two additions and two multiplications), while the Intel processor can do only two. With 256-bit vectors, AMD and Intel can both do two instructions per clock.

Т.е. тоже самое.
Но FMA быстрее на Интел.

Intel beats AMD on 256-bit fused multiply-and-add instructions, where AMD can do one while Intel can do two per clock. Intel is also better than AMD on 256-bit memory writes, where Intel has one 256-bit write port while the AMD processor has one 128-bit write port.

Очень интересная статья! Кстати, кроме как тем, что нужен свой формат хранения, матрички из GLM не такие же шустрые? Кажется, там была какая-то поддержка SIMD.

А можно такой вопрос — от неспециалиста в SIMD инструкциях.
По поводу разрядности. Я так понимаю что float в данном случае 32 бита?
А если точность нужна выше? 64 бита на Double. Этот же код сможет работать с double — данные влезут в регистры? Или надо будет модифицировать?
И возможен ли вариант на 80 бит (extended)? Или SIMD не поддерживает 80-битный float сопроцессора x86?
1 регистр SSE — это 4 float или 2 double
1 регистр AVX — это 8 float или 4 double
Код просто будет аналогичным.

80-битные числа SSE/AVX не поддерживаются
Благодарю!
А вот насчет 80 бит жаль :-( Странно что Intel не дала возможности использовать его в своих же расширениях. Пусть бы и не так быстро как float32.
Странно что Intel не дала возможности использовать его в своих же расширениях.
Нет в этом ничего странного. Операции с 80битами сильно дороже, чем с 64мя. А используются редко. Тратить драгоценные транзисторы на то, что не используется никто не будет.

Собственно изначально MMX работал с целыми числами, SSE — с 32-битными float'ами, SSE2 — с 64битными double'ами. И каждый раз собиралась статистика: кто будет использовать эти команды и зачем. Очевидно, что когда «выпекали» SSE3/4/etc толп людей, желающих воспользоваться 90битными флоатами у маркетиногово отдела Intel — не стояло…
Код просто будет аналогичным.

Если бы! AVX — это не в два раза больший SSE модуль, это два SSE модуля, работающих в паре. Если для сложения и умножения это ничего не меняет, то например _mm256_shuffle_ps ведёт себя совсем иначе.

На AVX2 есть _mm256_permute4x64_pd, на AVX действительно придётся допиливать напильником. Но в целом, признаюсь, я далёк от double тематики.
это не в два раза больший SSE модуль

В 2 раза больший.
это два SSE модуля, работающих в паре.

Амд-проблемы.

Я бы не назвал это прям уж проблемой. Иногда ymm нужен именно как спаренный sse (для float индекс будет в диапазоне 0-3), а иногда как один большой регистр с независимыми числами (в этом случае для float уже нужен индекс в диапазоне 0-7). Проблема, что они вторую часть реализуют в следующем наборе, или не реализуют вообще. Как пример, выше упомянутая _mm256_permute4x64_pd.

Вы вообще не поняли о чем я, видимо ни разу не использовали AVX. Я не про какие-то реализации каких-то процессоров, я про весь набор команд и их архитектуру. Посмотрите как работает например _mm256_shuffle_ps — в dst[127:0] могут попасть только биты из a[127:0] и b[127:0], но не из a[255:128] и b[255:128]. То есть нижние 128 бит работают только с нижними 128 битами, старшие со старшими — это больше похоже на два модуля SSE, чем на один модуль SSE, но в два раза больше. И это касается большинства команд. То же самое можно сказать про AVX512 — это скорее 4 модуля SSE, работающих параллельно, чем один очень длинный модуль SSE.

Посмотрите как работает например _mm256_shuffle_ps — в dst[127:0] могут попасть только биты из a[127:0] и b[127:0], но не из a[255:128] и b[255:128]. То есть нижние 128 бит работают только с нижними 128 битами, старшие со старшими — это больше похоже на два модуля SSE, чем на один модуль SSE, но в два раза больше. И это касается большинства команд. То же самое можно сказать про AVX512 — это скорее 4 модуля SSE, работающих параллельно, чем один очень длинный модуль SSE.

Похоже вы ничего не знаете об микроархитектуре, странно как-то вы используете avx и ничего о ней не знаете.

Во-первых, что значит «например» — если это единственный пример(вернее там много похожих инструкций, но суть одна). Зачем вы делаете вид, будто бы у вас много примеров?

Во-вторых ваша аналитика ничего не имеет общего с реальность. И если бы вы действительно использовали тот самых avx, то вы бы знали в чём именно заключается проблема, а заключается проблема в 8битном селекторе.

Опять же, ваши рассуждения имели бы смысл, если данные операции работали в рамках 128 бит. Т.е. если бы из А можно было перемещать в рамках первой части dst, а из Б в рамках второй. Это было бы так, если бы они попросту склеили два 128битных шафла, но это не так. Они почему-то шафлят aa bb aa bb. Неужели там mmx?

И да, никаких «модулей» не существует. Это обывательская терминология, ничего общего с реальностью не имеющая.

Нет. Видимо, вы путаете AVX и AVX2. В AVX все команды работают с регистрами как с двумя 128-битными независимыми половинками (за исключением пары команд, которые обращаются с 128-битными частями как единым целым).


И если бы вы действительно использовали тот самых avx, то вы бы знали в чём именно заключается проблема, а заключается проблема в 8битном селекторе.

Когда 8-битного селектора не хватает, то в качестве аргумента команды может использоваться другой регистр, например, в команде AVX2 _mm256_permutevar8x32_ps.


Они почему-то шафлят aa bb aa bb

Потому что lower(result) = op(lower(aa), lower(bb)), higher(result) = op(higher(aa), higher(bb)).


SELECT4(src, control)
{
  CASE(control[1:0])
    0: tmp[31:0] := src[31:0]
    1: tmp[31:0] := src[63:32]
    2: tmp[31:0] := src[95:64]
    3: tmp[31:0] := src[127:96]
  ESAC RETURN tmp[31:0]
}

dst[31:0] := SELECT4(a[127:0], imm8[1:0])
dst[63:32] := SELECT4(a[127:0], imm8[3:2])
dst[95:64] := SELECT4(b[127:0], imm8[5:4])
dst[127:96] := SELECT4(b[127:0], imm8[7:6])
dst[159:128] := SELECT4(a[255:128], imm8[1:0])
dst[191:160] := SELECT4(a[255:128], imm8[3:2])
dst[223:192] := SELECT4(b[255:128], imm8[5:4])
dst[255:224] := SELECT4(b[255:128], imm8[7:6])
dst[MAX:256] := 0
Во-первых, что значит «например» — если это единственный пример(вернее там много похожих инструкций, но суть одна). Зачем вы делаете вид, будто бы у вас много примеров?

Ооох, да пожалуйста:


__m256 _mm256_unpacklo_ps (__m256 a, __m256 b)

INTERLEAVE_DWORDS(src1[127:0], src2[127:0]){
    dst[31:0] := src1[31:0] 
    dst[63:32] := src2[31:0] 
    dst[95:64] := src1[63:32] 
    dst[127:96] := src2[63:32] 
    RETURN dst[127:0]
}   

dst[127:0] := INTERLEAVE_DWORDS(a[127:0], b[127:0])
dst[255:128] := INTERLEAVE_DWORDS(a[255:128], b[255:128])
dst[MAX:256] := 0

__m256i _mm256_packs_epi32 (__m256i a, __m256i b)

dst[15:0] := Saturate_Int32_To_Int16 (a[31:0])
dst[31:16] := Saturate_Int32_To_Int16 (a[63:32])
dst[47:32] := Saturate_Int32_To_Int16 (a[95:64])
dst[63:48] := Saturate_Int32_To_Int16 (a[127:96])
dst[79:64] := Saturate_Int32_To_Int16 (b[31:0])
dst[95:80] := Saturate_Int32_To_Int16 (b[63:32])
dst[111:96] := Saturate_Int32_To_Int16 (b[95:64])
dst[127:112] := Saturate_Int32_To_Int16 (b[127:96])
dst[143:128] := Saturate_Int32_To_Int16 (a[159:128])
dst[159:144] := Saturate_Int32_To_Int16 (a[191:160])
dst[175:160] := Saturate_Int32_To_Int16 (a[223:192])
dst[191:176] := Saturate_Int32_To_Int16 (a[255:224])
dst[207:192] := Saturate_Int32_To_Int16 (b[159:128])
dst[223:208] := Saturate_Int32_To_Int16 (b[191:160])
dst[239:224] := Saturate_Int32_To_Int16 (b[223:192])
dst[255:240] := Saturate_Int32_To_Int16 (b[255:224])
dst[MAX:256] := 0

__m256 _mm256_dp_ps (__m256 a, __m256 b, const int imm8)

DP(a[127:0], b[127:0], imm8[7:0]) {
    FOR j := 0 to 3
        i := j*32
        IF imm8[(4+j)%8]
            temp[i+31:i] := a[i+31:i] * b[i+31:i]
        ELSE
            temp[i+31:i] := 0
        FI
    ENDFOR

    sum[31:0] := (temp[127:96] + temp[95:64]) + (temp[63:32] + temp[31:0])

    FOR j := 0 to 3
        i := j*32
        IF imm8[j%8]
            tmpdst[i+31:i] := sum[31:0]
        ELSE
            tmpdst[i+31:i] := 0
        FI
    ENDFOR
    RETURN tmpdst[127:0]
}

dst[127:0] := DP(a[127:0], b[127:0], imm8[7:0])
dst[255:128] := DP(a[255:128], b[255:128], imm8[7:0])
dst[MAX:256] := 0

И в все остальные команды арифметических и логических действий соблюдают это правило: результат работы с нижними 128 битами записывается в нижние 128 битов. А теперь ваша очередь показать все остальные инструкции спокойно перекидываются битами из нижних и верхних частей? Такие команды конечно есть, но работают они со всеми 128 битами одновременно, то есть являются инструкциями копирования, а никак не SIMD-операциями.

В основном да. Но исключения всё-таки есть. Та же _mm256_permute4x64_pd.

SELECT4(src, control){
	CASE(control[1:0])
	0:	tmp[63:0] := src[63:0]
	1:	tmp[63:0] := src[127:64]
	2:	tmp[63:0] := src[191:128]
	3:	tmp[63:0] := src[255:192]
	ESAC
	RETURN tmp[63:0]
}

dst[63:0] := SELECT4(a[255:0], imm8[1:0])
dst[127:64] := SELECT4(a[255:0], imm8[3:2])
dst[191:128] := SELECT4(a[255:0], imm8[5:4])
dst[255:192] := SELECT4(a[255:0], imm8[7:6])
dst[MAX:256] := 0
Вот только это команда AVX2, а не AVX.
Увы, да. Однако это проблема тех, кому нужен double.
Почему же? Есть аналогичная команда для float _mm256_permutevar8x32_ps.

Это проблема тех, кто до сих пор пользуется Sandy Bridge, в котором нет AVX2 (например, я).
При миграции с SSE на AVX (повышение разрядности инструкций), сможет с минимальными изменениями. В SSE регистр вмещалось 4 float, в AVX 4 double. Если на той же разрядности регистров нужен double, тогда придётся переписывать.

80 только FPU. SSE: 32*4=64*2=128bit, AVX: 64*4=32*8=256bit.
По сути точность несколько меньше, плата за производительность.
Благодарю, а такой вопрос — если код написан под AVX — а процессор поддерживает только SSE, то в этом случае надо писать два варианта кода для поддержки обоих архитектур или есть возможность рассчитывать на эмуляцию с какой-нибудь стороны?
Для эмуляции есть способ. Где то находил в интернетах SDK для эмуляции. С практической точки зрения SDK интересно лишь для отладки.

Здесь вскользь упомянут вариант, как делать для разных платформ. Пишем разные реализации, и потом после анализа процессора через cpuid подставляем нужный указателю на функцию. И потом работаем с данным функционалом через указатель. Данный подход полезен для больших, ёмких операций. Если использовать его для одной инструкции — точно будет только хуже.

Как пример: функция для умножения матрицы на вектор — станет медленнее, функция для умножения матрицы на массив векторов — имеет смысл.

Гибкий вариант есть здесь optimizing_cpp.pdf. Искать по: CriticalFunction_Dispatch.

Хорошо, буду думать, а то тоже хотелось ускорить умножение матриц 4*4, но немного — не более 2-5 десятков. Но в идеале extended на 80 бит. Значит или многопоточный вариант с FPU или переходить на double для float64.
Разные руководства по оптимизации рекомендуют лишний раз не дёргать кэш для массовых операций сохранения.

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


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

Про это было отмечено в статье, что большие массивы данных нужны скорее на умножении матрицы на массив вершин.

Не понял к чему вы это. Финальный код содержит стриминговые инструкции, которые в данном случае сильно вредны. Ваше нежелание или невозможность тестировать на реальных кейсах подводит вас.

В чём то вы оказались правы. Но как то странно получилось. Stream никак не сказалось на SSE версии, зато сильно портила жизнь AVX1 (ускорение 6-7 против 8-9). Архитектура IvyBridge. В целом корректно было бы IVB такты измерить. А это лишь IACA 2 может. Но они как то по разному мерят вторая с третьей.
А теперь нет. Видимо над тестом придётся ещё поработать.

У вас данные выровнены по линиям кеша? 4x4x4 — это как раз одна линия. Возможно, если данные пересекают линии, результат будет отличаться.

Sign up to leave a comment.

Articles