Приемы программирования SSE/AVX [Обмен опытом]

[GAA]

В широко известных книгах: Юров Assembler. 2-е изд., 2003; Магда Ассемблер для процессоров Intel., 2006 использование SSE описано практически без примеров.
В руководствах Intel® 64 and IA-32 Architectures Optimization Reference Manual, 2013 и Intel® 64 and IA-32 Architectures Software Developer’s Manual, на мой взгляд, изложение чрезвычайно сжатое.
Предлагаю в этой теме обмениваться опытом. Будет здорово, если другие участники приведут свои примеры.

Начну с довольно тривиального примера сравнения двух чисел с плавающей точкой.

1. [Floating Point] Использование comiss (comisd) и cmpps (cmppd)

В SSE для (скалярного) сравнения двух single (double) есть инструкция, модифицирующая флаги ZF, PF и CF:
comiss xmm1, xmm2/r32 (comisd xmm1, xmm2/r64, соответственно).

Также есть инструкции cmpss (cmpsd) для (скалярного) сравнения двух single (double), которая заносит в приемник маску: единицы, если условие выполнено и нули в противном случае:
cmpss xmm1, xmm2/r32, imm8 (cmpsd xmm1, xmm2/r64 imm8).
и для сравнения упакованных (packed) single (double):
cmpps xmm1, xmm2/m128, imm8 (cmppd xmm1, xmm2/m128, imm8)

Значения масок обычно используется двумя способами:

A. Старший бит каждой маски переносится в регистр CPU
movmskps reg, xmm (movmskpd reg, xmm)
и побитно анализируется (командами, модифицирующими флаги, например bt) с последующим переходом на нужную ветвь. Либо значение reg используется в качестве индекса для перехода по адресу (хранимому в массиве адресов) или для загрузки констант из массива констант.

B. При помощи масок и логических инструкций (pand, pandn, por,…) модифицируются данные непосредственно в регистрах xmm.

Поскольку условные переходы и зависимость по данным замедляют скорость выполнения программы, то по возможности надо использовать второй способ (способ B).

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

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

1. Скалярный 1. Загружать в регистр xmm по очереди элементы массива и сравнивать при помощи comisd. Затем в случае выполнения условия — накапливать сумму, а в случае не выполнения — обходить сложение при помощи инструкции условного перехода.

2. Скалярный 2. Загружать в регистр xmm по очереди элементы массива, копировать в другой регистр и сравнивать копию с заданной величиной при помощи cmppd. Затем выполнять pand маски со значением элемента массива и безусловно выполнять сложение (в результате pand c содержащей нули маской будет получен ноль и прибавление нуля не приведёт к изменению суммы).

Прим. В AVX2 многие инструкции SSE стали не «разрушающими» (с тремя операндами; результат сравнения заносится в третий операнд). В этом наборе инструкций копировать элемент массива для сравнения не нужно.

3. Упакованный (packed). Выполняется накопление сумм четных и нечетных элементов массива одновременно. Для сравнения используется cmppd, как в варианте 2. После вычисления сумм четных и нечетных элементов они складываются (инструкция горизонтального сложения). Если массив содержит нечетное число элементов, то перед началом цикла подсчета сумм четных и нечетных элементов, выполняется скалярное накопление суммы как в варианте 2.

Ниже вставлен исходный текст 64-битного проекта Embarcadero Delphi XE7.
Прим. 1. Код легко может быть модифицирован под Borland Delphi 7 (BD7). Потребуется только выровнять вручную массивы на границу 128 бит и заменить использование 64-битных регистров общего назначения на 32-битные.
Upd 17.11.2015 BD7 не поддерживает SSE3.
movddup придется заменить на movsd + копирование в старшую часть xmm3;
haddpd — на перенос double из старшей половины xmm0 в младшую половину xmm1 и обычное сложение.[/Upd]

Если элементы массива не упорядочены, то на Sandy Bridge (SB) c медленной памятью:

Вариант 3 (sum3) оказывается быстрее варианта 1 для всех тестированных длин массива: от 2 до 100000. Для малых длин время sum3 составляет 0.8 и менее от времени sum1. Для больших длин — менее 0.2

Вариант 2 (sum2), как и ожидается, быстрее варианта 3 для очень коротких массивов (2, 3 элемента). Для размеров массива приблизительно от 100 и до 100 000 вариант 2 медленнее 3 чуть более чем в два раза.

Несмотря на то, что в варианте 2 присутствует дополнительное копирование, pand и всегда выполняется сложение — вариант 2 почти для всех размеров массива, если и медленнее, то совсем незначительно, а для массивов большой длины — составляет менее 0.5 от времени выполнения варианта 1.

Пример несколько искусственный, но он мне показал насколько вредно использовать условные переходы (даже на SB).

Используется синтаксис Delphi

 program Sum64;
uses sysutils;
{$APPTYPE CONSOLE}
const
 N  =  5000;
 Rep = 500000;

type
 TDArray = Array[1..N] of Double;
var
 DArray : TDArray;
 f: Double;

procedure Sum1;
{rdx - @DArray; rcx - N; xmm0 - sum(a), xmm1 - a,  xmm3 - f}
asm
 mov     rax,  Rep
@rep:
{Вычисление суммы}
 mov     rcx,  N
 pxor    xmm0, xmm0
 pxor    xmm1, xmm1
 pxor    xmm3, xmm3
 lea     rdx,  DArray
 movsd   xmm3, [f]
@next:
 movsd   xmm1, [rdx+ rcx*8-8]
 comisd  xmm1, xmm3
 jbe     @continue
 addsd   xmm0, xmm1
@continue:
 sub     rcx,  1
 jnz     @next
{Вычисление суммы end}
 sub     rax,  1
 jnz     @rep
end;

procedure Sum2;
{rdx - @DArray; rcx - N; xmm0 - sum(a), xmm1 - a, xmm2 - copy a, xmm3 - f}
asm
 mov     rax,  Rep
@rep:
{Вычисление суммы}
 mov     rcx,  N
 pxor    xmm0, xmm0
 pxor    xmm1, xmm1
 pxor    xmm2, xmm2
 pxor    xmm3, xmm3
 lea     rdx,  DArray
 movsd   xmm3, [f]
@next:
 movsd   xmm1, [rdx+ rcx*8-8]
 movaps  xmm2, xmm1
 cmppd   xmm1, xmm3,6          {Сравниваем с f?}
 pand    xmm1, xmm2            {Обнуляем если меньше f}
 addpd   xmm0, xmm1
 sub     rcx,  1
 jnz     @next
{Вычисление суммы end}
 sub     rax,  1
 jnz     @rep
end;

Procedure Sum3;
{rdx - @DArray; rcx - N; xmm0 - sum(a), xmm1 - a, xmm2 - copy a, xmm3 - f}
asm
 mov     rax,  Rep
@Rep:
{Вычисление суммы}
 lea     rdx,  DArray
 pxor    xmm0, xmm0
 pxor    xmm1, xmm1
 pxor    xmm2, xmm2
 mov     rcx,  N
 movddup xmm3, [f]             {Переносим с дублироваием f}
 btr     rcx,  0               {Кол-во элеменов массива нечётноё?}
 jnc     @next                 {Нет - переходим на суммирование упакованных}
 movsd   xmm1, [rdx+ rcx*8]    {Да - загружаем нечетный элемент}
 movaps  xmm2, xmm1            {Делаем копию для срвнения}
 cmpsd   xmm1, xmm3,6          {Сравниваем с f?}
 pand    xmm1, xmm2            {Обнуляем если меньше или равно f}
 addsd   xmm0, xmm1            {Накапливаем сумму}
 jrcxz   @res                  {Если больше элементов нет, то обходим цикл}
@next:
 movaps  xmm1, [rdx+ rcx*8-16] {Загружаем два очередных элемента}
 movaps  xmm2, xmm1            {Делаем копию для срвнения}
 cmppd   xmm1, xmm3,6          {Сравниваем с f?}
 pand    xmm1, xmm2            {Обнуляем если меньше f}
 addpd   xmm0, xmm1
 sub     rcx,  2
 jnz     @next
@res:
 haddpd  xmm0, xmm0
{Вычисление суммы end}
 sub     rax,  1
 jnz     @Rep
end;

var
 i:                        integer;
 Time1, Time2, d1, d2, d3: TDateTime;
 t: text;
begin
  f:= 0.5;
  randomize; for i:= 1 to N do DArray[i]:= random;

  Time1:= Time;
   Sum1;
  Time2:= Time;
  d1:= Time2-Time1;

  Time1:= Time;
   Sum2;
  Time2:= Time;
  d2:= Time2-Time1;

  Time1:= Time;
   Sum3;
  Time2:= Time;
  d3:= Time2-Time1;

  Assign(t, 'sumtst.txt');
{$I-}
  Append(t);
  if IOResult <> 0
   then Rewrite(t);
{$I+}
  Writeln(t, 'N= ', N, '   Rep=', Rep);
  if d3 > 0
   then writeln(t,  'd2/d1=', d2/d1:6:3, '  d3/d1=', d3/d1:6:3, '  d3/d2=', d3/d2:6:3)
   else writeln(t, 'd1=', d1, '   d2=', d2,', ',  '  d3=', d3);
  close(t);
end.

[bondkim137]

GAA в сообщении #1074098 писал(а):

Поскольку условные переходы и зависимость по данным замедляют скорость выполнения программы, то по возможности надо использовать второй способ (способ B)

Это не всегда так. Например, правильно предсказанный одиночный неосуществленный переход на многих конвеерах вообще не ест ни одного такта. Если данные не очень рандомные, то часто одновременный расчет по двум веткам и условное присвоение бывает неоправданным, несмотря на планарность кода.
Другое дело, что SSE, будучи SIMD, позволяет одновременно обрабатывать несколько значений - и тут уже не до выбора естественно.

Большое внимание заслуживает SSE4.2, (особенно целочисленный) - в нем появились операции блочного вычисления функции ошибки для сравнения двух паттернов - невероятный прирост производительности по сравнению с C++ реализацией, почти на уровне специализированных DSP-чипов

[GAA]

Да, это не всегда так. И на менее «рандомных» данных будут получены другие результаты. Поэтому и исходник программы тестирования вставлен в сообщение. Можно посмотреть, что будет при разных данных. Сформулировать правило лучшее, чем рекомендации Intel, вряд ли получится.

bondkim137, инструкции SSE мы видели, в крайнем случае, можем в документацию Intel заглянуть. Поучительный пример — в студию! [С исходником.]

Для любителей fpu.
Код, демонстрирующий близкие с SSE характерные особенности, связанные с увеличением длины массива: при малых длинах скорости выполнения обеих функций приблизительно одинаковы, при больших длинах скорость первой функции становится более чем в 3.5 раза выше скорости второй. Возрастание не монотонное. Думаю, для разных объёмов кеша третьего уровня зависимость отношения времен от длины массива будет разной.

Используется синтаксис Delphi

 program Sum32;
uses sysutils;
(* Тест двух вариантов функции вычисления суммы элементов массива,            *)
(*   больших некоторой величины (хранимой в переменной f)                     *)
(* В Sum2 используется условный переход.                                      *)
(* В Sum1 вместо условного перехода используется условное перемещение (fcmov) *) 

{$APPTYPE CONSOLE}
 const  N = 1000000;

 type   TDArray = Array[1..N] of Double;
 var    DArray : TDArray;
        f:       Double;
        Na:      LongWord;

 function Sum1: Double;
  asm
   mov      ecx, [Na]
   fld      [f]
   fldz
   fld      st(0)
 @Next:
   fld      qword ptr [DArray+ecx*8-8]
   fcomi    st(0), st(3)
   fcmovbe  st(0), st(2)
   faddp    st(1), st(0)
   sub      ecx, 1
   jnz      @Next
   fstp     st(1)
   fstp     st(1)
  end;

  function Sum2: Double;
  asm
   mov     ecx, [Na]
   fld     [f]
   fldz
   fld     st(0)
 @Next:
   fld     qword ptr [DArray+ecx*8-8]
   fcomi   st(0), st(3)
   jbe     @Continue
   fadd    st(1), st(0)
 @Continue:
   fstp    st(0)
   sub     ecx, 1
   jnz     @Next
   fstp    st(1)
   fstp    st(1)
  end;

const
 Rep = 10000;
var
 i:                    integer;
 res1, res2:           Double;
 Time1, Time2, d1, d2: TDateTime;
 t: text;
begin
  Na := N;
  f:= N div 2;
    f:= 0.5;
  randomize; for i:= 1 to N do DArray[i]:= random;
  Time1:= Time;
   for i:= 1 to Rep do res1:= Sum1;
  Time2:= Time;
  d1:= Time2-Time1;

  Time1:= Time;
   for i:= 1 to Rep do res2:= Sum2;
  Time2:= Time;
  d2:= Time2-Time1;

 Assign(t, 'sumtst.txt');
{$I-}
  Append(t);
  if IOResult <> 0
   then Rewrite(t);
{$I+}
  Writeln(t, 'N= ', N, '   Rep=', Rep);
  Writeln(t, 'd1=', d1/Rep/N*1E+15, '   d2=', d2/Rep/N*1E+15,'  d2/d1=', d2/d1);
  close(t);
end.

[bondkim137]

GAA в сообщении #1074389 писал(а):

bondkim137, инструкции SSE мы видели, в крайнем случае, можем в документацию Intel заглянуть. Поучительный пример — в студию! [С исходником.]

какого рода пример Вы хотите?

GAA в сообщении #1074389 писал(а):

Для любителей fpu...

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

for (int i = 0; i < n; i++)
  if (a[i] > 0) {
    a[i] = a[i]*b[i] + с[i];
  } else {
    a[i] = a[i] + d[i];
  }

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

a[i]*b[i] + c[i]

и

a[i] + d[i]

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

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

что касается конкретно функций FPU, то вроде как они больше не развиваются, и многие компиляторы умеют использовать SSE вместо FPU (хотя бы не в SIMD виде). выше я говорю об использовании условного присвоения vs условного перехода, не SSE vs FPU

[GAA]

bondkim137 в сообщении #1074400 писал(а):

какого рода пример Вы хотите?

Тот, который Вы считаете поучительным, и не близкий к разобранным в “Intel® 64 and IA-32 Architectures Optimization Reference Manual”, 2013, 2015.

-- Вт 17.11.2015 21:13:53 --

Текст на ассемблере. На ассемблере достаточно функций. Саму программу тестирования на Си, Delphi, Фортране, что Вам удобней. Каждый, думаю, для себя переведёт.

[bondkim137]

GAA в сообщении #1074414 писал(а):

Тот, который Вы считаете поучительным

Пример из моего предыдущего поста Вам не нравится?

[GAA]

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

Пример подсчёт суммы элементов массива (больших некоторой величины) был выбран по причине простоты. Но если погонять программу, то будут обнаружены, на мой взгляд, интересные особенности. И эти особенности не связаны непосредственно с SSE. Это должен был Вам продемонстрировать код с fpu.

Вообще говоря, я ожидал поучительный пример на тему SSE4.2.

Но если Вы сможете, то перед приведением примера с SSE4.2, приведите асм-текст для примера из сообщения post1074400.html#p1074400. Я с ходу не вижу ничего поучительного в этом примере. Вы меня заинтриговали. Тут действительно нужен текст Вашего цикла на асме с SSE (FPU) и комментариями, а также примерами данных, демонстрирующими. И давайте не спешить постить сообщения, а то ветка быстро замусориться.

Upd. Интересно, где такой цикл может быть использован? Есть очень отдалённое сходство с прямым ходом метода немонотонной прогонки. О методе немонотонной прогонке см. в
Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений, 1978.

В реализации прямого хода метода немонотонной прогонки есть что-то поучительное?

[bondkim137]

GAA в сообщении #1074439 писал(а):

Пример подсчёт суммы элементов массива (больших некоторой величины) был выбран по причине простоты. Но если погонять программу, то будут обнаружены, на мой взгляд, интересные особенности. И эти особенности не связаны непосредственно с SSE. Это должен был Вам продемонстрировать код с fpu.

Если Вы о сравнении Sum1 и Sum2 в вашем коде, то значительный прирост производительности, связан с отсутствием условного перехода. Условные
переходы обрабатывается предсказателем и в случае ошибок наказываются значительными задержками. Потому по возможности при обработке данных с
плохой предсказуемостью их следует избегать. Деградация будет нелинейная - там вообще сложная и малопонятная зависимость.

Потому да, условные присвоения (либо через AND и OR с маской) часто позволяют ускориться за счет отказа от условных переходов. Плюс ко всему,
при SIMD-обработке (packed, как вы его называли) это единственный способ строить условные конструкции - ведь операции выполняются один раз, одни

и те же для нескольких элементов сразу. Потому такие фокусы очень популярны при SIMD-обработке, со времен появления MMX.

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

int32 AbsBranchless(int32 x) {
  int32 mask = (x>>31);
  return (x^mask) - (mask);
}

Подобным образом раньше строили функцию ошибки между 8-битными grayscale-изображениями с помощью MMX. Потом MMX расширили (как раз в SSE) и
добавили туда как абсолютное значение, так и горизонтальную сумму, а также разрешили MMX операции над 128-ми разрядными регистрами SSE и пошло
веселье.

GAA в сообщении #1074439 писал(а):

Вообще говоря, я ожидал поучительный пример на тему SSE4.2.

А, понял. Чуть позже постараюсь отписать. Но смысл этого примера уже ясен наверное - реализация вычисления функции ошибки, как суммы абсолютных
значений попиксельной разницы двух блоков 16x16 на SSE4.2 работает почти в восемь раз быстрее, чем на MMX и в несколько десятков раз быстрее,
чем на голом языке.

GAA в сообщении #1074439 писал(а):

Но если Вы сможете, то перед приведением примера с SSE4.2, приведите асм-текст для примера из сообщения
post1074400.html#p1074400. Я сходу не вижу ничего поучительного в этом примере. Вы меня заинтриговали.

как-то так (мог опечататься, проверить сейчас негде):

  if (a[i] > 0) {
    a[i] = a[i]*b[i] + с[i];
  } else {
    a[i] = a[i] + d[i];
  }

 pxor    xmm4, xmm4
 shl      rcx, 2

@next:
 sub    rcx, 16
 movaps  xmm0, [rax + rcx] // a
 movaps  xmm1, [rbx + rcx] // b
 movaps  xmm2, [rsx + rcx] // c
 movaps  xmm3, [rdx + rcx] // d

 mulps   xmm1, xmm0
 addps   xmm1, xmm2     // a[i] = a[i]*b[i] + с[i]
 addps   xmm3, xmm0     // a[i] + d[i]

 cmpps   xmm0, xmm4, 6  // a[i] > 0
 pand    xmm1, xmm0     
 pandn    xmm0, xmm3
 por     xmm0, xmm1

 movaps  [rax + rcx], xmm0
 jnz     @next

GAA в сообщении #1074439 писал(а):

Upd. Интересно, где такой цикл может быть использован? Есть очень отдалённое сходство с прямым ходом метода немонотонной прогонки.

Недавно что-то похожее писалось, но я точно формулы не помню - так что пример выше в таком виде пока можно считать абстрактным.

[bondkim137]

bondkim137 в сообщении #1074464 писал(а):

  if (a[i] > 0) {
    a[i] = a[i]*b[i] + с[i];
  } else {
    a[i] = a[i] + d[i];
  }

Помимо фокусов с предсказателем переходов, тут еще присутсвуют условные обращения к памяти. В результате чего можно подобрать такие данные, что SSE-SIMD реализация будет работать медленее, чем на голом языке с условным переходом, несмотря на то, что за раз в первом случае обрабатывается 4 элемента и нет проблем с предсказателем! Возникнуть такое может при очень больших массивах отрицательных a[i] - тогда не будет происходить чтения из массивов b и c - и дополнительных кеш-промахов.

[Dmitriy40]

Не знаю насколько допустимо AVX2 (как расширении/продолжении SSE) в данной теме, но предлагаю ознакомиться с моим вопросом. Применены несколько специфичные именно для AVX моменты (впрочем, многие применимы и для SSE на относительно старых процессорах), как то:
1. Обнуление старшей половины регистра при записи в младшую.
2. Более редкое использование команд более медленных перестановок сквозь весь регистр в противовес перестановкам внутри половин регистров.
3. Более раняя перестановка отдельных чисел в более удобную позицию для дальнейших операций (вместо перестановки результатов).
4. Использование сдвига для выделения нужного бита, перенос его в нужную позицию для дальнейших вычислений и обнуление остальных битов регистра.
5. Использование совместно с предыдущим пунктом команд перестановок для обнуления неиспользуемой части регистра.
PS. Извините за отсутствие полного примера программы, но считаю его тривиальным для тех, кого это реально заинтересует.

[GAA]

bondkim137, я не знаю термин «условное присвоение» (но и не читаю ассемблер, поэтому за терминологией не слежу). Если знаете книгу, статью, приведите, пожалуйста, ссылку. Мне просто интересно.

bondkim137 в сообщении #1074364 писал(а):

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

bondkim137 в сообщении #1074400 писал(а):

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

a[i]*b[i] + c[i]

и

a[i] + d[i]

а затем использовать условное присвоение для выбора ответа.

Вот это я сначала не понял. [Я даже об одновременном исполнении двух веток на процессоре Itanium вспомнил. Эти ассоциации не так абсурдны, как кажется на первый взгляд. В новостях пробегало, что часть подходов, реализованных в Itanium, Intel планирует внедрить в x86.] На форуме встречаются люди, использующие разную терминологию, и кто-то, вполне возможно, неправильную или неточную. Поэтому во избежание переспросов и конфликтов крайне желательно сразу приводить исходный текст или, по крайней мере, его фрагменты.
По поводу Вашего примера. Я насочинял почти тот же текст (ну немного хуже), что и привели Вы. А вот попытки выполнить параллельное вычисление веток у меня не удались: накладные расходы становятся велики, пользы от такого параллельного вычисления нет. [Я играл с Double и пытался в старшей половине регистров вычислять одну ветку, а в младшей — другую.] Но они здесь и не уместны.

Я перечитал начальное сообщение темы. Написанное на скорую руку, оно плохо читаемо. То, что мне казалось очевидным, Вы уточнили и расширили.
В общем, спасибо за обсуждение!

По поводу поучительного примера на тему SSE4.2. Хорошо бы сразу подробно, и в начале новой подтемы (в теле сообщения) полужирным указать раздел, например [Strings], и заголовок самой подтемы.

Dmitriy40, я расширил тему (добавил в заголовке AVX). Но лучше обсуждать Вашу задачу в соответствующей ветке.

Данная ветка о приемах программирования и программистских трюках, а не о разных архитектурах. Обсуждение АVX вне контекста программирования выделено в отдельную ветку:
«Core X-series LGA-2066 (платформа Basin Falls)»

[Dmitriy40]

Может кому будет полезно. AVX. Модулярная арифметика.

Понадобилось тут много раз повторять $d=a+b \pmod c, \; a,b<c$ (условие критично!) для длинных массивов, разумеется прикрутил AVX, уместилось буквально в 4 команды:

Используется синтаксис ASM

vmovdqu         ymm0,[a]
vmovdqu         ymm1,[b]
vmovdqu         ymm2,[c]
                                ;Порты     Такты
vpaddb          ymm0,ymm0,ymm1  ;15     1/0.5   ymm0 - новые суммы
vpcmpgtb        ymm3,ymm2,ymm0  ;15     1/0.5   ymm3 - флаги превышения модулей
vpandn          ymm3,ymm3,ymm2  ;015    1/0.33  ymm3 - вычитаемое из сумм
vpsubb          ymm0,ymm0,ymm3  ;15     1/0.5   ymm0 - скорректированные суммы

vmovdqu         [d],ymm0

Мне было достаточно диапазона чисел $[0..63]$, потому обошёлся байтовыми операциями, зато одновременно обрабатываются 32 числа.
Загрузка и сохранение приведены лишь для полноты картины, в реальности обращений к памяти может и не быть, а работать только с регистрами.
Несмотря на полную зависимость по данным при выполнении в цикле итерации цикла разворачиваются процессором и реальная скорость составила именно 2 такта на 32 числа (их же кстати достаточно и для всех 4-х обращений в память вплоть до размера кэша L2, начиная с Haswell).

[Dmitriy40]

Не знаю полезно ли это кому, но пусть будет.
Понадобилось тут миллиарды раз выполнить цикл (с разными a)

Используется синтаксис C

uint64 offset[];//Большой массив
uint mask[128];//Флаги допустимости остатка от деления на простое число
uint p;//Простое число 3..127
uint128 a;//Начало диапазона проверки
for(i=0; i<20000000; i++) if(mask[(a + offset[i]) % p] != 0) {...}

т.е. для разных значений a проверить некий список смещений на допустимость остатков от деления на простое число. При решении в тупую основным тормозом будет операция целочисленного деления, в x86 режиме их будет 4 штуки на каждое число из массива offset[], остальные операции займут считанные такты и общее время будет порядка 100-130 тактов на каждое число.

Первое что напрашивается - преобразование (a+offset[i])%p в a%p+offset[i]%p (но надо проверять чтобы сумма не превысила p). Это позволяет выполнять не 4 деления, а лишь два (т.к. offset[] всего лишь uint64), т.е. ускоряет вдвое.
Второе что напрашивается - предрасчёт a%p и offsetmod[i]=offset[i]%p, т.к. это не меняется в цикле и является фактически константами. Это исключает вообще операции деления из цикла и оставляет лишь доступ к массиву mask[] и суммирование по модулю. Суммирование можно выполнить без условных переходов, а вот косвенный доступ к массиву сильно тормозит. Скорость работы составила примерно 3 такта на число.

И тут на сцену выходит AVX2 с его широкими регистрами и командами сдвига. Грузим в регистр сразу 32 числа offsetmod[] (они каждое влезают в байт) и пытаемся считать сразу их все. Превращаем массив mask[] в битовый, грузим его в регистр AVX и оперируем вообще без обращений к памяти (ну кроме конечно offsetmod[]). К сожалению в AVX2 нет команд сдвига 128 бит, максимум только 64, потому приходится сначала выделять нужную половинку маски в 64 бита, потом сдвигать, и это для каждого из 32-х байтов. Кроме того, нет и операций косвенного сдвига для счётчика сдвига в байте, потому приходится сдвиг 64-х битной маски повторять 8 раз чтобы обработать все 32 числа в регистре, да ещё и разбирать/собирать величины сдвигов и результаты. Итого, обработка 32-х чисел занимает примерно 25 тактов, т.е. даже менее такта на число. В этот момент конвейеры загружены примерно на 80% и казалось заметно это уже не ускорить.

Но потом вспомнил про супермощную команду косвенной перестановки байтов vpshufb в трёхрегистровом формате и оказалось что можно не делить 128 битную маску на кусочки и не делить 32 числа на группы, а выполнить сразу перестановку байтов маски по байтам чисел. Правда она переставляет байты, а нужно биты, но это обходится повтором этой команды для перестановки байтов и масок битов в байтах и трёх дополнительных команд vpand (две делят число на номер байта и номер бита в байте, третья объединяет результаты перестановок в один регистр). В итоге время обработки 32-х чисел сократилось до примерно 8 тактов, т.е. 4 числа на такт.

(Финальный исходник)

Используется синтаксис ASM

                                                        ;Порты     Такты
                vmovdqa         ymm7,[const1]           ;23     3/0.5   два повтора последовательности байтов 0x80,0x40,0x20,0x10,0x08,0x04,0x02,0x01
                vmovdqa         ymm6,[const2]           ;23     3/0.5   две 128 битных маски mask[]
                vmovdqa         ymm5,[a_p]              ;23     3/0.5   32 копии величины a%p
                vpxor           ymm4,ymm4,ymm4          ;015    1/0.33  32 байта 0x00
                vmovdqa         ymm3,[const3]           ;23     3/0.5   32 байта 0x78
                vmovdqa         ymm2,[const4]           ;23     3/0.5   32 байта 0x07
;Начало тела цикла
loop_32:        vpaddb          ymm0,ymm5,[ESI]         ;15     1/0.5   32 разных числа offsetmod[i]
                vpsubb          ymm1,ymm0,[p]           ;15     1/0.5   32 копии числа p
                vpminub         ymm0,ymm0,ymm1          ;15     1/0.5   32 суммы по модулю a+offsetmod[i] %p
                vpand           ymm1,ymm0,ymm3          ;015    1/0.33  выделяем номер байта маски
                vpand           ymm0,ymm0,ymm2          ;015    1/0.33  выделяем номер бита маски
                vpsrlw          ymm1,ymm1,3             ;0      1/1     преобразуем номер байта к диапазону 0..15
                vpshufb         ymm0,ymm7,ymm0          ;5      1/1     получаем 32 правильных масок битов
                vpshufb         ymm1,ymm6,ymm1          ;5      1/1     получаем 32 правильных байта из маски
                vpand           ymm0,ymm0,ymm1          ;015    1/0.33  перемножаем и получаем ноль или не ноль в зависимости от бита маски
                vpcmpeqb        ymm0,ymm0,ymm4          ;15     1/0.5   перемещаем бит в произвольной позиции байта в старший бит в байте
                vpmovmskb       EAX,ymm0                ;0      3/1     EAX - 32 бита для 32-х чисел
;Дальнейшая обработка битового массива
loop_bit:       bsf             EDX,EAX
                jz              next
                btr             EAX,EDX
                ;...
                jmp             loop_bit
;Конец тела цикла
next:           add             ESI,32
                sub             ECX,32
                jnz             loop_32

Вместо 8-ми повторов 4-х команд (vpshufb, vpsllvq, vpshufb, vpor) осталось лишь 7 команд (vpand, vpand, vpsrlw, vpshufb, vpshufb, vpand, vpcmpeqb).

Итого, число алгоритмические методы дали ускорение в 30-40 раз, применение AVX2 ускорило дополнительно ещё в 10-12 раз. Пришлось писать два десятка команд ассемблера.

[bondkim137]

(Оффтоп)

Элегантно. А на ARM с последними NEON-командами можно то же самое эффективно сделать?
Не знаю, может еще более новое что-то появилось в последних мобильных процессорах. Я давно оптимизациями на SIMD уже не занимался

[bondkim137]

GAA в сообщении #1074699 писал(а):

я не знаю термин «условное присвоение»

Извиняюсь за некропостинг. Только сейчас заметил, что упустил вопрос из виду. Если вопрос еще актуален, то это тоже самое, что упомянутое Вами условное перемещение.

GAA в сообщении #1074389 писал(а):

В Sum1 вместо условного перехода используется условное перемещение (fcmov)

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

bondkim137 в сообщении #1074464 писал(а):

cmpps xmm0, xmm4, 6 // a[i] > 0
pand xmm1, xmm0
pandn xmm0, xmm3
por xmm0, xmm1

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

-- 11.10.2017, 02:00 --

GAA в сообщении #1074699 писал(а):

По поводу поучительного примера на тему SSE4.2. Хорошо бы сразу подробно, и в начале новой подтемы (в теле сообщения) полужирным указать раздел, например [Strings], и заголовок самой подтемы

Постараюсь найти время и сделать полноценный пост. В двух словах - в видео-сжатии и анализе движения часто используется такой простой паттерн, как расчет функции ошибки для прямоугольного макроблока (прямоугольного куска изображения) для двух кадров. На вход функция помимо двух изображений и позиции макроблока берет двухмерное смещение (вектор) макроблока между двумя кадрами (собственно ищется лучший вектор).
На C расчет подобной ошибки означает вычисление суммы попиксельного модуля разницы (или квадрата разницы) между всеми элементами макроблока. На SSE4 можно считать всю горизонтальную строку двух макроблоков, например 16x16 или 8x8. За раз вычислить разницу между всеми элементами строки, за раз вычислить модуль или квадрат. За раз просуммировать все элементы строки с другой строкой. И за раз сложить элементы строки между собой, получив результат функции ошибки. Характерная сложность из квадратичной (например, 16x16~256) превращается в линейную (16), с учетом конвееров и стоимости всякой всячины, вроде расчета адресаци элемента при использовании C, производительность будет отличаться в десятки раз.
К условным переходам/присвоениям это отношения не имеет, просто речь о том, когда SIMD безоговорочно эффективнее C. Хотя возможно, я утратил нить беседы.