How to efficiently add two vectors in C++

Question

Suppose I have two vectors a and b, stored as a vector. I want to make a += b or a +=b * k, where k is a number.

I can for sure do the following,

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

But what are the possible ways to easily leverage SIMD instructions such as SSE2?

Answer 1

The only thing you should need is to enable auto-vectorization with your compiler.

For example, compiling your code (assuming float) with GCC (5.2.0) -O3 produces this main loop

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 also vectorizes the loop but also unrolls four times. Unrolling may help on some processors even though there is no dependency chain especially on Haswell. In fact, you can get GCC to unroll by adding -funroll-loops. GCC will unroll to eight independent operations in this case unlike in the case when there is a dependency chain.

One problem you may encounter is that your compiler may need to add some code to determine if the arrays overlap and make two branches one without vectorization for when they do overlap and one with vectorization for when they don't overlap. GCC and Clang both do this. But ICC does not vectorize the loop.

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

To fix this you need to tell the compiler the arrays don't overlap using the __restrict keyword.

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

In this case ICC produces two branches. One for when the arrays are 16 byte aligned and one for when they are not. Here is the aligned branch

..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 unrolls twice in both cases. Even though GCC and Clang produce a vectorized and unvectorize branch without __restrict you may want to use __restrict anyway to remove the overhead of the code to determine which branch to use.

The last thing you can try is to to tell the compiler the arrays are aligned. This will work with GCC and 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 produces in this case

.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

Lastly if you compiler supports OpenMP 4.0 you can use OpenMP like this

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 produces the same code in this case as when using __builtin_assume_aligned. This should work for a more recent version of ICC (which I don't have).

I did not check MSVC. I expect it vectorizes this loop as well.

For more details about restrict and the compiler producing different branches with and without overlap and for aligned and not aligned see sum-of-overlapping-arrays-auto-vectorization-and-restrict.

---

Here is one more suggestion to consider. If you know that the range of the loop is a multiple of the the SIMD width the compiler will not have to use cleanup code. The following code

// 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];
    }
}

produces the simplest assembly so far.

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

I used a multiple8 and 32 byte alignment because then just by using the compiler switch -mavx the compiler produces nice AVX vectorization.

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

I am not sure how the preamble could be made simpler but the only improvement I see left is to remove one of the iterators and a compare. Namely the addl $1, %ecx instruction should not be necessary. Niether should the cmpl %edx, %ecx be necessary. I'm not sure how to get GCC to fix this. I had a problem like before with GCC (Produce loops without cmp instruction in GCC).

Answer 2

The functions SAXPY (single-precision), DAXPY (double-precision), CAXPY (complex single-precision), and ZAXPY (complex double-precision) compute exactly the expression you want:

Y = a * X + Y

where a is a scalar constant, and X and Y are vectors.

These functions are provided by BLAS libraries and optimized for all practical platforms: for CPUs the best BLAS implementations are OpenBLAS, Intel MKL (optimized for Intel x86 processors and Xeon Phi co-processors only), BLIS, and Apple Accelerate (OS X only); for nVidia GPUs look at cuBLAS (part of CUDA SDK), for any GPUs - ArrayFire.

These libraries are well-optimized and deliver better performance than whatever implementation you can quickly hack up.