从两个数组的点积测量内存带宽

Question

两个数组的点积

for(int i=0; i<n; i++) {
    sum += x[i]*y[i];
}

不重用数据,因此它应该是一个内存绑定操作.因此,我应该能够从点积测量内存带宽.

使用代码 为什么 - 矢量化 - 循环 - 没有性能改进 我的系统带宽为9.3 GB/s.但是,当我尝试使用点积计算带宽时,我获得单个线程的速率的两倍以及使用多个线程的速率超过三倍(我的系统有四个核心/八个超线程).这对我没有意义,因为内存绑定操作不应该受益于多个线程.以下代码的输出如下:

Xeon E5-1620, GCC 4.9.0, Linux kernel 3.13
dot 1 thread:      1.0 GB, sum 191054.81, time 4.98 s, 21.56 GB/s, 5.39 GFLOPS
dot_avx 1 thread   1.0 GB, sum 191043.33, time 5.16 s, 20.79 GB/s, 5.20 GFLOPS
dot_avx 2 threads: 1.0 GB, sum 191045.34, time 3.44 s, 31.24 GB/s, 7.81 GFLOPS
dot_avx 8 threads: 1.0 GB, sum 191043.34, time 3.26 s, 32.91 GB/s, 8.23 GFLOPS

有人可以向我解释为什么我使用多个线程获得一个线程带宽的两倍以及带宽的三倍以上？

这是我使用的代码:

//g++ -O3 -fopenmp -mavx -ffast-math dot.cpp
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdint.h>
#include <x86intrin.h>
#include <omp.h>

extern "C" inline float horizontal_add(__m256 a) {
    __m256 t1 = _mm256_hadd_ps(a,a);
    __m256 t2 = _mm256_hadd_ps(t1,t1);
    __m128 t3 = _mm256_extractf128_ps(t2,1);
    __m128 t4 = _mm_add_ss(_mm256_castps256_ps128(t2),t3);
    return _mm_cvtss_f32(t4);
}

extern "C" float dot_avx(float * __restrict x, float * __restrict y, const int n) {
    x = (float*)__builtin_assume_aligned (x, 32);
    y = (float*)__builtin_assume_aligned (y, 32);
    float sum = 0;
    #pragma omp parallel reduction(+:sum)
    {
        __m256 sum1 = _mm256_setzero_ps();
        __m256 sum2 = _mm256_setzero_ps();
        __m256 sum3 = _mm256_setzero_ps();
        __m256 sum4 = _mm256_setzero_ps();
        __m256 x8, y8;
        #pragma omp for
        for(int i=0; i<n; i+=32) {
            x8 = _mm256_loadu_ps(&x[i]);
            y8 = _mm256_loadu_ps(&y[i]);
            sum1 = _mm256_add_ps(_mm256_mul_ps(x8,y8),sum1);
            x8 = _mm256_loadu_ps(&x[i+8]);
            y8 = _mm256_loadu_ps(&y[i+8]);
            sum2 = _mm256_add_ps(_mm256_mul_ps(x8,y8),sum2);
            x8 = _mm256_loadu_ps(&x[i+16]);
            y8 = _mm256_loadu_ps(&y[i+16]);
            sum3 = _mm256_add_ps(_mm256_mul_ps(x8,y8),sum3);
            x8 = _mm256_loadu_ps(&x[i+24]);
            y8 = _mm256_loadu_ps(&y[i+24]);
            sum4 = _mm256_add_ps(_mm256_mul_ps(x8,y8),sum4);
        }
        sum += horizontal_add(_mm256_add_ps(_mm256_add_ps(sum1,sum2),_mm256_add_ps(sum3,sum4)));
    }
    return sum; 
}

extern "C" float dot(float * __restrict x, float * __restrict y, const int n) {
    x = (float*)__builtin_assume_aligned (x, 32);
    y = (float*)__builtin_assume_aligned (y, 32);
    float sum = 0;
    for(int i=0; i<n; i++) {
        sum += x[i]*y[i];
    }
    return sum;
}

int main(){
    uint64_t LEN = 1 << 27;
    float *x = (float*)_mm_malloc(sizeof(float)*LEN,64);
    float *y = (float*)_mm_malloc(sizeof(float)*LEN,64);
    for(uint64_t i=0; i<LEN; i++) { x[i] = 1.0*rand()/RAND_MAX - 0.5; y[i] = 1.0*rand()/RAND_MAX - 0.5;}

    uint64_t size = 2*sizeof(float)*LEN;

    volatile float sum = 0;
    double dtime, rate, flops;  
    int repeat = 100;

    dtime = omp_get_wtime();
    for(int i=0; i<repeat; i++) sum += dot(x,y,LEN);
    dtime = omp_get_wtime() - dtime;
    rate = 1.0*repeat*size/dtime*1E-9;
    flops = 2.0*repeat*LEN/dtime*1E-9;
    printf("%f GB, sum %f, time %f s, %.2f GB/s, %.2f GFLOPS\n", 1.0*size/1024/1024/1024, sum, dtime, rate,flops);

    sum = 0;
    dtime = omp_get_wtime();
    for(int i=0; i<repeat; i++) sum += dot_avx(x,y,LEN);
    dtime = omp_get_wtime() - dtime;
    rate = 1.0*repeat*size/dtime*1E-9;
    flops = 2.0*repeat*LEN/dtime*1E-9;

    printf("%f GB, sum %f, time %f s, %.2f GB/s, %.2f GFLOPS\n", 1.0*size/1024/1024/1024, sum, dtime, rate,flops);
}

我刚刚按照Jonathan Dursi的建议下载,编译并运行了STREAM,结果如下:

一个线程

Function      Rate (MB/s)   Avg time     Min time     Max time
Copy:       14292.1657       0.0023       0.0022       0.0023
Scale:      14286.0807       0.0023       0.0022       0.0023
Add:        14724.3906       0.0033       0.0033       0.0033
Triad:      15224.3339       0.0032       0.0032       0.0032

八个线程

Function      Rate (MB/s)   Avg time     Min time     Max time
Copy:       24501.2282       0.0014       0.0013       0.0021
Scale:      23121.0556       0.0014       0.0014       0.0015
Add:        25263.7209       0.0024       0.0019       0.0056
Triad:      25817.7215       0.0020       0.0019       0.0027

Answer 1

这里有一些事情可以归结为:

• 你必须努力工作才能从内存子系统中获得最后一点性能; 和
• 不同的基准衡量不同的东西.

第一个有助于解释为什么需要多个线程来满足可用内存带宽.内存系统中存在大量并发性,并且利用它通常需要CPU代码中的一些并发性.多个执行线程帮助的一个重要原因是延迟隐藏 - 当一个线程停止等待数据到达时,另一个线程可能能够利用刚才可用的其他一些数据.

硬件可以帮助你在这种情况下,单个线程很多-因为内存访问是这样预测的,硬件可以提前预取的,当你需要它,给你一些延迟,即使一个线程隐藏的优势的数据; 但是预取可以做什么是有限的.例如,预取器不会将其自身跨越页面边界.对于大部分内容的规范参考是Ulrich Drepper对每个程序员应该了解的记忆,现在已经足够大,一些差距开始显现(英特尔的Sandy Bridge处理器的热芯片概述在这里 - 特别注意更紧密的集成与CPU的内存管理硬件).

至于与memset,mbw或STREAM进行比较的问题,在基准测试中进行比较总会引起令人头疼的问题,即使是那些声称测量相同的基准测试.特别是,"内存带宽"不是单个数字 - 性能根据操作的不同而有很大差异.mbw和Stream都执行某种版本的复制操作,此处拼写了STREAMs操作(直接从网页中获取,所有操作数都是双精度浮点数):

------------------------------------------------------------------
name        kernel                  bytes/iter      FLOPS/iter
------------------------------------------------------------------
COPY:       a(i) = b(i)                 16              0
SCALE:      a(i) = q*b(i)               16              1
SUM:        a(i) = b(i) + c(i)          24              1
TRIAD:      a(i) = b(i) + q*c(i)        24              2
------------------------------------------------------------------

因此在这些情况下大约1/2-1/3的内存操作是写入(在memset的情况下,所有内容都是写入).虽然单个写入可能比读取慢一点,但更大的问题是使用写入使内存子系统饱和要困难得多,因为当然你不能完全相同地预取写入.交叉读取和写入有所帮助,但是基本上所有读取的点积示例将是针对内存带宽的最佳可能情况.

此外,STREAM基准测试是(有意)完全可移植地编写的,只有一些编译器编译指示建议向量化,因此击败STREAM基准测试不一定是警告信号,特别是当您正在做的是两个流式读取时.