c++ - 如何在 C++ 中有效地添加两个 vector

Question

假设我有两个 vector a 和 b，存储为一个 vector 。我想生成 a += b 或 a +=b * k，其中 k 是一个数字。

我肯定可以做到以下几点，

while (size--) {
    (*a++) += (*b++) * k;
}

但有哪些可能的方法可以轻松利用 SSE2 等 SIMD 指令？

Answer 1

您应该唯一需要的是使用您的编译器启用自动矢量化。

例如，使用 GCC (5.2.0) -O3 编译您的代码(假设为 float)会生成此主循环

L8:
    movups  (%rsi,%rax), %xmm1
    addl    $1, %r11d
    mulps   %xmm2, %xmm1
    addps   (%rdi,%rax), %xmm1
    movaps  %xmm1, (%rdi,%rax)
    addq    $16, %rax
    cmpl    %r11d, %r10d
    ja  .L8

Clang 也对循环进行矢量化，但也展开了四次。即使没有依赖链，展开也可能对某些处理器有所帮助 especially on Haswell .事实上，您可以通过添加 -funroll-loops 让 GCC 展开。在这种情况下，GCC 将展开到八个独立操作 unlike in the case when there is a dependency chain .

您可能会遇到的一个问题是，您的编译器可能需要添加一些代码来确定数组是否重叠，并使两个分支在重叠时不进行矢量化，在不重叠时进行矢量化。 GCC 和 Clang 都这样做。但是 ICC 不会对循环进行矢量化。

带有 -O3 的 ICC 13.0.01

..B1.4:                         # Preds ..B1.2 ..B1.4
        movss     (%rsi), %xmm1                                 #3.21
        incl      %ecx                                          #2.5
        mulss     %xmm0, %xmm1                                  #3.28
        addss     (%rdi), %xmm1                                 #3.11
        movss     %xmm1, (%rdi)                                 #3.11
        movss     4(%rsi), %xmm2                                #3.21
        addq      $8, %rsi                                      #3.21
        mulss     %xmm0, %xmm2                                  #3.28
        addss     4(%rdi), %xmm2                                #3.11
        movss     %xmm2, 4(%rdi)                                #3.11
        addq      $8, %rdi                                      #3.11
        cmpl      %eax, %ecx                                    #2.5
        jb        ..B1.4        # Prob 63%                      #2.5

要解决此问题，您需要使用 __restrict 关键字告诉编译器数组不重叠。

void foo(float * __restrict a, float * __restrict b, float k, int size) {
    while (size--) {
        (*a++) += (*b++) * k;
    }
}

在这种情况下，ICC 产生两个分支。一种用于数组为 16 字节对齐时，另一种用于不对齐时。这是对齐的分支

..B1.16:                        # Preds ..B1.16 ..B1.15
        movaps    (%rsi), %xmm2                                 #3.21
        addl      $8, %r8d                                      #2.5
        movaps    16(%rsi), %xmm3                               #3.21
        addq      $32, %rsi                                     #1.6
        mulps     %xmm1, %xmm2                                  #3.28
        mulps     %xmm1, %xmm3                                  #3.28
        addps     (%rdi), %xmm2                                 #3.11
        addps     16(%rdi), %xmm3                               #3.11
        movaps    %xmm2, (%rdi)                                 #3.11
        movaps    %xmm3, 16(%rdi)                               #3.11
        addq      $32, %rdi                                     #1.6
        cmpl      %ecx, %r8d                                    #2.5
        jb        ..B1.16       # Prob 82%                      #2.5

在这两种情况下，ICC 都会展开两次。即使 GCC 和 Clang 生成了一个没有 __restrict 的矢量化和非矢量化分支，您可能仍然希望使用 __restrict 来消除代码的开销以确定要使用哪个分支。

您可以尝试的最后一件事是告诉编译器数组已对齐。这将适用于 GCC 和 Clang (3.6)

void foo(float * __restrict a, float * __restrict b, float k, int size) {
    a = (float*)__builtin_assume_aligned (a, 32);
    b = (float*)__builtin_assume_aligned (b, 32);
    while (size--) {
        (*a++) += (*b++) * k;
    }
}

GCC 在这种情况下产生

.L4:
    movaps  (%rsi,%r8), %xmm1
    addl    $1, %r10d
    mulps   %xmm2, %xmm1
    addps   (%rdi,%r8), %xmm1
    movaps  %xmm1, (%rdi,%r8)
    addq    $16, %r8
    cmpl    %r10d, %eax
    ja  .L4

最后，如果您的编译器支持 OpenMP 4.0，您可以像这样使用 OpenMP

void foo(float * __restrict a, float * __restrict b, float k, int size) {
    #pragma omp simd aligned(a:32) aligned(b:32)
    for(int i=0; i<size; i++) {
        a[i] += k*b[i];
    }
}

GCC 在这种情况下生成与使用 __builtin_assume_aligned 时相同的代码。这应该适用于较新版本的 ICC(我没有)。

我没有检查 MSVC。我希望它也矢量化这个循环。

有关 restrict 和编译器生成有重叠和无重叠的不同分支以及对齐和不对齐的更多详细信息，请参阅 sum-of-overlapping-arrays-auto-vectorization-and-restrict .

---

这里还有一个要考虑的建议。如果您知道循环范围是 SIMD 宽度的倍数，则编译器将不必使用清理代码。以下代码

// gcc -O3
// n = size/8
void foo(float * __restrict a, float * __restrict b, float k, int n) {
    a = (float*)__builtin_assume_aligned (a, 32);
    b = (float*)__builtin_assume_aligned (b, 32);
    //#pragma omp simd aligned(a:32) aligned(b:32)
    for(int i=0; i<n*8; i++) {
        a[i] += k*b[i];
    }
}

生成迄今为止最简单的程序集。

foo(float*, float*, float, int):
    sall    $2, %edx
    testl   %edx, %edx
    jle .L1
    subl    $4, %edx
    shufps  $0, %xmm0, %xmm0
    shrl    $2, %edx
    xorl    %eax, %eax
    xorl    %ecx, %ecx
    addl    $1, %edx
.L4:
    movaps  (%rsi,%rax), %xmm1
    addl    $1, %ecx
    mulps   %xmm0, %xmm1
    addps   (%rdi,%rax), %xmm1
    movaps  %xmm1, (%rdi,%rax)
    addq    $16, %rax
    cmpl    %edx, %ecx
    jb  .L4
.L1:
    rep ret

我使用了多个8 和 32 字节对齐，因为这样一来，只需使用编译器开关 -mavx，编译器就会生成出色的 AVX 向量化。

foo(float*, float*, float, int):
    sall    $3, %edx
    testl   %edx, %edx
    jle .L5
    vshufps $0, %xmm0, %xmm0, %xmm0
    subl    $8, %edx
    xorl    %eax, %eax
    shrl    $3, %edx
    xorl    %ecx, %ecx
    addl    $1, %edx
    vinsertf128 $1, %xmm0, %ymm0, %ymm0
.L4:
    vmulps  (%rsi,%rax), %ymm0, %ymm1
    addl    $1, %ecx
    vaddps  (%rdi,%rax), %ymm1, %ymm1
    vmovaps %ymm1, (%rdi,%rax)
    addq    $32, %rax
    cmpl    %edx, %ecx
    jb  .L4
    vzeroupper
.L5:
    rep ret

我不确定如何使序言更简单，但我看到的唯一改进是删除了一个迭代器和一个比较。即 addl $1, %ecx 指令应该不是必需的。 cmpl %edx, %ecx 是否必要。我不确定如何让 GCC 解决这个问题。我在使用 GCC ( Produce loops without cmp instruction in GCC ) 时遇到了一个问题。