用 SSE 计算 4d 向量平均值
Calculate 4d vectors average with SSE
我尝试加快计算放置在数组中的 4d 向量的平均值。这是我的代码：
#include <sys/time.h>
#include <sys/param.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <xmmintrin.h>
typedef float dot[4];
#define N 1000000
double gettime ()
{
struct timeval tv;
gettimeofday (&tv, 0);
return (double)tv.tv_sec + (0.000001 * (double)tv.tv_usec);
}
void calc_avg1 (dot res, const dot array[], int n)
{
int i,j;
memset (res, 0, sizeof (dot));
for (i = 0; i < n; i++)
{
for (j = 0; j<4; j++) res[j] += array[i][j];
}
for (j = 0; j<4; j++) res[j] /= n;
}
void calc_avg2 (dot res, const dot array[], int n)
{
int i;
__v4sf r = _mm_set1_ps (0.0);
for (i=0; i<n; i++) r += _mm_load_ps (array[i]);
r /= _mm_set1_ps ((float)n);
_mm_store_ps (res, r);
}
int main ()
{
void *space = malloc (N*sizeof(dot)+15);
dot *array = (dot*)(((unsigned long)space+15) & ~(unsigned long)15);
dot avg __attribute__((aligned(16)));
int i;
double time;
for (i = 0; i < N; i++)
{
array[i][0] = 1.0*random();
array[i][1] = 1.0*random();
array[i][2] = 1.0*random();
}
time = gettime();
calc_avg1 (avg, array, N);
time = gettime() - time;
printf ("%f\n%f %f %f\n", time, avg[0], avg[1], avg[2]);
time = gettime();
calc_avg2 (avg, array, N);
time = gettime() - time;
printf ("%f\n%f %f %f\n", time, avg[0], avg[1], avg[2]);
return 0;
}
正如您所见，calc_avg1 使用从 0 到 4 的简单循环，calc_avg2 将它们替换为 SSE 指令。我用 clang 3.4 编译这段代码：
cc -O2 -o test test.c
这里是 calc_avgX 函数的反汇编：
0000000000400860 <calc_avg1>:
400860: 55 push %rbp
400861: 48 89 e5 mov %rsp,%rbp
400864: 85 d2 test %edx,%edx
400866: 0f 57 c0 xorps %xmm0,%xmm0
400869: 0f 11 07 movups %xmm0,(%rdi)
40086c: 7e 42 jle 4008b0 <calc_avg1+0x50>
40086e: 48 83 c6 0c add [=12=]xc,%rsi
400872: 0f 57 c0 xorps %xmm0,%xmm0
400875: 89 d0 mov %edx,%eax
400877: 0f 57 c9 xorps %xmm1,%xmm1
40087a: 0f 57 d2 xorps %xmm2,%xmm2
40087d: 0f 57 db xorps %xmm3,%xmm3
400880: f3 0f 58 5e f4 addss -0xc(%rsi),%xmm3
400885: f3 0f 11 1f movss %xmm3,(%rdi)
400889: f3 0f 58 56 f8 addss -0x8(%rsi),%xmm2
40088e: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
400893: f3 0f 58 4e fc addss -0x4(%rsi),%xmm1
400898: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
40089d: f3 0f 58 06 addss (%rsi),%xmm0
4008a1: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008a6: 48 83 c6 10 add [=12=]x10,%rsi
4008aa: ff c8 dec %eax
4008ac: 75 d2 jne 400880 <calc_avg1+0x20>
4008ae: eb 0c jmp 4008bc <calc_avg1+0x5c>
4008b0: 0f 57 c0 xorps %xmm0,%xmm0
4008b3: 0f 57 c9 xorps %xmm1,%xmm1
4008b6: 0f 57 d2 xorps %xmm2,%xmm2
4008b9: 0f 57 db xorps %xmm3,%xmm3
4008bc: f3 0f 2a e2 cvtsi2ss %edx,%xmm4
4008c0: f3 0f 5e dc divss %xmm4,%xmm3
4008c4: f3 0f 11 1f movss %xmm3,(%rdi)
4008c8: f3 0f 5e d4 divss %xmm4,%xmm2
4008cc: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
4008d1: f3 0f 5e cc divss %xmm4,%xmm1
4008d5: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
4008da: f3 0f 5e c4 divss %xmm4,%xmm0
4008de: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008e3: 5d pop %rbp
4008e4: c3 retq
4008e5: 66 66 2e 0f 1f 84 00 nopw %cs:0x0(%rax,%rax,1)
4008ec: 00 00 00 00
00000000004008f0 <calc_avg2>:
4008f0: 55 push %rbp
4008f1: 48 89 e5 mov %rsp,%rbp
4008f4: 85 d2 test %edx,%edx
4008f6: 0f 57 c0 xorps %xmm0,%xmm0
4008f9: 7e 10 jle 40090b <calc_avg2+0x1b>
4008fb: 89 d0 mov %edx,%eax
4008fd: 0f 1f 00 nopl (%rax)
400900: 0f 58 06 addps (%rsi),%xmm0
400903: 48 83 c6 10 add [=12=]x10,%rsi
400907: ff c8 dec %eax
400909: 75 f5 jne 400900 <calc_avg2+0x10>
40090b: 66 0f 6e ca movd %edx,%xmm1
40090f: 66 0f 70 c9 00 pshufd [=12=]x0,%xmm1,%xmm1
400914: 0f 5b c9 cvtdq2ps %xmm1,%xmm1
400917: 0f 5e c1 divps %xmm1,%xmm0
40091a: 0f 29 07 movaps %xmm0,(%rdi)
40091d: 5d pop %rbp
40091e: c3 retq
40091f: 90 nop
结果如下：
> ./test
0.004287
1073864320.000000 1074018048.000000 1073044224.000000
0.003661
1073864320.000000 1074018048.000000 1073044224.000000
所以 SSE 版本快了 1.17 倍。但是当我尝试做看似相同的工作时，即计算数组中单精度标量的平均值（例如此处 SSE reduction of float vector 所述），SSE 版本运行速度快 3.32 倍。这是 calc_avgX 函数的代码：
float calc_avg1 (const float array[], int n)
{
int i;
float avg = 0;
for (i = 0; i < n; i++) avg += array[i];
return avg / n;
}
float calc_avg3 (const float array[], int n)
{
int i;
__v4sf r = _mm_set1_ps (0.0);
for (i=0; i<n; i+=4) r += _mm_load_ps (&(array[i]));
r = _mm_hadd_ps (r, r);
r = _mm_hadd_ps (r, r);
return r[0] / n;
}
所以我的问题是：为什么我在最后一个示例（计算单个浮点标量的平均值）中从 SSE 中获益如此之多，而在第一个示例（计算 4d 向量的平均值）中却没有？对我来说，这些工作几乎是一样的。如果我做错了，在第一个例子中加速计算的正确方法是什么？
更新：
如果您认为相关，我还提供了第二个示例的反汇编，其中计算了标量的平均值（也使用 clang3.4 -O2 编译）。
0000000000400860 <calc_avg1>:
400860: 55 push %rbp
400861: 48 89 e5 mov %rsp,%rbp
400864: 85 d2 test %edx,%edx
400866: 0f 57 c0 xorps %xmm0,%xmm0
400869: 0f 11 07 movups %xmm0,(%rdi)
40086c: 7e 42 jle 4008b0 <calc_avg1+0x50>
40086e: 48 83 c6 0c add [=15=]xc,%rsi
400872: 0f 57 c0 xorps %xmm0,%xmm0
400875: 89 d0 mov %edx,%eax
400877: 0f 57 c9 xorps %xmm1,%xmm1
40087a: 0f 57 d2 xorps %xmm2,%xmm2
40087d: 0f 57 db xorps %xmm3,%xmm3
400880: f3 0f 58 5e f4 addss -0xc(%rsi),%xmm3
400885: f3 0f 11 1f movss %xmm3,(%rdi)
400889: f3 0f 58 56 f8 addss -0x8(%rsi),%xmm2
40088e: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
400893: f3 0f 58 4e fc addss -0x4(%rsi),%xmm1
400898: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
40089d: f3 0f 58 06 addss (%rsi),%xmm0
4008a1: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008a6: 48 83 c6 10 add [=15=]x10,%rsi
4008aa: ff c8 dec %eax
4008ac: 75 d2 jne 400880 <calc_avg1+0x20>
4008ae: eb 0c jmp 4008bc <calc_avg1+0x5c>
4008b0: 0f 57 c0 xorps %xmm0,%xmm0
4008b3: 0f 57 c9 xorps %xmm1,%xmm1
4008b6: 0f 57 d2 xorps %xmm2,%xmm2
4008b9: 0f 57 db xorps %xmm3,%xmm3
4008bc: f3 0f 2a e2 cvtsi2ss %edx,%xmm4
4008c0: f3 0f 5e dc divss %xmm4,%xmm3
4008c4: f3 0f 11 1f movss %xmm3,(%rdi)
4008c8: f3 0f 5e d4 divss %xmm4,%xmm2
4008cc: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
4008d1: f3 0f 5e cc divss %xmm4,%xmm1
4008d5: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
4008da: f3 0f 5e c4 divss %xmm4,%xmm0
4008de: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008e3: 5d pop %rbp
4008e4: c3 retq
4008e5: 66 66 2e 0f 1f 84 00 nopw %cs:0x0(%rax,%rax,1)
4008ec: 00 00 00 00
00000000004008d0 <calc_avg3>:
4008d0: 55 push %rbp
4008d1: 48 89 e5 mov %rsp,%rbp
4008d4: 31 c0 xor %eax,%eax
4008d6: 85 f6 test %esi,%esi
4008d8: 0f 57 c0 xorps %xmm0,%xmm0
4008db: 7e 0f jle 4008ec <calc_avg3+0x1c>
4008dd: 0f 1f 00 nopl (%rax)
4008e0: 0f 58 04 87 addps (%rdi,%rax,4),%xmm0
4008e4: 48 83 c0 04 add [=15=]x4,%rax
4008e8: 39 f0 cmp %esi,%eax
4008ea: 7c f4 jl 4008e0 <calc_avg3+0x10>
4008ec: 66 0f 70 c8 01 pshufd [=15=]x1,%xmm0,%xmm1
4008f1: f3 0f 58 c8 addss %xmm0,%xmm1
4008f5: 66 0f 70 d0 03 pshufd [=15=]x3,%xmm0,%xmm2
4008fa: 0f 12 c0 movhlps %xmm0,%xmm0
4008fd: f3 0f 58 c1 addss %xmm1,%xmm0
400901: f3 0f 58 c2 addss %xmm2,%xmm0
400905: 0f 57 c9 xorps %xmm1,%xmm1
400908: f3 0f 2a ce cvtsi2ss %esi,%xmm1
40090c: f3 0f 5e c1 divss %xmm1,%xmm0
400910: 5d pop %rbp
400911: c3 retq
400912: 66 66 66 66 66 2e 0f nopw %cs:0x0(%rax,%rax,1)
400919: 1f 84 00 00 00 00 00
抱歉，这个回答有点冗长和杂乱无章。我 运行 一些基准测试，但在考虑其他尝试后，我没有花很长时间编辑早期的东西。
您的工作集为 15.25MiB (16MB)。通常要对这样的例程进行基准测试，您会多次平均一个较小的缓冲区，因此它适合缓存。您看不到慢速版本和快速版本之间的太大差异，因为差异被内存瓶颈隐藏了。
calc_avg1 根本不会自动向量化（注意 addss。ss 表示标量，单精度，与 addps 相对（打包单精度-精确））。我认为即使内联到 main 中它也不能自动向量化，因为它不能确定第 4 个向量位置中没有 NaNs，这会导致标量代码没有的 FP 异常。我尝试使用 gcc 4.9.2 -O3 -march=native -ffast-math 和 clang-3.5 为 Sandybridge 编译它，但都没有成功。
即便如此，内联到main的版本也只是运行速度稍慢，因为内存是瓶颈。当访问主内存时，32 位加载几乎可以跟上 128b 加载。 （但是，非内联版本会很糟糕：每个 += 结果都存储到 res 数组中，因为循环直接累积到可能有其他引用的内存中。所以它必须商店可见的每个操作。顺便说一句，这是您发布反汇编的版本。通过 -S -fverbose-asm.)
编译协助找出 main 的哪一部分
有点令人失望的是，clang 和 gcc 无法将 __v4sf 从 4 宽 AVX 自动矢量化为 8 宽 AVX。
从头开始，在对 calc_avgX 的调用包装 for (int i=0; i<4000 ; i++) 并将 N 减少到 10k 之后，gcc -O3 将 avg1 的内部循环进入：
400690: c5 f8 10 08 vmovups (%rax),%xmm1
400694: 48 83 c0 20 add [=10=]x20,%rax
400698: c4 e3 75 18 48 f0 01 vinsertf128 [=10=]x1,-0x10(%rax),%ymm1,%ymm1
40069f: c5 fc 58 c1 vaddps %ymm1,%ymm0,%ymm0
4006a3: 48 39 d8 cmp %rbx,%rax
4006a6: 75 e8 jne 400690 <main+0xe0>
$ (get CPU to max-turbo frequency) && time ./a.out
0.016515
1071570752.000000 1066917696.000000 1073897344.000000
0.032875
1071570944.000000 1066916416.000000 1073895680.000000
这很奇怪；我不知道为什么它不只使用 32B 负载。它确实使用了 32B vaddps，这是处理适合 L2 缓存的数据集时的瓶颈。
IDK 为什么它在另一个循环内时设法自动矢量化内部循环。请注意，此 仅 适用于内联到 main 的版本。可调用版本仍然是纯标量的。还要注意只有 gcc 管理这个。 clang 3.5 没有。也许 gcc 知道它将以 return 清零缓冲区的方式使用 malloc（因此它不必担心第 4 个元素中的 NaNs）？
我也很惊讶 clang 的非矢量化 avg1 并不慢，当一切都适合缓存时。 N=10000，重复计数 = 40k。
3.3GHz SNB i5 2500k, max turbo = 3.8GHz.
avg1: 0.350422s: clang -O3 -march=native (not vectorized. loop of 6 scalar addss with memory operands)
avg2: 0.320173s: clang -O3 -march=native
avg1: 0.497040s: clang -O3 -march=native -ffast-math (haven't looked at asm to see what happened)
avg1: 0.160374s: gcc -O3 -march=native (256b addps, with 2 128b loads)
avg2: 0.321028s: gcc -O3 -march=native (128b addps with a memory operand)
avg2: ~0.16: clang, unrolled with 2 dependency chains to hide latency (see below).
avg2: ~0.08: unrolled with 4 dep chains
avg2: ~0.04: in theory unrolled-by-4 with 256b AVX. I didn't try unrolling the one gcc auto-vectorized with 256b addps
所以最大的惊喜是纯标量 clang avg1 代码跟上了 avg2。也许循环携带的依赖链是更大的瓶颈？
perf 显示 clang 的非矢量化 avg1 每个周期有 1.47 个 insns，这可能使端口 1 上的 FP 添加单元饱和。（大多数循环指令都是添加）。
然而，avg2，使用 128b addps 和内存 ope运行d，每个周期仅获得 0.58 insns。将数组大小再减少 10 倍，达到 N=1000，每个周期获得 0.60 insns，可能是因为在 prologue/epilogue 中花费了更多时间。所以我认为循环携带的依赖链存在一个严重的问题。 clang 将循环展开 4，但只使用一个累加器。该循环有 7 条指令，解码为 10 微指令。 （每个 vaddps 是 2，因为它与内存 ope运行d 一起使用，具有 2 寄存器寻址模式，防止微融合。 cmp 和 jne 宏-保险丝）。 http://www.brendangregg.com/perf.html 说 UOPS_DISPATCHED.CORE 的 perf 事件是 r2b1，所以：
$ perf stat -d -e cycles,instructions,r2b1 ./a.out
0.031793
1053298112.000000 1052673664.000000 1116960256.000000
Performance counter stats for './a.out':
118,453,541 cycles
71,181,299 instructions # 0.60 insns per cycle
102,025,443 r2b1 # this is uops, but perf doesn't have a nice name for it
40,256,019 L1-dcache-loads
21,254 L1-dcache-load-misses # 0.05% of all L1-dcache hits
9,588 LLC-loads
0 LLC-load-misses:HG # 0.00% of all LL-cache hits
0.032276233 seconds time elapsed
这证实了我对 uops 分析的 7:10 说明。这实际上与此处的性能问题无关：循环 运行 way 比每个循环上限 4 微指令慢。更改内部循环以获得两个独立的 dep 链使吞吐量翻倍（60M 周期而不是 117M，但 81M insns 而不是 71M）：
for (i=0; i<n-1; i+=2) { // TODO: make sure the loop end condition is correct
r0 += _mm_load_ps (array[i]);
r1 += _mm_load_ps (array[i+1]);
}
r0 += r1;
展开 4（在循环结束时合并 4 个累加器）再次使性能翻倍。 （低至 42M 周期，81M insns，112M uops。）内循环有 4x vaddps -0x30(%rcx),%xmm4,%xmm4（和类似）、2x add、cmp、jl。这种形式的 vaddps 应该微融合，但我仍然看到比指令多得多的微指令，所以我猜 r2b1 算作未融合域的微指令。 （Linux perf 没有针对特定平台硬件事件的良好文档）。 C运行 再次提高 N，以确保它是完全控制所有计数的最内层循环，我看到 uop:insn 比率为 1.39，这与 8 insns、11 uops 非常匹配(1.375)（将 vaddps 算作 2，但将 cmp + jl 算作 1）。我找到了 http://www.bnikolic.co.uk/blog/hpc-prof-events.html，其中包含支持的性能事件的完整列表，包括它们用于 Sandybridge 的代码。 （以及如何为任何其他 CPU 转储 table 的说明）。 （在每个块中查找 Code: 行。您需要一个 umask 字节，然后是代码，作为 perf 的参数。）
# a.out does only avg2, as an unrolled-by-4 version.
$ perf stat -d -e cycles,instructions,r14a1,r2b1,r10e,r2c2,r1c2 ./a.out
0.011331
1053298752.000000 1052674496.000000 1116959488.000000
Performance counter stats for './a.out':
42,250,312 cycles [34.11%]
56,103,429 instructions # 1.33 insns per cycle
20,864,416 r14a1 # UOPS_DISPATCHED_PORT: 0x14=port2&3 loads
111,943,380 r2b1 # UOPS_DISPATCHED: (2->umask 00 -> this core, any thread).
72,208,772 r10e # UOPS_ISSUED: fused-domain
71,422,907 r2c2 # UOPS_RETIRED: retirement slots used (fused-domain)
111,597,049 r1c2 # UOPS_RETIRED: ALL (unfused-domain)
0 L1-dcache-loads
18,470 L1-dcache-load-misses # 0.00% of all L1-dcache hits
5,717 LLC-loads [66.05%]
0 LLC-load-misses:HG # 0.00% of all LL-cache hits
0.011920301 seconds time elapsed
是的，看起来这可以计算融合域和未融合域的 uops！
按 8 展开根本没有帮助：仍然是 4200 万个周期。 （但减少到 61M insns 和 97M uops，这要归功于更少的循环开销）。整洁，clang 使用 sub $-128, %rsi 而不是 add，因为 -128 适合 imm8，但 +128 不适合。所以我想展开 4 足以使 FP 添加端口饱和。
至于你的 1avg 函数 return 一个浮点数，而不是一个向量，clang 没有自动向量化第一个，但 gcc 可以。它发出一个巨大的序言和结尾来进行标量求和，直到它到达一个对齐的地址，然后在一个小循环中执行 32B AVX vaddps。您说您发现它们的速度差异更大，但您是否可能使用更小的缓冲区进行测试？这将说明矢量代码与非矢量代码相比有很大的加速。
我尝试加快计算放置在数组中的 4d 向量的平均值。这是我的代码：
#include <sys/time.h>
#include <sys/param.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <xmmintrin.h>
typedef float dot[4];
#define N 1000000
double gettime ()
{
struct timeval tv;
gettimeofday (&tv, 0);
return (double)tv.tv_sec + (0.000001 * (double)tv.tv_usec);
}
void calc_avg1 (dot res, const dot array[], int n)
{
int i,j;
memset (res, 0, sizeof (dot));
for (i = 0; i < n; i++)
{
for (j = 0; j<4; j++) res[j] += array[i][j];
}
for (j = 0; j<4; j++) res[j] /= n;
}
void calc_avg2 (dot res, const dot array[], int n)
{
int i;
__v4sf r = _mm_set1_ps (0.0);
for (i=0; i<n; i++) r += _mm_load_ps (array[i]);
r /= _mm_set1_ps ((float)n);
_mm_store_ps (res, r);
}
int main ()
{
void *space = malloc (N*sizeof(dot)+15);
dot *array = (dot*)(((unsigned long)space+15) & ~(unsigned long)15);
dot avg __attribute__((aligned(16)));
int i;
double time;
for (i = 0; i < N; i++)
{
array[i][0] = 1.0*random();
array[i][1] = 1.0*random();
array[i][2] = 1.0*random();
}
time = gettime();
calc_avg1 (avg, array, N);
time = gettime() - time;
printf ("%f\n%f %f %f\n", time, avg[0], avg[1], avg[2]);
time = gettime();
calc_avg2 (avg, array, N);
time = gettime() - time;
printf ("%f\n%f %f %f\n", time, avg[0], avg[1], avg[2]);
return 0;
}
正如您所见，calc_avg1 使用从 0 到 4 的简单循环，calc_avg2 将它们替换为 SSE 指令。我用 clang 3.4 编译这段代码：
cc -O2 -o test test.c
这里是 calc_avgX 函数的反汇编：
0000000000400860 <calc_avg1>:
400860: 55 push %rbp
400861: 48 89 e5 mov %rsp,%rbp
400864: 85 d2 test %edx,%edx
400866: 0f 57 c0 xorps %xmm0,%xmm0
400869: 0f 11 07 movups %xmm0,(%rdi)
40086c: 7e 42 jle 4008b0 <calc_avg1+0x50>
40086e: 48 83 c6 0c add [=12=]xc,%rsi
400872: 0f 57 c0 xorps %xmm0,%xmm0
400875: 89 d0 mov %edx,%eax
400877: 0f 57 c9 xorps %xmm1,%xmm1
40087a: 0f 57 d2 xorps %xmm2,%xmm2
40087d: 0f 57 db xorps %xmm3,%xmm3
400880: f3 0f 58 5e f4 addss -0xc(%rsi),%xmm3
400885: f3 0f 11 1f movss %xmm3,(%rdi)
400889: f3 0f 58 56 f8 addss -0x8(%rsi),%xmm2
40088e: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
400893: f3 0f 58 4e fc addss -0x4(%rsi),%xmm1
400898: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
40089d: f3 0f 58 06 addss (%rsi),%xmm0
4008a1: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008a6: 48 83 c6 10 add [=12=]x10,%rsi
4008aa: ff c8 dec %eax
4008ac: 75 d2 jne 400880 <calc_avg1+0x20>
4008ae: eb 0c jmp 4008bc <calc_avg1+0x5c>
4008b0: 0f 57 c0 xorps %xmm0,%xmm0
4008b3: 0f 57 c9 xorps %xmm1,%xmm1
4008b6: 0f 57 d2 xorps %xmm2,%xmm2
4008b9: 0f 57 db xorps %xmm3,%xmm3
4008bc: f3 0f 2a e2 cvtsi2ss %edx,%xmm4
4008c0: f3 0f 5e dc divss %xmm4,%xmm3
4008c4: f3 0f 11 1f movss %xmm3,(%rdi)
4008c8: f3 0f 5e d4 divss %xmm4,%xmm2
4008cc: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
4008d1: f3 0f 5e cc divss %xmm4,%xmm1
4008d5: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
4008da: f3 0f 5e c4 divss %xmm4,%xmm0
4008de: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008e3: 5d pop %rbp
4008e4: c3 retq
4008e5: 66 66 2e 0f 1f 84 00 nopw %cs:0x0(%rax,%rax,1)
4008ec: 00 00 00 00
00000000004008f0 <calc_avg2>:
4008f0: 55 push %rbp
4008f1: 48 89 e5 mov %rsp,%rbp
4008f4: 85 d2 test %edx,%edx
4008f6: 0f 57 c0 xorps %xmm0,%xmm0
4008f9: 7e 10 jle 40090b <calc_avg2+0x1b>
4008fb: 89 d0 mov %edx,%eax
4008fd: 0f 1f 00 nopl (%rax)
400900: 0f 58 06 addps (%rsi),%xmm0
400903: 48 83 c6 10 add [=12=]x10,%rsi
400907: ff c8 dec %eax
400909: 75 f5 jne 400900 <calc_avg2+0x10>
40090b: 66 0f 6e ca movd %edx,%xmm1
40090f: 66 0f 70 c9 00 pshufd [=12=]x0,%xmm1,%xmm1
400914: 0f 5b c9 cvtdq2ps %xmm1,%xmm1
400917: 0f 5e c1 divps %xmm1,%xmm0
40091a: 0f 29 07 movaps %xmm0,(%rdi)
40091d: 5d pop %rbp
40091e: c3 retq
40091f: 90 nop
结果如下：
> ./test
0.004287
1073864320.000000 1074018048.000000 1073044224.000000
0.003661
1073864320.000000 1074018048.000000 1073044224.000000
所以 SSE 版本快了 1.17 倍。但是当我尝试做看似相同的工作时，即计算数组中单精度标量的平均值（例如此处 SSE reduction of float vector 所述），SSE 版本运行速度快 3.32 倍。这是 calc_avgX 函数的代码：
float calc_avg1 (const float array[], int n)
{
int i;
float avg = 0;
for (i = 0; i < n; i++) avg += array[i];
return avg / n;
}
float calc_avg3 (const float array[], int n)
{
int i;
__v4sf r = _mm_set1_ps (0.0);
for (i=0; i<n; i+=4) r += _mm_load_ps (&(array[i]));
r = _mm_hadd_ps (r, r);
r = _mm_hadd_ps (r, r);
return r[0] / n;
}
所以我的问题是：为什么我在最后一个示例（计算单个浮点标量的平均值）中从 SSE 中获益如此之多，而在第一个示例（计算 4d 向量的平均值）中却没有？对我来说，这些工作几乎是一样的。如果我做错了，在第一个例子中加速计算的正确方法是什么？
更新： 如果您认为相关，我还提供了第二个示例的反汇编，其中计算了标量的平均值（也使用 clang3.4 -O2 编译）。
0000000000400860 <calc_avg1>:
400860: 55 push %rbp
400861: 48 89 e5 mov %rsp,%rbp
400864: 85 d2 test %edx,%edx
400866: 0f 57 c0 xorps %xmm0,%xmm0
400869: 0f 11 07 movups %xmm0,(%rdi)
40086c: 7e 42 jle 4008b0 <calc_avg1+0x50>
40086e: 48 83 c6 0c add [=15=]xc,%rsi
400872: 0f 57 c0 xorps %xmm0,%xmm0
400875: 89 d0 mov %edx,%eax
400877: 0f 57 c9 xorps %xmm1,%xmm1
40087a: 0f 57 d2 xorps %xmm2,%xmm2
40087d: 0f 57 db xorps %xmm3,%xmm3
400880: f3 0f 58 5e f4 addss -0xc(%rsi),%xmm3
400885: f3 0f 11 1f movss %xmm3,(%rdi)
400889: f3 0f 58 56 f8 addss -0x8(%rsi),%xmm2
40088e: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
400893: f3 0f 58 4e fc addss -0x4(%rsi),%xmm1
400898: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
40089d: f3 0f 58 06 addss (%rsi),%xmm0
4008a1: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008a6: 48 83 c6 10 add [=15=]x10,%rsi
4008aa: ff c8 dec %eax
4008ac: 75 d2 jne 400880 <calc_avg1+0x20>
4008ae: eb 0c jmp 4008bc <calc_avg1+0x5c>
4008b0: 0f 57 c0 xorps %xmm0,%xmm0
4008b3: 0f 57 c9 xorps %xmm1,%xmm1
4008b6: 0f 57 d2 xorps %xmm2,%xmm2
4008b9: 0f 57 db xorps %xmm3,%xmm3
4008bc: f3 0f 2a e2 cvtsi2ss %edx,%xmm4
4008c0: f3 0f 5e dc divss %xmm4,%xmm3
4008c4: f3 0f 11 1f movss %xmm3,(%rdi)
4008c8: f3 0f 5e d4 divss %xmm4,%xmm2
4008cc: f3 0f 11 57 04 movss %xmm2,0x4(%rdi)
4008d1: f3 0f 5e cc divss %xmm4,%xmm1
4008d5: f3 0f 11 4f 08 movss %xmm1,0x8(%rdi)
4008da: f3 0f 5e c4 divss %xmm4,%xmm0
4008de: f3 0f 11 47 0c movss %xmm0,0xc(%rdi)
4008e3: 5d pop %rbp
4008e4: c3 retq
4008e5: 66 66 2e 0f 1f 84 00 nopw %cs:0x0(%rax,%rax,1)
4008ec: 00 00 00 00
00000000004008d0 <calc_avg3>:
4008d0: 55 push %rbp
4008d1: 48 89 e5 mov %rsp,%rbp
4008d4: 31 c0 xor %eax,%eax
4008d6: 85 f6 test %esi,%esi
4008d8: 0f 57 c0 xorps %xmm0,%xmm0
4008db: 7e 0f jle 4008ec <calc_avg3+0x1c>
4008dd: 0f 1f 00 nopl (%rax)
4008e0: 0f 58 04 87 addps (%rdi,%rax,4),%xmm0
4008e4: 48 83 c0 04 add [=15=]x4,%rax
4008e8: 39 f0 cmp %esi,%eax
4008ea: 7c f4 jl 4008e0 <calc_avg3+0x10>
4008ec: 66 0f 70 c8 01 pshufd [=15=]x1,%xmm0,%xmm1
4008f1: f3 0f 58 c8 addss %xmm0,%xmm1
4008f5: 66 0f 70 d0 03 pshufd [=15=]x3,%xmm0,%xmm2
4008fa: 0f 12 c0 movhlps %xmm0,%xmm0
4008fd: f3 0f 58 c1 addss %xmm1,%xmm0
400901: f3 0f 58 c2 addss %xmm2,%xmm0
400905: 0f 57 c9 xorps %xmm1,%xmm1
400908: f3 0f 2a ce cvtsi2ss %esi,%xmm1
40090c: f3 0f 5e c1 divss %xmm1,%xmm0
400910: 5d pop %rbp
400911: c3 retq
400912: 66 66 66 66 66 2e 0f nopw %cs:0x0(%rax,%rax,1)
400919: 1f 84 00 00 00 00 00
抱歉，这个回答有点冗长和杂乱无章。我 运行 一些基准测试，但在考虑其他尝试后，我没有花很长时间编辑早期的东西。
您的工作集为 15.25MiB (16MB)。通常要对这样的例程进行基准测试，您会多次平均一个较小的缓冲区，因此它适合缓存。您看不到慢速版本和快速版本之间的太大差异，因为差异被内存瓶颈隐藏了。
calc_avg1 根本不会自动向量化（注意 addss。ss 表示标量，单精度，与 addps 相对（打包单精度-精确））。我认为即使内联到 main 中它也不能自动向量化，因为它不能确定第 4 个向量位置中没有 NaNs，这会导致标量代码没有的 FP 异常。我尝试使用 gcc 4.9.2 -O3 -march=native -ffast-math 和 clang-3.5 为 Sandybridge 编译它，但都没有成功。
即便如此，内联到main的版本也只是运行速度稍慢，因为内存是瓶颈。当访问主内存时，32 位加载几乎可以跟上 128b 加载。 （但是，非内联版本会很糟糕：每个 += 结果都存储到 res 数组中，因为循环直接累积到可能有其他引用的内存中。所以它必须商店可见的每个操作。顺便说一句，这是您发布反汇编的版本。通过 -S -fverbose-asm.)
有点令人失望的是，clang 和 gcc 无法将 __v4sf 从 4 宽 AVX 自动矢量化为 8 宽 AVX。
从头开始，在对 calc_avgX 的调用包装 for (int i=0; i<4000 ; i++) 并将 N 减少到 10k 之后，gcc -O3 将 avg1 的内部循环进入：
400690: c5 f8 10 08 vmovups (%rax),%xmm1
400694: 48 83 c0 20 add [=10=]x20,%rax
400698: c4 e3 75 18 48 f0 01 vinsertf128 [=10=]x1,-0x10(%rax),%ymm1,%ymm1
40069f: c5 fc 58 c1 vaddps %ymm1,%ymm0,%ymm0
4006a3: 48 39 d8 cmp %rbx,%rax
4006a6: 75 e8 jne 400690 <main+0xe0>
$ (get CPU to max-turbo frequency) && time ./a.out
0.016515
1071570752.000000 1066917696.000000 1073897344.000000
0.032875
1071570944.000000 1066916416.000000 1073895680.000000
这很奇怪；我不知道为什么它不只使用 32B 负载。它确实使用了 32B vaddps，这是处理适合 L2 缓存的数据集时的瓶颈。
IDK 为什么它在另一个循环内时设法自动矢量化内部循环。请注意，此 仅 适用于内联到 main 的版本。可调用版本仍然是纯标量的。还要注意只有 gcc 管理这个。 clang 3.5 没有。也许 gcc 知道它将以 return 清零缓冲区的方式使用 malloc（因此它不必担心第 4 个元素中的 NaNs）？
我也很惊讶 clang 的非矢量化 avg1 并不慢，当一切都适合缓存时。 N=10000，重复计数 = 40k。
3.3GHz SNB i5 2500k, max turbo = 3.8GHz.
avg1: 0.350422s: clang -O3 -march=native (not vectorized. loop of 6 scalar addss with memory operands)
avg2: 0.320173s: clang -O3 -march=native
avg1: 0.497040s: clang -O3 -march=native -ffast-math (haven't looked at asm to see what happened)
avg1: 0.160374s: gcc -O3 -march=native (256b addps, with 2 128b loads)
avg2: 0.321028s: gcc -O3 -march=native (128b addps with a memory operand)
avg2: ~0.16: clang, unrolled with 2 dependency chains to hide latency (see below).
avg2: ~0.08: unrolled with 4 dep chains
avg2: ~0.04: in theory unrolled-by-4 with 256b AVX. I didn't try unrolling the one gcc auto-vectorized with 256b addps
所以最大的惊喜是纯标量 clang avg1 代码跟上了 avg2。也许循环携带的依赖链是更大的瓶颈？
perf 显示 clang 的非矢量化 avg1 每个周期有 1.47 个 insns，这可能使端口 1 上的 FP 添加单元饱和。（大多数循环指令都是添加）。
然而，avg2，使用 128b addps 和内存 ope运行d，每个周期仅获得 0.58 insns。将数组大小再减少 10 倍，达到 N=1000，每个周期获得 0.60 insns，可能是因为在 prologue/epilogue 中花费了更多时间。所以我认为循环携带的依赖链存在一个严重的问题。 clang 将循环展开 4，但只使用一个累加器。该循环有 7 条指令，解码为 10 微指令。 （每个 vaddps 是 2，因为它与内存 ope运行d 一起使用，具有 2 寄存器寻址模式，防止微融合。 cmp 和 jne 宏-保险丝）。 http://www.brendangregg.com/perf.html 说 UOPS_DISPATCHED.CORE 的 perf 事件是 r2b1，所以：
$ perf stat -d -e cycles,instructions,r2b1 ./a.out
0.031793
1053298112.000000 1052673664.000000 1116960256.000000
Performance counter stats for './a.out':
118,453,541 cycles
71,181,299 instructions # 0.60 insns per cycle
102,025,443 r2b1 # this is uops, but perf doesn't have a nice name for it
40,256,019 L1-dcache-loads
21,254 L1-dcache-load-misses # 0.05% of all L1-dcache hits
9,588 LLC-loads
0 LLC-load-misses:HG # 0.00% of all LL-cache hits
0.032276233 seconds time elapsed
这证实了我对 uops 分析的 7:10 说明。这实际上与此处的性能问题无关：循环 运行 way 比每个循环上限 4 微指令慢。更改内部循环以获得两个独立的 dep 链使吞吐量翻倍（60M 周期而不是 117M，但 81M insns 而不是 71M）：
for (i=0; i<n-1; i+=2) { // TODO: make sure the loop end condition is correct
r0 += _mm_load_ps (array[i]);
r1 += _mm_load_ps (array[i+1]);
}
r0 += r1;
展开 4（在循环结束时合并 4 个累加器）再次使性能翻倍。 （低至 42M 周期，81M insns，112M uops。）内循环有 4x vaddps -0x30(%rcx),%xmm4,%xmm4（和类似）、2x add、cmp、jl。这种形式的 vaddps 应该微融合，但我仍然看到比指令多得多的微指令，所以我猜 r2b1 算作未融合域的微指令。 （Linux perf 没有针对特定平台硬件事件的良好文档）。 C运行 再次提高 N，以确保它是完全控制所有计数的最内层循环，我看到 uop:insn 比率为 1.39，这与 8 insns、11 uops 非常匹配(1.375)（将 vaddps 算作 2，但将 cmp + jl 算作 1）。我找到了 http://www.bnikolic.co.uk/blog/hpc-prof-events.html，其中包含支持的性能事件的完整列表，包括它们用于 Sandybridge 的代码。 （以及如何为任何其他 CPU 转储 table 的说明）。 （在每个块中查找 Code: 行。您需要一个 umask 字节，然后是代码，作为 perf 的参数。）
# a.out does only avg2, as an unrolled-by-4 version.
$ perf stat -d -e cycles,instructions,r14a1,r2b1,r10e,r2c2,r1c2 ./a.out
0.011331
1053298752.000000 1052674496.000000 1116959488.000000
Performance counter stats for './a.out':
42,250,312 cycles [34.11%]
56,103,429 instructions # 1.33 insns per cycle
20,864,416 r14a1 # UOPS_DISPATCHED_PORT: 0x14=port2&3 loads
111,943,380 r2b1 # UOPS_DISPATCHED: (2->umask 00 -> this core, any thread).
72,208,772 r10e # UOPS_ISSUED: fused-domain
71,422,907 r2c2 # UOPS_RETIRED: retirement slots used (fused-domain)
111,597,049 r1c2 # UOPS_RETIRED: ALL (unfused-domain)
0 L1-dcache-loads
18,470 L1-dcache-load-misses # 0.00% of all L1-dcache hits
5,717 LLC-loads [66.05%]
0 LLC-load-misses:HG # 0.00% of all LL-cache hits
0.011920301 seconds time elapsed
是的，看起来这可以计算融合域和未融合域的 uops！
按 8 展开根本没有帮助：仍然是 4200 万个周期。 （但减少到 61M insns 和 97M uops，这要归功于更少的循环开销）。整洁，clang 使用 sub $-128, %rsi 而不是 add，因为 -128 适合 imm8，但 +128 不适合。所以我想展开 4 足以使 FP 添加端口饱和。
至于你的 1avg 函数 return 一个浮点数，而不是一个向量，clang 没有自动向量化第一个，但 gcc 可以。它发出一个巨大的序言和结尾来进行标量求和，直到它到达一个对齐的地址，然后在一个小循环中执行 32B AVX vaddps。您说您发现它们的速度差异更大，但您是否可能使用更小的缓冲区进行测试？这将说明矢量代码与非矢量代码相比有很大的加速。