我正在尝试在 Intel 处理器上为以下功能获取 L1 缓存中的全部带宽
float triad(float *x, float *y, float *z, const int n) {
float k = 3.14159f;
for(int i=0; i<n; i++) {
z[i] = x[i] + k*y[i];
}
}
这是来自STREAM的三元组函数。
使用具有此功能的 SandyBridge/IvyBridge 处理器（使用带有 NASM 的汇编），我得到了大约 95% 的峰值。但是，除非我展开循环，否则使用 Haswell 只能达到峰值的 62%。如果我展开 16 次，我会得到 92%。我不明白这一点。
我决定使用 NASM 在汇编中编写我的函数。汇编中的主循环如下所示。
.L2:
vmovaps ymm1, [rdi+rax]
vfmadd231ps ymm1, ymm2, [rsi+rax]
vmovaps [rdx+rax], ymm1
add rax, 32
jne .L2
事实证明，在Agner Fog 的 Optimizing Assembly 手册y[i] = y[i] +k*x[i]中的示例 12.7-12.11 中，他对 Pentium M、Core 2、Sandy Bridge、FMA4 和 FMA3做了几乎相同的事情（但针对）。我设法或多或少地自己重现了他的代码（实际上他在广播时在 FMA3 示例中有一个小错误）。他给出了除 FMA4 和 FMA3 之外的每个处理器的指令大小计数、融合操作、执行端口。我试图自己为 FMA3 制作这张桌子。
ports
size μops-fused 0 1 2 3 4 5 6 7
vmovaps 5 1 ½ ½
vfmadd231ps 6 1 ½ ½ ½ ½
vmovaps 5 1 1 1
add 4 ½ ½
jne 2 ½ ½
--------------------------------------------------------------
total 22 4 ½ ½ 1 1 1 0 1 1
大小是指以字节为单位的指令长度。add和指令有半个微操作的原因jne是它们被融合到一个宏操作中（不要与仍然使用多个端口的微操作融合混淆）并且只需要端口6和一个微操作。 指令可以使用。。为了与 Agner Fog 的表格保持一致，并且因为我认为说一条可以到达不同端口的指令在 1/2 的时间里平等地到达每个端口更有意义，所以我为端口分配了 1/2vfmadd231ps端口 0 或端口 1。我选择了端口 0。负载vmovaps可以使用端口 2 或 3。我选择了 2 并vfmadd231ps使用了端口 3vmovaps并且vmadd231ps可以去至。
基于此表以及所有 Core2 处理器每个时钟周期都可以执行 4 μops 的事实，看起来这个循环应该每个时钟周期都可以，但我还没有设法获得它。有人可以向我解释为什么不展开就无法在 Haswell 上接近此功能的峰值带宽吗？如果不展开，这可能吗？如果可以，怎么做？让我明确一点，我真的在尝试最大化此功能的 ILP（我不仅想要最大带宽），所以这就是我不想展开的原因。
编辑： 这是一个更新，因为 Iwillnotexist Idonotexist 使用 IACA 显示商店从不使用端口 7。我设法在不展开的情况下打破 66% 的障碍，并且在每次迭代的一个时钟周期内完成此操作而无需展开（理论上）。让我们首先解决商店问题。
Stephen Canon 在评论中提到，端口 7 中的地址生成单元 (AGU) 只能处理简单的操作，例如[base + offset]and not [base + index]。在英特尔优化参考手册中，我发现的唯一内容是对 port7 的评论，上面写着“Simple_AGU”，但没有定义简单的含义。但后来 Iwillnotexist Idonotexist 在IACA的评论中发现，这个问题在 6 个月前已经提到过，英特尔的一名员工在 2014 年 3 月 11 日写道：
Port7 AGU 只能在具有简单内存地址（无索引寄存器）的存储上工作。
Stephen Canon 建议“使用存储地址作为加载操作数的偏移量”。我试过这样
vmovaps ymm1, [rdi + r9 + 32*i]
vfmadd231ps ymm1, ymm2, [rsi + r9 + 32*i]
vmovaps [r9 + 32*i], ymm1
add r9, 32*unroll
cmp r9, rcx
jne .L2
这确实导致商店使用端口7。但是，它还有另一个问题，vmadd231ps即您可以从 IACA 看到的负载不熔断。它还需要cmp我原来的功能没有的指令。所以商店少用了一个微操作，但是cmp（或者更确切地说，add因为cmp宏与 融合jne）需要一个。IACA 报告块吞吐量为 1.5。在实践中，这只得到大约 57% 的峰值。
但我找到了一种让vmadd231ps指令也与负载融合的方法。这只能使用像这样寻址 [绝对 32 位地址 + 索引] 的静态数组来完成。Evgeny Kluev 最初提出了这一点。
vmovaps ymm1, [src1_end + rax]
vfmadd231ps ymm1, ymm2, [src2_end + rax]
vmovaps [dst_end + rax], ymm1
add rax, 32
jl .L2
其中src1_end、src2_end和dst_end是静态数组的结束地址。
这用我期望的四个融合微操作重现了我的问题中的表格。如果将其放入 IACA，它会报告 1.0 的块吞吐量。从理论上讲，这应该与 SSE 和 AVX 版本一样好。在实践中，它得到了峰值的 72% 左右。这打破了 66% 的障碍，但距离我展开 16 次的 92% 仍有很长的路要走。因此，在 Haswell 上，接近顶峰的唯一选择就是展开。通过 Ivy Bridge 在 Core2 上这不是必需的，但在 Haswell 上是必需的。
结束编辑：
这是用于测试它的 C/C++ Linux 代码。NASM 代码发布在 C/C++ 代码之后。您唯一需要更改的是频率编号。在该行中double frequency = 1.3;，将 1.3 替换为处理器的工作（非标称）频率（如果 i5-4250U 在 BIOS 中禁用涡轮增压，则为 1.3 GHz）。
编译
nasm -f elf64 triad_sse_asm.asm
nasm -f elf64 triad_avx_asm.asm
nasm -f elf64 triad_fma_asm.asm
g++ -m64 -lrt -O3 -mfma tests.cpp triad_fma_asm.o -o tests_fma
g++ -m64 -lrt -O3 -mavx tests.cpp triad_avx_asm.o -o tests_avx
g++ -m64 -lrt -O3 -msse2 tests.cpp triad_sse_asm.o -o tests_sse
C/C++ 代码
#include <x86intrin.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#define TIMER_TYPE CLOCK_REALTIME
extern "C" float triad_sse_asm_repeat(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_sse_asm_repeat_unroll16(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_avx_asm_repeat(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_avx_asm_repeat_unroll16(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_fma_asm_repeat(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_fma_asm_repeat_unroll16(float *x, float *y, float *z, const int n, int repeat);
#if (defined(__FMA__))
float triad_fma_repeat(float *x, float *y, float *z, const int n, int repeat) {
float k = 3.14159f;
int r;
for(r=0; r<repeat; r++) {
int i;
__m256 k4 = _mm256_set1_ps(k);
for(i=0; i<n; i+=8) {
_mm256_store_ps(&z[i], _mm256_fmadd_ps(k4, _mm256_load_ps(&y[i]), _mm256_load_ps(&x[i])));
}
}
}
#elif (defined(__AVX__))
float triad_avx_repeat(float *x, float *y, float *z, const int n, int repeat) {
float k = 3.14159f;
int r;
for(r=0; r<repeat; r++) {
int i;
__m256 k4 = _mm256_set1_ps(k);
for(i=0; i<n; i+=8) {
_mm256_store_ps(&z[i], _mm256_add_ps(_mm256_load_ps(&x[i]), _mm256_mul_ps(k4, _mm256_load_ps(&y[i]))));
}
}
}
#else
float triad_sse_repeat(float *x, float *y, float *z, const int n, int repeat) {
float k = 3.14159f;
int r;
for(r=0; r<repeat; r++) {
int i;
__m128 k4 = _mm_set1_ps(k);
for(i=0; i<n; i+=4) {
_mm_store_ps(&z[i], _mm_add_ps(_mm_load_ps(&x[i]), _mm_mul_ps(k4, _mm_load_ps(&y[i]))));
}
}
}
#endif
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;
}
int main () {
int bytes_per_cycle = 0;
double frequency = 1.3; //Haswell
//double frequency = 3.6; //IB
//double frequency = 2.66; //Core2
#if (defined(__FMA__))
bytes_per_cycle = 96;
#elif (defined(__AVX__))
bytes_per_cycle = 48;
#else
bytes_per_cycle = 24;
#endif
double peak = frequency*bytes_per_cycle;
const int n =2048;
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);
char *c = b+n*sizeof(float);
float *x = (float*)a;
float *y = (float*)b;
float *z = (float*)c;
for(int i=0; i<n; i++) {
x[i] = 1.0f*i;
y[i] = 1.0f*i;
z[i] = 0;
}
int repeat = 1000000;
timespec time1, time2;
#if (defined(__FMA__))
triad_fma_repeat(x,y,z2,n,repeat);
#elif (defined(__AVX__))
triad_avx_repeat(x,y,z2,n,repeat);
#else
triad_sse_repeat(x,y,z2,n,repeat);
#endif
while(1) {
double dtime, rate;
clock_gettime(TIMER_TYPE, &time1);
#if (defined(__FMA__))
triad_fma_asm_repeat(x,y,z,n,repeat);
#elif (defined(__AVX__))
triad_avx_asm_repeat(x,y,z,n,repeat);
#else
triad_sse_asm_repeat(x,y,z,n,repeat);
#endif
clock_gettime(TIMER_TYPE, &time2);
dtime = time_diff(time1,time2);
rate = 3.0*1E-9*sizeof(float)*n*repeat/dtime;
printf("unroll1 rate %6.2f GB/s, efficency %6.2f%%, error %d\n", rate, 100*rate/peak, memcmp(z,z2, sizeof(float)*n));
clock_gettime(TIMER_TYPE, &time1);
#if (defined(__FMA__))
triad_fma_repeat(x,y,z,n,repeat);
#elif (defined(__AVX__))
triad_avx_repeat(x,y,z,n,repeat);
#else
triad_sse_repeat(x,y,z,n,repeat);
#endif
clock_gettime(TIMER_TYPE, &time2);
dtime = time_diff(time1,time2);
rate = 3.0*1E-9*sizeof(float)*n*repeat/dtime;
printf("intrinsic rate %6.2f GB/s, efficency %6.2f%%, error %d\n", rate, 100*rate/peak, memcmp(z,z2, sizeof(float)*n));
clock_gettime(TIMER_TYPE, &time1);
#if (defined(__FMA__))
triad_fma_asm_repeat_unroll16(x,y,z,n,repeat);
#elif (defined(__AVX__))
triad_avx_asm_repeat_unroll16(x,y,z,n,repeat);
#else
triad_sse_asm_repeat_unroll16(x,y,z,n,repeat);
#endif
clock_gettime(TIMER_TYPE, &time2);
dtime = time_diff(time1,time2);
rate = 3.0*1E-9*sizeof(float)*n*repeat/dtime;
printf("unroll16 rate %6.2f GB/s, efficency %6.2f%%, error %d\n", rate, 100*rate/peak, memcmp(z,z2, sizeof(float)*n));
}
}
使用 System V AMD64 ABI 的 NASM 代码。
triad_fma_asm.asm：
global triad_fma_asm_repeat
;RDI x, RSI y, RDX z, RCX n, R8 repeat
;z[i] = y[i] + 3.14159*x[i]
pi: dd 3.14159
;align 16
section .text
triad_fma_asm_repeat:
shl rcx, 2
add rdi, rcx
add rsi, rcx
add rdx, rcx
vbroadcastss ymm2, [rel pi]
;neg rcx
align 16
.L1:
mov rax, rcx
neg rax
align 16
.L2:
vmovaps ymm1, [rdi+rax]
vfmadd231ps ymm1, ymm2, [rsi+rax]
vmovaps [rdx+rax], ymm1
add rax, 32
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
global triad_fma_asm_repeat_unroll16
section .text
triad_fma_asm_repeat_unroll16:
shl rcx, 2
add rcx, rdi
vbroadcastss ymm2, [rel pi]
.L1:
xor rax, rax
mov r9, rdi
mov r10, rsi
mov r11, rdx
.L2:
%assign unroll 32
%assign i 0
%rep unroll
vmovaps ymm1, [r9 + 32*i]
vfmadd231ps ymm1, ymm2, [r10 + 32*i]
vmovaps [r11 + 32*i], ymm1
%assign i i+1
%endrep
add r9, 32*unroll
add r10, 32*unroll
add r11, 32*unroll
cmp r9, rcx
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
triad_ava_asm.asm：
global triad_avx_asm_repeat
;RDI x, RSI y, RDX z, RCX n, R8 repeat
pi: dd 3.14159
align 16
section .text
triad_avx_asm_repeat:
shl rcx, 2
add rdi, rcx
add rsi, rcx
add rdx, rcx
vbroadcastss ymm2, [rel pi]
;neg rcx
align 16
.L1:
mov rax, rcx
neg rax
align 16
.L2:
vmulps ymm1, ymm2, [rdi+rax]
vaddps ymm1, ymm1, [rsi+rax]
vmovaps [rdx+rax], ymm1
add rax, 32
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
global triad_avx_asm_repeat2
;RDI x, RSI y, RDX z, RCX n, R8 repeat
;pi: dd 3.14159
align 16
section .text
triad_avx_asm_repeat2:
shl rcx, 2
vbroadcastss ymm2, [rel pi]
align 16
.L1:
xor rax, rax
align 16
.L2:
vmulps ymm1, ymm2, [rdi+rax]
vaddps ymm1, ymm1, [rsi+rax]
vmovaps [rdx+rax], ymm1
add eax, 32
cmp eax, ecx
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
global triad_avx_asm_repeat_unroll16
align 16
section .text
triad_avx_asm_repeat_unroll16:
shl rcx, 2
add rcx, rdi
vbroadcastss ymm2, [rel pi]
align 16
.L1:
xor rax, rax
mov r9, rdi
mov r10, rsi
mov r11, rdx
align 16
.L2:
%assign unroll 16
%assign i 0
%rep unroll
vmulps ymm1, ymm2, [r9 + 32*i]
vaddps ymm1, ymm1, [r10 + 32*i]
vmovaps [r11 + 32*i], ymm1
%assign i i+1
%endrep
add r9, 32*unroll
add r10, 32*unroll
add r11, 32*unroll
cmp r9, rcx
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
triad_sse_asm.asm：
global triad_sse_asm_repeat
;RDI x, RSI y, RDX z, RCX n, R8 repeat
pi: dd 3.14159
;align 16
section .text
triad_sse_asm_repeat:
shl rcx, 2
add rdi, rcx
add rsi, rcx
add rdx, rcx
movss xmm2, [rel pi]
shufps xmm2, xmm2, 0
;neg rcx
align 16
.L1:
mov rax, rcx
neg rax
align 16
.L2:
movaps xmm1, [rdi+rax]
mulps xmm1, xmm2
addps xmm1, [rsi+rax]
movaps [rdx+rax], xmm1
add rax, 16
jne .L2
sub r8d, 1
jnz .L1
ret
global triad_sse_asm_repeat2
;RDI x, RSI y, RDX z, RCX n, R8 repeat
;pi: dd 3.14159
;align 16
section .text
triad_sse_asm_repeat2:
shl rcx, 2
movss xmm2, [rel pi]
shufps xmm2, xmm2, 0
align 16
.L1:
xor rax, rax
align 16
.L2:
movaps xmm1, [rdi+rax]
mulps xmm1, xmm2
addps xmm1, [rsi+rax]
movaps [rdx+rax], xmm1
add eax, 16
cmp eax, ecx
jne .L2
sub r8d, 1
jnz .L1
ret
global triad_sse_asm_repeat_unroll16
section .text
triad_sse_asm_repeat_unroll16:
shl rcx, 2
add rcx, rdi
movss xmm2, [rel pi]
shufps xmm2, xmm2, 0
.L1:
xor rax, rax
mov r9, rdi
mov r10, rsi
mov r11, rdx
.L2:
%assign unroll 8
%assign i 0
%rep unroll
movaps xmm1, [r9 + 16*i]
mulps xmm1, xmm2,
addps xmm1, [r10 + 16*i]
movaps [r11 + 16*i], xmm1
%assign i i+1
%endrep
add r9, 16*unroll
add r10, 16*unroll
add r11, 16*unroll
cmp r9, rcx
jne .L2
sub r8d, 1
jnz .L1
ret