Самый быстрый способ сделать горизонтальную векторную сумму с плавающей точкой на x86
У вас есть вектор из трех (или четырех) поплавков. Какой самый быстрый способ их суммировать?
SSE (movaps, shuffle, add, movd) всегда быстрее, чем x87? Стоят ли инструкции горизонтального добавления в SSE4.2? Сколько стоит перейти на FPU, затем на faddp, faddp? Какая самая быстрая конкретная последовательность команд?
"Постарайтесь упорядочить вещи так, чтобы вы могли суммировать четыре вектора за раз", не будут приняты в качестве ответа:-)
5 ответов
Вот некоторые версии, настроенные на основе руководства по микроархам руководства Agner Fog и таблиц инструкций. Смотрите также x86 tag wiki. Они должны быть эффективными на любом процессоре, без каких-либо серьезных узких мест. (например, я избегал вещей, которые могли бы немного помочь одному уарху, но были бы медленными на другом уарше). Размер кода также минимизирован.
Общее 2x hadd идиома хороша только для размера кода, а не для скорости на любых существующих процессорах. Есть варианты использования для этого (см. Ниже), но это не один из них.
Я также включил версию AVX. Любое горизонтальное уменьшение с AVX / AVX2 должно начинаться с vextractf128 и "вертикальная" операция для уменьшения до одного XMM (__m128) вектор.
Смотрите вывод asm из всего этого кода в Godbolt Compiler Explorer. Смотрите также мои улучшения в Библиотеке векторного класса C++ Agner Fog horizontal_add функции. ( ветка доски сообщений и код на github). Я использовал макросы CPP для выбора оптимальных тасов для размера кода для SSE2, SSE4 и AVX, а также для избежания movdqa когда AVX недоступен.
Есть компромиссы для рассмотрения:
- Размер кода: чем меньше, тем лучше по причинам I-кэша L1 и для выборки кода с диска (меньшие двоичные файлы). Общий размер двоичного файла в основном имеет значение для решений компилятора, принимаемых неоднократно по всей программе. Если вы потрудитесь написать что-то вручную с помощью встроенных функций, стоит потратить несколько байтов кода, если это даст какое-либо ускорение для всей программы (будьте осторожны с микробенчмарками, которые делают развертывание хорошо выглядящим).
- Размер uop-кеша: часто более ценный, чем L1 I$. 4 однопользовательских инструкции могут занимать меньше места, чем 2
haddpsтак что это очень актуально здесь. - задержка: иногда актуально
- пропускная способность: обычно не имеет значения, горизонтальные суммы не должны быть в самой внутренней петле.
- total uops fused-domain: Если окружающий код не является узким местом на том же порту, который использует hsum, это прокси для влияния hsum на пропускную способность всего этого.
Когда горизонтальное добавление нечасто:
Процессоры без Uop-кеша могли бы отдать предпочтение 2x haddps: Это медленно, когда он работает, но это не часто. Быть только 2 инструкциями минимизирует влияние на окружающий код (размер I $).
Процессоры с кешем uop, вероятно, предпочтут что-то, что займет меньше мопов, даже если это больше инструкций / больше размер кода x86. Общее количество используемых строк кэша мопов - это то, что мы хотим минимизировать, что не так просто, как минимизация общего количества мопов (взятые ветви и границы 32B всегда начинают новую строку кэша мопов).
В любом случае, с учетом вышесказанного, горизонтальные суммы встречаются очень часто, поэтому я попытаюсь тщательно составить несколько версий, которые хорошо компилируются. Не тестируется на каком-либо реальном оборудовании или даже тщательно проверяется. Там могут быть ошибки в случайных константах или что-то в этом роде.
Если вы делаете запасную / базовую версию своего кода, помните, что запускать его будут только старые процессоры; более новые процессоры будут работать с вашей версией AVX или SSE4.1 или чем-то еще.
Старые процессоры, такие как K8 и Core2(merom) и более ранние, имеют только 64-битные тасовые устройства. Core2 имеет 128-битные исполнительные блоки для большинства команд, но не для случайных. (Pentium M и K8 обрабатывают все 128-битные векторные инструкции как две 64-битные половины).
Тасует как movhlps эти данные также перемещаются в 64-битных блоках (без перестановок в 64-битных половинах).
На старых процессорах с медленным перемешиванием:
movhlps(Merom: 1uop) значительно быстрее, чемshufps(Мером:3 мопс). На Пентиуме-М дешевле чемmovaps, Кроме того, он работает в домене FP на Core2, избегая задержек обхода из-за других перемешиваний.unpcklpdбыстрее чемunpcklps,pshufdмедленный,pshuflw/pshufhwбыстрые (потому что они тасуют только 64-битную половину)pshufb mm0(MMX) быстро,pshufb xmm0медленный.haddpsочень медленный (6 моп на Merom и Pentium M)movshdup(Merom: 1uop) интересно: это единственный insop 1uop, который тасует в элементах 64b.
shufps на Core2 (включая Penryn) приносит данные в целочисленную область, вызывая задержку обхода, чтобы вернуть их к исполнительным блокам FP для addps, но movhlps полностью в области FP. shufpd также работает в домене float.
movshdup работает в целочисленной области, но только один моп.
AMD K10, Intel Core2 (Penryn / Wolfdale) и все последующие процессоры запускают все xmm-тасовки как один моп. (Но обратите внимание на задержку обхода с shufps на Пенрин, избегать с movhlps)
Без AVX, избегая впустую movaps / movdqa Инструкция требует тщательного выбора шаффлов. Только несколько перемешиваний работают как копирование и перемешивание, а не изменяют место назначения. Тасования, которые объединяют данные из двух входов (например, unpck* или же movhlps) может использоваться с переменной tmp, которая больше не нужна, вместо _mm_movehl_ps(same,same),
Некоторые из них могут быть сделаны быстрее (за исключением MOVAPS), но более уродливыми / менее "чистыми", если взять фиктивный аргумент для использования в качестве пункта назначения для первоначального перемешивания. Например:
// Use dummy = a recently-dead variable that vec depends on,
// so it doesn't introduce a false dependency,
// and the compiler probably still has it in a register
__m128d highhalf_pd(__m128d dummy, __m128d vec) {
#ifdef __AVX__
// With 3-operand AVX instructions, don't create an extra dependency on something we don't need anymore.
(void)dummy;
return _mm_unpackhi_pd(vec, vec);
#else
// Without AVX, we can save a MOVAPS with MOVHLPS into a dead register
__m128 tmp = _mm_castpd_ps(dummy);
__m128d high = _mm_castps_pd(_mm_movehl_ps(tmp, _mm_castpd_ps(vec)));
return high;
#endif
}
SSE1 (он же SSE):
float hsum_ps_sse1(__m128 v) { // v = [ D C | B A ]
__m128 shuf = _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 3, 0, 1)); // [ C D | A B ]
__m128 sums = _mm_add_ps(v, shuf); // sums = [ D+C C+D | B+A A+B ]
shuf = _mm_movehl_ps(shuf, sums); // [ C D | D+C C+D ] // let the compiler avoid a mov by reusing shuf
sums = _mm_add_ss(sums, shuf);
return _mm_cvtss_f32(sums);
}
# gcc 5.3 -O3: looks optimal
movaps xmm1, xmm0 # I think one movaps is unavoidable, unless we have a 2nd register with known-safe floats in the upper 2 elements
shufps xmm1, xmm0, 177
addps xmm0, xmm1
movhlps xmm1, xmm0 # note the reuse of shuf, avoiding a movaps
addss xmm0, xmm1
# clang 3.7.1 -O3:
movaps xmm1, xmm0
shufps xmm1, xmm1, 177
addps xmm1, xmm0
movaps xmm0, xmm1
shufpd xmm0, xmm0, 1
addss xmm0, xmm1
Я сообщил об ошибке лягушатника о пессимизации шаффлов. Он имеет свое собственное внутреннее представление для тасования и превращает его обратно в тасования. GCC чаще использует инструкции, которые непосредственно соответствуют встроенному вами.
Часто clang работает лучше, чем gcc, в коде, где выбор инструкций не настраивается вручную, или постоянное распространение может упростить вещи, даже если внутренние значения оптимальны для непостоянного случая. В целом, хорошо, что компиляторы работают как полноценный компилятор для встроенных функций, а не просто как ассемблер. Компиляторы часто могут генерировать хороший ассм из скалярного C, который даже не пытается работать так, как это делал бы хороший ассемблер. В конечном итоге компиляторы будут воспринимать встроенные функции как еще один оператор Си как входные данные для оптимизатора.
SSE3
float hsum_ps_sse3(__m128 v) {
__m128 shuf = _mm_movehdup_ps(v); // broadcast elements 3,1 to 2,0
__m128 sums = _mm_add_ps(v, shuf);
shuf = _mm_movehl_ps(shuf, sums); // high half -> low half
sums = _mm_add_ss(sums, shuf);
return _mm_cvtss_f32(sums);
}
# gcc 5.3 -O3: perfectly optimal code
movshdup xmm1, xmm0
addps xmm0, xmm1
movhlps xmm1, xmm0
addss xmm0, xmm1
Это имеет несколько преимуществ:
не требует никаких
movapsкопии для работы с деструктивными шаффлами (без AVX):movshdup xmm1, xmm2's назначение только для записи, поэтому он создаетtmpиз мертвого регистра для нас. Это также, почему я использовалmovehl_ps(tmp, sums)вместоmovehl_ps(sums, sums),маленький размер кода. Инструкции перетасовки малы:
movhlps3 байта,movshdupсоставляет 4 байта (так же, какshufps). Не требуется немедленный байт, поэтому с AVXvshufps5 байт ноvmovhlpsа такжеvmovshdupоба 4.
Я мог бы сохранить еще один байт с addps вместо addss, Поскольку это не будет использоваться во внутренних контурах, дополнительная энергия для переключения дополнительных транзисторов, вероятно, незначительна. Исключения FP из верхних 3 элементов не являются риском, потому что все элементы содержат действительные данные FP. Тем не менее, clang/LLVM на самом деле "понимает" тасования векторов и генерирует лучший код, если знает, что важен только младший элемент.
Как и в версии SSE1, добавление нечетных элементов к самим себе может вызвать исключения FP (например, переполнение), которые не произошли бы в противном случае, но это не должно быть проблемой. Денормалы медленные, но IIRC, дающий результат +Inf, не на большинстве уарчей.
SSE3 оптимизирует под размер кода
Если размер кода является вашей главной заботой, два haddps (_mm_hadd_ps) инструкции сделают свое дело (ответ Пола Р.). Это также самый простой для ввода и запоминания. Это не быстро, хотя. Даже Intel Skylake по-прежнему декодирует каждый haddps до 3 моп, с 6 циклами. Таким образом, несмотря на то, что он сохраняет байты машинного кода (I-кэш L1), он занимает больше места в более ценном uop-кэше. Реальные варианты использования для haddps: проблема транспонирования и суммирования или выполнение некоторого масштабирования на промежуточном этапе в этом SSE atoi() реализация.
AVX:
Эта версия сохраняет байт кода против ответа Марата на вопрос AVX.
#ifdef __AVX__
float hsum256_ps_avx(__m256 v) {
__m128 vlow = _mm256_castps256_ps128(v);
__m128 vhigh = _mm256_extractf128_ps(v, 1); // high 128
vlow = _mm_add_ps(vlow, vhigh); // add the low 128
return hsum_ps_sse3(vlow); // and inline the sse3 version, which is optimal for AVX
// (no wasted instructions, and all of them are the 4B minimum)
}
#endif
vmovaps xmm1,xmm0 # huh, what the heck gcc? Just extract to xmm1
vextractf128 xmm0,ymm0,0x1
vaddps xmm0,xmm1,xmm0
vmovshdup xmm1,xmm0
vaddps xmm0,xmm1,xmm0
vmovhlps xmm1,xmm1,xmm0
vaddss xmm0,xmm0,xmm1
vzeroupper
ret
Двойная точность:
double hsum_pd_sse2(__m128d vd) { // v = [ B | A ]
__m128 undef = _mm_undefined_ps(); // don't worry, we only use addSD, never touching the garbage bits with an FP add
__m128 shuftmp= _mm_movehl_ps(undef, _mm_castpd_ps(vd)); // there is no movhlpd
__m128d shuf = _mm_castps_pd(shuftmp);
return _mm_cvtsd_f64(_mm_add_sd(vd, shuf));
}
# gcc 5.3.0 -O3
pxor xmm1, xmm1 # hopefully when inlined, gcc could pick a register it knew wouldn't cause a false dep problem, and avoid the zeroing
movhlps xmm1, xmm0
addsd xmm0, xmm1
# clang 3.7.1 -O3 again doesn't use movhlps:
xorpd xmm2, xmm2 # with #define _mm_undefined_ps _mm_setzero_ps
movapd xmm1, xmm0
unpckhpd xmm1, xmm2
addsd xmm1, xmm0
movapd xmm0, xmm1 # another clang bug: wrong choice of operand order
// This doesn't compile the way it's written
double hsum_pd_scalar_sse2(__m128d vd) {
double tmp;
_mm_storeh_pd(&tmp, vd); // store the high half
double lo = _mm_cvtsd_f64(vd); // cast the low half
return lo+tmp;
}
# gcc 5.3 -O3
haddpd xmm0, xmm0 # Lower latency but less throughput than storing to memory
# ICC13
movhpd QWORD PTR [-8+rsp], xmm0 # only needs the store port, not the shuffle unit
addsd xmm0, QWORD PTR [-8+rsp]
Хранение в памяти и обратно позволяет избежать ALU UOP. Это хорошо, если давление порта перетасовки или ALU Uops в целом являются узким местом. (Обратите внимание, что это не нужно sub rsp, 8 или что-то еще, потому что x86-64 SysV ABI предоставляет красную зону, на которую обработчики сигналов не будут наступать.)
Некоторые люди хранят в массиве и суммируют все элементы, но компиляторы обычно не понимают, что младший элемент массива все еще находится в регистре перед хранением.
Integer:
pshufd это удобное копирование и перемешивание. Сдвиги битов и байтов, к сожалению, на месте, и punpckhqdq ставит верхнюю половину цели в нижнюю половину результата, противоположную movhlps может извлечь верхнюю половину в другой регистр.
С помощью movhlps для первого шага может быть хорошо на некоторых процессорах, но только если у нас есть чистая рег. pshufd это безопасный выбор, и быстро на все после Merom.
int hsum_epi32_sse2(__m128i x) {
#ifdef __AVX__
__m128i hi64 = _mm_unpackhi_epi64(x, x); // 3-operand non-destructive AVX lets us save a byte without needing a mov
#else
__m128i hi64 = _mm_shuffle_epi32(x, _MM_SHUFFLE(1, 0, 3, 2));
#endif
__m128i sum64 = _mm_add_epi32(hi64, x);
__m128i hi32 = _mm_shufflelo_epi16(sum64, _MM_SHUFFLE(1, 0, 3, 2)); // Swap the low two elements
__m128i sum32 = _mm_add_epi32(sum64, hi32);
return _mm_cvtsi128_si32(sum32); // SSE2 movd
//return _mm_extract_epi32(hl, 0); // SSE4, even though it compiles to movd instead of a literal pextrd r32,xmm,0
}
# gcc 5.3 -O3
pshufd xmm1,xmm0,0x4e
paddd xmm0,xmm1
pshuflw xmm1,xmm0,0x4e
paddd xmm0,xmm1
movd eax,xmm0
int hsum_epi32_ssse3_slow_smallcode(__m128i x){
x = _mm_hadd_epi32(x, x);
x = _mm_hadd_epi32(x, x);
return _mm_cvtsi128_si32(x);
}
На некоторых процессорах безопасно использовать FP-тасовки для целочисленных данных. Я этого не делал, поскольку на современных процессорах, которые максимально сохраняют 1 или 2 байта кода, без увеличения скорости (кроме размера кода / эффектов выравнивания).
SSE2
Все четыре:
const __m128 t = _mm_add_ps(v, _mm_movehl_ps(v, v));
const __m128 sum = _mm_add_ss(t, _mm_shuffle_ps(t, t, 1));
r1 + r2 + r3:
const __m128 t1 = _mm_movehl_ps(v, v);
const __m128 t2 = _mm_add_ps(v, t1);
const __m128 sum = _mm_add_ss(t1, _mm_shuffle_ps(t2, t2, 1));
Я обнаружил, что скорость примерно равна двойной HADDPS (но я не измерял слишком близко).
Вы можете сделать это в два HADDPS инструкции в SSE3:
v = _mm_hadd_ps(v, v);
v = _mm_hadd_ps(v, v);
Это помещает сумму во все элементы.
Я бы определенно попробовал SSE 4.2. Если вы делаете это несколько раз (я полагаю, что это так, если производительность является проблемой), вы можете предварительно загрузить регистр с помощью (1,1,1,1), а затем сделать несколько точек 4(my_vec(s), one_vec) в теме. Да, это делает избыточное умножение, но в наши дни это довольно дешево, и в такой операции, вероятно, будут преобладать горизонтальные зависимости, которые могут быть более оптимизированы в новой функции точечного продукта SSE. Вы должны проверить, превосходит ли это двойное горизонтальное добавление Пола Р.
Я также предлагаю сравнить его с прямым скалярным (или скалярным SSE) кодом - как ни странно, он часто быстрее (обычно потому, что внутренне он сериализован, но плотно конвейеризован с использованием обхода регистра, где специальные горизонтальные инструкции могут быть не быстрыми (пока)), если только вы выполняются SIMT-подобный код, который звучит так, как будто вы не являетесь (в противном случае вы бы сделали четыре продукта с точками).
Часто вопрос о самом быстром пути предполагает задачу, которую необходимо выполнить несколько раз в критическом по времени цикле.
Тогда возможно, что самым быстрым методом может быть итерационный метод, работающий попарно, который амортизирует часть работы между итерациями.
Общая стоимость сокращения путем разделения вектора на младшие/высокие части составляет O(log2(N)), в то время как амортизированная стоимость разделения вектора на четные/нечетные последовательности составляет O(1).
inline vec update(vec context, vec data) {
vec even = get_evens(context, data);
vec odd = get_odds(context, data);
return vertical_operation(even, odd);
}
void my_algo(vec *data, int N, vec_element_type *out) {
vec4 context{0,0,0,0};
context = update(context, data[0]);
int i;
for (int i = 0; i < N-1; i++) {
context = update(context, data[i+1]);
output[i] = extract_lane(context, 1);
}
context = update(context, anything);
output[N-1] = extract_lane(context, 1);
}
Искомая сумма будет найдена из второго элемента (индекс 1) аккумулятора (после 1 итерации), в то время как первый элемент будет содержать общее сокращение всех элементов на данный момент.
Reduct = [ -- ][ -- ][ -- ][ -- ]
New input = [i0 ][ i1 ][ i2 ][ i3 ]
evens = [ -- ][ -- ][ i0 ][ i2 ]
odds = [ -- ][ -- ][ i1 ][ i3 ]
------- vertical arithmetic reduction ----
Reduct = [ -- ][ -- ][ 01 ][ 23 ]
input = [ 4 ][ 5 ][ 6 ][ 7 ]
evens = [ -- ][ 01 ][ 4 ][ 6 ]
odds = [ -- ][ 23 ][ 5 ][ 7 ]
Reduct = [ -- ][ 0123 ][ 45 ][ 67 ]
New input: [ 8 ] [ 9 ] [ a ] [ b ]
evens = [ -- ][ 45 ][ 8 ][ a ]
odds = [0123][ 67 ][ 9 ][ b ]
------------------------------
Reduct = [0123][4567][ 89 ][ ab ]
У меня есть сомнения, окажется ли это быстрее для длины вектора 3 или 4, чем представлено г-ном Кордесом, однако для 16 или 8-битных данных этот метод должен оказаться полезным. Затем, конечно, нужно выполнить 3 или 4 раунда соответственно, прежде чем можно будет получить результат.
Если горизонтальная операция является суммой, то на самом деле можно использовать только одну операцию.