Почему SSE скалярный sqrt (x) медленнее, чем rsqrt (x) * x?

106

Я профилировал некоторые из наших основных математических вычислений на Intel Core Duo, и, глядя на различные подходы к вычислению квадратного корня, я заметил кое-что странное: используя скалярные операции SSE, быстрее получить обратный квадратный корень и умножить его. чтобы получить sqrt, чем использовать собственный код операции sqrt!

Я тестирую это с помощью цикла, например:

inline float TestSqrtFunction( float in );

void TestFunc()
{
  #define ARRAYSIZE 4096
  #define NUMITERS 16386
  float flIn[ ARRAYSIZE ]; // filled with random numbers ( 0 .. 2^22 )
  float flOut [ ARRAYSIZE ]; // filled with 0 to force fetch into L1 cache

  cyclecounter.Start();
  for ( int i = 0 ; i < NUMITERS ; ++i )
    for ( int j = 0 ; j < ARRAYSIZE ; ++j )
    {
       flOut[j] = TestSqrtFunction( flIn[j] );
       // unrolling this loop makes no difference -- I tested it.
    }
  cyclecounter.Stop();
  printf( "%d loops over %d floats took %.3f milliseconds",
          NUMITERS, ARRAYSIZE, cyclecounter.Milliseconds() );
}

Я пробовал это с несколькими разными телами для TestSqrtFunction, и у меня есть некоторые тайминги, которые действительно ломают мне голову. Хуже всего было использование встроенной функции sqrt () и возможность «оптимизировать» «умный» компилятор. При 24ns / float с использованием FPU x87 это было ужасно плохо:

inline float TestSqrtFunction( float in )
{  return sqrt(in); }

Следующее, что я пробовал, - это использовать встроенную функцию, чтобы заставить компилятор использовать скалярный код операции sqrt SSE:

inline void SSESqrt( float * restrict pOut, float * restrict pIn )
{
   _mm_store_ss( pOut, _mm_sqrt_ss( _mm_load_ss( pIn ) ) );
   // compiles to movss, sqrtss, movss
}

Это было лучше, 11,9 нс / плавучесть. Я также попробовал дурацкую технику аппроксимации Ньютона-Рафсона Кармака , которая работала даже лучше, чем оборудование, при 4,3 нс / число с плавающей запятой, хотя и с ошибкой 1 из 2 10 (что слишком много для моих целей).

Неприятность была, когда я попробовал операцию SSE для получения обратного квадратного корня, а затем использовал умножение, чтобы получить квадратный корень (x * 1 / √x = √x). Даже если это занимает две зависимые операции, он был самым быстрым решением на сегодняшний день, в 1.24ns / поплавком и с точностью до 2 -14 :

inline void SSESqrt_Recip_Times_X( float * restrict pOut, float * restrict pIn )
{
   __m128 in = _mm_load_ss( pIn );
   _mm_store_ss( pOut, _mm_mul_ss( in, _mm_rsqrt_ss( in ) ) );
   // compiles to movss, movaps, rsqrtss, mulss, movss
}

Мой вопрос в основном в том, что дает ? Почему встроенный аппаратный код операции извлечения квадратного корня в SSE медленнее, чем синтез его из двух других математических операций?

Я уверен, что это действительно стоимость самой операции, потому что я подтвердил:

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

( edit : stephentyrone правильно указывает, что операции с длинными строками чисел должны использовать векторизованные SIMD-упакованные операции, например, rsqrtps- но структура данных массива здесь предназначена только для целей тестирования: то, что я действительно пытаюсь измерить, - это скалярная производительность для использования в коде которые нельзя векторизовать.)

Crashworks
источник
13
х / sqrt (х) = sqrt (х). Или, другими словами: x ^ 1 * x ^ (- 1/2) = x ^ (1 - 1/2) = x ^ (1/2) = sqrt (x)
Crashworks
6
конечно inline float SSESqrt( float restrict fIn ) { float fOut; _mm_store_ss( &fOut, _mm_sqrt_ss( _mm_load_ss( &fIn ) ) ); return fOut; }. Но это плохая идея, потому что это может легко вызвать остановку загрузки-попадания-сохранения, если ЦП записывает числа с плавающей запятой в стек, а затем сразу же считывает их обратно - в частности, перестановка из векторного регистра в регистр с плавающей точкой для возвращаемого значения плохие новости. Кроме того, коды операций базовой машины, которые встроенные функции SSE представляют, в любом случае принимают адресные операнды.
Crashworks
4
Насколько важна LHS, зависит от конкретного поколения и степпинга данного x86: мой опыт показывает, что на любых устройствах вплоть до i7 перемещение данных между наборами регистров (например, с FPU на SSE на eax) очень плохое, в то время как круговой обход между xmm0 и стеком и обратно нет из-за переадресации магазина Intel. Вы можете рассчитать время сами, чтобы убедиться в этом. Как правило, самый простой способ увидеть потенциальную LHS - это посмотреть на выпущенную сборку и увидеть, где данные перебираются между наборами регистров; ваш компилятор может сделать умную вещь, а может и нет. Что касается нормализации векторов, я записал свои результаты здесь: bit.ly/9W5zoU
Crashworks
2
Для PowerPC - да: у IBM есть симулятор ЦП, который может прогнозировать LHS и многие другие пузыри конвейера с помощью статического анализа. Некоторые PPC также имеют аппаратный счетчик LHS, который вы можете опросить. Для x86 сложнее; хороших инструментов для профилирования меньше (VTune в наши дни несколько сломан), а переупорядоченные конвейеры менее детерминированы. Вы можете попытаться измерить его эмпирически, измерив количество инструкций за цикл, что можно сделать точно с помощью аппаратных счетчиков производительности. Регистры « инструкции удалены » и «общее количество циклов» могут быть прочитаны, например, с помощью PAPI или PerfSuite ( bit.ly/an6cMt ).
Crashworks
2
Вы также можете просто написать несколько перестановок в функции и рассчитать их время, чтобы увидеть, не пострадает ли какая-либо из них особенно от задержек. Intel не публикует много подробностей о том, как работают их конвейеры (то, что они LHS вообще - это своего рода грязный секрет), поэтому многое из того, что я узнал, я узнал, глядя на сценарий, который вызывает остановку на других арках (например, PPC ), а затем построить управляемый эксперимент, чтобы увидеть, есть ли он у x86.
Crashworks

Ответы:

216

sqrtssдает правильно округленный результат. rsqrtssдает приближение к обратной величине с точностью примерно до 11 бит.

sqrtssгенерирует гораздо более точный результат, когда требуется точность. rsqrtssсуществует для случаев, когда достаточно приближения, но требуется скорость. Если вы прочитаете документацию Intel, вы также найдете последовательность инструкций (аппроксимация обратного квадратного корня, за которой следует один шаг Ньютона-Рафсона), которая дает почти полную точность (~ 23 бита точности, если я правильно помню), и все же несколько быстрее чем sqrtss.

edit: Если скорость критична, и вы действительно вызываете это в цикле для многих значений, вы должны использовать векторизованные версии этих инструкций, rsqrtpsили sqrtpsоба из которых обрабатывают четыре числа с плавающей запятой на инструкцию.

Стивен Кэнон
источник
3
Шаг n / r дает точность 22 бита (удваивает ее); 23 бита - это полная точность.
Джаспер Беккерс
7
@ Джаспер Беккерс: Нет, не будет. Во-первых, float имеет точность 24 бита. Во- вторых, sqrtssэто правильно закругленные , который требует ~ 50 бит до округления, и не может быть достигнуто с помощью простой N / R итерации в одинарной точности.
Стивен Кэнон
1
Это определенно причина. Чтобы расширить этот результат: проект Intel Embree ( software.intel.com/en-us/articles/… ) использует векторизацию для своей математики. Вы можете скачать исходный код по этой ссылке и посмотреть, как они делают свои 3/4 D. Векторы. Их векторная нормализация использует rsqrt, за которым следует итерация newton-raphson, которая тогда очень точна и все еще быстрее, чем 1 / ssqrt!
Brandon Pelfrey
7
Небольшое предостережение: x rsqrt (x) приводит к NaN, если x равен нулю или бесконечности. 0 * rsqrt (0) = 0 * INF = NaN. INF rsqrt (INF) = INF * 0 = NaN. По этой причине CUDA на графических процессорах NVIDIA вычисляет приблизительные квадратные корни одинарной точности как обратное значение (rsqrt (x)), при этом оборудование обеспечивает как быстрое приближение к обратному, так и обратному квадратному корню. Очевидно, что явные проверки, обрабатывающие два особых случая, также возможны (но будут медленнее на GPU).
njuffa 04
@BrandonPelfrey В каком файле вы нашли шаг Newton Rhapson?
fredoverflow
7

Это верно и для деления. MULSS (a, RCPSS (b)) намного быстрее, чем DIVSS (a, b). Фактически, он все еще быстрее, даже если вы увеличите его точность с помощью итерации Ньютона-Рафсона.

Intel и AMD рекомендуют этот метод в своих руководствах по оптимизации. В приложениях, которые не требуют соответствия IEEE-754, единственной причиной использования div / sqrt является удобочитаемость кода.

Плевать
источник
1
Broadwell и более поздние версии имеют лучшую производительность деления FP, поэтому компиляторы, такие как clang, предпочитают не использовать обратный + Newton для скаляра на последних процессорах, потому что обычно он не быстрее. В большинстве циклов divэто не единственная операция, поэтому общая пропускная способность uop часто оказывается узким местом, даже если есть divpsили divss. См. Раздел Деление с плавающей запятой против умножения с плавающей запятой , где в моем ответе есть раздел о том, почему rcppsбольше не выигрывает пропускная способность. (Или выигрыш в задержке) и цифры по разделению пропускной способности / задержки.
Питер Кордес
Если ваши требования к точности настолько низки, что вы можете пропустить итерацию Ньютона, тогда да, a * rcpss(b)может быть быстрее, но это все равно больше, чем a/b!
Питер Кордес
5

Вместо того, чтобы давать ответ, который на самом деле может быть неверным (я также не собираюсь проверять или спорить о кеше и других вещах, скажем, они идентичны), я попытаюсь указать вам источник, который может ответить на ваш вопрос.
Разница может заключаться в том, как вычисляются sqrt и rsqrt. Вы можете прочитать больше здесь http://www.intel.com/products/processor/manuals/ . Я предлагаю начать с чтения о функциях процессора, которые вы используете, есть некоторая информация, особенно о rsqrt (процессор использует внутреннюю таблицу поиска с огромным приближением, что значительно упрощает получение результата). Может показаться, что rsqrt настолько быстрее, чем sqrt, что 1 дополнительная операция mul (что не слишком дорого) может не изменить ситуацию.

Изменить: несколько фактов, которые стоит упомянуть:
1. Однажды я делал некоторые микрооптимизации для своей графической библиотеки и использовал rsqrt для вычисления длины векторов. (вместо sqrt я умножил свою сумму квадратов на rsqrt, что и было сделано в ваших тестах), и он работал лучше.
2. Вычисление rsqrt с использованием простой таблицы поиска может быть проще, как и для rsqrt, когда x переходит в бесконечность, 1 / sqrt (x) переходит в 0, поэтому для малых x значения функции не меняются (сильно), тогда как для sqrt - он уходит в бесконечность, так что это тот простой случай;).

Кроме того, пояснение: я не уверен, где я нашел его в книгах, которые я связал, но я почти уверен, что читал, что rsqrt использует некоторую таблицу поиска, и ее следует использовать только тогда, когда результат не обязательно быть точным, хотя - я тоже могу ошибаться, как это было некоторое время назад :).

Марцин Дептула
источник
4

Ньютон-Рафсон сходится к нулю, f(x)используя приращения, равные -f/f' где f'- производная.

Для x=sqrt(y), вы можете попробовать решить f(x) = 0для xиспользования f(x) = x^2 - y;

Тогда приращение: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x с медленным делением.

Вы можете попробовать другие функции (например f(x) = 1/y - 1/x^2), но они будут столь же сложными.

Посмотрим 1/sqrt(y)сейчас. Вы можете попробовать f(x) = x^2 - 1/y, но это будет не менее сложно: dx = 2xy / (y*x^2 - 1)например. Один неочевидный альтернативный вариант f(x):f(x) = y - 1/x^2

Затем: dx = -f/f' = (y - 1/x^2) / (2/x^3) = 1/2 * x * (1 - y * x^2)

Ах! Это нетривиальное выражение, но в нем есть только умножение, а не деление. => Быстрее!

И: new_x = x + dxзатем полный этап обновления гласит:

x *= 3/2 - y/2 * x * x что тоже легко.

скаль
источник
2

Уже несколько лет назад на этот вопрос есть ряд других ответов. Вот что было правильным при консенсусе:

  • Инструкции rsqrt * вычисляют приближение к обратному квадратному корню, что составляет примерно 11–12 бит.
  • Это реализовано с помощью таблицы поиска (то есть ROM), индексированной мантиссой. (Фактически, это сжатая справочная таблица, похожая на старые математические таблицы, в которой используются корректировки младших битов для экономии на транзисторах.)
  • Причина, по которой он доступен, заключается в том, что это начальная оценка, используемая FPU для «реального» алгоритма извлечения квадратного корня.
  • Также есть примерная обратная инструкция rcp. Обе эти инструкции являются ключом к пониманию того, как FPU реализует извлечение квадратного корня и деление.

Вот что ошиблось в консенсусе:

  • FPU эпохи SSE не используют метод Ньютона-Рафсона для вычисления квадратных корней. Это отличный программный метод, но было бы ошибкой реализовывать его таким образом на оборудовании.

Алгоритм NR для вычисления обратного квадратного корня имеет этот шаг обновления, как отмечали другие:

x' = 0.5 * x * (3 - n*x*x);

Это много умножений, зависящих от данных, и одного вычитания.

Далее следует алгоритм, который фактически используют современные FPU.

Дано b[0] = n, предположим, что мы можем найти ряд чисел Y[i], b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2приближающихся к 1. Затем рассмотрим:

x[n] = b[0] * Y[0] * Y[1] * ... * Y[n]
y[n] = Y[0] * Y[1] * ... * Y[n]

Ясно x[n]подходит sqrt(n)и y[n]подходит 1/sqrt(n).

Мы можем использовать шаг обновления Ньютона-Рафсона для получения обратного квадратного корня, чтобы получить хорошее Y[i]:

b[i] = b[i-1] * Y[i-1]^2
Y[i] = 0.5 * (3 - b[i])

Затем:

x[0] = n Y[0]
x[i] = x[i-1] * Y[i]

и:

y[0] = Y[0]
y[i] = y[i-1] * Y[i]

Следующее ключевое наблюдение таково b[i] = x[i-1] * y[i-1]. Так:

Y[i] = 0.5 * (3 - x[i-1] * y[i-1])
     = 1 + 0.5 * (1 - x[i-1] * y[i-1])

Затем:

x[i] = x[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = x[i-1] + x[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))
y[i] = y[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = y[i-1] + y[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))

То есть, учитывая начальные x и y, мы можем использовать следующий шаг обновления:

r = 0.5 * (1 - x * y)
x' = x + x * r
y' = y + y * r

Или, что еще интереснее, мы можем установить h = 0.5 * y. Это инициализация:

Y = approx_rsqrt(n)
x = Y * n
h = Y * 0.5

И это шаг обновления:

r = 0.5 - x * h
x' = x + x * r
h' = h + h * r

Это алгоритм Гольдшмидта, и он имеет огромное преимущество, если вы реализуете его на аппаратном уровне: «внутренний цикл» состоит из трех операций умножения и сложения и ничего больше, а два из них независимы и могут быть конвейеризованы.

В 1999 году FPU уже нуждались в конвейерной схеме сложения / вычитания и конвейерной схеме умножения, иначе SSE не был бы очень «потоковым». В 1999 году потребовалась только одна из каждой схемы, чтобы реализовать этот внутренний цикл полностью конвейерным способом, не тратя много оборудования только на извлечение квадратного корня.

Сегодня, конечно, мы представили программисту слияние умножения-сложения. Опять же, внутренний цикл - это три конвейерных FMA, которые (опять же) обычно полезны, даже если вы не вычисляете квадратные корни.

Псевдоним
источник
1
Связано: Как sqrt () GCC работает после компиляции? Какой метод рута используется? Ньютон-Рафсон? содержит ссылки на проекты аппаратных модулей исполнения div / sqrt. Быстрый векторизованный rsqrt и обратный с SSE / AVX в зависимости от точности - одна итерация Ньютона в программном обеспечении, с FMA или без него, для использования с _mm256_rsqrt_psанализом производительности Haswell. Обычно это хорошая идея, если у вас нет другой работы в цикле, и это сильно затруднит пропускную способность делителя. HW sqrt - это один uop, поэтому его можно смешивать с другой работой.
Питер Кордес
-2

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

Витек
источник
Очевидно неверно. FMA зависит от текущего режима округления, но имеет пропускную способность два за такт на Haswell и более поздних версиях. Имея два полностью конвейерных блока FMA, Haswell может иметь в полете до 10 FMA одновременно. Правильный ответ rsqrt«s намного ниже точность, что означает гораздо меньше работы (или вообще?) После табличной получить начальное предположение.
Питер Кордес