Рассчитать среднее 4d векторов с SSE
Я пытаюсь ускорить вычисление среднего из 4d векторов, помещенных в массив. Вот мой код:
#include <sys/time.h>
#include <sys/param.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <xmmintrin.h>
typedef float dot[4];
#define N 1000000
double gettime ()
{
struct timeval tv;
gettimeofday (&tv, 0);
return (double)tv.tv_sec + (0.000001 * (double)tv.tv_usec);
}
void calc_avg1 (dot res, const dot array[], int n)
{
int i,j;
memset (res, 0, sizeof (dot));
for (i = 0; i < n; i++)
{
for (j = 0; j<4; j++) res[j] += array[i][j];
}
for (j = 0; j<4; j++) res[j] /= n;
}
void calc_avg2 (dot res, const dot array[], int n)
{
int i;
__v4sf r = _mm_set1_ps (0.0);
for (i=0; i<n; i++) r += _mm_load_ps (array[i]);
r /= _mm_set1_ps ((float)n);
_mm_store_ps (res, r);
}
int main ()
{
void *space = malloc (N*sizeof(dot)+15);
dot *array = (dot*)(((unsigned long)space+15) & ~(unsigned long)15);
dot avg __attribute__((aligned(16)));
int i;
double time;
for (i = 0; i < N; i++)
{
array[i][0] = 1.0*random();
array[i][1] = 1.0*random();
array[i][2] = 1.0*random();
}
time = gettime();
calc_avg1 (avg, array, N);
time = gettime() - time;
printf ("%f\n%f %f %f\n", time, avg[0], avg[1], avg[2]);
time = gettime();
calc_avg2 (avg, array, N);
time = gettime() - time;
printf ("%f\n%f %f %f\n", time, avg[0], avg[1], avg[2]);
return 0;
}
Итак, как вы можете видеть calc_avg1 использует наивные циклы от 0 до 4 и calc_avg2 заменяет их инструкциями SSE. Я компилирую этот код с помощью Clang 3.4:
cc -O2 -o test test.c
Вот разборка функций calc_avgX:
0000000000400860 <calc_avg1>:
400860: 55 push %rbp
400861: 48 89 e5 mov %rsp,%rbp
400864: 85 d2 test %edx,%edx
400866: 0f 57 c0 xorps %xmm0,%xmm0
400869: 0f 11 07 movups %xmm0,(%rdi)
40086c: 7e 42 jle 4008b0 <calc_avg1+0x50>
40086e: 48 83 c6 0c add $0xc,%rsi
400872: 0f 57 c0 xorps %xmm0,%xmm0
400875: 89 d0 mov %edx,%eax
400877: 0f 57 c9 xorps %xmm1,%xmm1
40087a: 0f 57 d2 xorps %xmm2,%xmm2
40087d: 0f 57 db xorps %xmm3,%xmm3
400880: f3 0f 58 5e f4 addss -0xc(%rsi),%xmm3
400885: f3 0f 11 1f movss %xmm3,(%rdi)
400889: f3 0f 58 56 f8 addss -0x8(%rsi),%xmm2
40088e: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
400893: f3 0f 58 4e fc addss -0x4(%rsi),%xmm1
400898: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
40089d: f3 0f 58 06 addss (%rsi),%xmm0
4008a1: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008a6: 48 83 c6 10 add $0x10,%rsi
4008aa: ff c8 dec %eax
4008ac: 75 d2 jne 400880 <calc_avg1+0x20>
4008ae: eb 0c jmp 4008bc <calc_avg1+0x5c>
4008b0: 0f 57 c0 xorps %xmm0,%xmm0
4008b3: 0f 57 c9 xorps %xmm1,%xmm1
4008b6: 0f 57 d2 xorps %xmm2,%xmm2
4008b9: 0f 57 db xorps %xmm3,%xmm3
4008bc: f3 0f 2a e2 cvtsi2ss %edx,%xmm4
4008c0: f3 0f 5e dc divss %xmm4,%xmm3
4008c4: f3 0f 11 1f movss %xmm3,(%rdi)
4008c8: f3 0f 5e d4 divss %xmm4,%xmm2
4008cc: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
4008d1: f3 0f 5e cc divss %xmm4,%xmm1
4008d5: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
4008da: f3 0f 5e c4 divss %xmm4,%xmm0
4008de: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008e3: 5d pop %rbp
4008e4: c3 retq
4008e5: 66 66 2e 0f 1f 84 00 nopw %cs:0x0(%rax,%rax,1)
4008ec: 00 00 00 00
00000000004008f0 <calc_avg2>:
4008f0: 55 push %rbp
4008f1: 48 89 e5 mov %rsp,%rbp
4008f4: 85 d2 test %edx,%edx
4008f6: 0f 57 c0 xorps %xmm0,%xmm0
4008f9: 7e 10 jle 40090b <calc_avg2+0x1b>
4008fb: 89 d0 mov %edx,%eax
4008fd: 0f 1f 00 nopl (%rax)
400900: 0f 58 06 addps (%rsi),%xmm0
400903: 48 83 c6 10 add $0x10,%rsi
400907: ff c8 dec %eax
400909: 75 f5 jne 400900 <calc_avg2+0x10>
40090b: 66 0f 6e ca movd %edx,%xmm1
40090f: 66 0f 70 c9 00 pshufd $0x0,%xmm1,%xmm1
400914: 0f 5b c9 cvtdq2ps %xmm1,%xmm1
400917: 0f 5e c1 divps %xmm1,%xmm0
40091a: 0f 29 07 movaps %xmm0,(%rdi)
40091d: 5d pop %rbp
40091e: c3 retq
40091f: 90 nop
И вот результат:
> ./test
0.004287
1073864320.000000 1074018048.000000 1073044224.000000
0.003661
1073864320.000000 1074018048.000000 1073044224.000000
Так что версия SSE в 1,17 раза быстрее. Но когда я пытаюсь выполнить, казалось бы, ту же самую работу, которая заключается в вычислении среднего значения скаляров одинарной точности в массиве (как описано здесь, например, уменьшение вектора с плавающей запятой в SSE), версия SSE работает в 3,32 раза быстрее. Вот код функций calc_avgX:
float calc_avg1 (const float array[], int n)
{
int i;
float avg = 0;
for (i = 0; i < n; i++) avg += array[i];
return avg / n;
}
float calc_avg3 (const float array[], int n)
{
int i;
__v4sf r = _mm_set1_ps (0.0);
for (i=0; i<n; i+=4) r += _mm_load_ps (&(array[i]));
r = _mm_hadd_ps (r, r);
r = _mm_hadd_ps (r, r);
return r[0] / n;
}
Поэтому мой вопрос таков: почему я так выигрываю от SSE в последнем примере (вычисление среднего для скаляров с плавающей запятой), а не в первом (вычисление среднего для 4d векторов)? Для меня эти работы практически одинаковы. Как правильно ускорить вычисления в первом примере, если я делаю это неправильно?
UPD: Если вы считаете это уместным, я также приведу разборку второго примера, где рассчитывается среднее значение скаляров (также скомпилировано с помощью clang3.4 -O2).
0000000000400860 <calc_avg1>:
400860: 55 push %rbp
400861: 48 89 e5 mov %rsp,%rbp
400864: 85 d2 test %edx,%edx
400866: 0f 57 c0 xorps %xmm0,%xmm0
400869: 0f 11 07 movups %xmm0,(%rdi)
40086c: 7e 42 jle 4008b0 <calc_avg1+0x50>
40086e: 48 83 c6 0c add $0xc,%rsi
400872: 0f 57 c0 xorps %xmm0,%xmm0
400875: 89 d0 mov %edx,%eax
400877: 0f 57 c9 xorps %xmm1,%xmm1
40087a: 0f 57 d2 xorps %xmm2,%xmm2
40087d: 0f 57 db xorps %xmm3,%xmm3
400880: f3 0f 58 5e f4 addss -0xc(%rsi),%xmm3
400885: f3 0f 11 1f movss %xmm3,(%rdi)
400889: f3 0f 58 56 f8 addss -0x8(%rsi),%xmm2
40088e: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
400893: f3 0f 58 4e fc addss -0x4(%rsi),%xmm1
400898: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
40089d: f3 0f 58 06 addss (%rsi),%xmm0
4008a1: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008a6: 48 83 c6 10 add $0x10,%rsi
4008aa: ff c8 dec %eax
4008ac: 75 d2 jne 400880 <calc_avg1+0x20>
4008ae: eb 0c jmp 4008bc <calc_avg1+0x5c>
4008b0: 0f 57 c0 xorps %xmm0,%xmm0
4008b3: 0f 57 c9 xorps %xmm1,%xmm1
4008b6: 0f 57 d2 xorps %xmm2,%xmm2
4008b9: 0f 57 db xorps %xmm3,%xmm3
4008bc: f3 0f 2a e2 cvtsi2ss %edx,%xmm4
4008c0: f3 0f 5e dc divss %xmm4,%xmm3
4008c4: f3 0f 11 1f movss %xmm3,(%rdi)
4008c8: f3 0f 5e d4 divss %xmm4,%xmm2
4008cc: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
4008d1: f3 0f 5e cc divss %xmm4,%xmm1
4008d5: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
4008da: f3 0f 5e c4 divss %xmm4,%xmm0
4008de: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008e3: 5d pop %rbp
4008e4: c3 retq
4008e5: 66 66 2e 0f 1f 84 00 nopw %cs:0x0(%rax,%rax,1)
4008ec: 00 00 00 00
00000000004008d0 <calc_avg3>:
4008d0: 55 push %rbp
4008d1: 48 89 e5 mov %rsp,%rbp
4008d4: 31 c0 xor %eax,%eax
4008d6: 85 f6 test %esi,%esi
4008d8: 0f 57 c0 xorps %xmm0,%xmm0
4008db: 7e 0f jle 4008ec <calc_avg3+0x1c>
4008dd: 0f 1f 00 nopl (%rax)
4008e0: 0f 58 04 87 addps (%rdi,%rax,4),%xmm0
4008e4: 48 83 c0 04 add $0x4,%rax
4008e8: 39 f0 cmp %esi,%eax
4008ea: 7c f4 jl 4008e0 <calc_avg3+0x10>
4008ec: 66 0f 70 c8 01 pshufd $0x1,%xmm0,%xmm1
4008f1: f3 0f 58 c8 addss %xmm0,%xmm1
4008f5: 66 0f 70 d0 03 pshufd $0x3,%xmm0,%xmm2
4008fa: 0f 12 c0 movhlps %xmm0,%xmm0
4008fd: f3 0f 58 c1 addss %xmm1,%xmm0
400901: f3 0f 58 c2 addss %xmm2,%xmm0
400905: 0f 57 c9 xorps %xmm1,%xmm1
400908: f3 0f 2a ce cvtsi2ss %esi,%xmm1
40090c: f3 0f 5e c1 divss %xmm1,%xmm0
400910: 5d pop %rbp
400911: c3 retq
400912: 66 66 66 66 66 2e 0f nopw %cs:0x0(%rax,%rax,1)
400919: 1f 84 00 00 00 00 00
1 ответ
Извините, этот ответ получился немного длинным и бессвязным. Я провел несколько тестов, но я не стал долго редактировать предыдущие вещи, подумав о чем-то другом, чтобы попробовать.
Ваш рабочий набор составляет 15,25 МБ (16 МБ). Обычно для тестирования такой процедуры вы должны усреднять меньший буфер несколько раз, чтобы он помещался в кэш. Вы не видите большой разницы между медленной и быстрой версиями, потому что разница скрыта узким местом в памяти.
calc_avg1 не векторизация вообще (обратите внимание на addss, ss означает скаляр, с одинарной точностью, в отличие от addps (упакован с одинарной точностью)). Я думаю, что он не может автоматически векторизоваться, даже если встроен в main, потому что не может быть уверен, что нет NaNs в положении 4-го вектора, что вызовет исключение FP, которого не будет скалярный код. Я попытался скомпилировать его для Sandybridge с помощью gcc 4.9.2 -O3 -march=native -ffast-mathи с clang-3.5, но не повезло ни с одним.
Несмотря на это, версия встроена в main только работает немного медленнее, потому что память является узким местом. 32-разрядные нагрузки могут почти не отставать от 128-разрядных при загрузке основной памяти. (Не встроенная версия будет плохой, хотя: каждый += результат сохраняется в res массив, потому что цикл накапливается непосредственно в памяти, которая может иметь другие ссылки на него. Таким образом, он должен сделать каждую операцию видимой в магазине. Это версия, для которой вы разместили разборку, кстати. Разбираемся, какая часть main была получена при помощи -S -fverbose-asm.)
К сожалению, Clang и GCC не могут автоматически векторизовать __v4sf от 4 до 8 AVX.
Поцарапайте это, после упаковки for (int i=0; i<4000 ; i++) вокруг звонков calc_avgXи сокращение N до 10 Кб, гкц -O3 превращает внутренний внутренний цикл avg1 в:
400690: c5 f8 10 08 vmovups (%rax),%xmm1
400694: 48 83 c0 20 add $0x20,%rax
400698: c4 e3 75 18 48 f0 01 vinsertf128 $0x1,-0x10(%rax),%ymm1,%ymm1
40069f: c5 fc 58 c1 vaddps %ymm1,%ymm0,%ymm0
4006a3: 48 39 d8 cmp %rbx,%rax
4006a6: 75 e8 jne 400690 <main+0xe0>
$ (get CPU to max-turbo frequency) && time ./a.out
0.016515
1071570752.000000 1066917696.000000 1073897344.000000
0.032875
1071570944.000000 1066916416.000000 1073895680.000000
Это странно; Я понятия не имею, почему он не просто использует 32B нагрузки. Это действительно использует 32B vaddps, что является узким местом при работе с набором данных, который помещается в кэш L2.
IDK, почему ему удалось автоматически векторизовать внутренний цикл, когда он был внутри другого цикла. Обратите внимание, что это относится только к версии, встроенной в main, Вызываемая версия все еще только для скаляров. Также обратите внимание, что только gcc справился с этим. лязг 3.5 не сделал. Может быть, GCC знал, что он будет использовать malloc таким образом, чтобы вернуть обнуленный буфер (так что не нужно беспокоиться о NaNс 4-го элемента)?
Я также удивлен, что Clang не векторизован avg1 не медленнее, когда все помещается в кеш. N=10000, счетчик повторений = 40к.
3.3GHz SNB i5 2500k, max turbo = 3.8GHz.
avg1: 0.350422s: clang -O3 -march=native (not vectorized. loop of 6 scalar addss with memory operands)
avg2: 0.320173s: clang -O3 -march=native
avg1: 0.497040s: clang -O3 -march=native -ffast-math (haven't looked at asm to see what happened)
avg1: 0.160374s: gcc -O3 -march=native (256b addps, with 2 128b loads)
avg2: 0.321028s: gcc -O3 -march=native (128b addps with a memory operand)
avg2: ~0.16: clang, unrolled with 2 dependency chains to hide latency (see below).
avg2: ~0.08: unrolled with 4 dep chains
avg2: ~0.04: in theory unrolled-by-4 with 256b AVX. I didn't try unrolling the one gcc auto-vectorized with 256b addps
Так что большой сюрприз в том, что скалярный лязг avg1 код идет в ногу с avg2, Возможно, цепочка зависимостей, переносимая циклами, является более узким местом?
perf показывает 1,47 insns за цикл для не векторизованного Clang avg1, что, вероятно, насыщает модуль добавления FP на порту 1. (Большинство инструкций цикла - это добавления).
Тем не мение, avg2с использованием 128b addps с операндом памяти, получает только 0,58 insns за цикл. Уменьшение размера массива еще в 10 раз до N=1000, получает 0,60 insns за цикл, вероятно, из-за большего количества времени, проводимого в прологе / эпилоге. Так что я думаю, что есть серьезная проблема с цепочкой зависимостей, переносимых циклами. clang развертывает цикл на 4, но использует только один аккумулятор. Цикл имеет 7 инструкций, которые декодируют до 10 моп. (Каждый vaddps 2, так как он используется с операндом памяти с режимом адресации с 2-мя регистрами, что предотвращает микроплавление. cmp а также jne макро-предохранитель). http://www.brendangregg.com/perf.html говорит, что perf событие для UOPS_DISPATCHED.CORE является r2b1, так:
$ perf stat -d -e cycles,instructions,r2b1 ./a.out
0.031793
1053298112.000000 1052673664.000000 1116960256.000000
Performance counter stats for './a.out':
118,453,541 cycles
71,181,299 instructions # 0.60 insns per cycle
102,025,443 r2b1 # this is uops, but perf doesn't have a nice name for it
40,256,019 L1-dcache-loads
21,254 L1-dcache-load-misses # 0.05% of all L1-dcache hits
9,588 LLC-loads
0 LLC-load-misses:HG # 0.00% of all LL-cache hits
0.032276233 seconds time elapsed
Это подтверждает мои 7:10 инструкции к анализу мопов. Это на самом деле не относится к проблеме производительности: цикл работает намного медленнее, чем верхний предел 4 моп / цикл. Изменение внутреннего цикла для получения двух отдельных цепочек развертывания удваивает пропускную способность (60M циклов вместо 117M, но 81M insns вместо 71M):
for (i=0; i<n-1; i+=2) { // TODO: make sure the loop end condition is correct
r0 += _mm_load_ps (array[i]);
r1 += _mm_load_ps (array[i+1]);
}
r0 += r1;
Развертывание на 4 (с 4 аккумуляторами, которые вы объединяете в конце цикла) снова удваивает производительность. (до 42M циклов, 81M Insns, 112M Uops.) Внутренний цикл имеет 4x vaddps -0x30(%rcx),%xmm4,%xmm4 (и аналогичные), 2x add, cmp, jl, Эта форма vaddps должен микро-предохранитель, но я все еще вижу гораздо больше мопов, чем инструкции, поэтому я думаю, r2b1 считает неиспользованные домены мопов (Linux perf не имеет хороших документов для платформных событий HW). Сгибать N Опять же, чтобы убедиться, что это самый внутренний цикл, который полностью доминирует во всех подсчетах, я вижу соотношение uop:insn, равное 1,39, которое хорошо соответствует 8 insns, 11 uops (1,375) (считая vaddps как 2, но считая cmp + jl как один). Я нашел http://www.bnikolic.co.uk/blog/hpc-prof-events.html, в котором есть полный список поддерживаемых событий perf, включая их коды для Sandybridge. (И инструкции о том, как сбросить таблицу для любого другого процессора). (Ищите Code: строка в каждом блоке. Вам нужен байт umask, а затем код, как аргумент perf.)
# a.out does only avg2, as an unrolled-by-4 version.
$ perf stat -d -e cycles,instructions,r14a1,r2b1,r10e,r2c2,r1c2 ./a.out
0.011331
1053298752.000000 1052674496.000000 1116959488.000000
Performance counter stats for './a.out':
42,250,312 cycles [34.11%]
56,103,429 instructions # 1.33 insns per cycle
20,864,416 r14a1 # UOPS_DISPATCHED_PORT: 0x14=port2&3 loads
111,943,380 r2b1 # UOPS_DISPATCHED: (2->umask 00 -> this core, any thread).
72,208,772 r10e # UOPS_ISSUED: fused-domain
71,422,907 r2c2 # UOPS_RETIRED: retirement slots used (fused-domain)
111,597,049 r1c2 # UOPS_RETIRED: ALL (unfused-domain)
0 L1-dcache-loads
18,470 L1-dcache-load-misses # 0.00% of all L1-dcache hits
5,717 LLC-loads [66.05%]
0 LLC-load-misses:HG # 0.00% of all LL-cache hits
0.011920301 seconds time elapsed
Так что да, похоже, что это может подсчитать число мопов в fused-domain и unused-domain!
Развертывание на 8 не помогает совсем: все еще 42M циклов. (но до 61 млн. insns и 97 млн. мопов благодаря меньшим накладным расходам петли). Аккуратно, лязг использует sub $-128, %rsi вместо добавления, потому что -128 вписывается в imm8+128 нет. Так что я думаю, что разворачивания на 4 достаточно, чтобы насытить порт добавления FP.
Что касается ваших 1avg-функций, которые возвращают один float, а не вектор, хорошо, clang не автоматически векторизовал первую, а gcc делает. Он выдает гигантский пролог и эпилог для скалярных сумм, пока не доберется до выровненного адреса, а затем делает 32B AVX vaddps в маленькой петле. Вы говорите, что нашли с ними гораздо больший разность скоростей, но вы, возможно, тестировали с меньшим буфером? Это объясняло бы ускорение векторного кода по сравнению с не векторным.