如何在c++中有效地添加两个向量

How to efficiently add two vectors in C++

本文关键字：两个 向量 添加 c++ 有效地      更新时间：2023-10-16

假设我有两个向量a和b，它们存储为一个向量。我想让a += b或a +=b * k，其中k是一个数字。

我肯定能做到以下几点，

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

但是，有什么可能的方法可以轻松地利用SIMD指令，如SSE2?

您应该唯一需要的是在编译器中启用自动向量化。

例如，使用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也对循环进行矢量化，但也展开了四次。展开在某些处理器上可能会有帮助，即使没有依赖链，特别是在Haswell上。实际上，您可以通过添加-funroll-loops来使GCC展开。在这种情况下，GCC将展开为8个独立的操作，这与存在依赖链的情况不同。

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

ICC 13.0.01 with -O3 

..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中遇到了一个问题(在GCC中生成没有cmp指令的循环)。

函数SAXPY(单精度)、DAXPY(双精度)、CAXPY(复杂单精度)和ZAXPY(复杂双精度)精确计算您想要的表达式:

Y = a * X + Y

其中a为标量常数，X和Y为向量。

这些功能由BLAS库提供并针对所有实际平台进行了优化:对于cpu，最佳BLAS实现是OpenBLAS, Intel MKL(仅针对Intel x86处理器和Xeon Phi协处理器进行了优化)，BLIS和Apple Accelerate(仅限OS X);对于nVidia gpu，请查看cuBLAS (CUDA SDK的一部分)，对于任何gpu - ArrayFire。

这些库经过了很好的优化，提供了比任何实现都更好的性能。

相关文章：