Быстрая векторизация rsqrt и обратная с SSE/AVX в зависимости от точности

Question

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

__m128 recip_float4_ieee(__m128 x) { return _mm_div_ps(_mm_set1_ps(1.0f), x); }
__m128 rsqrt_float4_ieee(__m128 x) { return _mm_div_ps(_mm_set1_ps(1.0f), _mm_sqrt_ps(x)); }

Это работает отлично, но медленно: согласно руководству, они занимают 14 и 28 циклов на Sandy Bridge (пропускная способность). Соответствующие версии AVX занимают почти столько же времени на Haswell.

С другой стороны, вместо этого можно использовать следующие версии:

__m128 recip_float4_half(__m128 x) { return _mm_rcp_ps(x); }
__m128 rsqrt_float4_half(__m128 x) { return _mm_rsqrt_ps(x); }

Они занимают всего один или два цикла времени (пропускная способность), обеспечивая значительное повышение производительности. Однако они ОЧЕНЬ приблизительны: они дают результат с относительной погрешностью менее 1,5 * 2^-12. Учитывая, что машинный эпсилон чисел с плавающей точкой одинарной точности равен 2 ^ -24, мы можем сказать, что это приближение имеет приблизительно половину точности.

Кажется, что итерация Ньютона-Рафсона может быть добавлена ​​для получения результата с одинарной точностью (возможно, не настолько точной, как того требует стандарт IEEE), см. GCC, ICC, обсуждения в LLVM. Теоретически, тот же метод может использоваться для значений двойной точности, получая половину или одинарную или двойную точность.

Мне интересны рабочие реализации этого подхода для типов данных с плавающей запятой и двойных значений, а также для всех (половинная, одинарная, двойная) точность. Обработка особых случаев (деление на ноль, sqrt(-1), inf/nan и т. П.) Не требуется. Кроме того, мне не ясно, какая из этих подпрограмм будет быстрее, чем тривиальные IEEE-совместимые решения, а какие будет медленнее.

Вот несколько небольших ограничений на ответы, пожалуйста:

1. Используйте встроенные функции в своих примерах кода. Сборка зависит от компилятора, поэтому менее полезна.
2. Используйте аналогичное соглашение об именах для функций.
3. Реализуйте подпрограммы, используя один регистр SSE/AVX, содержащий плотно упакованные значения типа float/double в качестве входных данных. Если есть значительное повышение производительности, вы также можете опубликовать подпрограммы, принимающие в качестве входных данных несколько регистров (два регистра могут быть жизнеспособными).
4. Не публикуйте обе версии SSE/AVX, если они абсолютно равны до изменения _mm на _mm256 и наоборот.

Любые оценки производительности, измерения, обсуждения приветствуются.

РЕЗЮМЕ

Вот версии для чисел с плавающей запятой одинарной точности с одной итерацией NR:

__m128 recip_float4_single(__m128 x) {
  __m128 res = _mm_rcp_ps(x);
  __m128 muls = _mm_mul_ps(x, _mm_mul_ps(res, res));
  return res =  _mm_sub_ps(_mm_add_ps(res, res), muls);
}
__m128 rsqrt_float4_single(__m128 x) {
  __m128 three = _mm_set1_ps(3.0f), half = _mm_set1_ps(0.5f);
  __m128 res = _mm_rsqrt_ps(x); 
  __m128 muls = _mm_mul_ps(_mm_mul_ps(x, res), res); 
  return res = _mm_mul_ps(_mm_mul_ps(half, res), _mm_sub_ps(three, muls)); 
}

Ответ, данный Питером Кордесом, объясняет, как создавать другие версии, и содержит тщательный теоретический анализ производительности.

Все реализованные решения с бенчмарком вы можете найти здесь: rece_rsqrt_benchmark.

Полученные результаты по пропускной способности на Ivy Bridge представлены ниже. Были протестированы только реализации SSE с одним регистром. Затраченное время дается в циклах за звонок. Первое число для половинной точности (без NR), второе для одинарной точности (1 NR итерация), третье для 2 NR итераций.

1. Возврат на плавание занимает 1, 4 цикла против 7 циклов.
2. rsqrt on float занимает 1, 6 цикла против 14 циклов.
3. Возврат на удвоение занимает 3, 6, 9 циклов против 14 циклов.
4. rsqrt on double занимает 3, 8, 13 циклов против 28 циклов.

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

Answer 1

Есть много примеров алгоритма на практике. Например:

Ньютон Рафсон с SSE2 - может кто-нибудь объяснить мне, что у этих трех строк есть ответ, объясняющий итерацию, использованную одним из примеров Intel.

Для анализа производительности, скажем, Haswell (который может выполнять FP mul на двух портах исполнения, в отличие от предыдущих проектов), я воспроизведу код здесь (по одному оператору на строку). См. http://agner.org/optimize/ для таблиц пропускной способности инструкций и задержек, а также для документации о том, как понять больше предыстории.

// execute (aka dispatch) on cycle 1, results ready on cycle 6
nr = _mm_rsqrt_ps( x );

// both of these execute on cycle 6, results ready on cycle 11
xnr = _mm_mul_ps( x, nr );         // dep on nr
half_nr = _mm_mul_ps( half, nr );  // dep on nr

// can execute on cycle 11, result ready on cycle 16
muls = _mm_mul_ps( xnr , nr );     // dep on xnr

// can execute on cycle 16, result ready on cycle 19
three_minus_muls = _mm_sub_ps( three, muls );  // dep on muls

// can execute on cycle 19, result ready on cycle 24
result = _mm_mul_ps( half_nr, three_minus_muls ); // dep on three_minus_muls
// result is an approximation of 1/sqrt(x), with ~22 to 23 bits of precision in the mantissa.

Здесь очень много места для других вычислений, которые могут перекрываться, если они не являются частью цепочки зависимостей. Однако, если данные для каждой итерации вашего кода зависят от данных из предыдущего, вам может быть лучше с задержкой в ​​11 циклов: sqrtps, Или даже если каждая итерация цикла достаточно длинная, чтобы неупорядоченное выполнение не могло скрыть все это путем наложения независимых итераций.

Получить sqrt(x) вместо 1/sqrt(x), умножить на x (и исправить, если x может быть ноль, например, путем маскировки (_mm_andn_ps) с результатом CMPPS против 0,0). Оптимальным способом является замена half_nr с half_xnr = _mm_mul_ps( half, xnr );, Это не удлиняет цепь dep, потому что half_xnr может начинаться с цикла 11, но не требуется до конца (цикл 19). То же самое с доступным FMA: без увеличения общей задержки.

Если достаточно 11 бит точности (без итераций Ньютона), руководство по оптимизации Intel (раздел 11.12.3) предлагает использовать rcpps(rsqrt(x)), что хуже, чем умножение на исходный x, по крайней мере с AVX. Возможно, он сохраняет инструкцию movdqa с 128-битным SSE, но 256b rcpps медленнее, чем 256b mul или fma. (И это позволяет вам добавить результат sqrt к чему-то бесплатно с FMA вместо mul для последнего шага).

---

Версия этого цикла для SSE, без учета каких-либо инструкций по перемещению, составляет 6 моп. Это означает, что пропускная способность Haswell должна составлять один на 3 цикла (учитывая, что два исполнительных порта могут обрабатывать FP mul, а rsqrt находится на противоположном порту от add / sub FP). На SnB/IvB (и, вероятно, Nehalem) его пропускная способность должна составлять один на 5 циклов, поскольку mulps и rsqrtps конкурируют за порт 0. subps находится на port1, что не является узким местом.

Для Haswell мы можем использовать FMA, чтобы объединить вычитание с мул. Однако модуль divrs / sqrt не имеет ширины 256b, поэтому, в отличие от всего остального, divps / sqrtps / rsqrtps / rcpps для регистров ymm требует дополнительных мопов и дополнительных циклов.

// vrsqrtps ymm has higher latency
// execute on cycle 1, results ready on cycle 8
nr = _mm256_rsqrt_ps( x );

// both of can execute on cycle 8, results ready on cycle 13
xnr = _mm256_mul_ps( x, nr );         // dep on nr
half_nr = _mm256_mul_ps( half, nr );  // dep on nr

// can execute on cycle 13, result ready on cycle 18
three_minus_muls = _mm256_fnmadd_ps( xnr, nr, three );  // -(xnr*nr) + 3

// can execute on cycle 18, result ready on cycle 23
result = _mm256_mul_ps( half_nr, three_minus_muls ); // dep on three_minus_muls

Мы экономим 3 цикла с FMA. Это компенсируется использованием 2-cyle-slower 256b rsqrt, для чистого усиления на 1 такт меньше задержки (довольно хорошо для вдвое большей ширины). SnB/IvB AVX будет иметь задержку 24c для 128b, задержку 26c для 256b.

Версия 256b FMA использует всего 7 моп. (VRSQRTPS равен 3 моп, 2 для p0 и один для p1/5.) 256b mulps и fma являются инструкциями с одним мопом, и оба могут работать на порте 0 или порте 1. (p0 только на pre-Haswell). Таким образом, пропускная способность должна составлять: один на 3c, если движок OOO отправляет мопы на порты оптимального исполнения. (т. е. случайное перемещение из rsqrt всегда переходит к p5, а не к p1, где он будет занимать полосу пропускания mul/fma.) Что касается перекрытия с другими вычислениями (а не только с независимым выполнением самого себя), опять же это довольно легко. Только 7 мопов с цепочкой деп из 23 циклов оставляют много места для других вещей, пока эти мопы находятся в буфере переупорядочения.

Если это шаг в гигантской цепочке депов, и больше ничего не происходит (даже независимой следующей итерации), то 256b vsqrtps задержка 19 циклов с пропускной способностью один на 14 циклов. (Хасуэллы). Если вам все еще действительно нужно взаимное, то 256b vdivps также имеет задержку 18-21c, с пропускной способностью один на 14 c. Так что для обычного sqrt, инструкция имеет меньшую задержку. Для кванта приема итерация приближения имеет меньшую задержку. (И FAR больше пропускной способности, если он перекрывается с самим собой. Если перекрывается с другими вещами, которые не делят единицу, sqrtps не проблема.)

sqrtps может быть выигрыш пропускной способности против rsqrt + Ньютоновская итерация, если она является частью тела цикла с достаточным количеством других делений, не связанных с делением и чтением, что единица деления не насыщена.

Это особенно верно, если вам нужно sqrt(x) не 1/sqrt(x), например, в Haswell с AVX2 цикл copy+arcsinh над массивом с плавающей точкой, который помещается в кэш L3, реализован как fastlog(v + sqrt(v*v + 1)), работает примерно с той же пропускной способностью с реальным VSQRTPS или с VRSQRTPS + итерацией Ньютона-Рафсона. (Даже с очень быстрой аппроксимацией для log (), поэтому общее тело цикла составляет около 9 операций FMA/add/mul/convert и 2 логических значения, плюс VSQRTPS ymm. Существует ускорение от использования только fastlog(v2_plus_1 * rsqrt(v2_plus_1) + v2_plus_1) таким образом, это не узкое место в пропускной способности памяти, но это может быть узким местом в задержке (так что неупорядоченное выполнение не может использовать весь параллелизм независимых итераций).

Другие точности

Для обеспечения полуточности нет никаких инструкций для выполнения вычислений на половинных числах. Вы должны конвертировать на лету при загрузке / хранении, используя инструкции по конвертации.

Для двойной точности нет rsqrtpd, Предположительно, мысль заключается в том, что если вам не нужна полная точность, вы должны использовать float. Таким образом, вы можете преобразовать в float и обратно, а затем сделать точно такой же алгоритм, но с pd вместо ps. встроенные функции Или вы можете хранить ваши данные как плавающие некоторое время. например, преобразовать два ymm-регистра двойников в один ymm-регистр синглов.

К сожалению, нет ни одной инструкции, которая бы принимала два миллиметровых регистра двойных и выводила миллиметровый сингл. Вы должны идти ymm->xmm дважды, тогда _mm256_insertf128_ps один хмм до высокого 128 другого. Но тогда вы можете передать этот вектор 256b ммм по тому же алгоритму.

Если вы собираетесь преобразовать обратно в удвоение сразу после этого, возможно, имеет смысл выполнить итерацию rsqrt + Ньютона-Рафсона для двух 128b регистров синглов по отдельности. Дополнительные 2 мопа для вставки / извлечения и 2 дополнительных мопа для 256b rsqrt начинают складываться, не говоря уже о 3-тактной задержке vinsertf128 / vextractf128, Вычисления будут перекрываться, и оба результата будут готовы через пару циклов. (Или 1 цикл, если у вас есть специальная версия вашего кода для чередования операций на 2 входах одновременно).

Помните, что одинарная точность имеет меньший диапазон показателей, чем двойная, поэтому преобразование может переполниться до +Inf или понизиться до 0,0. Если вы не уверены, просто используйте обычный _mm_sqrt_pd,

---

С AVX512F есть _mm512_rsqrt14_pd( __m512d a), С AVX512ER (KNL, но не SKX или Cannonlake), есть _mm512_rsqrt28_pd, _ps версии конечно тоже существуют. 14 бит точности мантиссы может быть достаточно для использования без итерации Ньютона в большинстве случаев.

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