用SSE计算平均4d向量
我尝试加速计算放置在数组中的4d向量的平均值。 这是我的代码：
#include #include #include #include #include #include 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 : 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 40086e: 48 83 c6 0c add $0xc,%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 $0x10,%rsi 4008aa: ff c8 dec %eax 4008ac: 75 d2 jne 400880 4008ae: eb 0c jmp 4008bc 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 : 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 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 $0x10,%rsi 400907: ff c8 dec %eax 400909: 75 f5 jne 400900 40090b: 66 0f 6e ca movd %edx,%xmm1 40090f: 66 0f 70 c9 00 pshufd $0x0,%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减少 ），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向量的平均值）？ 对我来说，这些工作几乎是一样的。 如果我做错了，在第一个例子中加快计算的正确方法是什么？
UPD：如果您认为它是相关的，我还提供了第二个示例的反汇编，其中计算了标量的平均值（也使用clang3.4 -O2编译）。
0000000000400860 : 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 40086e: 48 83 c6 0c add $0xc,%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 $0x10,%rsi 4008aa: ff c8 dec %eax 4008ac: 75 d2 jne 400880 4008ae: eb 0c jmp 4008bc 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 : 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 4008dd: 0f 1f 00 nopl (%rax) 4008e0: 0f 58 04 87 addps (%rdi,%rax,4),%xmm0 4008e4: 48 83 c0 04 add $0x4,%rax 4008e8: 39 f0 cmp %esi,%eax 4008ea: 7c f4 jl 4008e0 4008ec: 66 0f 70 c8 01 pshufd $0x1,%xmm0,%xmm1 4008f1: f3 0f 58 c8 addss %xmm0,%xmm1 4008f5: 66 0f 70 d0 03 pshufd $0x3,%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表示标量，单精度，而不是addps （压缩单精度））。 我认为即使内联到main也无法自动向量化，因为它无法确定第4个向量位置中是否存在NaN ，这会导致标量代码不具备的FPexception。 我尝试使用gcc 4.9.2 -O3 -march=native -ffast-math和clang-3.5为Sandybridge编译它，但是没有运气。
即便如此，内联到main版本的版本运行速度稍慢，因为内存是瓶颈。 在击中主存储器时，32位负载几乎可以跟上128b负载。 （但非内联版本会很糟糕：每个+=结果都存储到res数组中，因为循环直接累积到可能有其他引用的内存中。因此它必须使每个操作都可以通过存储看到。这是您发布反汇编的版本，BTW。通过使用-S -fverbose-asm进行编译来排序主要部分是哪个部分。）
有点令人失望的是，clang和gcc无法将__v4sf从4宽到8宽的AVX自动矢量化。
在包含for (int i=0; i<4000 ; i++)围绕对calc_avgX的调用并将N减少到10k之后，gcc -O3将avg1的内部内循环转换为：
400690: c5 f8 10 08 vmovups (%rax),%xmm1 400694: 48 83 c0 20 add $0x20,%rax 400698: c4 e3 75 18 48 f0 01 vinsertf128 $0x1,-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 $ (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管理这个。 铿锵3.5没有。 也许gcc知道它会以一种返回零化缓冲区的方式使用malloc （所以它不必担心第四元素中的NaN ）？
当一切都适合缓存时，我也很惊讶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
所以最令人惊讶的是，标量avg1代码与avg2保持avg2 。 也许循环携带的依赖链是更大的瓶颈？
对于clang的非矢量化avg1 ， perf每个周期显示1.47个insn，这可能会使端口1上的FP add单元饱和。（大多数循环指令都是加法）。
但是，使用128b addps和内存操作数的addps每个周期只能获得0.58个insn。 将数组大小减少10倍， N=1000 ，每个周期得到0.60个insn，这可能是因为在序言/尾声中花费的时间更多。 因此，我认为循环携带的依赖链存在严重问题。 clang将循环展开4，但只使用一个累加器。 该循环有7条指令，解码为10微秒。 （每个vaddps为2，因为它与具有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指令。 这实际上与此处的性能问题无关：循环运行速度比每循环上限4个uop慢。 更改内部循环以获得两个单独的dep链，使吞吐量增加一倍（60M周期而不是117M，但是81M insns而不是71M）：
for (i=0; i
展开4（在循环结束时合并4个累加器）再次使性能翻倍。 （低至42M周期，81M insn，112M vaddps -0x30(%rcx),%xmm4,%xmm4 。）内环具有4x vaddps -0x30(%rcx),%xmm4,%xmm4 （和类似），2x add ， cmp ， jl 。 这种forms的vaddps应该微融合，但我仍然看到比指令更多的r2b1 ，所以我猜r2b1计算未融合的域r2b1 。 （Linux perf没有针对平台特定的HW事件的良好文档）。 再次开始N ，以确保它是完全控制所有计数的最内层循环，我看到uop：insn比率为1.39，与8个insn匹配，11 vaddps （1.375）（将vaddps计为2，但计算cmp + jl作为一个）。 我找到了http://www.bnikolic.co.uk/blog/hpc-prof-events.html ，其中包含支持的perf事件的完整列表，包括他们的Sandybridge代码。 （以及有关如何为任何其他CPU转储表的说明）。 （在每个块中查找Code:行。您需要一个umask字节，然后是代码，作为一个arg来执行。）
# 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展开根本没有帮助：仍然是42M周期。 （但由于循环开销较少，因此低至61M insn和97M uops）。 整洁，clang使用sub $-128, %rsi而不是add，因为-128适合imm8 ，但+128不适合。 所以我猜4展开就足以让FP添加端口饱和了。
对于返回单个float而不是vector的1avg函数，well clang没有自动向量化第一个，但是gcc可以。 它发出一个巨大的序幕和结尾来做标量和直到它到达一个对齐的地址，然后在一个vaddps中做32B AVX vaddps 。 你说你发现它们有更大的速度差异，但你是否可以使用较小的缓冲区进行测试？ 这将导致矢量代码与非矢量代码相比出现大幅加速。