Производить циклы без команды cmp в GCC

У меня есть ряд жестких циклов, которые я пытаюсь оптимизировать с помощью GCC и intrinsics. Рассмотрим, например, следующую функцию.

void triad(float *x, float *y, float *z, const int n) {
    float k = 3.14159f;
    int i;
    __m256 k4 = _mm256_set1_ps(k);
    for(i=0; i<n; i+=8) {
        _mm256_store_ps(&z[i], _mm256_add_ps(_mm256_load_ps(&x[i]), _mm256_mul_ps(k4, _mm256_load_ps(&y[i]))));
    }
}

Это создает основной цикл, подобный этому

20: vmulps ymm0,ymm1,[rsi+rax*1]
25: vaddps ymm0,ymm0,[rdi+rax*1]
2a: vmovaps [rdx+rax*1],ymm0
2f: add    rax,0x20
33: cmp    rax,rcx
36: jne    20 

Но инструкция cmp не нужна. Вместо того, чтобы rax начинать с нуля и заканчивать в sizeof(float)*n, мы можем установить базовые указатели (rsi, rdi и rdx) в конец массива и установить rax в -sizeof(float)*n а затем проверить на нуль. Я могу сделать это с помощью моего собственного кода сборки, подобного этому

.L2  vmulps          ymm1, ymm2, [rdi+rax]
     vaddps          ymm0, ymm1, [rsi+rax]
     vmovaps         [rdx+rax], ymm0
     add             rax, 32
     jne             .L2

но я не могу заставить GCC сделать это. Сейчас у меня несколько тестов, где это имеет существенное значение. До недавнего времени GCC и intrinsics меня очень сильно ранили, поэтому мне интересно, есть ли компилятор или способ изменить порядок/изменить код, поэтому команда cmp не создается с помощью GCC.

Я попробовал следующее, но он все еще производит cmp. Все варианты, которые я попытался, все еще создают cmp.

void triad2(float *x, float *y, float *z, const int n) {
    float k = 3.14159f;
    float *x2 = x+n;
    float *y2 = y+n;
    float *z2 = z+n;    
    int i;
    __m256 k4 = _mm256_set1_ps(k);
    for(i=-n; i<0; i+=8) {
        _mm256_store_ps(&z2[i], _mm256_add_ps(_mm256_load_ps(&x2[i]), _mm256_mul_ps(k4, _mm256_load_ps(&y2[i]))));
    }
}

Изменить: Меня интересует максимизация уровня команд parallelism (ILP) для этих функций для массивов, которые вписываются в кеш L1 (фактически для n=2048). Хотя разворачивание может использоваться для улучшения полосы пропускания, оно может уменьшить ILP (при условии, что полная пропускная способность может быть достигнута без разворота).

Изменить: Вот таблица результатов для Core2 (pre Nehalem), IvyBridge и системы Haswell. Intrinsics - это результат использования intrinsics, unroll1 - мой код сборки, не использующий cmp, а unroll16 - мой код сборки, разворачивающий 16 раз. Процентное соотношение - это процентная доля максимальной производительности (частота * num_bytes_cycle, где num_bytes_cycle - 24 для SSE, 48 для AVX и 96 для FMA).

                 SSE         AVX         FMA
intrinsic      71.3%       90.9%       53.6%      
unroll1        97.0%       96.1%       63.5%
unroll16       98.6%       90.4%       93.6%
ScottD         96.5%
32B code align             95.5%

Для SSE я получаю почти такой же результат без разворачивания, как при разворачивании, но только если я не использую cmp. В AVX я получаю лучший результат без разворачивания и без использования cmp. Интересно, что на IB разворачивание на самом деле хуже. На Хасуэле я получаю лучший результат при разворачивании. Вот почему я спросил об этом question. Исходный код для проверки этого может быть найден в этом вопросе.

Edit:

Основываясь на ответе ScottD, теперь я получаю почти 97% с внутренними функциями моей системы Core2 (64-разрядный режим до Nehalem). Я не уверен, почему вопрос cmp имеет значение, поскольку он должен так или иначе, за два такта. Для Sandy Bridge получается, что потеря эффективности связана с выравниванием кода не с дополнительным cmp. В любом случае Haswell работает только разворачивание.

Ответы

Ответ 1

Как насчет этого. Компилятор - gcc 4.9.0 mingw x64:

void triad(float *x, float *y, float *z, const int n) {
    float k = 3.14159f;
    intptr_t i;
    __m256 k4 = _mm256_set1_ps(k);

    for(i = -n; i < 0; i += 8) {
        _mm256_store_ps(&z[i+n], _mm256_add_ps(_mm256_load_ps(&x[i+n]), _mm256_mul_ps(k4, _mm256_load_ps(&y[i+n]))));
    }
}

gcc -c-O3 -march = corei7 -mavx2 triad.c

0000000000000000 <triad>:
   0:   44 89 c8                mov    eax,r9d
   3:   f7 d8                   neg    eax
   5:   48 98                   cdqe
   7:   48 85 c0                test   rax,rax
   a:   79 31                   jns    3d <triad+0x3d>
   c:   c5 fc 28 0d 00 00 00 00 vmovaps ymm1,YMMWORD PTR [rip+0x0]
  14:   4d 63 c9                movsxd r9,r9d
  17:   49 c1 e1 02             shl    r9,0x2
  1b:   4c 01 ca                add    rdx,r9
  1e:   4c 01 c9                add    rcx,r9
  21:   4d 01 c8                add    r8,r9

  24:   c5 f4 59 04 82          vmulps ymm0,ymm1,YMMWORD PTR [rdx+rax*4]
  29:   c5 fc 58 04 81          vaddps ymm0,ymm0,YMMWORD PTR [rcx+rax*4]
  2e:   c4 c1 7c 29 04 80       vmovaps YMMWORD PTR [r8+rax*4],ymm0
  34:   48 83 c0 08             add    rax,0x8
  38:   78 ea                   js     24 <triad+0x24>

  3a:   c5 f8 77                vzeroupper
  3d:   c3                      ret

Как и ваш ручной код, gcc использует 5 инструкций для цикла. В коде gcc используется масштаб = 4, в котором используется scale = 1. Мне удалось заставить gcc использовать scale = 1 с циклом 5 команд, но код C неудобен, а 2 из инструкций AVX в цикле растут от 5 до 6 байтов.

Ответ 2

Декодер команд на мосту Intel Ivy или позже может слить cmp и jne в одну операцию в конвейере (так называемый макро-op-fusion), поэтому на этих последних процессорах cmp все равно должен исчезнуть.

Ответ 3

Конечный код:

#define SF sizeof(float)
#ifndef NO                   //floats per vector, compile with -DNO = 1,2,4,8,...
#define NO 8                 //MUST be power of two
#endif

void triadfinaler(float const *restrict x, float const *restrict y,   \
                  float *restrict z, size_t n)
{
  float *restrict d = __builtin_assume_aligned(z, NO*SF);       //gcc builtin,
  float const *restrict m = __builtin_assume_aligned(y, NO*SF); //optional but produces
  float const *restrict a = __builtin_assume_aligned(x, NO*SF); //better code
  float const k = 3.14159f;
  n*=SF;
  while (n &= ~((size_t)(NO*SF)-1))    //this is why NO*SF must be power of two
    {
      size_t nl = n/SF;
      for (size_t i = 0; i<NO; i++)
        {
          d[nl-NO+i] = k * m[nl-NO+i] + a[nl-NO+i];
        }
      n -= (NO*SF);
    }
}

Я предпочитаю, чтобы компилятор выбирал инструкции, а не использовал intrinsics (не в последнюю очередь потому, что вы использовали intel-intrinsics, что gcc не очень нравится). В любом случае, следующий код создает хорошую сборку для меня на gcc 4.8:

void triad(float *restrict x, float *restrict y, float *restrict z, size_t n)
//I hope you weren't aliasing any function arguments... Oh, an it void, not float
{
  float *restrict d = __builtin_assume_aligned(z, 32);  // Uh, make sure your arrays
  float *restrict m = __builtin_assume_aligned(y, 32);  // are aligned? Faster that way
  float *restrict a = __builtin_assume_aligned(x, 32);  //
  float const k = 3.14159f;
  while (n &= ~((size_t)0x7))       //black magic, causes gcc to omit code for non-multiples of 8 floats
    {
      n -= 8;                       //You were always computing on 8 floats at a time, right?
      d[n+0] = k * m[n+0] + a[n+0]; //manual unrolling
      d[n+1] = k * m[n+1] + a[n+1];
      d[n+2] = k * m[n+2] + a[n+2];
      d[n+3] = k * m[n+3] + a[n+3];
      d[n+4] = k * m[n+4] + a[n+4];
      d[n+5] = k * m[n+5] + a[n+5];
      d[n+6] = k * m[n+6] + a[n+6];
      d[n+7] = k * m[n+7] + a[n+7];
    }
}

Это создает хороший код для моего corei7avx2, с -O3:

triad:
    andq    $-8, %rcx
    je  .L8
    vmovaps .LC0(%rip), %ymm1

.L4:
    subq    $8, %rcx
    vmovaps (%rsi,%rcx,4), %ymm0
    vfmadd213ps (%rdi,%rcx,4), %ymm1, %ymm0
    vmovaps %ymm0, (%rdx,%rcx,4)
    andq    $-8, %rcx
    jne .L4
    vzeroupper
.L8:
    rep ret
    .cfi_endproc

.LC0:
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000

Изменить: Я был немного разочарован тем, что компилятор не оптимизировал этот код до последней инструкции, поэтому я немного перепутал его. Простое изменение порядка вещей в цикле избавилось от AND, выпущенного компилятором, что привело меня к правильному пути. Затем я должен был заставить его не делать ненужного вычисления адресов в цикле. Вздох.

void triadtwo(float *restrict x, float *restrict y, float *restrict z, size_t n)
{
  float *restrict d = __builtin_assume_aligned(z, 32);
  float *restrict m = __builtin_assume_aligned(y, 32);
  float *restrict a = __builtin_assume_aligned(x, 32);
  float const k = 3.14159f;
  n<<=2;
  while (n &= -32)
    {
      d[(n>>2)-8] = k * m[(n>>2)-8] + a[(n>>2)-8];
      d[(n>>2)-7] = k * m[(n>>2)-7] + a[(n>>2)-7];
      d[(n>>2)-6] = k * m[(n>>2)-6] + a[(n>>2)-6];
      d[(n>>2)-5] = k * m[(n>>2)-5] + a[(n>>2)-5];
      d[(n>>2)-4] = k * m[(n>>2)-4] + a[(n>>2)-4];
      d[(n>>2)-3] = k * m[(n>>2)-3] + a[(n>>2)-3];
      d[(n>>2)-2] = k * m[(n>>2)-2] + a[(n>>2)-2];
      d[(n>>2)-1] = k * m[(n>>2)-1] + a[(n>>2)-1];
      n -= 32;
    }
}

Уродливый код? Да. Но сборка:

triadtwo:
    salq    $2, %rcx
    andq    $-32, %rcx
    je  .L54
    vmovaps .LC0(%rip), %ymm1

.L50:
    vmovaps -32(%rsi,%rcx), %ymm0
    vfmadd213ps -32(%rdi,%rcx), %ymm1, %ymm0
    vmovaps %ymm0, -32(%rdx,%rcx)
    subq    $32, %rcx
    jne .L50
    vzeroupper
.L54:
    rep ret
    .cfi_endproc
.LC0:
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000

Mmmmhhh, славные пять инструкций в цикле, макро-op с возможностью вычитания и ветвления...