{var%20f='http://v.t.sina.com.cn/share/share.php?appkey=1515056452',u=z||d.location,p=['&url=',e(u),'&title=',e(t||d.title),'&source=',e(r),'&sourceUrl=',e(l),'&content=',c||'gb2312','&pic=',e(p||'')].join('');function%20a(){if(!window.open([f,p].join(''),'mb',['toolbar=0,status=0,resizable=1,width=440,height=430,left=',(s.width-440)/2,',top=',(s.height-430)/2].join('')))u.href=[f,p].join('');};if(/Firefox/.test(navigator.userAgent))setTimeout(a,0);else%20a();})(screen,document,encodeURIComponent,'','','https://www.manongdao.com/data/attach/logo/logo.png', '推荐 萌系小妹纸 的问题《L1 memory bandwidth: 50% drop in efficiency using》','https://www.manongdao.com/q-196204.html','页面编码gb2312|utf-8默认gb2312'));){kind=link}
I want to achieve the maximum bandwidth of the following operations with Intel processors.
for(int i=0; i<n; i++) z[i] = x[i] + y[i]; //n=2048
where x, y, and z are float arrays. I am doing this on Haswell, Ivy Bridge , and Westmere systems.
I originally allocated the memory like this
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;
When I did this I got about 50% of the peak bandwidth I expected for each system.
The peak values is calculated as frequency * average bytes/clock_cycle. The average bytes/clock cycle for each system is:
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
This means that e.g. on Haswell I I only observe 48 bytes/clock cycle (could be two reads in one clock cycle and one write the next clock cycle).
I printed out the difference in the address of b-a and c-b and each are 8256 bytes. The value 8256 is 8192+64. So they are each larger than the array size (8192 bytes) by one cache-line.
On a whim I tried allocating the memory like this.
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;
This nearly doubled my peak bandwidth so that I now get around 90% of the peak bandwidth. However, when I tried k=1 it dropped back to 50%. I have tried other values of k and found that e.g. k=2, k=33, k=65 only gets 50% of the peak but e.g. k=10, k=32, k=63 gave the full speed. I don't understand this.
In Agner Fog's micrarchitecture manual he says that there is a false dependency with memory address with the same set and offset
It is not possible to read and write simultaneously from addresses that are spaced by a multiple of 4 Kbytes.
But that's exactly where I see the biggest benefit! When k=0 the memory address differ by exactly 2*4096 bytes. Agner also talks about Cache bank conflicts. But Haswell and Westmere are not suppose to have these bank conflicts so that should not explain what I am observing. What's going on!?
I understand that the OoO execution decides which address to read and write so even if the arrays' memory addresses differ by exactly 4096 bytes that does not necessarily mean the processor reads e.g. &x[0] and writes &z[0] at the same time but then why would being off by a single cache line cause it to choke?
Edit: Based on Evgeny Kluev's answer I now believe this is what Agner Fog calls a "bogus store forwarding stall". In his manual under the Pentium Pro, II and II he writes:
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
Edit: Here is table of the efficiencies on each system for k=0 and k=1.
k=0 k=1
Westmere: 99% 66%
Ivy Bridge: 98% 44%
Haswell: 90% 49%
I think I can explain these numbers if I assume that for k=1 that writes and reads cannot happen in the same clock cycle.
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%
This theory works out pretty well. Ivy bridge is a bit lower than I would expect but Ivy Bridge suffers from bank cache conflicts where the others don't so that may be another effect to consider.
Below is working code to test this yourself. On a system without AVX compile with g++ -O3 sum.cpp otherwise compile with g++ -O3 -mavx sum.cpp. Try varying the value 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);
}
}
I think the gap between
aandbdoes not really matter. After leaving only one gap betweenbandcI've got the following results on Haswell:Since Haswell is known to be free of bank conflicts, the only remaining explanation is false dependence between memory addresses (and you've found proper place in Agner Fog's microarchitecture manual explaining exactly this problem). The difference between bank conflict and false sharing is that bank conflict prevents accessing the same bank twice during the same clock cycle while false sharing prevents reading from some offset in 4K piece of memory just after you've written something to same offset (and not only during the same clock cycle but also for several clock cycles after the write).
Since your code (for
k=0) writes to any offset just after doing two reads from the same offset and would not read from it for a very long time, this case should be considered as "best", so I placedk=0at the end of the table. Fork=1you always read from offset that is very recently overwritten, which means false sharing and therefore performance degradation. With largerktime between write and read increases and CPU core has more chances to pass written data through all memory hierarchy (which means two address translations for read and write, updating cache data and tags and getting data from cache, data synchronization between cores, and probably many more stuff).k=12or 24 clocks (on my CPU) is enough for every written piece of data to be ready for subsequent read operations, so starting with this value performance gets back to usual. Looks not very different from 20+ clocks on AMD (as said by @Mysticial).First, observe that the three arrays have a total size of 24KB. In addition, since you're initializing the arrays before executing the main loop, most accesses in the main loop will hit into the L1D, which is 32KB in size and 8-way associative on modern Intel processors. So we don't have to worry about misses or hardware prefetching. The most important performance event in this case is
LD_BLOCKS_PARTIAL.ADDRESS_ALIAS, which occurs when a partial address comparison involving a later load results in a match with an earlier store and all of the conditions of store forwarding are satisfied, but the target locations are actually different. Intel refers to this situation as 4K aliasing or false store forwarding. The observable performance penalty of 4K aliasing depends on the surrounding code.By measuring
cycles,LD_BLOCKS_PARTIAL.ADDRESS_ALIASandMEM_UOPS_RETIRED.ALL_LOADS, we can see that for all values ofkwhere the achieved bandwidth is much smaller than the peak bandwidth,LD_BLOCKS_PARTIAL.ADDRESS_ALIASandMEM_UOPS_RETIRED.ALL_LOADSare almost equal. Also for all values ofkwhere the achieved bandwidth is close to the peak bandwidth,LD_BLOCKS_PARTIAL.ADDRESS_ALIASis very small compared toMEM_UOPS_RETIRED.ALL_LOADS. This confirms that bandwidth degradation is occurring due to most loads suffering from 4K aliasing.The Intel optimization manual Section 12.8 says the following:
That is, there are two conditions for a later load to alias with an earlier store:
Only then 4K aliasing may occur.
However, this appears to not be accurate. The following listing shows the AVX-based (32-byte) sequence of memory accesses and the least significant 12 bits of their addresses for different values of
k.On my Haswell processor, the efficiencies for different values of
kare the following:where 86% is the highest and 48% is the lowest.
Examine bits 5-11 of the accessed locations for the three different cases (different values of
k). Whenkis 0 or 10 (or larger I think), there is no load whose bits 5-11 of its address match that of any recent previous store. This violates the first condition of 4K aliasing. The results are consistent with the manual. However, whenkis between 0 and 10, it's difficult to see what's happening. We need a simulator.I wrote a program that basically emulates the addresses of the memory accesses and calculates the total number of loads that suffered 4K aliasing for different values of
k. One issue I faced was that we don't know, for any given load, the number of stores that are still in the store buffer (have not been committed yet). Therefore, to perform a comprehensive analysis, I tried with different store buffer occupancies. I assume that at any point in time, there are exactly a fixed number of stores in the store buffer. So the current load needs to be compared against all of these stores.You can find the raw results here. You can then create an interactive 3D plot of the results by simply copying them into the data box here and click on the 'Draw Graph" button.
By using the
LD_BLOCKS_PARTIAL.ADDRESS_ALIASperformance counter, we have determined that the bandwidth is sub-optimal forklarger than 0 and smaller than 10. Hence, we need to find the smallest store buffer occupancy that achieves this behavior. From the results, it turns out that it is 18. That is, the most recent 18 stores need to be in the store buffer so that when k=10 or k=0, the L1D bandwidth is optimal and sub-optimal otherwise.Note that the total number of 32-byte loads is 512.
The change in
aliasCountis consistent with the change inLD_BLOCKS_PARTIAL.ADDRESS_ALIAS, but it changes much slower. That is,aliasCountdecreases until it reaches the optimal value whenk= 10, but there is a steep drop fromk= 9.@BeeOnRope suggested to:
occupancyparameter.ADDRESS_ALIASrather than due to actual 4K aliasing.By implementing both (the code is here), the results are as follows:
Much better! Although there seems that there is still room for improvement in the simulator.
The following graph plots
ADDRESS_ALIAS / MEM_LOADSvs.kin three cases:_mm256_store_ps(z1+64*i+ 0,_mm256_add_ps(_mm256_load_ps(x1+64*i+ .), _mm256_load_ps(y1+64*i+ .)));._mm256_store_ps(z1+64*i+ 0,_mm256_add_ps(_mm256_load_ps(x1+64*i+ .), _mm256_set_ps(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)));. In this case,MEM_LOADSis halved._mm256_store_ps(z1+64*i+ 0,_mm256_add_ps(_mm256_set_ps(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0), _mm256_load_ps(y1+64*i+ .)));. In this case,MEM_LOADSis halved.The weird thing here is that the loads from
yshould never alias any store, butADDRESS_ALIASshows otherwise.I've improved the simulator so that it can simulate different store throughputs for different values of
k, which seems to reflect very well what's actually happening. An explanation of why the store throughputs might be so wildly different even though most of them hit in the L1 is still needed.The following store throughputs were used:
For example, a value of 0.83 means that the probability that a store will leave the store buffer (be committed) in any given cycle is 83%.
I've conducted additional experiments to determine whether the two conditions for 4K aliasing to occur as described in the Intel manual are accurate and whether
LD_BLOCKS_PARTIAL.ADDRESS_ALIASworks as expected. Consider the following loop body:LD_BLOCKS_PARTIAL.ADDRESS_ALIASin this case increases with the number of iterations, indicating that there is aliasing. To determine which of the two loads alias with a previous store, the same experiment can be repeat with one of the loads removed from the code. Without the first load,LD_BLOCKS_PARTIAL.ADDRESS_ALIASstill increase with the number of iterations. In addition, the number of events is a little less than the number of loads. Without the second load, however,LD_BLOCKS_PARTIAL.ADDRESS_ALIASis fixed at some small value (with high standard deviation) no matter what the number of iterations is. This indicates that the second load is responsible for most or all of the aliasing events. Let's examine the addresses:Remember that the size of each access is 8 bytes. Consider the first store-load-load group. The store targets the 8-byte range 0000 0000 1001 - 0000 0001 0000. The first load targets the 8-byte range 0000 0001 0000 - 0000 0001 0111. These two ranges overlap (bits 5-11 match and bits 0-4 overlap). The second load targets the 8-byte range 0000 0100 1000 - 0000 0100 1111. The range of this load does not overlap with the range of the previous store. In fact, it does not overlap with any previous store to the same page. This reasoning applies to all other values of bits 0-11 of
r13except when an access crosses a cache line or page boundary.TBC.