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

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

__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, можно сказать, что это приближение имеет примерно половину точности.

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

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

Вот несколько незначительных ограничений на ответы:

Используйте встроенные примеры кода. Сборка зависит от компилятора, поэтому менее полезна.
Используйте аналогичные соглашения об именах для функций.
Реализовать подпрограммы с использованием одного регистра SSE/AVX, содержащего плотно упакованные значения float/double в качестве входных данных. Если есть значительный прирост производительности, вы также можете отправлять процедуры, принимающие несколько регистров в качестве входных данных (два регистра могут быть жизнеспособными).
Не отправляйте обе версии 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)); 
}

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

Здесь вы можете найти все реализованные решения: recip_rsqrt_benchmark.

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

recip on float занимает циклы 1, 4 против циклов 7.
rsqrt на float занимает циклы 1, 6 против 14.
recip on double принимает циклы 3, 6, 9 по сравнению с 14.
rsqrt на двойных циклах цикла 3, 8, 13 против циклов 28.

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

Ответ 1

На практике существует множество примеров алгоритма. Например:

Ньютон Рафсон с SSE2 - может кто-нибудь объяснить мне эти 3 строки имеет ответ, объясняющий итерацию, используемую один из примеров 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 (и fixup, если x может быть нулевым, например, путем маскировки (_mm_andn_ps) с результатом CMPPS против 0,0). Оптимальным способом является замена half_nr на half_xnr = _mm_mul_ps( half, xnr );. Это не удлиняет цепочку отрезков, потому что half_xnr может начинаться в цикле 11, но не требуется до конца (цикл 19). То же самое с доступным FMA: нет увеличения общей задержки.

Если достаточно 11 бит точности (без итерации Newton), Руководство по оптимизации 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 находится на противоположном порту от FP add/sub). На SnB/IvB (и, вероятно, Nehalem) он должен иметь пропускную способность один на 5 циклов, так как mulps и rsqrtps конкурируют за порт 0. subps находится на порту1, что не является узким местом.

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

// 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 будет латентностью 24 с для задержки 128b, 26c для 256b.

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

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

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

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

Другие замечания

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

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

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

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

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

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

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