Обработка нулей в _mm256_rsqrt_ps()
При условии _mm256_sqrt_ps() является относительно медленным, и что значения, которые я генерирую, немедленно усекаются с _mm256_floor_ps()Оглядываясь вокруг, кажется, что делаешь:
_mm256_mul_ps(_mm256_rsqrt_ps(eightFloats),
eightFloats);
Это способ добиться дополнительной производительности и избежать остановки конвейера.
К сожалению, с нулевыми значениями, я, конечно, получаю сбой расчета 1/sqrt(0), Как лучше обойти это? Я пробовал это (что работает и быстрее), но есть ли лучший способ, или я столкнусь с проблемами при определенных условиях?
_mm256_mul_ps(_mm256_rsqrt_ps(_mm256_max_ps(eightFloats,
_mm256_set1_ps(0.1))),
eightFloats);
Мой код предназначен для вертикального приложения, поэтому я могу предположить, что оно будет работать на процессоре Haswell (i7-4810MQ), поэтому можно использовать FMA/AVX2. Оригинальный код примерно:
float vals[MAX];
int sum = 0;
for (int i = 0; i < MAX; i++)
{
int thisSqrt = (int) floor(sqrt(vals[i]));
sum += min(thisSqrt, 0x3F);
}
Все значения vals должны быть целочисленными значениями. (Почему все не просто int это другой вопрос...)
1 ответ
tl; dr: смотрите конец кода, который компилируется и должен работать.
Чтобы просто решить 0.0 проблема, вы могли бы также особые случаи 0.0 с FP сравнить источник с 0.0, Используйте результат сравнения в качестве маски для обнуления любых NaN, полученных в результате 0 * +Infinity в sqrt(x) = x * rsqrt(x)). Clang делает это при автовекторизации. (Но он использует blendps с нулевым вектором, вместо использования маски сравнения с andnps непосредственно к нулю или сохранить элементы.)
Также было бы возможно использовать sqrt(x) ~= recip(rsqrt(x)), как предложено njuffa. rsqrt(0) = +Inf, recip(+Inf) = 0, Однако использование двух приближений может составить относительную ошибку, что является проблемой.
То, что вам не хватает:
Усечение до целого числа (вместо округления) требует точного sqrt результат, когда вход является идеальным квадратом. Если результат для 25*rsqrt(25) равен 4.999999 или около того (вместо 5.00001), вы добавите 4 вместо 5,
Даже с итерацией Ньютона-Рафсона, rsqrtps не совсем точный способ sqrtps есть, так что еще может дать 5.0 - 1ulp, (1ulp = одна единица на последнем месте = самый младший бит мантиссы).
Также:
- Формула ньютона рафсона
- Производительность реализации Newton Raphson SSE (задержка / пропускная способность). Обратите внимание, что мы заботимся больше о пропускной способности, чем о задержке, поскольку мы используем ее в цикле, который больше ничего не делает. sqrt не является частью переносимой по цепочке цепи dep, поэтому разные итерации могут иметь свои sqrt calcs в полете сразу.
Может быть возможно убить 2 зайцев одним выстрелом, добавив небольшую константу перед выполнением (x+offset)*approx_rsqrt(x+offset) а затем обрезать до целого числа. Достаточно большой, чтобы преодолеть максимальную относительную погрешность 1,5 * 2 -12, но достаточно мал, чтобы не сталкиваться sqrt_approx(63*63-1+offset) вплоть до 63 (самый чувствительный случай).
63*1.5*2^(-12) == 0.023071...
approx_sqrt(63*63-1) == 62.99206... +/- 0.023068..
На самом деле, мы облажались без итераций Ньютона, даже не добавляя ничего. approx_sqrt(63*63-1) мог выйти выше 63.0 сам по себе. n=36 это наибольшее значение, где относительная ошибка в sqrt(n*n-1) + error меньше чем sqrt(n*n), GNU Calc:
define f(n) { local x=sqrt(n*n-1); local e=x*1.5*2^(-12); print x; print e, x+e; }
; f(36)
35.98610843089316319413
~0.01317850650545403926 ~35.99928693739861723339
; f(37)
36.9864840178138587015
~0.01354485498699237990 ~37.00002887280085108140
Имеют ли ваши исходные данные какие-либо свойства, которые означают, что вам не нужно беспокоиться о том, что они находятся чуть ниже большого идеального квадрата? например, это всегда идеальные квадраты?
Вы можете проверить все возможные входные значения, так как важный домен очень мал (целочисленные значения FP от 0,63*63), чтобы увидеть, достаточно ли мала практическая ошибка на Intel Haswell, но это будет хрупкая оптимизация, которая может сломайте ваш код на процессорах AMD или даже на будущих процессорах Intel. К сожалению, простое кодирование в ISA-спецификации гарантирует, что относительная ошибка до 1,5 * 2 -12 требует большего количества инструкций. Я не вижу никаких хитростей итераций NR.
Если бы ваш верхний предел был меньше (например, 20), вы могли бы просто сделать isqrt = static_cast<int> ((x+0.5)*approx_rsqrt(x+0.5)), Вы получите 20 за 20*20, но всегда 19 за 20*20-1,
; define test_approx_sqrt(x, off) { local s=x*x+off; local sq=s/sqrt(s); local sq_1=(s-1)/sqrt(s-1); local e=1.5*2^(-12); print sq, sq_1; print sq*e, sq_1*e; }
; test_approx_sqrt(20, 0.5)
~20.01249609618950056874 ~19.98749609130668473087 # (x+0.5)/sqrt(x+0.5)
~0.00732879495710064718 ~0.00731963968187500662 # relative error
Обратите внимание, что val * (x +/- err) = val*x +/- val*err, IEEE FP mul производит результаты, которые правильно округлены до 0.5ulp, поэтому это должно работать для относительных ошибок FP.
В любом случае, я думаю, вам нужна одна итерация Ньютона-Рафсона.
Лучше всего добавить 0.5 на ваш вход, прежде чем делать прибл. rsqrt, Это обходит проблему 0/0 = NaN и сдвигает диапазон ошибки +/- к одной стороне от точки отсечения целого числа (для чисел в диапазоне, который нас интересует).
Инструкции FP min/max имеют ту же производительность, что и FP add, и в любом случае будут на критическом пути. Использование add вместо max также решает проблему результатов для идеальных квадратов, потенциально находящихся на несколько ступеней ниже правильного результата.
Выход компилятора: достойная отправная точка
Я получаю довольно хорошие результаты автовекторизации из Clang 3.7.1 с sum_int, с -fno-math-errno -funsafe-math-optimizations, -ffinite-math-only не требуется (но даже при полной -ffast-math, лязг избегает sqrt(0) = NaN когда используешь rsqrtps).
sum_fp не авто-векторизация, даже с полным -ffast-math,
тем не мение clang Версия имеет ту же проблему, что и ваша идея: обрезать неточный результат из rsqrt + NR, потенциально давая неправильное целое число. IDK, если по этой причине gcc не выполняет авто-векторизацию, потому что он мог бы использовать sqrtps для большого ускорения без изменения результатов. (По крайней мере, пока все значения с плавающей точкой находятся между 0 и INT_MAX 2, в противном случае обратное преобразование в целое число даст "неопределенный" результат INT_MIN. (Установленный бит установлен, все остальные биты очищены). Это тот случай, когда -ffast-math ломает вашу программу, если вы не используете -mrecip=none или что-то.
Смотрите вывод asm на Godbolt от:
// autovectorizes with clang, but has rounding problems.
// Note the use of sqrtf, and that floorf before truncating to int is redundant. (removed because clang doesn't optimize away the roundps)
int sum_int(float vals[]){
int sum = 0;
for (int i = 0; i < MAX; i++) {
int thisSqrt = (int) sqrtf(vals[i]);
sum += std::min(thisSqrt, 0x3F);
}
return sum;
}
Чтобы вручную векторизовать с помощью встроенных функций, мы можем посмотреть на вывод asm из -fno-unroll-loops (чтобы все было просто). Я собирался включить это в ответ, но потом понял, что у него есть проблемы.
положить это вместе:
Я думаю, что преобразование в int внутри цикла лучше, чем использование floorf а потом addps, roundps является инструкцией по 2 мопа (задержка 6 с) в Haswell (1 моп в SnB/IvB). Хуже того, оба мопа требуют port1, поэтому они конкурируют с FP add / mul. cvttps2dq это инструкция 1-моп для port1 с задержкой 3c, а затем мы можем использовать целое число min и добавить для ограничения и накопления, так что port5 получает что-то делать. Использование целочисленного векторного аккумулятора также означает, что цепочка зависимостей, переносимых циклами, составляет 1 цикл, поэтому нам не нужно развертывать или использовать несколько аккумуляторов, чтобы сохранить несколько итераций в полете. Меньший код всегда лучше для общей картины (кэш Uop, I-кэш L1, предикторы ветвления).
Пока нам не грозит переполнение 32-битных аккумуляторов, это, кажется, лучший выбор. (Не проверяя ничего и даже не проверяя это).
Я не пользуюсь sqrt(x) ~= approx_recip(approx_sqrt(x)) метод, потому что я не знаю, как сделать итерацию Ньютона, чтобы уточнить его (вероятно, это будет включать в себя деление). И потому что сложная ошибка больше.
Горизонтальная сумма из этого ответа.
Полная, но не проверенная версия:
#include <immintrin.h>
#define MAX 4096
// 2*sqrt(x) ~= 2*x*approx_rsqrt(x), with a Newton-Raphson iteration
// dividing by 2 is faster in the integer domain, so we don't do it
__m256 approx_2sqrt_ps256(__m256 x) {
// clang / gcc usually use -3.0 and -0.5. We could do the same by using fnmsub_ps (add 3 = subtract -3), so we can share constants
__m256 three = _mm256_set1_ps(3.0f);
//__m256 half = _mm256_set1_ps(0.5f); // we omit the *0.5 step
__m256 nr = _mm256_rsqrt_ps( x ); // initial approximation for Newton-Raphson
// 1/sqrt(x) ~= nr * (3 - x*nr * nr) * 0.5 = nr*(1.5 - x*0.5*nr*nr)
// sqrt(x) = x/sqrt(x) ~= (x*nr) * (3 - x*nr * nr) * 0.5
// 2*sqrt(x) ~= (x*nr) * (3 - x*nr * nr)
__m256 xnr = _mm256_mul_ps( x, nr );
__m256 three_minus_muls = _mm256_fnmadd_ps( xnr, nr, three ); // -(xnr*nr) + 3
return _mm256_mul_ps( xnr, three_minus_muls );
}
// packed int32_t: correct results for inputs from 0 to well above 63*63
__m256i isqrt256_ps(__m256 x) {
__m256 offset = _mm256_set1_ps(0.5f); // or subtract -0.5, to maybe share constants with compiler-generated Newton iterations.
__m256 xoff = _mm256_add_ps(x, offset); // avoids 0*Inf = NaN, and rounding error before truncation
__m256 approx_2sqrt_xoff = approx_2sqrt_ps256(xoff);
__m256i i2sqrtx = _mm256_cvttps_epi32(approx_2sqrt_xoff);
return _mm256_srli_epi32(i2sqrtx, 1); // divide by 2 with truncation
// alternatively, we could mask the low bit to zero and divide by two outside the loop, but that has no advantage unless port0 turns out to be the bottleneck
}
__m256i isqrt256_ps_simple_exact(__m256 x) {
__m256 sqrt_x = _mm256_sqrt_ps(x);
__m256i isqrtx = _mm256_cvttps_epi32(sqrt_x);
return isqrtx;
}
int hsum_epi32_avx(__m256i x256){
__m128i xhi = _mm256_extracti128_si256(x256, 1);
__m128i xlo = _mm256_castsi256_si128(x256);
__m128i x = _mm_add_epi32(xlo, xhi);
__m128i hl = _mm_shuffle_epi32(x, _MM_SHUFFLE(1, 0, 3, 2));
hl = _mm_add_epi32(hl, x);
x = _mm_shuffle_epi32(hl, _MM_SHUFFLE(2, 3, 0, 1));
hl = _mm_add_epi32(hl, x);
return _mm_cvtsi128_si32(hl);
}
int sum_int_avx(float vals[]){
__m256i sum = _mm256_setzero_si256();
__m256i upperlimit = _mm256_set1_epi32(0x3F);
for (int i = 0; i < MAX; i+=8) {
__m256 v = _mm256_loadu_ps(vals+i);
__m256i visqrt = isqrt256_ps(v);
// assert visqrt == isqrt256_ps_simple_exact(v) or something
visqrt = _mm256_min_epi32(visqrt, upperlimit);
sum = _mm256_add_epi32(sum, visqrt);
}
return hsum_epi32_avx(sum);
}
Компилирует Godbolt в хороший код, но я его не проверял. clang делает немного более приятный код, чем gcc: clang использует широковещательную загрузку из 4B-местоположений для констант set1 вместо того, чтобы повторять их во время компиляции в 32B-константы. У gcc также есть причудливая movdqa для копирования регистра.
В любом случае, весь цикл получается только из 9 векторных инструкций, по сравнению с 12 для сгенерированного компилятором sum_int версия. Это, вероятно, не заметил x*initial_guess(x) общие подвыражения, которые встречаются в итерационной формуле Ньютона-Рафсона, когда вы умножаете результат на x, или что-то типа того. Он также делает дополнительные MULPS вместо psrld, потому что он делает *0.5 перед преобразованием в int. Так вот откуда берутся две дополнительные инструкции mulps, и есть cmpps / blendvps.
sum_int_avx(float*):
vpxor ymm3, ymm3, ymm3
xor eax, eax
vbroadcastss ymm0, dword ptr [rip + .LCPI4_0] ; set1(0.5)
vbroadcastss ymm1, dword ptr [rip + .LCPI4_1] ; set1(3.0)
vpbroadcastd ymm2, dword ptr [rip + .LCPI4_2] ; set1(63)
LBB4_1: ; latencies
vaddps ymm4, ymm0, ymmword ptr [rdi + 4*rax] ; 3c
vrsqrtps ymm5, ymm4 ; 7c
vmulps ymm4, ymm4, ymm5 ; x*nr ; 5c
vfnmadd213ps ymm5, ymm4, ymm1 ; 5c
vmulps ymm4, ymm4, ymm5 ; 5c
vcvttps2dq ymm4, ymm4 ; 3c
vpsrld ymm4, ymm4, 1 ; 1c this would be a mulps (but not on the critical path) if we did this in the FP domain
vpminsd ymm4, ymm4, ymm2 ; 1c
vpaddd ymm3, ymm4, ymm3 ; 1c
; ... (those 9 insns repeated: loop unrolling)
add rax, 16
cmp rax, 4096
jl .LBB4_1
;... horizontal sum
IACA считает, что без развертывания Haswell может поддерживать пропускную способность в одну итерацию на 4,15 цикла, что является узким местом на портах 0 и 1. Таким образом, потенциально вы можете сократить цикл, накапливая sqrt(x)*2 (с усечением до четных чисел с использованием _mm256_and_si256), и делим только на два вне цикла.
Также согласно IACA, задержка одной итерации составляет 38 циклов в Haswell. Я получаю только 31с, так что, вероятно, он включает задержку загрузки L1 или что-то в этом роде. В любом случае это означает, что для насыщения исполнительных блоков операции из 8 итераций должны выполняться сразу. Это 8 * ~14 unops-uops = 112 unsused-uops (или меньше с развертыванием clang), которые должны быть в полете сразу. Планировщик Haswell на самом деле только 60 записей, но ROB составляет 192 записи. Ранние мопы из ранних итераций уже будут выполнены, поэтому их нужно отслеживать только в ROB, а не в планировщике. Однако, многие из медленных мопов находятся в начале каждой итерации. Тем не менее, есть основания надеяться, что это приблизится к насыщению портов 0 и 1. Если данные не нагреваются в кеше L1, пропускная способность кеша / памяти, вероятно, станет узким местом.
Чередование операций из нескольких цепочек dep также будет лучше. Когда clang развертывается, он помещает все 9 инструкций для одной итерации перед всеми 9 инструкциями для другой итерации. В нем используется удивительно небольшое количество регистров, поэтому было бы возможно смешать инструкции для 2 или 4 итераций. Такого рода вещи должны быть хороши для компиляторов, но это неудобно для людей.:/
Также было бы немного более эффективно, если бы компилятор выбирал режим адресации с одним регистром, чтобы загрузка могла микросинхронизироваться с vaddps, GCC делает это.