Опять о современных 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:

__m128 dsum1=_mm_set1_ps(0.f),dsum2=_mm_set1_ps(0.f);
__m128 ones = _mm_set1_ps(1.f);
__m128 maxv = _mm_set1_ps(maximum-25);
for (row=0; row < bottom; row += 8)
   for (col=0; col < right; col += 8) {
      __m128 sum1=_mm_set1_ps(0.f),sum2=_mm_set1_ps(0.f);
      for (y=row; y < row+8 && y < bottom; y++)
        for (x=col; x < col+8 && x < right; x++)
        {
          __m128 pixel = _mm_load_ps(image[y*width+x]);
          __m128 mask = _mm_cmpge_ps(pixel,maxv);
          if(_mm_movemask_ps(mask)) goto skip_block;
          sum1 = _mm_add_ps(sum1,pixel);
          sum2 = _mm_add_ps(sum2,ones);
        }
       dsum1=_mm_add_ps(dsum1,sum1);
       dsum2=_mm_add_ps(dsum2,sum2);
skip_block: ;
   }
И там даже финальное вычисление (dsum[i+4]/dsum[i]) можно будет векторизовать, хотя это и не нужно c точки зрения перформанса.

И что бы вы думали:

  • 12-15ms, т.е. в 3-4 раза быстрее (а после распараллеливания - можно вообще об этом месте забыть)
  • работает на всех актуальных горшках
  • Я бы не сказал, что код стал сильно запутаннее.
Т.е. я опять к тому, что вера в современные компиляторы - сильно преувеличена.

Comments

А какие компиляторы и с какими опциями ты используешь?

2012-й Visual Studio, 2013-й интел и gcc 4.7 делают в этом месте примерно одинаковый код.
Он SSE-шный, но скалярный.

clang не пробовал?

Надо пробовать clang с Polly-оптимизатором, а у меня все руки не доходят это место вдумчиво изучить.

И, кстати, понятно почему они так делают: они не могут допереть, что попиксельное сравнение/сложение-или-выбрасывание-блока - эквивалентны тому же самому, но по 4 пикселя.

Потому что с точки зрения семантики C++ они не эквивалентны: после выполнения goto skip_block код никогда не загружает остальные пискелы/компоненты пикселов, и компилятор не может опеределить, что обращаться к остальным компонентам безопасно (т.е. не сгенерирует Access Violation).

Ну да.

Но чем вручную делать векторизуемый (конкретным) компилятором код - я вот лучше тупо нахреначу 5 команд _mm.

Да, это меня поимеет, когда я начну переносить на ARM. Но пусть лучше поимеет путем некомпилируемости, чем тем, что векторизуемый (одним компилятором) код окажется скалярным на другом компиляторе.

Интеловский компилятор векторизовал 3 из 7 циклов. Вижуал всего лишь один.

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

Внутренний цикл остался скалярным.

А вместо "переписать, чтобы смогли векторизовать" - я лучше уж сам векторизую. Быстрее же, никаких тебе проб и ошибок.

А если избавиться от _mm_movemask_ps у которой большая задержка, то совсем хорошо будет. Ну например использовать ptest.

Ну и http://agner.org/optimize/instruction_tables.pdf в качестве библии по скорости инструкций.

ptest - это же SSE4.1
У меня ограничение - SSE3 и не выше. Потому что AMD.

PTEST будет быстрее только на Sandy Bridge, и то только на один такт (по задержке)

Алексей, а можно Вам в привате несколько вопросов про интеллект компиляторов и паралеллизм задать?

Да, конечно, lexa@lexa.ru
Но обычно я даю ответы в духе "пробовать надо".

Спасибо. Сначала попробую то, что понял из презентации.

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

тупо на SSE2

Вместо этого можно использовать например Eigen:
http://eigen.tuxfamily.org/index.php?title=Main_Page

Или нечто-то подобное. Там эти intrinsics скрыты за слоем compile-time абстрации, и причём не только sse:

"Explicit vectorization is performed for SSE 2/3/4, ARM NEON, and AltiVec instruction sets, with graceful fallback to non-vectorized code."

Если же Eigen и подобные либы не подходят - можно сделать свою велосипедную абстракцию.

я скорее о legacy коде.

Ну вот есть, к примеру, dcraw. В ней есть буквально 5-8 мест, которые тормозят, ну 10.
Все эти 10 мест - без насилия над остальным кодом (9500 строк) - ускоряются в 3-4 раза ручной векторизацией и еще в 3-4 раза распараллеливанием. Получается отлично.

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

Ну если мест не так много, и не ломает делать везде ifdef'ы на разные архитектуры (либо вообще не нужно), то да - можно и вручную.

Я в GCC проверил два варианта
Первый вариант:

float in[]={0,0,0,0,0,0,0,0};
__m128 a, b,c;
scanf("%f%f%f%f%f%f%f%f", in,in+1,in+2,in+3,in+4,in+5,in+6,in+7);
a = _mm_loadu_ps(in);
b = _mm_loadu_ps(in+4);
c = _mm_add_ps(a , b );
float out[4];
_mm_store_ps(out, c );
printf("%f %f %f %f", out[0],out[1],out[2],out[3]);

Дает такой ассемблерный выхлоп в GCC:

...
call __isoc99_scanf
movups 0(%rbp), %xmm1
movups (%rbx), %xmm0
movl $.LC1, %edi
addps %xmm1, %xmm0
movl $4, %eax
movaps %xmm0, 64(%rsp)
cvtss2sd 64(%rsp), %xmm0
cvtss2sd 76(%rsp), %xmm3
cvtss2sd 72(%rsp), %xmm2
cvtss2sd 68(%rsp), %xmm1
call printf
...

И такой вариант с указателями вместо _mm_loadu_ps:

float in[]={0,0,0,0,0,0,0,0};
__m128 *a, *b,c;
scanf("%f%f%f%f%f%f%f%f", in,in+1,in+2,in+3,in+4,in+5,in+6,in+7);
a = (void *)in;
b = (void *)(in+4);
c = _mm_add_ps(*a , *b);
float out[4];
_mm_store_ps(out, c );
printf("%f %f %f %f", out[0],out[1],out[2],out[3]);

получаем:

...
call __isoc99_scanf
movaps 32(%rsp), %xmm0
movl $.LC1, %edi
addps 48(%rsp), %xmm0
movl $4, %eax
movaps %xmm0, 64(%rsp)
cvtss2sd 64(%rsp), %xmm0
cvtss2sd 76(%rsp), %xmm3
cvtss2sd 72(%rsp), %xmm2
cvtss2sd 68(%rsp), %xmm1
call printf
...

Если в первом случае вместо _mm_loadu_ps использовать _mm_load_ps, результаты будут одинаковы. Компилятор может определить ситуацию, когда в _mm_loadu_ps передаются выровненные данные и для этого случая использовать movaps вместо movups? И зачем тогда нужен _mm_load_ps, если можно через делать через указатели?

В чем ваш вопрос то?

Ну вопрос собственно вот:
Компилятор может определить ситуацию, когда в _mm_loadu_ps передаются выровненные данные и для этого случая использовать movaps вместо movups?
И зачем тогда нужен _mm_load_ps, если можно через делать через указатели?

Начнем с конца.
_mm_load_ps - читает 16 байт (в регистр). Инструкция нужна ровно для этого: если значение затем используется (несколько раз), то наверное его стоит поднять в регистр.

Далее. Компилятор, конечно, может породить код, который на рантайме проверяет, какой дали указатель и в зависимости от этого - выполняет movaps или movups. Теперь представим себе цикл, который выполняется, скажем, 80 млн. раз при обработке одного 80-mpix изображения и содержит в себе два load и 1 store.
У нас будет или 3 локальных проверки (на каждой итерации) или развосьмерение кода (и 3 проверки всего, на все итерации). Оба плохо, как мне кажется.

_mm_load_ps - читает 16 байт (в регистр).

Все немного сложнее.
Из кода

float in[]={0,0,0,0,0,0,0,0};
__m128 a, b,c;
scanf("%f%f%f%f%f%f%f%f", in,in+1,in+2,in+3,in+4,in+5,in+6,in+7);
a = _mm_load_ps(in);
b = _mm_load_ps(in+4);
c = _mm_add_ps(a , b );
float out[4];
_mm_store_ps(out, c );
printf("%f %f %f %f", out[0],out[1],out[2],out[3]);

получается только одно чтение в регистр. Операция сложения получается такой:
...
movaps 32(%rsp), %xmm0
movl $.LC1, %edi
addps 48(%rsp), %xmm0
...

Т.е. _mm_load_ps не обязательно поместит что-нибудь в регистр.

Компилятор, конечно, может породить код, который на рантайме проверяет, какой дали указатель

А может ли компилятор понять это через статический анализ? Например, если указатель изначально поставлен правильно и инкрементируется в цикле на правильное число?

> А может ли компилятор понять это через статический анализ? Например, если указатель изначально поставлен правильно и инкрементируется в цикле на правильное число?

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

А касаемо кода - воткните -O3 и код станет другой.

>А касаемо кода - воткните -O3 и код станет другой.
Никаких изменений
Linux gcc-4_5-branch revision 167585