c++ - 读/写一个步长远大于其宽度的矩阵会导致性能损失很大

Question

我正在对 1024x1024 矩阵进行密集矩阵乘法运算。我使用 64x64 block 使用循环阻塞/平铺来执行此操作。我创建了一个高度优化的 64x64 矩阵乘法函数(请参阅我的问题末尾的代码)。

gemm64(float *a, float *b, float *c, int stride).  

这是运行在图 block 上的代码。具有 16x16 block 的 1024x1204 矩阵。

for(int i=0; i<16; i++) {
    for(int j=0; j<16; j++) {
        for(int k=0; k<16; k++) {
            gemm64(&a[64*(i*1024 + k)], &b[64*(k*1024 + j)], &c[64*(i*1024 + j)], 1024);
        }
    }
}

但是，作为一个猜测，我尝试重新排列 b 的所有图 block 的内存(请参阅此问题的代码末尾)。新矩阵中的矩阵 b2因此每个图 block 的步幅为 64 而不是 1024。这有效地生成了一个步幅为 64 的 64x64 矩阵数组。

float *b2 = (float*)_mm_malloc(1024*1024*sizeof(float), 64);
reorder(b, b2, 1024);
for(int i=0; i<16; i++) {
    for(int j=0; j<16; j++) {
        for(int k=0; k<16; k++) {
            gemm64_v2(&a[64*(i*1024 + k)], &b2[64*64*(k*16 + j)], &c[64*(i*1024 + j)], 64);
        }
    }
}
_mm_free(b2);

注意 b 的偏移量是如何从 &b[64*(k*1024 + j)] 改变的至 &b2[64*64*(k*16 + j)]并且步幅传递给了gemm64已从 1024 变为 64。

在我的 Sandy Bridge 系统上，我的代码的性能从峰值触发器的不到 20% 跃升至 70%!

为什么以这种方式重新排列 b 矩阵中的图 block 会产生如此巨大的差异？

数组 a、b、b2 和 c 是 64 字节对齐的。

extern "C" void gemm64(float *a, float*b, float*c, int stride) {
    for(int i=0; i<64; i++) {
        row_m64x64(&a[1024*i], b, &c[1024*i], stride);
    }
}

void row_m64x64(const float *a, const float *b, float *c, int stride) {
    __m256 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
    tmp0 = _mm256_loadu_ps(&c[ 0]);
    tmp1 = _mm256_loadu_ps(&c[ 8]);
    tmp2 = _mm256_loadu_ps(&c[16]);
    tmp3 = _mm256_loadu_ps(&c[24]);
    tmp4 = _mm256_loadu_ps(&c[32]);
    tmp5 = _mm256_loadu_ps(&c[40]);
    tmp6 = _mm256_loadu_ps(&c[48]);
    tmp7 = _mm256_loadu_ps(&c[56]);

    for(int i=0; i<64; i++) {
        __m256 areg0 = _mm256_broadcast_ss(&a[i]);

        __m256 breg0 = _mm256_loadu_ps(&b[stride*i +  0]);
        tmp0 = _mm256_add_ps(_mm256_mul_ps(areg0,breg0), tmp0);    
        __m256 breg1 = _mm256_loadu_ps(&b[stride*i +  8]);
        tmp1 = _mm256_add_ps(_mm256_mul_ps(areg0,breg1), tmp1);
        __m256 breg2 = _mm256_loadu_ps(&b[stride*i + 16]);
        tmp2 = _mm256_add_ps(_mm256_mul_ps(areg0,breg2), tmp2);    
        __m256 breg3 = _mm256_loadu_ps(&b[stride*i + 24]);
        tmp3 = _mm256_add_ps(_mm256_mul_ps(areg0,breg3), tmp3);   
        __m256 breg4 = _mm256_loadu_ps(&b[stride*i + 32]);
        tmp4 = _mm256_add_ps(_mm256_mul_ps(areg0,breg4), tmp4);    
        __m256 breg5 = _mm256_loadu_ps(&b[stride*i + 40]);
        tmp5 = _mm256_add_ps(_mm256_mul_ps(areg0,breg5), tmp5);    
        __m256 breg6 = _mm256_loadu_ps(&b[stride*i + 48]);
        tmp6 = _mm256_add_ps(_mm256_mul_ps(areg0,breg6), tmp6);    
        __m256 breg7 = _mm256_loadu_ps(&b[stride*i + 56]);
        tmp7 = _mm256_add_ps(_mm256_mul_ps(areg0,breg7), tmp7);    
    }
    _mm256_storeu_ps(&c[ 0], tmp0);
    _mm256_storeu_ps(&c[ 8], tmp1);
    _mm256_storeu_ps(&c[16], tmp2);
    _mm256_storeu_ps(&c[24], tmp3);
    _mm256_storeu_ps(&c[32], tmp4);
    _mm256_storeu_ps(&c[40], tmp5);
    _mm256_storeu_ps(&c[48], tmp6);
    _mm256_storeu_ps(&c[56], tmp7);
}

重新排序 b 矩阵的代码。

reorder(float *b, float *b2, int stride) {
    //int k = 0;
    for(int i=0; i<16; i++) {
        for(int j=0; j<16; j++) {
            for(int i2=0; i2<64; i2++) {
                for(int j2=0; j2<64; j2++) {
                    //b2[k++] = b[1024*(64*i+i2) + 64*j + j2];
                    b2[64*64*(i*16 + j) + 64*i2+j2] = b[1024*(64*i+i2) + 64*j + j2];
                }
            }
        }
    }
}

Answer 1

我认为问题出在最内层的循环中，并且与预取有关。您正在执行以下操作:

1. 从c中加载256个连续字节(4个缓存行)
2. 从 b 中加载 256 个连续字节 64 次(总共 256 个缓存行)
3. 保存256个连续字节到c

在第 2 步中，当步幅为 64 时，您正在有效地读取一个长的连续内存块。当步幅为 1024 时，您正在以不连续的步骤读取内存，一次读取 256 个字节。

因此，当 stride 为 64 时，预取器会提前将连续的 block 读取到缓存中，并且每行最多有一个缓存未命中。当步幅为 1024 时，预取器会被不连续的读取弄糊涂(您交替读取连续的缓存行和广泛分离的缓存行)并且每行有 64 次缓存未命中。