c - 用SSE计算4d vector 平均值

Question

我尝试加速计算放置在数组中的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    $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 <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    $0x10,%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 $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 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矢量的平均值计算）中却没有如此受益？对我来说，这些工作几乎是相同的。如果我做错了，第一个示例中加快计算速度的正确方法是什么？

UPD：
如果您认为这是相关的，我还将提供第二个示例的反汇编，该示例计算标量的平均值（也使用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    $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 <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    $0x4,%rax
  4008e8:   39 f0                   cmp    %esi,%eax
  4008ea:   7c f4                   jl     4008e0 <calc_avg3+0x10>
  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 

Answer 1

抱歉，这个答案有点冗长而混乱。我运行了一些基准测试，但是在考虑了其他尝试之后，我花了很长时间来编辑早期的内容。

您的工作集为15.25MiB（16MB）。通常，为了对这样的例程进行基准测试，您将多次平均较小的缓冲区，因此它适合高速缓存。您不会在慢速版本和快速版本之间看到太多差异，因为差异被内存瓶颈所隐藏。

calc_avg1根本不会自动矢量化（请注意addss。ss表示标量，单精度，而不是addps（压缩的单精度））。我认为即使内联到main中也无法自动矢量化，因为无法确定在第4个矢量位置中没有NaN，这会导致标量代码没有的FP异常。我尝试使用gcc 4.9.2 -O3 -march=native -ffast-math和clang-3.5为Sandybridge编译它，但是两者都没有运气。

即使如此，内联到main的版本运行速度也稍慢，因为内存是瓶颈。当访问主内存时，32位负载几乎可以跟上128b负载。 （但是，非内联版本会很糟糕：每个+=结果都存储到res数组中，因为循环直接累积到可能有其他引用的内存中。因此，它必须使每个操作都可见商店。这是您为其发布反汇编的版本，BTW。整理出main的哪个部分是通过-S -fverbose-asm进行编译的。）

令人失望的是，clang和gcc无法自动将__v4sf从4宽AVX矢量化为8宽。

在将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 <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对此进行了管理。铛3.5没有。也许gcc知道它将以返回零缓冲区的方式使用 malloc（所以它不必担心第4个元素中的 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

因此，最大的惊喜是纯标量clang avg1代码与 avg2保持一致。也许循环携带的依赖链是更大的瓶颈？

perf显示的是clang的非矢量化 avg1每个周期1.47 insns，这很可能会使端口1上的FP加法器饱和（大多数循环指令都加了）。

但是，将128b avg2与内存操作数一起使用的 addps每个周期仅获得0.58 insns。将数组大小再减小10倍至 N=1000，每个周期可获得0.60 insns，这可能是因为在序言/结尾中花费了更多时间。我认为循环承载的依赖链存在一个严重的问题。 clang将循环展开4，但仅使用单个累加器。该循环有7条指令，可解码为10微指令。 （每个 vaddps为2，因为它与具有2寄存器寻址模式的内存操作数一起使用，防止了微融合。 cmp和 jne宏保险丝）。 http://www.brendangregg.com/perf.html表示 perf的 UOPS_DISPATCHED.CORE事件是 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

这确认了我的7:10指令分析。这实际上与这里的性能问题无关：循环的运行速度比每个周期上限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个）。内部循环具有4x vaddps -0x30(%rcx),%xmm4,%xmm4（和类似值），2x add， cmp， jl。这种形式的 vaddps应该是微熔丝，但是我仍然看到比指令更多的微指令，因此我想 r2b1会计数未融合的微指令。 （Linux perf没有针对特定平台的硬件事件的出色文档）。再次启动 N，以确保它是最内部的循环完全控制所有计数，我看到uop：insn的比值为1.39，与8 insns，11 uops（1.375）匹配得很好（将 vaddps计为2 ，但将 cmp + jl视为一）。我找到了 http://www.bnikolic.co.uk/blog/hpc-prof-events.html，其中包含受支持的perf事件的完整列表，包括它们对Sandybridge的代码。 （以及有关如何为任何其他CPU转储表的说明）。 （在每个块中查找 Code:行。您需要一个umask字节，然后是代码，作为 perf的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

是的，看起来这可以算出融合域和非融合域！

展开8完全没有帮助：仍为42M周期。 （但由于减少了循环开销，因此可以减少到61M insns和97M uops）。整洁，clang使用 sub $-128, %rsi而不是add，因为-128适合 imm8，但+128不适合。因此，我想展开4足以使FP添加端口饱和。

至于返回单个浮点而不是向量的1avg函数，clang不会自动向量化第一个浮点数，但是gcc会自动向量化。它发出一个巨大的序言和结语以进行标量求和，直到到达对齐的地址，然后在一个小循环中执行32B AVX vaddps。您说您发现它们的速度差异更大，但是您是否可能使用较小的缓冲区进行测试？这将导致矢量代码与非矢量的速度大大提高。