优化简单模板操作，将变量保存到寄存器中
Optimise Simple Stencil Operation keeping variables into registers
我试图使下面的代码更快，将两个变量（我们需要重用的变量）保留在寄存器中或任何比缓存更近的地方。该代码在位置 idx 处获取数组中的三个相邻元素并将它们相加。
void stencil(double * input, double * output){
unsigned int idx = 1;
output[0] = input[0] + input[1];
for(; idx < SIZE - 1; idx++){
output[idx] = input[idx-1] + input[idx] + input[idx+1];
}
output[idx] = input[idx-1] + input[idx];
}
我的实现是这样的：
void stencil(double * input, double * output){
unsigned int idx = 0;
double x , y = 0, z;
z = input[idx];
for(; idx < SIZE - 1; idx++){
x = y;
y = z;
z = input[idx + 1];
output[idx] = x + y + z;
}
output[idx] = y + z;
}
思路是复用前面操作的变量，让程序更快。
但是，该程序在速度和性能方面似乎并没有提高。我在 AMD Opteron(tm) Processor 6320 CPU 上使用 gcc，我正在使用以下标志编译代码：-march=native -O3 -Wall -std=c99.
我试过使用和不使用本机，生成的程序集不同，但我无法获得更好的性能。生成的没有 -march=native 标志的程序集如下所示：
stencil:
.LFB7:
.cfi_startproc
subl , %edx
movsd (%rdi), %xmm1
je .L4
movq %rsi, %rcx
xorpd %xmm0, %xmm0
xorl %eax, %eax
jmp .L3
.p2align 4,,10
.p2align 3
.L6:
movapd %xmm1, %xmm0
movapd %xmm2, %xmm1
.L3:
addl , %eax
addsd %xmm1, %xmm0
addq , %rcx
movl %eax, %r8d
movsd (%rdi,%r8,8), %xmm2
leaq 0(,%r8,8), %r9
addsd %xmm2, %xmm0
movsd %xmm0, -8(%rcx)
cmpl %edx, %eax
jne .L6
.L2:
addsd %xmm2, %xmm1
movsd %xmm1, (%rsi,%r9)
ret
.L4:
movapd %xmm1, %xmm2
xorl %r9d, %r9d
xorpd %xmm1, %xmm1
jmp .L2
加上 -march=native 标志看起来像这样:
stencil:
.LFB20:
.cfi_startproc
vmovsd (%rdi), %xmm1
vxorpd %xmm0, %xmm0, %xmm0
leaq 144(%rdi), %rdx
leaq 136(%rsi), %rax
xorl %ecx, %ecx
.p2align 4,,10
.p2align 3
.L2:
vaddsd %xmm1, %xmm0, %xmm0
vmovsd -136(%rdx), %xmm4
prefetcht0 (%rdx)
addl , %ecx
prefetchw (%rax)
addq , %rdx
addq , %rax
vaddsd %xmm1, %xmm4, %xmm1
vaddsd %xmm4, %xmm0, %xmm0
vmovsd %xmm0, -200(%rax)
vmovsd -192(%rdx), %xmm3
vaddsd %xmm3, %xmm1, %xmm1
vaddsd %xmm3, %xmm4, %xmm4
vmovsd %xmm1, -192(%rax)
vmovsd -184(%rdx), %xmm2
vaddsd %xmm2, %xmm4, %xmm4
vaddsd %xmm2, %xmm3, %xmm3
vmovsd %xmm4, -184(%rax)
vmovsd %xmm4, -184(%rax)
vmovsd -176(%rdx), %xmm0
vaddsd %xmm0, %xmm3, %xmm3
vaddsd %xmm0, %xmm2, %xmm2
vmovsd %xmm3, -176(%rax)
vmovsd -168(%rdx), %xmm1
vaddsd %xmm1, %xmm2, %xmm2
vaddsd %xmm1, %xmm0, %xmm0
vmovsd %xmm2, -168(%rax)
vmovsd -160(%rdx), %xmm2
vaddsd %xmm2, %xmm0, %xmm0
vaddsd %xmm2, %xmm1, %xmm1
vmovsd %xmm0, -160(%rax)
vmovsd -152(%rdx), %xmm0
vaddsd %xmm0, %xmm1, %xmm1
vaddsd %xmm0, %xmm2, %xmm2
vmovsd %xmm1, -152(%rax)
vmovsd -144(%rdx), %xmm1
vaddsd %xmm1, %xmm2, %xmm2
vmovsd %xmm2, -144(%rax)
cmpl 99999992, %ecx
jne .L2
movabsq 199999944, %rdx
movabsq 199999936, %rcx
addq %rdi, %rdx
addq %rsi, %rcx
xorl %eax, %eax
jmp .L3
.p2align 4,,7
.p2align 3
.L4:
vmovaps %xmm2, %xmm1
.L3:
vaddsd %xmm0, %xmm1, %xmm0
vmovsd (%rdx,%rax), %xmm2
vaddsd %xmm2, %xmm0, %xmm0
vmovsd %xmm0, (%rcx,%rax)
addq , %rax
vmovaps %xmm1, %xmm0
cmpq , %rax
jne .L4
vaddsd %xmm2, %xmm1, %xmm1
movabsq 199999992, %rax
vmovsd %xmm1, (%rsi,%rax)
ret
有没有人对如何让 GCC 将变量保存到寄存器中以提高代码速度有什么建议？或者任何其他方式让我的代码有效绕过缓存？
使用寄存器循环时，展开循环通常是个好主意。除非明确要求，否则 gcc 不会这样做。
这是一个 4 级循环展开的示例。
void stencil(double * input, double * output){
double x, y, z, w, u, v ;
x=0.0;
y=input[0];
int idx=0;
for(; idx < SIZE - 5; idx+=4){
z=input[idx+1];
w=input[idx+2];
u=input[idx+3];
v=input[idx+4];
output[idx] =x+y+z;
output[idx+1]=y+z+w;
output[idx+2]=z+w+u;
output[idx+3]=w+u+v;
x=u;
y=v;
}
z=input[idx+1];
w=input[idx+2];
u=input[idx+3];
output[idx] =x+y+z;
output[idx+1]=y+z+w;
output[idx+2]=z+w+u;
output[idx+3]=w+u;
}
idx值读写1个内存，每2个idx值1个寄存器拷贝。
可以尝试不同的展开级别，但每次迭代总是有 2 个寄存器副本，而 4 个似乎是一个很好的折衷方案。
如果 size 不是 4 的倍数，则需要序言。
void stencil(double * input, double * output){
double x, y, z, w, u, v ;
int idx=0;
int remain=SIZE%4;
x=0.0;y=input[0]
switch (remain) {
case 3: z=input[++idx]; output[idx-1]=x+y+z; x=y; y=z;
case 2: z=input[++idx]; output[idx-1]=x+y+z; x=y; y=z;
case 1: z=input[++idx]; output[idx-1]=x+y+z; x=y; y=z;
}
for(; idx < SIZE - 5; idx+=4){
z=input[idx+1];
....
不出所料，asm相当复杂，很难说会有什么收获。
您也可以尝试在您的原始代码上使用 -funroll-loops。编译器非常好，可能会提供更好的解决方案。
这是个好主意，但如果编译器知道它是安全的，它们已经为你做了。 使用 double *restrict output 和 const double *restrict input 来保证存储到 output[] 的编译器不会更改将从 input[].
读取的内容
但是使用 SIMD 的自动矢量化是一个更重要的优化，每条指令产生 2 或 4 个 double 结果。在检查重叠后，GCC 和 ICC 将在 -O3 执行此操作。 （但是 clang 无法自动矢量化它，只是用标量展开 [v]addsd 避免不必要的重新加载。
不幸的是，您的优化版本无法自动矢量化！（这是编译器的错误，即一个错过优化的错误，当它知道输出不重叠时，所以重新阅读来自内存的来源与否是等效的）。
看起来 gcc 在原始版本上做得很好，-O3 -march=native（特别是在针对 Intel 进行调整时，AVX 的更宽向量是值得的。）我计算了 4 double结果来自 3 个未对齐的负载和 2 个 vaddpd ymm.
它在使用矢量化循环之前检查重叠。您可以使用 double *restrict output 和 input 来保证指针不重叠，因此不需要回退循环。
L1d 缓存带宽在现代 CPU 上非常出色；重新加载相同的数据不是大问题（每个时钟加载 2 次）。指令吞吐量是一个更大的问题。内存源 addsd 并不比将数据保存在寄存器中花费更多。
如果使用 128 位向量进行向量化，保留 in[idx+1..2] 向量以用作下一次迭代的 in[idx+ -1..1] 向量是有意义的。 GCC 实际上确实这样做了。
但是当您每条指令生成 4 个结果时，一次迭代的 3 个输入向量中的 none 对下一次迭代直接有用。不过，通过随机播放来从加载结果中创建 3 个向量之一来节省一些加载端口带宽可能会很有用。如果我使用 __m256d 内在函数手动矢量化，我会尝试这样做。或者使用 float 和 128 位 __m128 向量。
#define SIZE 1000000
void stencil_restrict(double *restrict input, double *restrict output)
{
int idx = 1;
output[0] = input[0] + input[1];
for(; idx < SIZE - 1; idx++){
output[idx] = input[idx-1] + input[idx] + input[idx+1];
}
output[idx] = input[idx-1] + input[idx];
}
使用 gcc8.3 -O3 -Wall -std=c99 -march=broadwell -masm=intel、from the Godbolt compiler explorer 编译成此 asm（在这种情况下不需要 -ffast-math，并且不会生成与内循环的区别。）
stencil_restrict:
vmovsd xmm0, QWORD PTR [rdi]
vaddsd xmm0, xmm0, QWORD PTR [rdi+8]
xor eax, eax
vmovsd QWORD PTR [rsi], xmm0 # first iteration
### Main loop
.L12:
vmovupd ymm2, YMMWORD PTR [rdi+8+rax] # idx +0 .. +3
vaddpd ymm0, ymm2, YMMWORD PTR [rdi+rax] # idx -1 .. +2
vaddpd ymm0, ymm0, YMMWORD PTR [rdi+16+rax] # idx +1 .. +4
vmovupd YMMWORD PTR [rsi+8+rax], ymm0 # store idx +0 .. +3
add rax, 32 # byte offset += 32
cmp rax, 7999968
jne .L12
# cleanup of last few elements
vmovsd xmm1, QWORD PTR [rdi+7999976]
vaddsd xmm0, xmm1, QWORD PTR [rdi+7999968]
vaddsd xmm1, xmm1, QWORD PTR [rdi+7999984]
vunpcklpd xmm0, xmm0, xmm1
vaddpd xmm0, xmm0, XMMWORD PTR [rdi+7999984]
vmovups XMMWORD PTR [rsi+7999976], xmm0
vmovsd xmm0, QWORD PTR [rdi+7999984]
vaddsd xmm0, xmm0, QWORD PTR [rdi+7999992]
vmovsd QWORD PTR [rsi+7999992], xmm0
vzeroupper
ret
不幸的是，gcc 使用索引寻址模式，因此带有内存源的 vaddpd 指令在 SnB 系列（包括您的 Broadwell Xeon E5-2698 v4）上的前端未分层为 2 微指令。 Micro fusion and addressing modes
vmovupd ymm2, YMMWORD PTR [rdi+8+rax] # 1 uop, no micro-fusion
vaddpd ymm0, ymm2, YMMWORD PTR [rdi+rax] # 2 uops. (micro-fused in decoders/uop cache, unlaminates)
vaddpd ymm0, ymm0, YMMWORD PTR [rdi+16+rax] # 2 uops. (ditto)
vmovupd YMMWORD PTR [rsi+8+rax], ymm0 # 1 uop (stays micro-fused, but can't use the port 7 store AGU)
add rax, 32 # 1 uop
cmp rax, 7999968 # 0 uops, macro-fuses with JNE
jne .L12 # 1 uop
吞吐量分析，见https://agner.org/optimize/ and What considerations go into predicting latency for operations on modern superscalar processors and how can I calculate them by hand?
GCC 的循环是前端 issue/rename 阶段的 8 个融合域 uops 发送到无序后端。这意味着前端最大吞吐量是每 2 个周期 1 次迭代。
[v]addpd 在 Intel 之前，Skylake 只能在端口 1 上 运行，而 [v]mulpd 或 FMA 的吞吐量是其两倍。 （Skylake 放弃了专用的 FP 加法单元，运行s FP 加法与 mul 和 fma 相同。）所以这也是每个迭代瓶颈 2 个周期。
我们有 3 个加载 + 1 个存储，所有这些都需要端口 2 或 3 之一。（索引寻址模式存储不能在端口 7 上使用专用存储 AGU）。所以这是每个迭代瓶颈的另一个 2 个周期。但不是真的； 未对齐 跨高速缓存行边界的加载更昂贵。实验表明，Intel Skylake（可能还有 Broadwell）会重放被发现是缓存行拆分的加载微指令，因此它们会再次 运行 从第二个缓存行获取数据。 .
我们的数据是 8 字节对齐的，但我们将 32 字节负载均匀分布在 64 字节行内的所有 8 字节偏移量上。在这 8 个起始元素中的 5 个处，没有缓存行拆分。在另外 3 个，有。所以平均成本实际上是 3 * (8+3)/8 = 4.125 每次迭代分派的负载微指令。我不知道存储地址 uops 是否需要重播。可能不会;只是当数据从存储缓冲区提交到 L1d 时才重要，而不是存储地址或存储数据微指令。 （只要它不跨越 4k 边界分割，将 发生在输出未对齐的情况下。
假设除 output[1] 之外的任何输出对齐方式都是 32 字节对齐的。 asm 在循环外存储 output[0]，然后有效地执行 output[i*4 + 1]，因此每隔一个存储将是缓存行拆分。
在这种情况下，最好达到输出数组的对齐边界。 gcc7 和更早版本喜欢将其中一个指针与循环序言对齐，但不幸的是，他们选择了我们从所有对齐方式加载的输入。
无论如何，GCC 的实际瓶颈是端口 2 / 端口 3 吞吐量。 这 2 个端口每次迭代平均 5.125 uops = 理论最大平均吞吐量 1每 2.5625 次循环迭代（4 次加倍）.
使用非索引存储可以减少这个瓶颈。
但这忽略了 4k 拆分惩罚，这在 Broadwell 上约为 100 个周期，并假设完美的 HW 预取可以跟上 ~12.5 字节/周期的每一种方式（加载和存储）。 所以这更有可能会限制内存带宽，除非数据在 L2 缓存中已经很热。L1d 可以吸收相同字节的冗余负载，但仍然有大量的非冗余带宽。
一点展开会让乱序执行更进一步，并在硬件预取跟不上时帮助吸收缓存未命中的气泡。如果它为存储使用非索引寻址模式，它可以使用端口 7，从而减少端口 2/3 上的压力。这将使负载 运行 领先于添加，希望在穿过
时吸收气泡
具有 128 位向量的寄存器中的数据重用
来自gcc8.3 -O3 -Wall -std=c99 -march=broadwell -mno-avx
的内部循环
# prologue to reach an alignment boundary somewhere?
.L12:
movupd xmm2, XMMWORD PTR [rdi+rax]
movupd xmm1, XMMWORD PTR [rdi+8+rax]
addpd xmm0, xmm2
addpd xmm0, xmm1
movups XMMWORD PTR [rsi+rax], xmm0
add rax, 16
movapd xmm0, xmm1 # x = z
cmp rax, 7999992
jne .L12
这是对 gcc7.4 的回归，它避免了寄存器复制。 （但是 gcc7 将循环开销浪费在与数组索引分开的计数器上。）
# prologue to reach an alignment boundary so one load can be aligned.
# r10=input and r9=input+8 or something like that
# r8=output
.L18: # do {
movupd xmm0, XMMWORD PTR [r10+rdx]
add ecx, 1
addpd xmm0, xmm1 # x+y
movapd xmm1, XMMWORD PTR [r9+rdx] # z for this iteration, x for next
addpd xmm0, xmm1 # (x+y) + z
movups XMMWORD PTR [r8+rdx], xmm0
add rdx, 16
cmp ecx, r11d
jb .L18 # } while(i < max);
平均而言，这仍然可能比 AVX 256 位向量慢。
使用用于 128 位向量的 AVX（例如，为 Piledriver 调优），它可以避免单独的 movupd xmm0 加载，并使用 vaddpd xmm0, xmm1, [r10+rdx].
它们都无法使用对齐存储，而且在 input 中找到已知对齐后也无法利用将负载折叠到内存操作数中的优势 addpd :/
Skylake 上的实际性能实验表明，如果数据适合 L1d 缓存，实际性能与我的预测相当接近。
有趣的事实：对于像全局 double in[SIZE+10]; 这样的静态缓冲区，gcc 将创建一个使用非索引寻址模式的循环版本。对于 运行 在一个循环中多次使用 SIZE=1000，这将加速从 ~800ms 到 ~700ms。稍后会更新更多细节。
我试图使下面的代码更快，将两个变量（我们需要重用的变量）保留在寄存器中或任何比缓存更近的地方。该代码在位置 idx 处获取数组中的三个相邻元素并将它们相加。
void stencil(double * input, double * output){
unsigned int idx = 1;
output[0] = input[0] + input[1];
for(; idx < SIZE - 1; idx++){
output[idx] = input[idx-1] + input[idx] + input[idx+1];
}
output[idx] = input[idx-1] + input[idx];
}
我的实现是这样的：
void stencil(double * input, double * output){
unsigned int idx = 0;
double x , y = 0, z;
z = input[idx];
for(; idx < SIZE - 1; idx++){
x = y;
y = z;
z = input[idx + 1];
output[idx] = x + y + z;
}
output[idx] = y + z;
}
思路是复用前面操作的变量，让程序更快。
但是，该程序在速度和性能方面似乎并没有提高。我在 AMD Opteron(tm) Processor 6320 CPU 上使用 gcc，我正在使用以下标志编译代码：-march=native -O3 -Wall -std=c99.
我试过使用和不使用本机，生成的程序集不同，但我无法获得更好的性能。生成的没有 -march=native 标志的程序集如下所示：
stencil:
.LFB7:
.cfi_startproc
subl , %edx
movsd (%rdi), %xmm1
je .L4
movq %rsi, %rcx
xorpd %xmm0, %xmm0
xorl %eax, %eax
jmp .L3
.p2align 4,,10
.p2align 3
.L6:
movapd %xmm1, %xmm0
movapd %xmm2, %xmm1
.L3:
addl , %eax
addsd %xmm1, %xmm0
addq , %rcx
movl %eax, %r8d
movsd (%rdi,%r8,8), %xmm2
leaq 0(,%r8,8), %r9
addsd %xmm2, %xmm0
movsd %xmm0, -8(%rcx)
cmpl %edx, %eax
jne .L6
.L2:
addsd %xmm2, %xmm1
movsd %xmm1, (%rsi,%r9)
ret
.L4:
movapd %xmm1, %xmm2
xorl %r9d, %r9d
xorpd %xmm1, %xmm1
jmp .L2
加上 -march=native 标志看起来像这样:
stencil:
.LFB20:
.cfi_startproc
vmovsd (%rdi), %xmm1
vxorpd %xmm0, %xmm0, %xmm0
leaq 144(%rdi), %rdx
leaq 136(%rsi), %rax
xorl %ecx, %ecx
.p2align 4,,10
.p2align 3
.L2:
vaddsd %xmm1, %xmm0, %xmm0
vmovsd -136(%rdx), %xmm4
prefetcht0 (%rdx)
addl , %ecx
prefetchw (%rax)
addq , %rdx
addq , %rax
vaddsd %xmm1, %xmm4, %xmm1
vaddsd %xmm4, %xmm0, %xmm0
vmovsd %xmm0, -200(%rax)
vmovsd -192(%rdx), %xmm3
vaddsd %xmm3, %xmm1, %xmm1
vaddsd %xmm3, %xmm4, %xmm4
vmovsd %xmm1, -192(%rax)
vmovsd -184(%rdx), %xmm2
vaddsd %xmm2, %xmm4, %xmm4
vaddsd %xmm2, %xmm3, %xmm3
vmovsd %xmm4, -184(%rax)
vmovsd %xmm4, -184(%rax)
vmovsd -176(%rdx), %xmm0
vaddsd %xmm0, %xmm3, %xmm3
vaddsd %xmm0, %xmm2, %xmm2
vmovsd %xmm3, -176(%rax)
vmovsd -168(%rdx), %xmm1
vaddsd %xmm1, %xmm2, %xmm2
vaddsd %xmm1, %xmm0, %xmm0
vmovsd %xmm2, -168(%rax)
vmovsd -160(%rdx), %xmm2
vaddsd %xmm2, %xmm0, %xmm0
vaddsd %xmm2, %xmm1, %xmm1
vmovsd %xmm0, -160(%rax)
vmovsd -152(%rdx), %xmm0
vaddsd %xmm0, %xmm1, %xmm1
vaddsd %xmm0, %xmm2, %xmm2
vmovsd %xmm1, -152(%rax)
vmovsd -144(%rdx), %xmm1
vaddsd %xmm1, %xmm2, %xmm2
vmovsd %xmm2, -144(%rax)
cmpl 99999992, %ecx
jne .L2
movabsq 199999944, %rdx
movabsq 199999936, %rcx
addq %rdi, %rdx
addq %rsi, %rcx
xorl %eax, %eax
jmp .L3
.p2align 4,,7
.p2align 3
.L4:
vmovaps %xmm2, %xmm1
.L3:
vaddsd %xmm0, %xmm1, %xmm0
vmovsd (%rdx,%rax), %xmm2
vaddsd %xmm2, %xmm0, %xmm0
vmovsd %xmm0, (%rcx,%rax)
addq , %rax
vmovaps %xmm1, %xmm0
cmpq , %rax
jne .L4
vaddsd %xmm2, %xmm1, %xmm1
movabsq 199999992, %rax
vmovsd %xmm1, (%rsi,%rax)
ret
有没有人对如何让 GCC 将变量保存到寄存器中以提高代码速度有什么建议？或者任何其他方式让我的代码有效绕过缓存？
使用寄存器循环时，展开循环通常是个好主意。除非明确要求，否则 gcc 不会这样做。
这是一个 4 级循环展开的示例。
void stencil(double * input, double * output){
double x, y, z, w, u, v ;
x=0.0;
y=input[0];
int idx=0;
for(; idx < SIZE - 5; idx+=4){
z=input[idx+1];
w=input[idx+2];
u=input[idx+3];
v=input[idx+4];
output[idx] =x+y+z;
output[idx+1]=y+z+w;
output[idx+2]=z+w+u;
output[idx+3]=w+u+v;
x=u;
y=v;
}
z=input[idx+1];
w=input[idx+2];
u=input[idx+3];
output[idx] =x+y+z;
output[idx+1]=y+z+w;
output[idx+2]=z+w+u;
output[idx+3]=w+u;
}
idx值读写1个内存，每2个idx值1个寄存器拷贝。
可以尝试不同的展开级别，但每次迭代总是有 2 个寄存器副本，而 4 个似乎是一个很好的折衷方案。
如果 size 不是 4 的倍数，则需要序言。
void stencil(double * input, double * output){
double x, y, z, w, u, v ;
int idx=0;
int remain=SIZE%4;
x=0.0;y=input[0]
switch (remain) {
case 3: z=input[++idx]; output[idx-1]=x+y+z; x=y; y=z;
case 2: z=input[++idx]; output[idx-1]=x+y+z; x=y; y=z;
case 1: z=input[++idx]; output[idx-1]=x+y+z; x=y; y=z;
}
for(; idx < SIZE - 5; idx+=4){
z=input[idx+1];
....
不出所料，asm相当复杂，很难说会有什么收获。
您也可以尝试在您的原始代码上使用 -funroll-loops。编译器非常好，可能会提供更好的解决方案。
这是个好主意，但如果编译器知道它是安全的，它们已经为你做了。 使用 double *restrict output 和 const double *restrict input 来保证存储到 output[] 的编译器不会更改将从 input[].
但是使用 SIMD 的自动矢量化是一个更重要的优化，每条指令产生 2 或 4 个 double 结果。在检查重叠后，GCC 和 ICC 将在 -O3 执行此操作。 （但是 clang 无法自动矢量化它，只是用标量展开 [v]addsd 避免不必要的重新加载。
不幸的是，您的优化版本无法自动矢量化！（这是编译器的错误，即一个错过优化的错误，当它知道输出不重叠时，所以重新阅读来自内存的来源与否是等效的）。
看起来 gcc 在原始版本上做得很好，-O3 -march=native（特别是在针对 Intel 进行调整时，AVX 的更宽向量是值得的。）我计算了 4 double结果来自 3 个未对齐的负载和 2 个 vaddpd ymm.
它在使用矢量化循环之前检查重叠。您可以使用 double *restrict output 和 input 来保证指针不重叠，因此不需要回退循环。
L1d 缓存带宽在现代 CPU 上非常出色；重新加载相同的数据不是大问题（每个时钟加载 2 次）。指令吞吐量是一个更大的问题。内存源 addsd 并不比将数据保存在寄存器中花费更多。
如果使用 128 位向量进行向量化，保留 in[idx+1..2] 向量以用作下一次迭代的 in[idx+ -1..1] 向量是有意义的。 GCC 实际上确实这样做了。
但是当您每条指令生成 4 个结果时，一次迭代的 3 个输入向量中的 none 对下一次迭代直接有用。不过，通过随机播放来从加载结果中创建 3 个向量之一来节省一些加载端口带宽可能会很有用。如果我使用 __m256d 内在函数手动矢量化，我会尝试这样做。或者使用 float 和 128 位 __m128 向量。
#define SIZE 1000000
void stencil_restrict(double *restrict input, double *restrict output)
{
int idx = 1;
output[0] = input[0] + input[1];
for(; idx < SIZE - 1; idx++){
output[idx] = input[idx-1] + input[idx] + input[idx+1];
}
output[idx] = input[idx-1] + input[idx];
}
使用 gcc8.3 -O3 -Wall -std=c99 -march=broadwell -masm=intel、from the Godbolt compiler explorer 编译成此 asm（在这种情况下不需要 -ffast-math，并且不会生成与内循环的区别。）
stencil_restrict:
vmovsd xmm0, QWORD PTR [rdi]
vaddsd xmm0, xmm0, QWORD PTR [rdi+8]
xor eax, eax
vmovsd QWORD PTR [rsi], xmm0 # first iteration
### Main loop
.L12:
vmovupd ymm2, YMMWORD PTR [rdi+8+rax] # idx +0 .. +3
vaddpd ymm0, ymm2, YMMWORD PTR [rdi+rax] # idx -1 .. +2
vaddpd ymm0, ymm0, YMMWORD PTR [rdi+16+rax] # idx +1 .. +4
vmovupd YMMWORD PTR [rsi+8+rax], ymm0 # store idx +0 .. +3
add rax, 32 # byte offset += 32
cmp rax, 7999968
jne .L12
# cleanup of last few elements
vmovsd xmm1, QWORD PTR [rdi+7999976]
vaddsd xmm0, xmm1, QWORD PTR [rdi+7999968]
vaddsd xmm1, xmm1, QWORD PTR [rdi+7999984]
vunpcklpd xmm0, xmm0, xmm1
vaddpd xmm0, xmm0, XMMWORD PTR [rdi+7999984]
vmovups XMMWORD PTR [rsi+7999976], xmm0
vmovsd xmm0, QWORD PTR [rdi+7999984]
vaddsd xmm0, xmm0, QWORD PTR [rdi+7999992]
vmovsd QWORD PTR [rsi+7999992], xmm0
vzeroupper
ret
不幸的是，gcc 使用索引寻址模式，因此带有内存源的 vaddpd 指令在 SnB 系列（包括您的 Broadwell Xeon E5-2698 v4）上的前端未分层为 2 微指令。 Micro fusion and addressing modes
vmovupd ymm2, YMMWORD PTR [rdi+8+rax] # 1 uop, no micro-fusion
vaddpd ymm0, ymm2, YMMWORD PTR [rdi+rax] # 2 uops. (micro-fused in decoders/uop cache, unlaminates)
vaddpd ymm0, ymm0, YMMWORD PTR [rdi+16+rax] # 2 uops. (ditto)
vmovupd YMMWORD PTR [rsi+8+rax], ymm0 # 1 uop (stays micro-fused, but can't use the port 7 store AGU)
add rax, 32 # 1 uop
cmp rax, 7999968 # 0 uops, macro-fuses with JNE
jne .L12 # 1 uop
吞吐量分析，见https://agner.org/optimize/ and What considerations go into predicting latency for operations on modern superscalar processors and how can I calculate them by hand?
GCC 的循环是前端 issue/rename 阶段的 8 个融合域 uops 发送到无序后端。这意味着前端最大吞吐量是每 2 个周期 1 次迭代。
[v]addpd 在 Intel 之前，Skylake 只能在端口 1 上 运行，而 [v]mulpd 或 FMA 的吞吐量是其两倍。 （Skylake 放弃了专用的 FP 加法单元，运行s FP 加法与 mul 和 fma 相同。）所以这也是每个迭代瓶颈 2 个周期。
我们有 3 个加载 + 1 个存储，所有这些都需要端口 2 或 3 之一。（索引寻址模式存储不能在端口 7 上使用专用存储 AGU）。所以这是每个迭代瓶颈的另一个 2 个周期。但不是真的； 未对齐 跨高速缓存行边界的加载更昂贵。实验表明，Intel Skylake（可能还有 Broadwell）会重放被发现是缓存行拆分的加载微指令，因此它们会再次 运行 从第二个缓存行获取数据。
我们的数据是 8 字节对齐的，但我们将 32 字节负载均匀分布在 64 字节行内的所有 8 字节偏移量上。在这 8 个起始元素中的 5 个处，没有缓存行拆分。在另外 3 个，有。所以平均成本实际上是 3 * (8+3)/8 = 4.125 每次迭代分派的负载微指令。我不知道存储地址 uops 是否需要重播。可能不会;只是当数据从存储缓冲区提交到 L1d 时才重要，而不是存储地址或存储数据微指令。 （只要它不跨越 4k 边界分割，将 发生在输出未对齐的情况下。
假设除 output[1] 之外的任何输出对齐方式都是 32 字节对齐的。 asm 在循环外存储 output[0]，然后有效地执行 output[i*4 + 1]，因此每隔一个存储将是缓存行拆分。
在这种情况下，最好达到输出数组的对齐边界。 gcc7 和更早版本喜欢将其中一个指针与循环序言对齐，但不幸的是，他们选择了我们从所有对齐方式加载的输入。
无论如何，GCC 的实际瓶颈是端口 2 / 端口 3 吞吐量。 这 2 个端口每次迭代平均 5.125 uops = 理论最大平均吞吐量 1每 2.5625 次循环迭代（4 次加倍）.
使用非索引存储可以减少这个瓶颈。
但这忽略了 4k 拆分惩罚，这在 Broadwell 上约为 100 个周期，并假设完美的 HW 预取可以跟上 ~12.5 字节/周期的每一种方式（加载和存储）。 所以这更有可能会限制内存带宽，除非数据在 L2 缓存中已经很热。L1d 可以吸收相同字节的冗余负载，但仍然有大量的非冗余带宽。
一点展开会让乱序执行更进一步，并在硬件预取跟不上时帮助吸收缓存未命中的气泡。如果它为存储使用非索引寻址模式，它可以使用端口 7，从而减少端口 2/3 上的压力。这将使负载 运行 领先于添加，希望在穿过
时吸收气泡具有 128 位向量的寄存器中的数据重用
来自gcc8.3 -O3 -Wall -std=c99 -march=broadwell -mno-avx
# prologue to reach an alignment boundary somewhere?
.L12:
movupd xmm2, XMMWORD PTR [rdi+rax]
movupd xmm1, XMMWORD PTR [rdi+8+rax]
addpd xmm0, xmm2
addpd xmm0, xmm1
movups XMMWORD PTR [rsi+rax], xmm0
add rax, 16
movapd xmm0, xmm1 # x = z
cmp rax, 7999992
jne .L12
这是对 gcc7.4 的回归，它避免了寄存器复制。 （但是 gcc7 将循环开销浪费在与数组索引分开的计数器上。）
# prologue to reach an alignment boundary so one load can be aligned.
# r10=input and r9=input+8 or something like that
# r8=output
.L18: # do {
movupd xmm0, XMMWORD PTR [r10+rdx]
add ecx, 1
addpd xmm0, xmm1 # x+y
movapd xmm1, XMMWORD PTR [r9+rdx] # z for this iteration, x for next
addpd xmm0, xmm1 # (x+y) + z
movups XMMWORD PTR [r8+rdx], xmm0
add rdx, 16
cmp ecx, r11d
jb .L18 # } while(i < max);
平均而言，这仍然可能比 AVX 256 位向量慢。
使用用于 128 位向量的 AVX（例如，为 Piledriver 调优），它可以避免单独的 movupd xmm0 加载，并使用 vaddpd xmm0, xmm1, [r10+rdx].
它们都无法使用对齐存储，而且在 input 中找到已知对齐后也无法利用将负载折叠到内存操作数中的优势 addpd :/
Skylake 上的实际性能实验表明，如果数据适合 L1d 缓存，实际性能与我的预测相当接近。
有趣的事实：对于像全局 double in[SIZE+10]; 这样的静态缓冲区，gcc 将创建一个使用非索引寻址模式的循环版本。对于 运行 在一个循环中多次使用 SIZE=1000，这将加速从 ~800ms 到 ~700ms。稍后会更新更多细节。