c - L1内存带宽: 50% drop in efficiency using addresses which differ by 4096+64 bytes

Question

我想用英特尔处理器实现以下操作的最大带宽。

for(int i=0; i<n; i++) z[i] = x[i] + y[i]; //n=2048

其中 x、y 和 z 是 float 组。我在 Haswell、Ivy Bridge 和 Westmere 系统上执行此操作。

我原来是这样分配内存的

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% 的预期峰值带宽。

峰值计算为 频率 * 平均字节数/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 I 上，我只观察到 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 的微架构手册中，他说内存地址具有相同的集合和偏移量存在错误的依赖关系

It is not possible to read and write simultaneously from addresses that are spaced by a multiple of 4 Kbytes.

但这正是我看到最大好处的地方!当 k=0 时，内存地址正好相差 2*4096 个字节。 Agner 还谈到了缓存库冲突。但是 Haswell 和 Westmere 不应该有这些银行冲突，所以这不应该解释我所观察到的。 发生了什么事!？

我了解 OoO 执行决定了读取和写入哪个地址，因此即使数组的内存地址相差 4096 个字节，也不一定意味着处理器读取例如&x[0] 并同时写入 &z[0] 但是为什么会被单个缓存行关闭导致它阻塞呢？

编辑:根据 Evgeny Kluev 的回答，我现在相信这就是 Agner Fog 所说的“虚假商店转发摊位”。在他的 Pentium Pro、II 和 II 手册中，他写道:

Interestingly, you can get a get a bogus store forwarding stall when writing and reading completely different addresses if they happen to have the same set-value in different cache banks:

; 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 桥遭受银行缓存冲突，而其他人没有，因此这可能是另一个需要考虑的影响。

以下是您自己测试的工作代码。在没有 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);
    }
}

Answer 1

我认为 a 和 b 之间的差距并不重要。在 b 和 c 之间只留下一个差距后，我在 Haswell 上得到了以下结果:

k   %
-----
1  48
2  48
3  48
4  48
5  46
6  53
7  59
8  67
9  73
10 81
11 85
12 87
13 87
...
0  86

由于众所周知 Haswell 没有存储库冲突，因此唯一剩下的解释是内存地址之间的错误依赖(并且您在 Agner Fog 的微架构手册中找到了正确解释此问题的位置)。 bank 冲突和错误共享之间的区别在于，bank 冲突防止在同一时钟周期内访问同一个 bank 两次，而错误共享阻止在您将某些内容写入相同偏移量之后从 4K 内存中的某个偏移量读取(不仅在同一个时钟周期内，但也在写入后的几个时钟周期内)。

由于您的代码(对于 k=0)仅在 对同一偏移量进行两次读取之后写入任何偏移量，并且在很长一段时间内都不会从中读取，这种情况应该算是“最好的”，所以我把 k=0 放在了表格的最后。对于 k=1，您总是从最近被覆盖的偏移量中读取，这意味着错误共享，因此性能下降。随着写入和读取之间的 k 时间增加，CPU 内核有更多机会通过所有内存层次结构传递写入的数据(这意味着读取和写入的两个地址转换，更新缓存数据和标签以及从缓存，核心之间的数据同步，可能还有更多东西)。 k=12 或 24 个时钟(在我的 CPU 上)足以让每条写入的数据为后续读取操作做好准备，因此从这个值开始，性能恢复正常。看起来与 AMD 上的 20 多个时钟差别不大(如 @Mysticial 所说)。