SSE

Про AVX2 и размывание кэшей

Вот наконец удалось ощутить офигенную пользу от AVX2, причем двойную. Вот такой вот код:

_mm256_stream_si256((__m256i *)&drowp[col], _mm256_i32gather_epi32((const int*)table, _mm256_cvttps_epi32(_mm256_load_ps(&srowp[col])), 4));

в 4.5-5 раз быстрее, чем простой SSE2 аналог (в котором, понятно, нет gather) и в ~6 раз быстрее скалярного C-кода:

drowp[col] = table[(unsigned)srowp[col]];

Рассмотрение всего хозяйства под микроскопом показало, что основной взнос в результат дает _mm256_stream, а вовсе не gather. Стоит заменить stream на store, как все сразу портится. По достаточно очевидной причине: длина строки и drowp и srowp -...

Опять о современных CPU

Продолжаю набирать материал к вот этой презентации

Вот такой вот код (это автобаланс белого из dcraw в моем вольном изложении, слегка сокращенный, только для 4-компонентных а не байеровских изображений):

float sum[8],dsum[8]={0,0,0,0,0,0,0,0};
for (row=0; row < bottom; row += 8)
   for (col=0; col < right; col += 8) {
      for(i=0;i<8;i++) sum[i]=0.f;
      for (y=row; y < row+8 && y < bottom; y++)
        for (x=col; x < col+8 && x < right; x++)
          for(i=0; i<4; i++) {
            val = image[y*width+x][i];
            if (val > maximum-25) goto skip_block;
            sum[i] += val;
            sum[i+4]++;
           }
       for(i=0;i<8;i++) dsum[i] += sum[i];
skip_block: ;
   }
То есть, если простыми словами: берем блоки 8x8, если в этом блоке все пиксели имеют значения меньшие (maximum-25), то значения пикселов усредняем, если хоть один пиксел вылетает за указанный лимит - блок выбрасываем.

На моей тестовой 9-Mpix картинке (Nikon D800 "интерполированный в half") этот код работает 45-50 миллисекунд. Что для автобаланса белого как-то неприлично много.

Перепишем его тупо на SSE2:

Об эффективном использовании современных CPU

Практика вот к этой презентации:

Берем 36Mpix файлик с D800, распаковываем егонный RAW в 16-битный однокомпонентный битмеп, дальше начинаем процессить.

Процессим без "интерполяции", т.е. 4 пикселя исходного байера образуют один выходной пиксель (режим half_size у LibRaw/dcraw). Получаем такие вот времена:

  1. LibRaw::dcraw_process() плюс формирование RGBA-битмепа: 420ms.
  2. Перепишем этот самый dcraw_process() на SSE3, процессить будем в плавучке (с эмуляцией особенностей dcraw), выдаем такой же 8-битный RGBA: 110ms (и более-менее понятно где еще выиграть миллисекунд 20).
  3. Добавим в предыдущий суп еще: подсчет RAW-гистограммы и сохранение исходного float-битмепа нетронутым (и без дублирования по памяти), т.е. баланс белого и конверсию цвета делаем два раза, один раз для подсчета гистограммы результирующего файла, а второй раз - прямо на вывод. 180ms.
Это все был один поток. Распараллелим его:
  1. "распараллеленный dcraw_process": 130ms (так в RawDigger сделано, тамошний RGB render устроен именно так, но гистограммы и статистика в эти 130ms не входят, равно как и битмеп для показа там готовится иначе и потому дольше).
  2. "распаралеленный ассемблерный dcraw_process": не делал, ожидаю <50ms (потому что вариант с двойным вычислением, как следующий, но без гистограмм - 57ms).
  3. параллельная ассемблерная версия с гистограммами, сохранением float RAW: 75ms
Сравнивая последний вариант с последовательным C-шным, нужно понимать, что в C-шном варианте еще где-то 200ms придется на гистограммы и еще 200 на конверсию int16 - float. То есть реальное ускорение от SSE и параллельной обработки - раз в 10 (75ms против 800). И это оптимизированный C-шный, в LibRaw это место заметно пооптимизировано в сравнении с исходным dcraw.

Минусы: потрачено изрядно времени (несколько дней), окончательного кода 25кб, а промежуточных вариантов было в разы больше. Ну и две баги ловились достаточно мучительно.

Презентация с Highload

По просьбам трудящихся масс, моя презентация с Хайлоада:

Никаких откровений нет, задача была - показать что multicore/simd - это очень просто и стоит того. Читатели моего уютненького легко узнают примеры 1.5-2-летней давности.

Анимированность Slideshare порезала, но вроде накрывающих друг друга картинок у меня нет.

Про AVX и ISPC

Разработчики Intel SPMD Program Compiler, который в этом блоге уже несколько раз поминался, выпустили версию 1.0.9 в которой

=== v1.0.9 === (26 September 2011)

The binary release of v1.0.9 is the first that supports AVX code generation. Two targets are provided: "avx", which runs with a programCount of 8, and "avx-x2" which runs 16 program instances simultaneously.

Честь им за это и хвала, мои попытки (не слишком настойчивые) самостоятельно собрать ISPC с LLVM3 так и не увенчались успехом, то какие-то ошибки в LLVM-овских H-файлах, то не линкуется, не больно то и хотелось.

Так как предыдущие тесты никуда не делись, я их переделал с новым ISPC, как с SSE4, так и с AVX.

О векторном умножении: нет гигапикселя в секунду

Тема матричного цветопреобразования не отпускает нас.

Наш читатель maratyszcza намекает нам, что haddps - тоже хорошая инструкция.

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

О векторном умножении - финал

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

Ассемблерные упражнения на тему векторных расширений.

Опять про movntps

Читал комментарии к прошлому посту про movntps, много думал. Вспоминал, что на i7-920 выигрыш от movntps был и значительный.

Пришлось потратить час-полтора на очередной набор микробенчмарок. Картинка получилась куда сложнее, чем я сначала подумал.

О векторных расширениях gcc/clang (2)

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

К сожалению, предложенное там решение (офигенно быстрое!) считает неправильно, но направление движение указано верно и мы приходим к такому варианту:

  1. транспонируем матрицу, на которую умножаем, дополним нулями правую колонку, чтобы вышло 4x4
  2. Каждое из (четырех) входных значений - размножим на вектор.
  3. Нужный нам результат - это SIMD-сумма SIMD-произведений вышеупомянутых векторов на строки вышеупомянутой транспонированной матрицы.
Короче, проще кодом:

О компиляторах и процессорах

В комментариях к моему посту о целых числах и плавающей точке мне посоветовали обратить внимание на векторные типы данных и сравнить их по производительности со скалярными типами (ну и ручным SSE-ассемблером тоже).

Я обратил и сравнил, но в процессе получилась масса побочных результатов (разные архитектуры, разные компиляторы), которые жалко выкинуть, а хочется опубликовать.

О legacy и форматах данных

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

А я вот тут читаю ассемблер, порожденный компилятором из C-шного кода обработки изображения. Много думаю.

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

  • от малейшего чиха переполняется или превращается в тыкву обнуляется, нужно постоянно следить за диапазоном;
  • векторные (SSE,AVX) операции с этим типом - очень ограничены, да и не векторные - тоже.
В результате, сплошные мучения компилятору, а результат - медленно работает. Скажем, вот такой вот код (из dcraw), который делает преобразование по матричному профилю для 3-4 входных каналов (цветов) и трех выходных:
out[0] = out[1] = out[2] = 0;
for (сc=0;сc<colors;сc++) {
  out[0] += out_cam[0][сc] * img[сc];
  out[1] += out_cam[1][сc] * img[сc];
  out[2] += out_cam[2][сc] * img[сc];
}
for(cс=0;сc<3;сc++) img[сc] = CLIP((int) out[cс]);
img[] - unsigned short, out[] - int, out_cam[] - float.

Смотрю на код, который порождает интеловский компилятор. Ну код, да. (три-)четыре load по 2 байта (количество loads зависит от colors /количества цветов/), дальше "вручную" выписанный dot product (нормальный, насколько это возможно), ну и обрезание до диапазона 0-65к и store.

Скорость работы этого кода - 100 мегапикселей в секунду на Sandy Bridge 4.5Ghz (в один поток, понятно что параллелится это на ура). Как-то не очень....

Да, считаем в мегапикселях т.к. в unsigned short у нас 8 байт на (4-компонентный) пиксель, а для float/int - 16 байт.

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

Провел на поминавшемся вчера 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 и можно ли вообще что-то.

Об Intel ISPC

Интел выкатил в опенсорс такую вот игрушку: Intel SPMD Program Compiler.

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

void vector_sum(uniform float in[], uniform int count, reference uniform float out[]) {
    float sum = 0.0;
    for (uniform int i = 0; i < count; i += programCount) {
        int index = i + programIndex;
        sum+= in[index];
    }
    out[programIndex] = sum;
}
И оно у вас пойдет фигачить шириной programCount, ну там дальше надо будет элементы out просуммировать уже снаружи.

OpenMP: Intel vs Visual Studio

Непонятки с производительностью OpenMP в случае VisualStudio (см. еще обсуждение в каментах) заставило разобраться.

Выяснилось, что две переменные, которые read-only и вообще объявлены как const - считаются разделяемыми. Как следствие - локи при обращении и прочие удовольствия.

После устранения этого недоразумения (#pragma omp ... private(...)) и замены _mm_store_ps() на _mm_stream_ps(), чтобы в обеих компиляторах выходной массив писался через movntps, все встало на свои места: С-шный код с powf() остается сильно быстрее у Интела т.к. он векторизует powf в _svml_powf4, а SSE-код работает с практически одинаковой скоростью.

Его пример другим наука:

#pragma omp parallel for default(none)
и явное прописывание shared/private/whatever для всех переменных.

О степенной функции, OpenMP и прочих граблях

Есть вот такое вот известное выражение, переложеное из википедии:

#define To_sRGB(q) (((q)<=0.0031308f)? 12.92f*(q) : (1+0.055f)*powf(q,1/2.4f)-0.055f)
В том смысле, что у sRGB внизу - линейный участок, а дальше степенная функция.

Если его напрямую запрограммировать (вот прямо с powf()), то в зависимости от компилятора получается от 7 Mpix/sec (Visual C++ 2010, один поток) до 100Mpix/sec (Intel Parallel Studio XE, включен OpenMP). Пиксель - это 3 компонента по 4 байта (float), то бишь гигабайт с небольшим в секунду - это максимум.

Как-то медленно.

В-принципе, бывает быстрое приближенное возведение в степень для ограниченного диапазона значений. Более того, можно просто приблизить sRGB-кривую полиномом не очень высокой степени и получить вполне приемлемый для практики результат. Но хотелось заодно решить задачу быстрого наложения произвольной кривой.

Полна чудес могучая природа

А в SSE4.1 оказывается есть минимум-максимум для 16-битных целых. 8 штук зараз.

Хреначит 11 гигабайт в секунду, сдается мне что параллелить по ядрам дальше смысла нет, упрешься в память.

И dot-product есть, правда только для float/double.

Чешутся руки забить на владельцев AMD и всего младше Penryn, держите меня...

Subscribe to SSE