Menu
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);
}
}
997
I think the gap between a and b does not really matter. After leaving only one gap between b and c I've got the following results on 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
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 placed k=0 at the end of the table. For k=1 you always read from offset that is very recently overwritten, which means false sharing and therefore performance degradation. With larger k time 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=12 or 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).
147
TL;DR: For certain values of k, too many 4K aliasing conditions occur, which is the main cause for the bandwidth degradation. In 4K aliasing, a load is stalled unnecessarily, thereby increasing the effective load latency and stalling all later dependent instructions. This in turn results in reduced L1 bandwidth utilization. For these values of k, most 4K aliasing conditions can be eliminated by splitting the loop as follows:
for(int i=0; i<n/64; i++) {
_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)));
}
for(int i=0; i<n/64; i++) {
_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)));
}
for(int i=0; i<n/64; i++) {
_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)));
}
for(int i=0; i<n/64; i++) {
_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)));
}
This split eliminates most 4K aliasing for the cases when k is an odd positive integer (such as 1). The achieved L1 bandwidth is improved by about 50% on Haswell. There is still room for improvement, for example, by unrolling the loop and figuring out a way to not use the indexed-addressing mode for loads and stores.
However, this split doesn't eliminate 4K aliasing for even values of k. So a different split needs to be used for even values of k. However, when k is 0, optimal performance can be achieved without splitting the loop. In this case, performance is backend-bound on ports 1, 2, 3, 4, and 7 simultaneously.
There could be a penalty of a few cycles in certain cases when performing a load and store at the same time, but in this particular case, this penalty basically does not exist because there are basically no such conflicts (i.e., the addresses of concurrent loads and stores are sufficiently far apart). In addition, the total working set size fits in the L1 so there is no L1-L2 traffic beyond the first execution of the loop.
The rest of this answer includes a detailed explanation of this summary.
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_ALIAS and MEM_UOPS_RETIRED.ALL_LOADS, we can see that for all values of k where the achieved bandwidth is much smaller than the peak bandwidth, LD_BLOCKS_PARTIAL.ADDRESS_ALIAS and MEM_UOPS_RETIRED.ALL_LOADS are almost equal. Also for all values of k where the achieved bandwidth is close to the peak bandwidth, LD_BLOCKS_PARTIAL.ADDRESS_ALIAS is very small compared to MEM_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:
4-KByte memory aliasing occurs when the code stores to one memory location and shortly after that it loads from a different memory location with a 4-KByte offset between them. For example, a load to linear address 0x400020 follows a store to linear address 0x401020.
The load and store have the same value for bits 5 - 11 of their addresses and the accessed byte offsets should have partial or complete overlap.
That is, there are two necessary conditions for a later load to alias with an earlier store:
On processors that support AVX-512, it seems to me that a single load uop can load up to 64 bytes. So I think the range for the first condition should be 6-11 instead of 5-11.
The following listing shows the AVX-based (32-byte) sequence of memory accesses and the least significant 12 bits of their addresses for two different values of k.
======
k=0
======
load x+(0*64+0)*4 = x+0 where x is 4k aligned 0000 000|0 0000
load y+(0*64+0)*4 = y+0 where y is 4k aligned 0000 000|0 0000
store z+(0*64+0)*4 = z+0 where z is 4k aligned 0000 000|0 0000
load x+(0*64+8)*4 = x+32 where x is 4k aligned 0000 001|0 0000
load y+(0*64+8)*4 = y+32 where y is 4k aligned 0000 001|0 0000
store z+(0*64+8)*4 = z+32 where z is 4k aligned 0000 001|0 0000
load x+(0*64+16)*4 = x+64 where x is 4k aligned 0000 010|0 0000
load y+(0*64+16)*4 = y+64 where y is 4k aligned 0000 010|0 0000
store z+(0*64+16)*4= z+64 where z is 4k aligned 0000 010|0 0000
load x+(0*64+24)*4 = x+96 where x is 4k aligned 0000 011|0 0000
load y+(0*64+24)*4 = y+96 where y is 4k aligned 0000 011|0 0000
store z+(0*64+24)*4 = z+96 where z is 4k aligned 0000 011|0 0000
load x+(0*64+32)*4 = x+128 where x is 4k aligned 0000 100|0 0000
load y+(0*64+32)*4 = y+128 where y is 4k aligned 0000 100|0 0000
store z+(0*64+32)*4= z+128 where z is 4k aligned 0000 100|0 0000
.
.
.
======
k=1
======
load x+(0*64+0)*4 = x+0 where x is 4k aligned 0000 000|0 0000
load y+(0*64+0)*4 = y+0 where y is 4k+64 aligned 0000 010|0 0000
store z+(0*64+0)*4 = z+0 where z is 4k+128 aligned 0000 100|0 0000
load x+(0*64+8)*4 = x+32 where x is 4k aligned 0000 001|0 0000
load y+(0*64+8)*4 = y+32 where y is 4k+64 aligned 0000 011|0 0000
store z+(0*64+8)*4 = z+32 where z is 4k+128 aligned 0000 101|0 0000
load x+(0*64+16)*4 = x+64 where x is 4k aligned 0000 010|0 0000
load y+(0*64+16)*4 = y+64 where y is 4k+64 aligned 0000 100|0 0000
store z+(0*64+16)*4= z+64 where z is 4k+128 aligned 0000 110|0 0000
load x+(0*64+24)*4 = x+96 where x is 4k aligned 0000 011|0 0000
load y+(0*64+24)*4 = y+96 where y is 4k+64 aligned 0000 101|0 0000
store z+(0*64+24)*4 = z+96 where z is 4k+128 aligned 0000 111|0 0000
load x+(0*64+32)*4 = x+128 where x is 4k aligned 0000 100|0 0000
load y+(0*64+32)*4 = y+128 where y is 4k+64 aligned 0000 110|0 0000
store z+(0*64+32)*4= z+128 where z is 4k+128 aligned 0001 000|0 0000
.
.
.
Note that when k=0, no load seem satisfy the two conditions of 4K aliasing. On the other hand, when k=1, all loads seem satisfy the conditions. However, it's tedious to do this manually for all iterations and all values of k. So I wrote a program that basically generates 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, I've designed the simulator so that it can use different store throughputs for different values of k, which seems to better reflect what's actually happening on a real processor. The code can be found here.
The following figure show the number of 4K aliasing cases produced by the simulator compared to the measured number using LD_BLOCKS_PARTIAL.ADDRESS_ALIAS on Haswell. I've tuned the store throughput used in the simulator for each value of k to make the two curves as similar as possible. The second figure shows the inverse store throughput (total cycles divided by total number of stores) used in the simulator and measured on Haswell. Note that the store throughput when k=0 doesn't matter because there is no 4K aliasing anyway. Since there are two loads for each store, the inverse load throughput is half of the inverse store throughput.
Obviously the amount of time each store remains in the store buffer is different on Haswell and the simulator, so I needed to use different throughputs to make the two curves similar. The simulator can be used to show how the store throughput can impact the number of 4K aliases. If the store throughput is very close to 1c/store, then the number of 4K aliasing cases would have been much smaller. 4K aliasing conditions don't result in pipeline flushes, but they may result in uop replays from the RS. In this particular case, I didn't observe any replays though.
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.
There is actually a penalty of a few cycles when executing a load and store at the same time, but they can only happen when the addresses of the load and store are within 64 bytes (but not equal) on Haswell or 32 bytes on Ivy Bridge and Sandy Bridge. Weird performance effects from nearby dependent stores in a pointer-chasing loop on IvyBridge. Adding an extra load speeds it up?. In this case, the addresses of all accesses are 32-byte aligned, but, on IvB, the L1 ports are all 16-byte in size, so the penalty can be incurred on Haswell and IvB. In fact, since loads and stores may take more time to retire and since there are more load buffers than store buffers, it's more likely that a later load will false-alias an earlier store. This raises the question, though, how the 4K alias penalty and the L1 access penalty interact with each other and contribute to the overall performance. Using the CYCLE_ACTIVITY.STALLS_LDM_PENDING event and the load latency performance monitoring facility MEM_TRANS_RETIRED.LOAD_LATENCY_GT_*, it seems to me that there is no observable L1 access penalty. This implies that most of the time the addresses of concurrent loads and stores don't induce the penalty. Hence, the 4K aliasing penalty is the main cause for bandwidth degradation.
I've used the following code to make measurements on Haswell. This is essentially the same code emitted by g++ -O3 -mavx.
%define SIZE 64*64*2
%define K_ 10
BITS 64
DEFAULT REL
GLOBAL main
EXTERN printf
EXTERN exit
section .data
align 4096
bufsrc1: times (SIZE+(64*K_)) db 1
bufsrc2: times (SIZE+(64*K_)) db 1
bufdest: times SIZE db 1
section .text
global _start
_start:
mov rax, 1000000
.outer:
mov rbp, SIZE/256
lea rsi, [bufsrc1]
lea rdi, [bufsrc2]
lea r13, [bufdest]
.loop:
vmovaps ymm1, [rsi]
vaddps ymm0, ymm1, [rdi]
add rsi, 256
add rdi, 256
add r13, 256
vmovaps[r13-256], ymm0
vmovaps ymm2, [rsi-224]
vaddps ymm0, ymm2, [rdi-224]
vmovaps [r13-224], ymm0
vmovaps ymm3, [rsi-192]
vaddps ymm0, ymm3, [rdi-192]
vmovaps [r13-192], ymm0
vmovaps ymm4, [rsi-160]
vaddps ymm0, ymm4, [rdi-160]
vmovaps [r13-160], ymm0
vmovaps ymm5, [rsi-128]
vaddps ymm0, ymm5, [rdi-128]
vmovaps [r13-128], ymm0
vmovaps ymm6, [rsi-96]
vaddps ymm0, ymm6, [rdi-96]
vmovaps [r13-96], ymm0
vmovaps ymm7, [rsi-64]
vaddps ymm0, ymm7, [rdi-64]
vmovaps [r13-64], ymm0
vmovaps ymm1, [rsi-32]
vaddps ymm0, ymm1, [rdi-32]
vmovaps [r13-32], ymm0
dec rbp
jg .loop
dec rax
jg .outer
xor edi,edi
mov eax,231
syscall
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
