Как эффективно добавить два вектора в C++

Question

Предположим, у меня есть два вектора a и b, сохраненные в виде вектора. Я хочу сделать a += b или a +=b * k, где k - это номер.

я могу точно сделать следующее,

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

Но каковы возможные способы легко LEVERAGE инструкции SIMD, такие как SSE2?

Answer 1

Функция SAXPY (с одинарной точностью), DAXPY (двойная точность), CAXPY (комплекс с одинарной точностью) и ZAXPY (комплекс с двойной точностью) вычислить точно выражение, которое вы хотите:

Y = a * X + Y 

где a - скалярная константа, а X и Y - векторы.

Этих функций обеспечиваются BLAS библиотек и оптимизированы для всех практических платформ: для процессоров лучших реализации BLAS являются OpenBLAS, Intel MKL (оптимизированы для x86 процессоров Intel и Xeon Phi сопроцессоров только), BLIS и Apple Accelerate (OS X только); для графических процессоров nVidia смотрите cuBLAS (часть CUDA SDK), для любых графических процессоров - ArrayFire.

Эти библиотеки хорошо оптимизированы и обеспечивают лучшую производительность, чем любая реализация, которую вы можете быстро взломать.

Answer 2

Единственное, что вам нужно должно Нужно включить автоин теризацию с вашим компилятором.

Например, компиляции кода (предполагая, поплавок) с GCC (5.2.0) -O3 производит этот основной цикл

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. Фактически, вы можете получить GCC для разворачивания, добавив -funroll-loops. GCC будет разворачиваться до восьми независимых операций в этом случае unlike in the case when there is a dependency chain.

Одна из проблем, с которыми вы можете столкнуться, состоит в том, что вашему компилятору может потребоваться добавить код, чтобы определить, перекрываются ли массивы и делают две ветви без векторизации, когда они перекрываются, а другая - с векторией, если они не перекрываются. GCC и Clang оба делают это. Но ICC не векторизовать цикл.

ICC 13.0.01 с -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]; 
    } 
} 

НКУ производит такой же код в этом случае, как при использовании __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 не требуется. Niether должен использовать cmpl %edx, %ecx. Я не уверен, как заставить GCC исправить это. У меня была проблема, как раньше, с GCC (Produce loops without cmp instruction in GCC).