c++ 從兩個陣列的點積測量儲存器頻寬

兩個陣列的點積

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

不重用資料,因此應該是記憶體繫結操作.因此,我應該能夠從點積測量記憶體頻寬.

使用程式碼

ofollow,noindex" target="_blank">why-vectorizing-the-loop-does-not-have-performance-improvement 我的系統頻寬為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 thread1.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,結果如下：

一線

FunctionRate (MB/s)Avg timeMin timeMax time
Copy:14292.16570.00230.00220.0023
Scale:14286.08070.00230.00220.0023
Add:14724.39060.00330.00330.0033
Triad:15224.33390.00320.00320.0032

八個執行緒

FunctionRate (MB/s)Avg timeMin timeMax time
Copy:24501.22820.00140.00130.0021
Scale:23121.05560.00140.00140.0015
Add:25263.72090.00240.00190.0056
Triad:25817.72150.00200.00190.0027

這裡有一些事情可以歸結為：

>你必須非常努力地從記憶體子系統中獲得最後一點的效能;和

不同的基準測量不同的東西.

第一個有助於解釋為什麼需要多個執行緒來飽和可用的記憶體頻寬.在記憶體系統中存在大量的併發性,並且利用這一點,您的CPU程式碼通常需要一些併發.多個執行執行緒幫助的一個重要原因是latency hiding – 一個執行緒停止等待資料到達,另一個執行緒可能可以利用剛剛可用的其他一些資料.

在這種情況下,硬體可以幫助您在一個單一的執行緒上 – 因為記憶體訪問是如此可預測,硬體可以在需要時預取資料,為您提供即使使用一個執行緒的延遲隱藏的一些優勢;但是預取可以做什麼有限制.例如,預取器將不會將其自身跨越頁面邊界.許多這方面的規範參考是What Every Programmer Should Know About Memory by Ulrich Drepper ,現在已經足夠了,有一些差距開始顯現(英特爾酷睿晶片Sandy Bridge處理器的概述是here – 尤其注意記憶體管理硬體與CPU的更緊密的整合) .

關於與memset,mbw 或STREAM 比較的問題,比較基準將總是引起頭痛,即使是聲稱測量同樣的事情的基準.特別地,“儲存器頻寬”不是單個數字 – 根據操作而不同,效能變化很大. mbw和Stream都執行某些版本的複製操作,STREAMs操作在此被拼寫(直接從網頁中獲取,所有運算元都是雙精度浮點數)：

------------------------------------------------------------------
namekernelbytes/iterFLOPS/iter
------------------------------------------------------------------
COPY:a(i) = b(i)160
SCALE:a(i) = q*b(i)161
SUM:a(i) = b(i) + c(i)241
TRIAD:a(i) = b(i) + q*c(i)242
------------------------------------------------------------------

所以在這些情況下大概是1 / 2-1 / 3的記憶體操作是寫入的(而且在memset的情況下都是寫入的).雖然單獨的寫入可能比讀取慢一點,但更大的問題是,使用寫入使記憶體子系統飽和很難,因為您當然不能等同於預取寫入.交錯閱讀和寫入有助於,但是您的點陣產品例項基本上都是讀取將是關於針對記憶體頻寬的最佳可能情況.

此外,STREAM基準測試(有意)完全可移植,僅使用一些編譯器編譯指示來建議向量化,因此,打擊STREAM基準測試不一定是一個警告訊號,特別是當您正在做的是兩個流讀取.