Как помочь gcc векторизовать код C

У меня есть следующий код. Первая часть просто читается в матрице комплексных чисел из стандарта в матрицу, называемую M Интересная часть – вторая часть.

 #include  #include  #include  #include  #include  int main() { int n, m, c, d; float re, im; scanf("%d %d", &n, &m); assert(n==m); complex float M[n][n]; for(c=0; c<n; c++) { for(d=0; d<n; d++) { scanf("%f%fi", &re, &im); M[c][d] = re + im * I; } } for(c=0; c<n; c++) { for(d=0; d<n; d++) { printf("%.2f%+.2fi ", creal(M[c][d]), cimag(M[c][d])); } printf("\n"); } /* Example:input 2 3 1+2i 2+3i 74-4i 3+4i 4+5i -7-8i */ /* Part 2. M is now an n by n matrix of complex numbers */ int s=1, i, j; int *f = malloc(n * sizeof *f); complex float *delta = malloc(n * sizeof *delta); complex float *v = malloc(n * sizeof *v); complex float p = 1, prod; for (i = 0; i < n; i++) { v[i] = 0; for (j = 0; j <n; j++) { v[i] += M[j][i]; } p *= v[i]; f[i] = i; delta[i] = 1; } j = 0; while (j < n-1) { prod = 1.; for (i = 0; i < n; i++) { v[i] -= 2.*delta[j]*M[j][i]; prod *= v[i]; } delta[j] = -delta[j]; s = -s; p += s*prod; f[0] = 0; f[j] = f[j+1]; f[j+1] = j+1; j = f[0]; } free(delta); free(f); free(v); printf("%f + i%f\n", creal(p/pow(2.,(n-1))), cimag(p/pow(2.,(n-1)))); return 0; } 

Я компилирую с gcc -fopt-info-vec-all -O3 -ffast-math -march=bdver2 permanent-in-cc -lm . Это объясняет мне, почему почти ни одна из циклов не подвергается векторизации.

Важнейшей частью производительности являются линии 47–50, которые:

 for (i = 0; i < n; i++) { v[i] -= 2.*delta[j]*M[j][i]; prod *= v[i]; } 

gcc говорит мне:

 permanent-in-cc:47:7: note: reduction used in loop. permanent-in-cc:47:7: note: Unknown def-use cycle pattern. permanent-in-cc:47:7: note: reduction used in loop. permanent-in-cc:47:7: note: Unknown def-use cycle pattern. permanent-in-cc:47:7: note: Unsupported pattern. permanent-in-cc:47:7: note: not vectorized: unsupported use in stmt. permanent-in-cc:47:7: note: unexpected pattern. [...] permanent-in-cc:48:26: note: SLP: step doesn't divide the vector-size. permanent-in-cc:48:26: note: Unknown alignment for access: IMAGPART_EXPR  permanent-in-cc:48:26: note: SLP: step doesn't divide the vector-size. permanent-in-cc:48:26: note: Unknown alignment for access: REALPART_EXPR  [...] permanent-in-cc:48:26: note: Build SLP failed: unrolling required in basic block SLP permanent-in-cc:48:26: note: Failed to SLP the basic block. permanent-in-cc:48:26: note: not vectorized: failed to find SLP opportunities in basic block. 

Как я могу исправить проблемы, которые останавливают эту часть от векторизации?


Любопытно, что эта часть векторизована, но я не уверен, почему:

 for (j = 0; j <n; j++) { v[i] += M[j][i]; 

Полный вывод gcc -fopt-info-vec-all -O3 -ffast-math -march = bdver2 постоянный-в-cc -lm находится в https://bpaste.net/show/18ebc3d66a53 .

    Думаю, я мог бы это понять. После большого количества проб / ошибок стало ясно, что gcc, встроенный в оптимизацию векторизации, является своего рода жестко закодированным и не «понимает» сложные числа должным образом. Я внесла некоторые изменения в код и получил ваш внутренний чувствительный к производительности цикл для векторизации, подтвержденный выходом gcc (хотя я не уверен, что желаемый результат будет эквивалентен вычислительному эквиваленту тому, что вы хотите). Хотя мое понимание ограничено тем, что вы хотите, чтобы код делал, вывод состоит в том, что он будет работать нормально, если вы будете вычислять реальный и образ отдельно. Посмотри:

      float t_r = 0.0, t_im = 0.0; // two new temporaries while (j < n-1) { prod = 1.; for (i = 0; i < n; i++) { // fill the temps after subtraction from V to avoid stmt error t_r = creal (v[i]) - (2. * creal(delta[j]) * creal (M[j][i])); t_im = cimag(v[i]) - (2. * cimag(delta[j]) * cimag (M[j][i])) * I; //v[i] = 2.*delta[j]*M[j][i]; v[i] = t_r + t_im; // sum of real and img prod *= v[i]; } delta[j] = -delta[j]; s = -s; p += s*prod; f[0] = 0; f[j] = f[j+1]; f[j+1] = j+1; j = f[0]; } 

    Рассмотрим сначала код. У нас есть

     complex float M[rows][cols]; complex float v[cols]; float delta[rows]; complex float p = 1.0f; float s = 1.0f; 

    Хотя delta[] имеет complex float тип complex float в коде OP, он всегда содержит только -1.0f или +1.0f . (Более того, вычисления можно было бы упростить, если бы они были -2.0f или +2.0f .) По этой причине я определил как реальный, а не комплексный.

    Аналогично, OP определяет s как int , но эффективно использует его как -1.0f и +1.0f (только в вычислениях). Вот почему я определил его явно как float .

    Я опускаю массив f , потому что есть тривиальный способ полностью избежать этого.

    Первый цикл интересной части кода,

      for (i = 0; i < n; i++) { v[i] = 0; for (j = 0; j  

    выполняет несколько вещей. Он инициализирует все элементы массива delta[] равными 1; он может (и, вероятно, должен) быть разделен на отдельный цикл.

    Поскольку внешний цикл возрастает в i , p будет произведением элементов из v ; он также может быть разделен на отдельный цикл.

    Поскольку внутренний цикл суммирует все элементы в столбце i с v[i] , внешняя и внутренняя петли просто суммируют каждую строку как вектор вектору v .

    Таким образом, мы можем переписать выше, в псевдокоде, так как

     Copy first row of matrix M to vector v For r = 1 .. rows-1: Add complex values in row r of matrix M to vector v p = product of complex elements in vector v delta = 1.0f, 1.0f, 1.0f, .., 1.0f, 1.0f 

    Давайте посмотрим на второй вложенный цикл:

      j = 0; while (j < n-1) { prod = 1.; for (i = 0; i < n; i++) { v[i] -= 2.*delta[j]*M[j][i]; prod *= v[i]; } delta[j] = -delta[j]; s = -s; p += s*prod; f[0] = 0; f[j] = f[j+1]; f[j+1] = j+1; j = f[0]; } 

    Трудно видеть, если вы не изучите значения j мере продвижения цикла, но последние 4 строки в теле внешнего цикла реализуют целую последовательность OEIS A007814 в j (0,1,0,2,0,1, 0,3,0,1,0,2,0,1,0,4, ...). Число итераций в этом цикле - 2 строки-1 -1. Эта часть последовательности симметрична и реализует бинарное дерево высот rows-1:

      4 3 3 2 2 2 2 (Read horizontally) 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 

    Оказывается, что если мы перебираем i = 1 .. 2 строки-1 , то r - это число нулевых младших бит в i . GCC предоставляет встроенную функцию __builtin_ctz() , которая вычисляет именно это. (Обратите внимание, что __builtin_ctz(0) дает неопределенное значение, поэтому не делайте этого, даже если оно вызывает определенное значение на вашем компьютере .)

    Внутренний цикл вычитает комплексные значения в строке j матрицы, масштабированной на 2*delta[j] , из вектора v[] . Он также вычисляет произведение комплексных записей в векторе v[] (после вычитания) на переменную prod .

    После внутреннего цикла delta[j] отрицается, как и масштабный коэффициент s . Значение переменной prod , масштабированное по s , добавляется к p .

    После цикла конечный (сложный) результат равен p разделенному на 2 строки-1 . Это лучше сделать с использованием функции ldexp() C99 (отдельно на реальных и сложных частях).

    Поэтому мы можем переписать второй вложенный цикл в псевдокоде, так как

     s = 1.0f For k = 1 .. rows-1, inclusive: r = __builtin_ctz(k), ie number of least significant bits that are zero in k Subtract the complex values on row r of matrix M, scaled by delta[r], from vector v[] prod = the product of (complex) elements in vector v[] Negate scale factor s (changing its sign) Add prod times s to result p 

    По моему опыту, лучше разделить действительную и мнимую части на отдельные векторы и матрицы. Рассмотрим следующие определения:

     typedef struct { float *real; float *imag; size_t floats_per_row; /* Often 'stride' */ size_t rows; size_t cols; } complex_float_matrix; /* Set an array of floats to a specific value */ void float_set(float *, float, size_t); /* Copy an array of floats */ void float_copy(float *, const float *, size_t); /* malloc() vector-aligned memory; size in floats */ float *float_malloc(size_t); /* Elementwise addition of floats */ void float_add(float *, const float *, size_t); /* Elementwise addition of floats, scaled by a real scale factor */ void float_add_scaled(float *, const float *, float, size_t); /* Complex float product, separate real and imag arrays */ complex float complex_float_product(const float *, const float *, size_t); 

    Все вышесказанное легко векторизовать, пока float_malloc() дает достаточно выровненные указатели (и компилятору сообщается, что через, например, атрибут функции GCC __attribute__ ((__assume_aligned__ (BYTES_IN_VECTOR))); ) и член floats_per_row в матрице является кратным количеству поплавков в векторе.

    (Я не знаю, может ли GCC автоматически векторизовать вышеуказанные функции, но я знаю, что они могут быть векторизованы «вручную» с использованием расширений вектора GCC.)

    При этом выше вся интересная часть кода в псевдо-C становится

     complex float permanent(const complex_float_matrix *m) { float *v_real, *v_imag; float *scale; /* OP used 'delta' */ complex float result; /* OP used 'p' */ complex float product; /* OP used 'prod' */ float coeff = 1.0f; /* OP used 's' */ size_t i = 1 << (m->rows - 1); size_t r; if (!m || m->cols < 1 || m->rows < 1 || !i) { /* TODO: No input matrix, or too many rows in input matrix */ } v_real = float_malloc(m->cols); v_imag = float_malloc(m->cols); scale = float_malloc(m->rows - 1); if (!v_real || !v_imag || !scale) { free(scale); free(v_imag); free(v_real); /* TODO: Out of memory error */ } float_set(scale, 2.0f, m->rows - 1); /* Sum matrix rows to v. */ float_copy(v_real, m->real, m->cols); for (r = 1; r < m->rows; r++) float_add(v_real, m->real + r * m->floats_per_row, m->cols); float_copy(v_imag, m->imag, m->cols); for (r = 1; r < m->rows; r++) float_add(v_imag, m->imag + r * m->floats_per_row, m->cols); result = complex_float_product(v_real, v_imag, m->cols); while (--i) { r = __builtin_ctz(i); scale[r] = -scale[r]; float_add_scaled(v_real, m->real + r * m->floats_per_row, m->cols); float_add_scaled(v_imag, m->imag + r * m->floats_per_row, m->cols); product = complex_float_product(v_real, v_imag, m->cols); coeff = -coeff; result += coeff * product; } free(scale); free(v_imag); free(v_real); return result; } 

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

    Затем я рассмотрю вывод сборки GCC ( -S ), чтобы убедиться, что он может достаточно векторизовать отдельные операции (функции, перечисленные ранее).

    Ручная векторизация функций с использованием расширений вектора GCC довольно проста. Объявление вектора float тривиально:

     typedef float vec2f __attribute__((vector_size (8), aligned (8))); /* 64 bits; MMX, 3DNow! */ typedef float vec4f __attribute__((vector_size (16), aligned (16))); /* 128 bits; SSE */ typedef float vec8f __attribute__((vector_size (32), aligned (32))); /* 256 bits; AVX, AVX2 */ typedef float vec16f __attribute__((vector_size (64), aligned (64))); /* 512 bits; AVX512F */ 

    Отдельные компоненты в каждом векторе могут быть решены с использованием нотации массива ( v[0] и v[1] для vec2f v; ). GCC может выполнять базовые операции по целым векторам по элементам; нам просто нужно добавить и умножить здесь. Следует избегать горизонтальных операций (операций, которые применяются между элементами в одном и том же векторе), и вместо этого элементы переупорядочиваются.

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

    Память, выделенная для вектора, должна быть достаточно выровнена. malloc() не обеспечивает достаточную согласованную память во всех случаях; вы должны использовать posix_memalign() . Выравниваемый атрибут может использоваться для увеличения выравнивания GCC для векторного типа при распределении локально или статически. В матрице вам нужно убедиться, что строки начинаются с достаточно выровненной границы; вот почему у меня есть переменная floats_per_row в структуре.

    В тех случаях, когда число элементов в векторе (или строке) велико, но не кратно количеству поплавков в векторе, вы должны наложить вектор на «инертные» значения - значения, которые не влияют на результат, как 0.0f для сложения и вычитания, и 1.0f для умножения.

    По крайней мере, на x86 и x86-64, GCC будет генерировать лучший код для циклов с использованием указателей. Например, это

     void float_set(float *array, float value, size_t count) { float *const limit = array + count; while (array < limit) *(array++) = value; } 

    дает лучший код, чем

     void float_set(float *array, float value, size_t count) { size_t i; for (i = 0; i < count; i++) array[i] = value; } 

    или же

     void float_set(float *array, float value, size_t count) { while (count--) *(array++) = value; } 

    (которые, как правило, производят аналогичный код). Лично я реализую его как

     void float_set(float *array, float value, size_t count) { if (!((uintptr_t)array & 7) && !(count & 1)) { uint64_t *const end = (uint64_t *)__builtin_assume_aligned((void *)(array + count), 8); uint64_t *ptr = (uint64_t *)__builtin_assume_aligned((void *)array, 8); uint64_t val; __builtin_memcpy(&val, &value, 4); __builtin_memcpy(4 + (char *)&val, &value, 4); while (ptr < end) *(ptr++) = val; } else { uint32_t *const end = (uint32_t *)__builtin_assume_aligned((void *)(array + count), 4); uint32_t *ptr = (uint32_t *)__builtin_assume_aligned((void *)array, 4); uint32_t val; __builtin_memcpy(&val, &value, 4); while (ptr < end) *(ptr++) = val; } } 

    и float_copy() как

     void float_copy(float *target, const float *source, size_t count) { if (!((uintptr_t)source & 7) && !((uintptr_t)target & 7) && !(count & 1)) { uint64_t *const end = (uint64_t *)__builtin_assume_aligned((void *)(array + count), 8); uint64_t *ptr = (uint64_t *)__builtin_assume_aligned((void *)target, 8); uint64_t *src = (uint64_t *)__builtin_assume_aligned((void *)source, 8); while (ptr < end) *(ptr++) = *(src++); } else { uint32_t *const end = (uint32_t *)__builtin_assume_aligned((void *)(array + count), 4); uint32_t *ptr = (uint32_t *)__builtin_assume_aligned((void *)array, 4); uint32_t *src = (uint32_t *)__builtin_assume_aligned((void *)source, 4); while (ptr < end) *(ptr++) = *(src++); } } 

    или что-то близкое к этому.

    Наиболее сложным для векторизации является complex_float_product() . Если вы заполняете неиспользуемые элементы в конечном векторе с 1.0f для реальной части и 0.0f для мнимой части, вы можете легко вычислить сложный продукт для каждого вектора. Помните, что

    ( a + b i) × ( c + d i) = ( a c - b d ) + ( a d + b c ) i

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

    (Короче говоря, «твердая» часть состоит в том, чтобы найти способ переупорядочить элементы, чтобы максимально использовать упакованное векторное умножение, и не нужно так много перемешать / перемещать, чтобы замедлить его.)

    Для функции float_add_scaled() вы должны создать вектор, заполненный масштабным коэффициентом; что-то вроде следующего,

     void float_add_scaled(float *array, const float *source, float scale, size_t count) { const vec4f coeff = { scale, scale, scale, scale }; vec4f *ptr = (vec4f *)__builtin_assume_aligned((void *)array, 16); vec4f *const end = (vec4f *)__builtin_assume_aligned((void *)(array + count), 16); const vec4f *src = (vec4f *)__builtin_assume_aligned((void *)source, 16); while (ptr < end) *(ptr++) += *(src++) * coeff; } 

    если мы игнорируем проверки выравнивания и размера и резервную реализацию.

    Журналы оптимизатора четко указывают

    Неизвестное выравнивание для доступа: …

    при попытке векторизации

     printf("%.2f%+.2fi ", creal(M[c][d]), cimag(M[c][d])); //24 v[i] += M[j][i]; //38 p *= v[i]; //40 v[i] -= 2.*delta[j]*M[j][i]; //48 

    Кажется, действительно связано, что вам нужно принудительно выравнивать ваши массивы M , delta и v в памяти.

    Авто-векторизация в GCC

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

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

     complex float * restrict delta; posix_memalign(&delta, 64, n * sizeof *delta); //to adapt 

    Какова ваша целевая среда? (ОС, ЦП)

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