Я хочу добиться максимальной пропускной способности следующих операций с процессорами Intel.
for(int i=0; i<n; i++) z[i] = x[i] + y[i]; //n=2048
где x, y и z - массивы с плавающей точкой. Я делаю это на системах Хасуэлла, Айви-Бридж и Уэстмира.
Я изначально выделял такую память
char *a = (char*)_mm_malloc(sizeof(float)*n, 64);
char *b = (char*)_mm_malloc(sizeof(float)*n, 64);
char *c = (char*)_mm_malloc(sizeof(float)*n, 64);
float *x = (float*)a; float *y = (float*)b; float *z = (float*)c;
Когда я это сделал, я получил около 50% максимальной пропускной способности, ожидаемой для каждой системы.
Пиковые значения рассчитываются как frequency * average bytes/clock_cycle. Средний байт/тактовый цикл для каждой системы:
Core2: two 16 byte reads one 16 byte write per 2 clock cycles -> 24 bytes/clock cycle
SB/IB: two 32 byte reads and one 32 byte write per 2 clock cycles -> 48 bytes/clock cycle
Haswell: two 32 byte reads and one 32 byte write per clock cycle -> 96 bytes/clock cycle
Это означает, что, например, на Haswell я Я наблюдаю только 48 байтов/такт (может быть два чтения за один такт и один записать следующий такт).
Я распечатал разницу в адресе b-a и c-b и каждый из них - 8256 байт. Значение 8256 равно 8192 + 64. Таким образом, каждый из них больше размера массива (8192 байта) одной строкой кэша.
По прихоти я попытался выделить такую память.
const int k = 0;
char *mem = (char*)_mm_malloc(1<<18,4096);
char *a = mem;
char *b = a+n*sizeof(float)+k*64;
char *c = b+n*sizeof(float)+k*64;
float *x = (float*)a; float *y = (float*)b; float *z = (float*)c;
Это почти удвоило мою максимальную пропускную способность, так что теперь я получаю около 90% максимальной пропускной способности. Однако, когда я пробовал k=1, он упал до 50%. Я пробовал другие значения k и обнаружил, что, например, k=2, k=33, k=65 получает только 50% пика, но, например, k=10, k=32, k=63 дали полную скорость. Я не понимаю этого.
В руководстве по микроархитектуре Agner Fog он говорит, что существует ложная зависимость с адресом памяти с тем же набором и смещением
Невозможно читать и писать одновременно с адресов которые разнесены на 4 Кбайта.
Но именно там, где я вижу самую большую выгоду! Когда k=0 адрес памяти отличается точно 2*4096 байтами. Агнер также рассказывает о конфликтах банка-кэша. Но Хасуэлл и Уэстмир не предполагают иметь эти банковские конфликты, чтобы не объяснять, что я наблюдаю. Что происходит!?
Я понимаю, что исполнение OoO решает, какой адрес читать и писать так, даже если адреса памяти массивов отличаются примерно на 4096 байт, что не обязательно означает, что процессор читает, например. &x[0] и записывает &z[0] одновременно, но тогда почему бы отключить одну строку кэша, чтобы он задохнулся?
Редактировать: На основании ответа Евгения Клюева я теперь верю, что это то, что Агнер Фог называет "фальшивым магазином переадресации". В своем руководстве под Pentium Pro, II и II он пишет:
Интересно, что при письме и чтении вы можете получить фальшивый магазин для пересылки совершенно разные адреса, если они имеют одинаковое заданное значение в разных кешках банки:
; Example 5.28. Bogus store-to-load forwarding stall
mov byte ptr [esi], al
mov ebx, dword ptr [esi+4092]
; No stall
mov ecx, dword ptr [esi+4096]
; Bogus stall
Изменить: Здесь приведена таблица эффективности каждой системы для k=0 и k=1.
k=0 k=1
Westmere: 99% 66%
Ivy Bridge: 98% 44%
Haswell: 90% 49%
Я думаю, что могу объяснить эти числа, если предположить, что для k=1, что записи и чтения не могут произойти в одном такте.
cycle Westmere Ivy Bridge Haswell
1 read 16 read 16 read 16 read 32 read 32
2 write 16 read 16 read 16 write 32
3 write 16
4 write 16
k=1/k=0 peak 16/24=66% 24/48=50% 48/96=50%
Эта теория работает очень хорошо. Мост Ivy немного ниже, чем я ожидал, но Ivy Bridge страдает от конфликтов банковского кэша, когда другие не делают этого, что может быть другим эффектом для рассмотрения.
Ниже приведен рабочий код, чтобы проверить это самостоятельно. В системе без компиляции AVX с g++ -O3 sum.cpp иначе скомпилируйте с помощью g++ -O3 -mavx sum.cpp. Попробуйте изменить значение k.
//sum.cpp
#include <x86intrin.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#define TIMER_TYPE CLOCK_REALTIME
double time_diff(timespec start, timespec end)
{
timespec temp;
if ((end.tv_nsec-start.tv_nsec)<0) {
temp.tv_sec = end.tv_sec-start.tv_sec-1;
temp.tv_nsec = 1000000000+end.tv_nsec-start.tv_nsec;
} else {
temp.tv_sec = end.tv_sec-start.tv_sec;
temp.tv_nsec = end.tv_nsec-start.tv_nsec;
}
return (double)temp.tv_sec + (double)temp.tv_nsec*1E-9;
}
void sum(float * __restrict x, float * __restrict y, float * __restrict z, const int n) {
#if defined(__GNUC__)
x = (float*)__builtin_assume_aligned (x, 64);
y = (float*)__builtin_assume_aligned (y, 64);
z = (float*)__builtin_assume_aligned (z, 64);
#endif
for(int i=0; i<n; i++) {
z[i] = x[i] + y[i];
}
}
#if (defined(__AVX__))
void sum_avx(float *x, float *y, float *z, const int n) {
float *x1 = x;
float *y1 = y;
float *z1 = z;
for(int i=0; i<n/64; i++) { //unroll eight times
_mm256_store_ps(z1+64*i+ 0,_mm256_add_ps(_mm256_load_ps(x1+64*i+ 0), _mm256_load_ps(y1+64*i+ 0)));
_mm256_store_ps(z1+64*i+ 8,_mm256_add_ps(_mm256_load_ps(x1+64*i+ 8), _mm256_load_ps(y1+64*i+ 8)));
_mm256_store_ps(z1+64*i+ 16,_mm256_add_ps(_mm256_load_ps(x1+64*i+16), _mm256_load_ps(y1+64*i+ 16)));
_mm256_store_ps(z1+64*i+ 24,_mm256_add_ps(_mm256_load_ps(x1+64*i+24), _mm256_load_ps(y1+64*i+ 24)));
_mm256_store_ps(z1+64*i+ 32,_mm256_add_ps(_mm256_load_ps(x1+64*i+32), _mm256_load_ps(y1+64*i+ 32)));
_mm256_store_ps(z1+64*i+ 40,_mm256_add_ps(_mm256_load_ps(x1+64*i+40), _mm256_load_ps(y1+64*i+ 40)));
_mm256_store_ps(z1+64*i+ 48,_mm256_add_ps(_mm256_load_ps(x1+64*i+48), _mm256_load_ps(y1+64*i+ 48)));
_mm256_store_ps(z1+64*i+ 56,_mm256_add_ps(_mm256_load_ps(x1+64*i+56), _mm256_load_ps(y1+64*i+ 56)));
}
}
#else
void sum_sse(float *x, float *y, float *z, const int n) {
float *x1 = x;
float *y1 = y;
float *z1 = z;
for(int i=0; i<n/32; i++) { //unroll eight times
_mm_store_ps(z1+32*i+ 0,_mm_add_ps(_mm_load_ps(x1+32*i+ 0), _mm_load_ps(y1+32*i+ 0)));
_mm_store_ps(z1+32*i+ 4,_mm_add_ps(_mm_load_ps(x1+32*i+ 4), _mm_load_ps(y1+32*i+ 4)));
_mm_store_ps(z1+32*i+ 8,_mm_add_ps(_mm_load_ps(x1+32*i+ 8), _mm_load_ps(y1+32*i+ 8)));
_mm_store_ps(z1+32*i+ 12,_mm_add_ps(_mm_load_ps(x1+32*i+12), _mm_load_ps(y1+32*i+ 12)));
_mm_store_ps(z1+32*i+ 16,_mm_add_ps(_mm_load_ps(x1+32*i+16), _mm_load_ps(y1+32*i+ 16)));
_mm_store_ps(z1+32*i+ 20,_mm_add_ps(_mm_load_ps(x1+32*i+20), _mm_load_ps(y1+32*i+ 20)));
_mm_store_ps(z1+32*i+ 24,_mm_add_ps(_mm_load_ps(x1+32*i+24), _mm_load_ps(y1+32*i+ 24)));
_mm_store_ps(z1+32*i+ 28,_mm_add_ps(_mm_load_ps(x1+32*i+28), _mm_load_ps(y1+32*i+ 28)));
}
}
#endif
int main () {
const int n = 2048;
const int k = 0;
float *z2 = (float*)_mm_malloc(sizeof(float)*n, 64);
char *mem = (char*)_mm_malloc(1<<18,4096);
char *a = mem;
char *b = a+n*sizeof(float)+k*64;
char *c = b+n*sizeof(float)+k*64;
float *x = (float*)a;
float *y = (float*)b;
float *z = (float*)c;
printf("x %p, y %p, z %p, y-x %d, z-y %d\n", a, b, c, b-a, c-b);
for(int i=0; i<n; i++) {
x[i] = (1.0f*i+1.0f);
y[i] = (1.0f*i+1.0f);
z[i] = 0;
}
int repeat = 1000000;
timespec time1, time2;
sum(x,y,z,n);
#if (defined(__AVX__))
sum_avx(x,y,z2,n);
#else
sum_sse(x,y,z2,n);
#endif
printf("error: %d\n", memcmp(z,z2,sizeof(float)*n));
while(1) {
clock_gettime(TIMER_TYPE, &time1);
#if (defined(__AVX__))
for(int r=0; r<repeat; r++) sum_avx(x,y,z,n);
#else
for(int r=0; r<repeat; r++) sum_sse(x,y,z,n);
#endif
clock_gettime(TIMER_TYPE, &time2);
double dtime = time_diff(time1,time2);
double peak = 1.3*96; //haswell @1.3GHz
//double peak = 3.6*48; //Ivy Bridge @ 3.6Ghz
//double peak = 2.4*24; // Westmere @ 2.4GHz
double rate = 3.0*1E-9*sizeof(float)*n*repeat/dtime;
printf("dtime %f, %f GB/s, peak, %f, efficiency %f%%\n", dtime, rate, peak, 100*rate/peak);
}
}
