Развертывание цикла для достижения максимальной пропускной способности с Ivy Bridge и Haswell
Я рассчитываю восемь точечных продуктов одновременно с AVX. В моем текущем коде я делаю что-то вроде этого (перед развертыванием):
Ivy-Bridge / Sandy-Bridge
__m256 areg0 = _mm256_set1_ps(a[m]);
for(int i=0; i<n; i++) {
__m256 breg0 = _mm256_load_ps(&b[8*i]);
tmp0 = _mm256_add_ps(_mm256_mul_ps(arge0,breg0), tmp0);
}
Haswell
__m256 areg0 = _mm256_set1_ps(a[m]);
for(int i=0; i<n; i++) {
__m256 breg0 = _mm256_load_ps(&b[8*i]);
tmp0 = _mm256_fmadd_ps(arge0, breg0, tmp0);
}
Сколько раз мне нужно развернуть цикл для каждого случая, чтобы обеспечить максимальную пропускную способность?
Для Haswell, использующего FMA3, я думаю, что ответ здесь - FLOPS за цикл для Sandy-Bridge и haswell SSE2 / AVX / AVX2. Мне нужно развернуть петлю 10 раз.
Для Ivy Bridge я думаю, что это 8. Вот моя логика. Добавление AVX имеет задержку 3, а умножение - задержку 5. Ivy Bridge может выполнять одно умножение AVX и одно сложение AVX одновременно, используя разные порты. Используя обозначения m для умножения, a для сложения и x для отсутствия операций, а также число для обозначения частичной суммы (например, m5 означает умножение с 5-й частичной суммой), я могу написать:
port0: m1 m2 m3 m4 m5 m6 m7 m8 m1 m2 m3 m4 m5 ...
port1: x x x x x a1 a2 a3 a4 a5 a6 a7 a8 ...
Таким образом, используя 8 частичных сумм после девяти тактов (четыре от нагрузки и пять от умножения), я могу отправлять одну загрузку AVX, одно сложение AVX и одно умножение AVX за каждый такт.
Я предполагаю, что это означает, что невозможно достичь максимальной пропускной способности для этих задач в 32-битном режиме в Ivy Bridge и Haswell, поскольку 32-битный режим имеет только восемь регистров AVX?
Изменить: Что касается щедрости. Мои основные вопросы остаются в силе. Я хочу получить максимальную пропускную способность функций Ivy Bridge или Haswell выше, n может быть любым значением, большим или равным 64. Я думаю, что это может быть сделано только с использованием разворачивания (восемь раз для Ivy Bridge и 10 раз для Haswell). Если вы думаете, что это можно сделать другим методом, давайте посмотрим. В некотором смысле это вариация того, как мне достичь теоретического максимума 4 FLOP за цикл?, Но вместо только умножения и сложения я ищу одну 256-битную загрузку (или две 128-битные загрузки), одно AVX-умножение и одно AVX-добавление каждый такт с Ivy Bridge или две 256-битные загрузки и две инструкции FMA3 за такт.
Я также хотел бы знать, сколько регистров необходимо. Для Ivy Bridge я думаю, что это 10. Один для трансляции, один для загрузки (только один из-за переименования регистров) и восемь для восьми неполных сумм. Поэтому я не думаю, что это можно сделать в 32-битном режиме (и действительно, когда я работаю в 32-битном режиме, производительность значительно падает).
Следует отметить, что компилятор может давать неверные результаты. Разница в производительности между MSVC и GCC для высокооптимизированного кода умножения матриц.
Текущая функция, которую я использую для Ivy Bridge, ниже. Это в основном умножает одну строку матрицы 64x64 a со всей матрицей 64x64 b (Я запускаю эту функцию 64 раза в каждом ряду a чтобы получить полную матрицу, умноженную на матрицу c).
#include <immintrin.h>
extern "C" void row_m64x64(const float *a, const float *b, float *c) {
const int vec_size = 8;
const int n = 64;
__m256 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
tmp0 = _mm256_loadu_ps(&c[0*vec_size]);
tmp1 = _mm256_loadu_ps(&c[1*vec_size]);
tmp2 = _mm256_loadu_ps(&c[2*vec_size]);
tmp3 = _mm256_loadu_ps(&c[3*vec_size]);
tmp4 = _mm256_loadu_ps(&c[4*vec_size]);
tmp5 = _mm256_loadu_ps(&c[5*vec_size]);
tmp6 = _mm256_loadu_ps(&c[6*vec_size]);
tmp7 = _mm256_loadu_ps(&c[7*vec_size]);
for(int i=0; i<n; i++) {
__m256 areg0 = _mm256_set1_ps(a[i]);
__m256 breg0 = _mm256_loadu_ps(&b[vec_size*(8*i + 0)]);
tmp0 = _mm256_add_ps(_mm256_mul_ps(areg0,breg0), tmp0);
__m256 breg1 = _mm256_loadu_ps(&b[vec_size*(8*i + 1)]);
tmp1 = _mm256_add_ps(_mm256_mul_ps(areg0,breg1), tmp1);
__m256 breg2 = _mm256_loadu_ps(&b[vec_size*(8*i + 2)]);
tmp2 = _mm256_add_ps(_mm256_mul_ps(areg0,breg2), tmp2);
__m256 breg3 = _mm256_loadu_ps(&b[vec_size*(8*i + 3)]);
tmp3 = _mm256_add_ps(_mm256_mul_ps(areg0,breg3), tmp3);
__m256 breg4 = _mm256_loadu_ps(&b[vec_size*(8*i + 4)]);
tmp4 = _mm256_add_ps(_mm256_mul_ps(areg0,breg4), tmp4);
__m256 breg5 = _mm256_loadu_ps(&b[vec_size*(8*i + 5)]);
tmp5 = _mm256_add_ps(_mm256_mul_ps(areg0,breg5), tmp5);
__m256 breg6 = _mm256_loadu_ps(&b[vec_size*(8*i + 6)]);
tmp6 = _mm256_add_ps(_mm256_mul_ps(areg0,breg6), tmp6);
__m256 breg7 = _mm256_loadu_ps(&b[vec_size*(8*i + 7)]);
tmp7 = _mm256_add_ps(_mm256_mul_ps(areg0,breg7), tmp7);
}
_mm256_storeu_ps(&c[0*vec_size], tmp0);
_mm256_storeu_ps(&c[1*vec_size], tmp1);
_mm256_storeu_ps(&c[2*vec_size], tmp2);
_mm256_storeu_ps(&c[3*vec_size], tmp3);
_mm256_storeu_ps(&c[4*vec_size], tmp4);
_mm256_storeu_ps(&c[5*vec_size], tmp5);
_mm256_storeu_ps(&c[6*vec_size], tmp6);
_mm256_storeu_ps(&c[7*vec_size], tmp7);
}
2 ответа
Для Sandy/Ivy Bridge вам нужно развернуть на 3:
- Только FP Add имеет зависимость от предыдущей итерации цикла
- FP Add может выдавать каждый цикл
- FP Add занимает три цикла для завершения
- Таким образом, развертывание на 3/1 = 3 полностью скрывает задержку
- FP Mul и FP Load не зависят от предыдущей итерации, и вы можете положиться на ядро OoO для их выдачи в почти оптимальном порядке. Эти инструкции могут влиять на коэффициент развертывания только в том случае, если они снижают пропускную способность FP Add (в данном случае это не так, FP Load + FP Add + FP Mul может выдавать каждый цикл).
Для Haswell вам нужно развернуть на 10:
- Только FMA зависит от предыдущей итерации цикла
- FMA может выпускать дважды каждый цикл (т.е. в среднем независимые инструкции занимают 0,5 цикла)
- FMA имеет задержку 5
- Таким образом, развертывание на 5/0,5 = 10 полностью скрывает задержку FMA
- Две микрооперации FP Load не зависят от предыдущей итерации и могут совместно использоваться с 2x FMA, поэтому они не влияют на коэффициент развертывания.
Я только отвечаю на свой вопрос здесь, чтобы добавить информацию.
Я пошел дальше и профилировал код Ivy Bridge. Когда я впервые проверил это в MSVC2012, развертывание более чем на два не сильно помогло. Однако я подозревал, что MSVC не реализовал встроенные функции оптимально, основываясь на моих наблюдениях за разницей в производительности между MSVC и GCC для высокооптимизированного кода матричного умножения. Поэтому я скомпилировал ядро в GCC с g++ -c -mavx -O3 -mabi=ms, преобразовал объект в COFF64 и поместил его в MSVC, и теперь я получаю, что развертывание на три дает лучший результат, подтверждающий ответ Марата Дункана.
Вот время в секундах, Xeon E5 1620 @3,6 ГГц MSVC2012
unroll time default time with GCC kernel
1 3.7 3.2
2 1.8 (2.0x faster) 1.6 (2.0x faster)
3 1.6 (2.3x faster) 1.2 (2.7x faster)
4 1.6 (2.3x faster) 1.2 (2.7x faster)
Вот времена на i5-4250U с использованием FMA с GCC в Linux (g++ -mavx -mfma -fopenmp -O3 main.cpp kernel_fma.cpp -o sum_fma)
unroll time
1 20.3
2 10.2 (2.0x faster)
3 6.7 (3.0x faster)
4 5.2 (4.0x faster)
8 2.9 (7.0x faster)
10 2.6 (7.8x faster)
Код ниже для Sandy-Bridge/Ivy Bridge. Для использования Haswell, например tmp0 = _mm256_fmadd_ps(a8,b8_1,tmp0) вместо.
kernel.cpp
#include <immintrin.h>
extern "C" void foo_unroll1(const int n, const float *b, float *c) {
__m256 tmp0 = _mm256_set1_ps(0.0f);
__m256 a8 = _mm256_set1_ps(1.0f);
for(int i=0; i<n; i+=8) {
__m256 b8 = _mm256_loadu_ps(&b[i + 0]);
tmp0 = _mm256_add_ps(_mm256_mul_ps(a8,b8), tmp0);
}
_mm256_storeu_ps(c, tmp0);
}
extern "C" void foo_unroll2(const int n, const float *b, float *c) {
__m256 tmp0 = _mm256_set1_ps(0.0f);
__m256 tmp1 = _mm256_set1_ps(0.0f);
__m256 a8 = _mm256_set1_ps(1.0f);
for(int i=0; i<n; i+=16) {
__m256 b8_1 = _mm256_loadu_ps(&b[i + 0]);
tmp0 = _mm256_add_ps(_mm256_mul_ps(a8,b8_1), tmp0);
__m256 b8_2 = _mm256_loadu_ps(&b[i + 8]);
tmp1 = _mm256_add_ps(_mm256_mul_ps(a8,b8_2), tmp1);
}
tmp0 = _mm256_add_ps(tmp0,tmp1);
_mm256_storeu_ps(c, tmp0);
}
extern "C" void foo_unroll3(const int n, const float *b, float *c) {
__m256 tmp0 = _mm256_set1_ps(0.0f);
__m256 tmp1 = _mm256_set1_ps(0.0f);
__m256 tmp2 = _mm256_set1_ps(0.0f);
__m256 a8 = _mm256_set1_ps(1.0f);
for(int i=0; i<n; i+=24) {
__m256 b8_1 = _mm256_loadu_ps(&b[i + 0]);
tmp0 = _mm256_add_ps(_mm256_mul_ps(a8,b8_1), tmp0);
__m256 b8_2 = _mm256_loadu_ps(&b[i + 8]);
tmp1 = _mm256_add_ps(_mm256_mul_ps(a8,b8_2), tmp1);
__m256 b8_3 = _mm256_loadu_ps(&b[i + 16]);
tmp2 = _mm256_add_ps(_mm256_mul_ps(a8,b8_3), tmp2);
}
tmp0 = _mm256_add_ps(tmp0,_mm256_add_ps(tmp1,tmp2));
_mm256_storeu_ps(c, tmp0);
}
extern "C" void foo_unroll4(const int n, const float *b, float *c) {
__m256 tmp0 = _mm256_set1_ps(0.0f);
__m256 tmp1 = _mm256_set1_ps(0.0f);
__m256 tmp2 = _mm256_set1_ps(0.0f);
__m256 tmp3 = _mm256_set1_ps(0.0f);
__m256 a8 = _mm256_set1_ps(1.0f);
for(int i=0; i<n; i+=32) {
__m256 b8_1 = _mm256_loadu_ps(&b[i + 0]);
tmp0 = _mm256_add_ps(_mm256_mul_ps(a8,b8_1), tmp0);
__m256 b8_2 = _mm256_loadu_ps(&b[i + 8]);
tmp1 = _mm256_add_ps(_mm256_mul_ps(a8,b8_2), tmp1);
__m256 b8_3 = _mm256_loadu_ps(&b[i + 16]);
tmp2 = _mm256_add_ps(_mm256_mul_ps(a8,b8_3), tmp2);
__m256 b8_4 = _mm256_loadu_ps(&b[i + 24]);
tmp3 = _mm256_add_ps(_mm256_mul_ps(a8,b8_4), tmp3);
}
tmp0 = _mm256_add_ps(_mm256_add_ps(tmp0,tmp1),_mm256_add_ps(tmp2,tmp3));
_mm256_storeu_ps(c, tmp0);
}
main.cpp
#include <stdio.h>
#include <omp.h>
#include <immintrin.h>
extern "C" void foo_unroll1(const int n, const float *b, float *c);
extern "C" void foo_unroll2(const int n, const float *b, float *c);
extern "C" void foo_unroll3(const int n, const float *b, float *c);
extern "C" void foo_unroll4(const int n, const float *b, float *c);
int main() {
const int n = 3*1<<10;
const int r = 10000000;
double dtime;
float *b = (float*)_mm_malloc(sizeof(float)*n, 64);
float *c = (float*)_mm_malloc(8, 64);
for(int i=0; i<n; i++) b[i] = 1.0f;
__m256 out;
dtime = omp_get_wtime();
for(int i=0; i<r; i++) foo_unroll1(n, b, c);
dtime = omp_get_wtime() - dtime;
printf("%f, ", dtime); for(int i=0; i<8; i++) printf("%f ", c[i]); printf("\n");
dtime = omp_get_wtime();
for(int i=0; i<r; i++) foo_unroll2(n, b, c);
dtime = omp_get_wtime() - dtime;
printf("%f, ", dtime); for(int i=0; i<8; i++) printf("%f ", c[i]); printf("\n");
dtime = omp_get_wtime();
for(int i=0; i<r; i++) foo_unroll3(n, b, c);
dtime = omp_get_wtime() - dtime;
printf("%f, ", dtime); for(int i=0; i<8; i++) printf("%f ", c[i]); printf("\n");
dtime = omp_get_wtime();
for(int i=0; i<r; i++) foo_unroll4(n, b, c);
dtime = omp_get_wtime() - dtime;
printf("%f, ", dtime); for(int i=0; i<8; i++) printf("%f ", c[i]); printf("\n");
}