Об автоматической векторизации

Провел на поминавшемся вчера ISPC еще один тест, на применимость ровно в том месте, куда он лучше всего приспособлен.

Есть такое ужасное место в обработке изображений (уже поминавшееся в этом блоге): преобразование из линейной гаммы в sRGB-гамму или в Lab. Там в формуле сначала линейный участок, а потом степенной. Вот как это выглядит, если делать в лоб для плавающей точки:

void linear2srgb(float *in, float *out)
{
   for(int i = 0; i< DATA_SIZE; i++)
        out[i] = ((in[i]<=0.0031308f)? 12.92f*in[i] : (1+0.055f)*powf(in[i],1/2.4f)-0.055f);
}

Ну, понятно, для (целочисленных) 8-16-битных данных составляют лукап-таблицу, а дальше фигачат по ней, а вот для большей битности таблица получается изрядно большой. Вычислять эту штуку в лоб мешают две вещи:

  • ветвление на каждое входное данное;
  • возведение в степень, которое тоже очень медленное: на SSE/AVX такой функции нет, на FP87 есть, но безобразно медленная.
К примеру, при обработке плавающих данных LCMS (преобразование в Lab, преобразование по матричному профилю в sRGB) процентов 90 времени уходит именно на вышепоказанную операцию (правда в LCMS это место еще сделано потрясающе неудачно с точки зрения производительности).

Как я уже писал, правильное решение заключается в замене вышепоказанной функции чем-то приличным, скажем для кубических сплайнов с таблицей в 4к строк максимальная ошибка по всему диапазону не превышает 10-6, что для всех применений достаточно, при скорости порядка 1.2-1.5Gb/sec на одно процессорное ядро. Но одна строчка кода превращается в несколько десятков, таблицу коэффициентов сплайнов надо еще построить, что мучительно.

Посмотрим, что можно сделать с помощью ISPC и можно ли вообще что-то. В этот раз мне оказалось проще считать производительность в мегабайтах в секунду, мегапиксели из прошлого текста какие-то непрозрачные (поди упомни, что надо на 12 умножать).

На i7-2600K, разогнанном до 4.5Ghz, в один поток получается:

1.Visual C++ 2010:

  • Параметры по умолчанию: 116 Mb/sec. Весь пар уходит в библиотечную Cipow_pentium4
  • /arch:SSE2 /fp:fast: 137Mb/sec, пар уходит в libm_sse2_pow

2. Берем тот же код (те же две строчки) и компилируем ISPC. Получаем 410Mb/sec (компиляция под SSE4, математическая библиотека встроенная). Втрое быстрее, чем Visual C++.

3. Читаем мануал на ISPC, находим там про all() и any() и модифицируем код:

export void srgb_ispc(uniform float in[], reference uniform float out[])
{
   for(uniform int i = 0; i < DATA_SIZE; i+=programCount)
   {
        float val = in[i+programIndex];
        bool b = val <=0.0031308;
        if(all(b)) // все на линейном участке
           out[i+programIndex] = 12.92*val;
        else if(!any(b)) // все на степенном участке
           out[i+programIndex] = (1.055)*pow(val,1.0/2.4)-0.055;
        else
           out[i+programIndex] = b? 12.92*val:(1.055)*pow(val,1.0/2.4)-0.055;  
   }
}
Вместо одной строчки в теле цикла получается 8 (терпимо), правда формула удвоилась и теперь если она поменяется, то придется в двух местах править. Правильно было бы подефайнить, конечно.

Получаем: 523Mb/sec, еще на четверть больше.

Переходим в 64-битный режим, получаем еще больше: 730Mb/sec.

4. А что же интеловский компилятор? А вот что:

  • /arch:SSE2 /fp:fast: 313Mb/sec, в 2.5 раза быстрее VS2010. За счет того, что зовет не libm_sse_pow, а libm_sse_powf, которая должна быть минимум вдвое быстрее.
  • /arch:AVX /fp:fast: те же ~315Mb/sec.
5. Автовекторизация

Интеловский компилятор может сделать автоматическую векторизацию цикла одним из двух способов

  1. Добавим #pragma ivdep и #pragma vector always к телу цикла и включим векторизацию явно.
  2. Или включим ее неявно таким способом:
    • Включим автоинлайн (/Ob2)
    • Код используется в том же файле, где и находится
Во втором случае компилятор видит контекст использования и видит, что векторизовать можно (в частности, входной и выходной массивы не перекрываются), в первом - мы сами это ему сказали.

При автовекторизации вместо libm используется SVML (short vector math library) и получается 530Mb/sec для /arch:SSE2 (степень вычисляется svml_powf4) и 720Mb/sec для /arch:AVX (svml_powf8). Судя по профайлу, маскирование результатов в зависимости от значений - тоже происходит, svml_pow.._mask() там фигурирует.

Мораль: В таком простом случае интеловский С++ компилятор справился с кодом лучше, чем ISPC. Для ISPC пришлось руками расписывать случаи срабатывания/не срабатывания всех условий, а C++ компилятор дал тот же перформанс на исходном однострочном варианте.

В случае использования Visual C++, использование ISPC для критических мест, очевидно, имеет смысл, если лень руками писать. В частности, встроенная математическая библиотека (pow(), exp(), log(), sin()/cos()/tan()) у ISPC, очевидно, сильно быстрее системной, даже SSE-шной.

На других компиляторах эксперименты пока не проводились, но будут.

Update: вот какой важный слон остался непримеченным. Intel C++ свои 720Mb/sec достигает при использовании AVX-библиотеки, по 8 элементов за раз. С SSE-библиотекой только 530. В то же время, ISPC пока AVX не умеет и свои 720Mb/sec получил на 4-way коде. По аналогии с тем что происходит у Intel, на AVX мы вправе ожидать еще процентов 40 прироста, что делает ISPC потенциально более интересным.

Comments

Я несколько не в тему, но вот что привлекло мое внимание:

> возведение в степень, которое тоже очень медленное: на SSE/AVX такой функции нет, на FP87 есть, но безобразно медленная

Я в своей работе использую (вернее, использовал) причесанный мною код некого Julien Pommier, который написал очень быстрые (и "достаточно" точные) SSE инлайны для sin, cos, tan, exp и log.

А та же powf(in[i],1/2.4f) - это ведь exp(1/2.4f * log(in[i]))

В диапазоне (0.0001, 1] применение этих приближенных функций дает ошибку до 6,0e-8 (по сравнению с powf)

А скорость такая: На моем i5-2500, разогнанном до 3,8 ГГц при всех 4 ядрах я получил 1,4e+11 обработок 4-компонентных векторов, или 5,6e+11 обработанных float. OpenMP использовал для загрузки всех 4 ядер.

А нет, вру. Ошибся с расчетами. 5,6e+8 обработанных float в секунду.

1 float - 4 байта.
Т.е. около 2Gb/sec если бы мы только экспоненту считали.

А так как надо и экспоненту и логарифм, то значит надо пополам делить?

Плюс еще ветвление и если оно случилось, то всякий blendps (а проверка всяко есть). Ну то есть интеловские/ispc-шные 700Mb/sec получаются вполне адекватными.

Нет, это вместе и экспонента, и логарифм. И blendps - это небольшая добавка.

Я по той ссылке, что Вы дали ниже, вот что увидел: else if (__math_lib == __math_lib_ispc_fast)

Посмотрите, там наверняка есть настройка, какую математическую библиотеку использовать. Наверняка если fast подкрутить, то будет побыстрее 700МБ/с.

Увы, но для экспоненты там только одна версия, а экспонента жрет чуть меньше половины времени (проверил, заменив pow() на log())

Ну то есть есть еще куда расти: LLVM генерирует код хуже, чем написано ручками, неудивительно :)

А, нет, тоже две. Но отчего-то не помогает :)

О, нашел.

--math-lib=fast

Экспонента и логарифм жрут одинаково примерно.

Ага. Но там начинаются чудеса out-of-order:

out=1.055*in; - 145msec (на полугигабайте данных)
out=1.055*log(in); - 381 msec, с exp - столько же. Т.е. добавка +236
out=1.055*pow(in,1/2.4); - 800 (а ожидал бы я 381+236 = 617)
out=1.055*pow(in,1/2.4)-0.055; - 730 (т.е. меньше, хотя добавили вычитание)

Это все при одних и тех же параметрах компилятора, только я забыл посмотреть каких :)

Надо код читать, но его реально дохрена и лень.

fast отчего-то банально медленнее у меня, чем default. Разбираться лень :)

До меня доперло.

Вы же на четырех ядрах упражнялись. А я - на одном.
Если урезать осетра разделить вчетверо, то как раз нужные ~500-600Mb/sec у вас и получилось. У меня доходит до over 9000 720, но и мегагерцев поболее будет.

Собственно, в регистре 16 байт, 720/16 = 45 млн. обработок вектора на ядре, на четырех ядрах было бы 1.8e+8. Разница в тактовых частотах: 4500/3800 = 1.18, разница в быстродействии 1.8/1.4=1.28, получается ISPC ручной код маленько уделал :).

Правда не факт, что встроенная реализация exp(log()) у ISPC такая же точная, как обсуждаемая. Ну и скорость памяти тоже важна, просто на чтение-запись должно процентов 15 уходить.

А, даже так. Удивительно.

А, точно :)

Не, у меня все в L1 крутится. Я ж не пропускную способность памяти мерял, в конце концов. Но для того, чтобы понять, кто кого уделал, нужные все-таки более точные измерения и расчеты :)

Э, фетч, запись и ветвление - это процентов 20 времени.

А так, конечно да, кто тут сильнее, слон или кит, это еще посмотреть.

Я наоборот рад буду, если ISPC быстрее. Значит, сравнительно неплохие функции используются.

Да там же tradeoff точность/скорость.

Про pow() вчера накопал вот это, кстати: http://www.hxa.name/articles/content/fast-pow-adjustable_hxa7241_2007.html

Там очень такой аккуратненький компактный полином в сорцах, надо при случае попробовать.

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

Если показатель степени известен на времени компиляции (как в рассматриваемом моем примере), то логарифм тоже на времени компиляции считается.

Нет, если основание степени известно на время компиляции, то можно было посчитать логарифм во время компиляции.

Там параметр есть в функции, ilog2 называется. Он равен:
for pow(r, val): ilog2 = log(r) / log(2) = log(r) * 1.44269504088896

И точно.
Мне показалось, что у них просто параметры в pow переставлены местами в тексте по ошибке.
Однако ж нет, из исходников следует, что ошибки вроде нет.

Тогда выигрыш не такой большой, как я думал.

Ну да, все так и делают. Либо по месту дотачивают (для нужной области входных значений), либо в общем виде.
Помянутые libm_sse2_pow(f) и svml_powf4(8) - это почти оно и есть, разве только не инлайны, а честные функции с call (инлайны тоже бывают).

В исходниках ispc, кстати, лежат готовые полиномы (которые потом попадают на вход компилятору и превращаются в инлайны). Вот тут, во второй половине файла: https://github.com/ispc/ispc/blob/master/stdlib.ispc

Только вот насчет e+11 - вы нолики не перепутали?
Ну собственно частота 3.8*1e+9, на 16 инструкций/такт (AVX, dual issue) на 4 ядра = 2.4e+11 flop/s (при очень либеральном способе оценки), это уже в пересчете на скаляры.

Да, про нолики ошибся на 3 штуки, выше написал :)

попробуйте ещё IntelMKL, там есть функции аля SVML, только для длинных векторов.

По моим прикидкам, большого выигрыша не будет. Да, длинные версии pow(log+exp) получаются на 10-20 процентов быстрее (не помню точных цифр), но ведь ветвление никуда не девается, поэтому придется честно делать второй проход (он будет раз в 5 быстрее) и для маленьких входных значений пересчитывать линейной функцией.
По прикидкам, выигрыш если и будет, то совсем уж копеечный. Ну разве что в MKL все уже многоядерное внутри и не придется #pragma omp писать.

По постановке:
У меня, собственно, задача была - пощупать за вымя ISPC на реальной, простой модельной задаче.

А правильное решение выглядит совсем иначе, кубический сплайн на 4k-таблице фигачит в пару раз быстрее при более чем приемлемой точности. И по ядрам масштабируется на OpenMP очень хорошо.

А почему бы не использовать приближение minimax полиномом? Например, полиномиальная дробь -0.0282646 + 24.8353 x + 11724.6 x^2 + 625030. x^3 +
6.38388*10^6 x^4 + 1.2464*10^7 x^5 + 2.45062*10^6 x^6 -
187084. x^7)/(1 + 819.465 x + 93987.2 x^2 + 2.15735*10^6 x^3 +
1.05524*10^7 x^4 + 8.94367*10^6 x^5) даёт для вашей функции относительную ошибку 4.79681*10^-7, и скорее всего будет вычисляться быстрее, чем pow

Там в тексте есть ссылка, правильное на мой вкус решение - это кубический сплайн (http://blog.lexa.ru/2010/12/27/o_stepennoi_funktsii_openmp_i_prochikh_gr...)

Но полиномом тоже можно, просто таблицу для сплайна я умею генерировать, а для minimax - нет.

Сплайн, если фигачить по 4 числа за раз, это
- 1 фетч данных (4 числа)
- 1 округление
- 4 фетча из таблицы сплайнов (она влезает в L2 с запасом)
- транспонирование
- ну и собственно расчет кубического полинома.

Для AVX я думаю что надо не кубический сплайн, а 7-й степени, но таблицу меньше. Еще не пробовал.

ИМХО, сплайн и прочие табличные методы плохо сочетаются с SIMD, особенно с out-of-order SIMD. А коэффициенты MiniMax полинома умеет считать Maple и (с некоторыми ограничениями) Mathematica. Да и вручную их несложно посчитать.

С точки зрения моей практической задачи (imaging), sRGB gamma и Lab - есть частый, но частный случай "наложения тоновой кривой" на изображение. На практике же это может быть тоновая кривая из профиля (произвольная, заданная точками или кусочным полиномом), может быть модифицированный профиль "a-la Lab", по смыслу похожий, а по коэффициентам - чуть другой и т.п.

То есть речь идет о том, что коэффициенты в любом случае табличные. Просто их может быть десяток-полтора (общий полином на всю область) или полтора десятка тысяч (сплайн с 4k точек).

Правда есть еще OpenCL (и вообще LLVM), то есть код для общего полинома можно на ходу сгенерировать. Тоже интересно.

А ваш полином я тоже гляну. Единственное что я не понял - это где должна стоять открывающая скобка (у вас скобки несбалансированы), но из общих соображений - в самом начале.

Я не поленился, формулу скопипейстил и посчитал дельты на диапазоне 0-1 с шагом 0.001

Хочу уточнить - вы всю формулу приблизили полиномом или только степенную часть? А то у меня получается, что на степенной части все отлично, а вот на стартовом линейном участке - полная катастрофа.

Возможно, конечно, что дело в точности в этом месте: если у меня x 1E-3, то в седьмой степени это будет минус 21-я.

Я считал minimax полином для степенной части функции на диапазоне от 0.0031308 до 1. Минимизировалась относительная ошибка.

А.

Для степенной части отлично работает, правда скорость я не мерял. И там, конечно, седьмая степень еще не больно бъется.
А вот в линейной части может быть неприятно, все-таки более-менее реальные линейные значения - это порядка 10^-4 (12 стопов в фотографическом смысле), а малоинтересные, но на практике имеющиеся - это 10^-5 (14-битный АЦП т.е. 2^-14). В седьмой степени - уже на грани фола для float, не говоря о накоплении ошибок.

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

А что мешает посчитать линейные значения отдельно + сделать blend? При грамотной реализации это будет почти бесплатно.

Да ничего, если мы говорим о конкретной формуле.

Но конкретная формула - это частный случай. А хочется, на самом деле, общего и достаточно быстрого решения для накладывания произвольной "кривой" (в терминах фотошопа) на пикселы. Сплайны, вот, работают отлично - если вылеплены руками. А полиномы еще не попробовал.

По алгоритму Ремеза стоит посмотреть вот эту реализацию для Matlab. Хотя эта реализация у меня так и не заработала, к ней прилагается pdf, которого мне хватило для реализации алгоритма Ремеза самому.

А вообще по этой теме рекомендую почитать:
Robert Green "Faster Math Functions" (статья, презентация (часть 1), презентация (часть 2))
Jean-Michel Muller "Elementary Functions: Algorithms and Implementation"

Ой какие клевые презентации, спасибо!

Книжку я видел и она у меня даже есть в PDF, но очень не хочется разбираться еще и в этом, я же "прикладной" программист :)

А какой нибудь _mm_maskmoveu_si128 и совсем без ветвлений это очень плохо будет?
Маску сделать вычитанием из 0 как целых. Мы вроде на x86. Должно работать.
Правда AMD еще есть. Там maskmovdqu вроде хуже работает.
(С оптимизациями сталкиваться приходилось давно и не много.)

Ну насколько я вижу, ISPC работает так:
- вычисляет две ветки (если all/any убрать)
- вычисляет маску через cmpltps
- смешивает результаты двух веток через blendvps
- пишет результат
По смыслу то же самое, но наверное смешение в регистре, а потом одна запись - выгоднее, чем две записи.

А у AMD blendvps вроде нету?
Хотя если для проверки что ISPC умеет...

А для SSE2 оно делает через and и or, маска та же. разница не особо принципиальная.