如何在C++中有效地添加兩個向量

Question

假設我有兩個向量a和b，存儲爲一個向量。我想製作a += b或a +=b * k，其中k是一個數字。

我可以肯定地做到以下幾點，

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

但什麼是可能的方式來輕鬆地利用SIMD指令，如SSE2？

Answer 1

功能SAXPY（單精度），DAXPY（雙精度），CAXPY（複雜的單精度）和ZAXPY（複雜的雙精度）計算正是你想要表達的：

Y = a * X + Y 

其中a是標量常數，而X和Y是向量。

這些功能由BLAS庫提供和優化對於所有的實際平臺：對於CPU最好BLAS實現是OpenBLAS，Intel MKL，BLIS（僅供的Intel x86處理器和Xeon披協處理器優化），和Apple Accelerate（OS X只要）;對於nVidia GPU來說，任何GPU都可以使用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 

鏘也向量化環路而且解開四次。即使沒有依賴鏈especially on Haswell，展開可能有助於某些處理器。實際上，您可以通過添加-funroll-loops來讓GCC展開。在這種情況下，GCC將展開8個獨立的操作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刪除代碼開銷以確定要使用哪個分支。

您可以嘗試的最後一件事是告訴編譯器數組已對齊。這將GCC和鏘（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指令不應該是必要的。 Niether應該cmpl %edx, %ecx是必要的。我不確定如何讓GCC來解決這個問題。我遇到過像GCC一樣的問題（Produce loops without cmp instruction in GCC）。