深入代碼優化 (二) 使用SIMD優化程序

lichwoo 2025-01-08 發布于中國香港

展開全文

在現代 CPU 中，并行性操作大致分為三種類型：

（1）指令級并行，主要由 cpu 流水線技術，亂序執行技術等技術完成。

（2）線程級并行，主要依靠多核多線程技術實現。

（3）數據級并行，主要依靠 SIMD (單指令多數據) 來實現。

指令級并行和線程級并行這兩種技術不在本文進行討論，本文將詳細介紹 SIMD 及其使用方法。

SIMD 介紹

SIMD 是 CPU 硬件設計的一部分，并且是可以通過指令集架構 (ISA) 直接訪問。SIMD 描述了具有多個處理元件的計算機同時對多個數據執行相同操作的過程。

以加法指令為例，單指令單數據 (SISD) 的 CPU 對加法指令譯碼后，執行部件先訪問內存，取得第一個操作數，之后再一次訪問內存，取得第二個操作數，隨后才能進行求和運算。而在支持 SIMD 的 CPU 中，指令譯碼后幾個執行部件同時訪問內存，一次性獲得所有操作數進行運算。這個特點使 SIMD 特別適合于多媒體應用等數據密集型運算。

英特爾的第一個 IA-32 SIMD 指令集是 MMX 指令集。它沿用了 x87 時代的浮點寄存器，使 CPU 無法對浮點數進行 SIMD 操作，只能處理整數。SIMD 流指令擴展 (SSE) 是 x86 架構的 SIMD 指令集擴展，最早是由 Intel 設計并在1999年推出的奔騰3系列 CPU 中引入。SSE 浮點指令在新的獨立寄存器 XMM 上運行，并擴展了一些在 MMX 寄存器上運行的整數指令。SSE 包含了70條新指令，其中大部分都適用于單精度浮點數據類型 (float) 。

SSE最初添加了8個新的128位寄存器，XMM0 - XMM7，在 AMD64 拓展里面，又額外添加了8個寄存器，XMM8 - XMM15，這個拓展也被引入到了 Intel64 位處理器 (IA-64) 架構中。不過寄存器XMM8 - XMM15 只能用來處理 64bit 的操作數。SSE2 進一步支持雙精度浮點數，由于寄存器長度沒有變長，所以只能支持2個雙精度浮點計算或是4個單精度浮點計算。另外，它在這組寄存器上實現了整型計算，從而代替了 MMX。

在 Intel 的 AVX 指令集中，將 SSE 的128位數據通道拓寬到256位，并由此產生的寄存器稱之為YMM。并且 AVX 全面兼容 SSE/SSE2/SSE3/SSE4，也就是 YMM 寄存器的低128位就是 XMM 寄存器。

現代編譯器有三種方式來支持 SIMD：

（1）編譯器能夠在沒有用戶干預的情況下生成 SIMD 代碼，稱之為自動矢量化。

（2）用戶可以插入 Intrinsics 函數實現 SIMD。

（3）用戶可以使用矢量 C++ 類 (僅限ICC編譯器) 來實現 SIMD。

1. 自動矢量化

在程序編寫過程中，可能會經常遇到以下循環的方式，一次執行許多數字的加法。

for (i = 0; i < n; i ++)
    c[i] = a[i] + b[i];

在 gcc 編譯器中，如果添加了 -ftree-vectorize 編譯選項 (O2已包含此優化選項)，那么編譯器可以自動將這類循環轉換成矢量操作序列。

下面以 vector.c 做實驗，看看編譯器怎么實現適量自動化。

[root@wuhan ~]# cat vector.c 
#include <stdio.h>
#define N 100000
#define M 1024
int a[M], b[M], c[M];
int main()
{
    int i, j;
    for (j = 0; j < N; j ++)
        for (i = 0; i < M; i ++)
            c[i] = a[i] + b[i];
    return 0;
}

對其進行編譯運行，第一次編譯不帶 -ftree-vectorize，即期望編譯器不對這段代碼進行自動矢量化。第二次帶上 -ftree-vectorize，期望編譯器能對其進行自動矢量化。由于 -ftree-vectorize 只在 -O1 優化下才生效，所以兩次編譯也都帶了 -O1 進行。

[root@wuhan ~]# gcc vector.c -O1 -o no_vector                   
[root@wuhan ~]# time ./no_vector 
real    0m0.110s
user    0m0.108s
sys     0m0.001s
[root@wuhan ~]# gcc vector.c -O1 -o with_vector -ftree-vectorize
[root@wuhan ~]# time ./with_vector 
real    0m0.029s
user    0m0.028s
sys     0m0.000s

從運行耗時可以看出，矢量化編譯 (-ftree-vectorize) 之后性能大大提升。再將這段程序的匯編碼打印出來，可以看出矢量化編譯之后，匯編碼里面使用了 XMM 寄存器。前面介紹過了，XMM 寄存器可以同時裝載4個 int 類型的數據，并對其進行相同的操作，這也就是性能提升的關鍵。

[root@wuhan ~]# gcc vector.c -O1 -S
[root@wuhan ~]# cat vector.s
.file "vector.c"
.text
.globl main
.type main, @function
main:
.LFB11:
.cfi_startproc
movl $100000, %ecx
jmp .L2
.L6:
subl $1, %ecx
je .L4
.L2:
movl $0, %eax
.L3:
movl b(%rax), %edx
addl a(%rax), %edx
movl %edx, c(%rax)
addq $4, %rax
cmpq $4096, %rax
jne .L3
jmp .L6
.L4:
movl $0, %eax
ret
.cfi_endproc
.LFE11:
.size main, .-main
.comm c,4096,32
.comm b,4096,32
.comm a,4096,32
.ident "GCC: (GNU) 8.4.1 20200928 (Red Hat 8.4.1-1)"
.section .note.GNU-stack,"",@progbits

[root@wuhan ~]# gcc vector.c -O1 -ftree-vectorize -S
[root@wuhan ~]# cat vector.s
.file "vector.c"
.text
.globl main
.type main, @function
main:
.LFB11:
.cfi_startproc
movl $100000, %edx
.L2:
movl $0, %eax
.L3:
movdqa a(%rax), %xmm0
paddd b(%rax), %xmm0
movaps %xmm0, c(%rax)
addq $16, %rax
cmpq $4096, %rax
jne .L3
subl $1, %edx
jne .L2
movl $0, %eax
ret
.cfi_endproc
.LFE11:
.size main, .-main
.comm c,4096,32
.comm b,4096,32
.comm a,4096,32
.ident "GCC: (GNU) 8.4.1 20200928 (Red Hat 8.4.1-1)"
.section .note.GNU-stack,"",@progbits

上面的代碼很容易自動矢量化，我們再來對比一下這兩段代碼：

#include <stdio.h>

#define N 100000
#define M 1024
int a[M+16], b[M], c[M];

int main()
{
int i, j;
for (j = 0; j < N; j ++)
for (i = 0; i < M; i ++)
{
a[i+1] = a[i];
}
return 0;
}

#include <stdio.h>

#define N 100000
#define M 1024
int a[M+16], b[M], c[M];

int main()
{
int i, j;
for (j = 0; j < N; j ++)
for (i = 0; i < M; i ++)
{
a[i+4] = a[i];
}
return 0;
}

區別之處已加粗顯示，分別將他們編成匯編碼，再進行對比：

[root@wuhan ~]# gcc vector2.c -O1 -ftree-vectorize -S

[root@wuhan ~]# cat vector2.s
.file "vector2.c"
.text
.globl main
.type main, @function
main:
.LFB11:
.cfi_startproc
movl $100000, %esi
movl $a+4096, %ecx
jmp .L2
.L6:
subl $1, %esi
je .L4
.L2:
movl $a, %eax
.L3:
movl (%rax), %edx
movl %edx, 4(%rax)
addq $4, %rax
cmpq %rax, %rcx
jne .L3
jmp .L6
.L4:
movl $0, %eax
ret
.cfi_endproc
.LFE11:
.size main, .-main
.comm c,4096,32
.comm b,4096,32
.comm a,4160,32
.ident "GCC: (GNU) 8.4.1 20200928 (Red Hat 8.4.1-1)"
.section .note.GNU-stack,"",@progbits

[root@wuhan ~]# gcc vector2.c -O1 -ftree-vectorize -S

[root@wuhan ~]# cat vector2.s
.file "vector2.c"
.text
.globl main
.type main, @function
main:
.LFB11:
.cfi_startproc
movl $100000, %ecx
movl $a+4096, %edx
.L2:
movl $a, %eax
.L3:
movdqa (%rax), %xmm0
movaps %xmm0, 16(%rax)
addq $16, %rax
cmpq %rax, %rdx
jne .L3
subl $1, %ecx
jne .L2
movl $0, %eax
ret
.cfi_endproc
.LFE11:
.size main, .-main
.comm c,4096,32
.comm b,4096,32
.comm a,4160,32
.ident "GCC: (GNU) 8.4.1 20200928 (Red Hat 8.4.1-1)"
.section .note.GNU-stack,"",@progbits

顯然右側程序的匯編碼使用了 XMM 寄存器，而左邊的程序卻沒有。

假如他們都會進行矢量化，那么以下4條操作是要同時進行的，假如 a[0] = 0, a[1] = 1, a[2] = 2...，那么左邊的程序運行完之后，得到的結果 a[1] = 0, a[2] = 1, a[3] = 2...。但實際上，左邊程序運行完之后應該得到的結果是 a[1] = 0, a[2] = 0, a[3] = 0...。所以，如果左邊的程序也矢量化，那么程序的結果就是錯誤的。而右邊的程序卻不受影響，雖然右邊程序 a[i+4] 的值也依賴于 a[i]，但是他們地址相差128位，而XMM寄存器剛好是128位寬，矢量化運行之后也不影響本來的結果。

a[1] = a[0];

a[2] = a[1];

a[3] = a[2];

a[4] = a[3];

a[4] = a[0];

a[5] = a[1];

a[6] = a[2];

a[7] = a[3];

自動矢量化，就像任何循環優化或其他編譯優化一樣，必須準確地保留程序行為。在執行期間必須遵守所有的依賴項，以防止出現錯誤結果。如果出現處理不了的循環依賴，那么循環依賴必須獨立于矢量化指令執行。

要矢量化一個程序，編譯器必須首先了解語句之間的依賴關系，并在必要時重新對齊它們。一旦映射了依賴關系，編譯器必須正確安排實現指令，將適當的候選者更改為矢量指令，并用這些指令對多個數據項進行操作。

編譯器進行矢量自動化優化通常會經歷一下三個步驟：

（1）建立依賴圖。識別哪些語句依賴于哪些其他語句。這包括檢查每個語句并識別語句訪問的每個數據項，將數組訪問修飾符映射到函數以及檢查每個訪問對所有語句中其他訪問的依賴關系。依賴圖包含了距離不大于矢量大小的所有局部依賴。如果矢量寄存器為 128 位，數組類型為 32 位，則矢量大小為 128/32 = 4。所有其他非循環依賴項都不應使其矢量化無效，因為他們不會調用相同的矢量指令。

（2）聚類。使用依賴圖，優化器可以對強連接組件 (SCC) 進行聚類，并將可矢量化語句與其余語句分開。例如，一個程序的循環內包含三個語句組：(SCC1+SCC2)、SCC3 和 SCC4，其中只有第二組 (SCC3) 可以矢量化。那么最終的程序將包含三個循環，每個循環一個語句組，只有中間的循環被矢量化。優化器不能在不違反語句執行順序的情況下將第一個與最后一個連接起來，因為這很可能會保證不了數據有效性。

（3）監測慣用語法，一些不明顯的依賴可以根據特定的習慣用法進一步優化。例如，下面的數據依賴項可以進行矢量化，因為可以獲取右側值然后將其存儲在左側值上，因此數據不會在賦值中發生變化。

a[i] = a[i+1];  /* 不同于a[i+1] = a[i]， a[i+1] = a[i]不能矢量化 */

2. 插入 Intrinsics 函數 (內在函數) 實現 SIMD

對于程序員來說，intrinsics 看起來就像普通的庫函數。只要包含了相關的頭文件，就可以使用內在函數。如果要將4個整數和另外4個整數相加，可以使用 _mm_add_epi32 內在函數。這個函數的聲明包含在 <emmintrin.h> 頭文件中。

intrinsics 與庫函數不同的是，intrinsics 是直接在編譯器中實現的。上面的 _mm_add_epi32 SSE2內在函數通常編譯成一條指令 paddd。__m128i 內置數據類型是四個整數型的矢量，每個 32 位，總共 128 位。編譯器將發出兩條指令：第一條將參數從內存加載到寄存器中，第二條將四個值相加。通常來說，CPU 調用一個庫函數所花費的時間，可能是調用 intrinsics 的數倍。

包含足夠數量的矢量 intrinsics 或嵌入等效匯編源代碼的過程稱為手動矢量化。現代編譯器和庫已經使用內在函數、匯編或兩者的組合實現了很多東西。例如，memset，memcpy 或 memmove 標準 C 庫函數的實現就使用 SSE2 指令以獲得更好的性能。然而，在高性能計算、游戲開發或編譯器開發等細分領域之外，即使是非常有經驗的 C 和 C++ 程序員在很大程度上也不熟悉 SIMD 內在函數。

下面函數是使用 intrinsics 函數來實現 c[i] = a[i] + b[i] 。

[root@wuhan ~]# cat vector_sse2.c
#include <stdio.h>
#include <emmintrin.h>
#define N 100000
#define M 1024
int a[M], b[M], c[M];
int add_sse2(int size, int *a, int *b, int *c) {
    int i = 0;
    for (; i + 4 <= size; i += 4) {
        /* 加載 a, b 數組的 128 位塊 */
        __m128i ma = _mm_loadu_si128((__m128i*) &a[i]);
        __m128i mb = _mm_loadu_si128((__m128i*) &b[i]);
        /* 128 位塊矢量相加 */
        ma = _mm_add_epi32(ma, mb);
        /* 將相加結果存儲到 c 數組的 128 位塊 */
        _mm_storeu_si128((__m128i*) &c[i], ma);
    }
}
int main()
{
    for (int i = 0; i < N; i ++)
        add_sse2(M, a, b, c);
    return 0;
}

將文件編成匯編碼，匯編碼包含了 paddd %xmm1, %xmm0

[root@wuhan ~]# gcc vector_sse2.c -O1 -S
[root@wuhan ~]# cat vector_sse2.s 
        .file   "vector_sse2.c"
        .text
        .comm   a,4160,32
        .comm   b,4096,32
        .comm   c,4096,32
        .globl  add_sse2
        .type   add_sse2, @function
add_sse2:
.LFB503:
        .cfi_startproc
        pushq   %rbp
        .cfi_def_cfa_offset 16
        .cfi_offset 6, -16
        movq    %rsp, %rbp
        .cfi_def_cfa_register 6
        subq    $40, %rsp
        movl    %edi, -132(%rbp)
        movq    %rsi, -144(%rbp)
        movq    %rdx, -152(%rbp)
        movq    %rcx, -160(%rbp)
        movl    $0, -4(%rbp)
        jmp     .L2
.L6:
        movl    -4(%rbp), %eax
        cltq
        leaq    0(,%rax,4), %rdx
        movq    -144(%rbp), %rax
        addq    %rdx, %rax
        movq    %rax, -128(%rbp)
        movq    -128(%rbp), %rax
        movdqu  (%rax), %xmm0
        movaps  %xmm0, -32(%rbp)
        movl    -4(%rbp), %eax
        cltq
        leaq    0(,%rax,4), %rdx
        movq    -152(%rbp), %rax
        addq    %rdx, %rax
        movq    %rax, -120(%rbp)
        movq    -120(%rbp), %rax
        movdqu  (%rax), %xmm0
        movaps  %xmm0, -48(%rbp)
        movdqa  -32(%rbp), %xmm0
        movaps  %xmm0, -96(%rbp)
        movdqa  -48(%rbp), %xmm0
        movaps  %xmm0, -112(%rbp)
        movdqa  -96(%rbp), %xmm1
        movdqa  -112(%rbp), %xmm0
        paddd   %xmm1, %xmm0
        movaps  %xmm0, -32(%rbp)
        movl    -4(%rbp), %eax
        cltq
        leaq    0(,%rax,4), %rdx
        movq    -160(%rbp), %rax
        addq    %rdx, %rax
        movq    %rax, -56(%rbp)
        movdqa  -32(%rbp), %xmm0
        movaps  %xmm0, -80(%rbp)
        movdqa  -80(%rbp), %xmm0
        movq    -56(%rbp), %rax
        movups  %xmm0, (%rax)
        addl    $4, -4(%rbp)
.L2:
        movl    -4(%rbp), %eax
        addl    $3, %eax
        cmpl    %eax, -132(%rbp)
        jg      .L6
        nop
        leave
        .cfi_def_cfa 7, 8
        ret
        .cfi_endproc
.LFE503:
        .size   add_sse2, .-add_sse2
        .globl  main
        .type   main, @function
main:
.LFB504:
        .cfi_startproc
        pushq   %rbp
        .cfi_def_cfa_offset 16
        .cfi_offset 6, -16
        movq    %rsp, %rbp
        .cfi_def_cfa_register 6
        subq    $16, %rsp
        movl    $0, -4(%rbp)
        jmp     .L8
.L9:
        movl    $c, %ecx
        movl    $b, %edx
        movl    $a, %esi
        movl    $1024, %edi
        call    add_sse2
        addl    $1, -4(%rbp)
.L8:
        cmpl    $99999, -4(%rbp)
        jle     .L9
        movl    $0, %eax
        leave
        .cfi_def_cfa 7, 8
        ret
        .cfi_endproc
.LFE504:
        .size   main, .-main
        .ident  "GCC: (GNU) 8.4.1 20200928 (Red Hat 8.4.1-1)"
        .section        .note.GNU-stack,"",@progbits

測試一下這個程序的性能，比前面介紹的直接使用矢量化優化的程序性能還要好一些。

[root@wuhan ~]# gcc vector.c -O1 -o with_vector -ftree-vectorize
[root@wuhan ~]# time ./with_vector

real 0m0.029s
user 0m0.028s
sys 0m0.000s

[root@wuhan ~]# gcc vector_sse2.c -o vector_sse2 -O1
[root@wuhan ~]# time ./vector_sse2

real 0m0.020s
user 0m0.020s
sys 0m0.000s

3. 利用 C++ 類 (ICC編譯器專用)

使用該方法依賴于環境里面安裝了 Intel ICC 編譯器。ICC 編譯器集成在 Intel oneAPI 開發套件里面，下載鏈接：Download the Intel? oneAPI Base Toolkit

直接在 linux 上用 yum 安裝：

（1）在 /temp 文件夾底下創建 YUM 的 repo 文件

tee > /tmp/oneAPI.repo << EOF
[oneAPI]
name=Intel? oneAPI repository
baseurl=https://yum.repos.intel.com/oneapi
enabled=1
gpgcheck=1
repo_gpgcheck=1
gpgkey=https://yum.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB
EOF

（2）將新建的 repo 文件移到 /etc/yum.repos.d

sudo mv /tmp/oneAPI.repo /etc/yum.repos.d

（3）使用 yum 命令安裝 Intel? oneAPI Base Toolkit

sudo yum install intel-basekit

（4）在使用編譯器之前運行setvars.sh設置環境變量

source /opt/intel/oneapi/setvars.sh intel64

（5）安裝結束，查看編譯器版本

[root@wuhan ~]# icpx --version
Intel(R) oneAPI DPC++/C++ Compiler 2021.4.0 (2021.4.0.20210924)
Target: x86_64-unknown-linux-gnu
Thread model: posix
InstalledDir: /opt/intel/oneapi/compiler/2021.4.0/linux/bin

使用 C++ 類進行 SIMD 操作允許在單個操作中對數組或數據向量進行操作。同樣是計算 c[i] = a[i] + b[i]，用傳統數組的方法表示如下：

int a[4], b[4], c[4]; 
for (i=0; i<4; i++)     /* 需要4次迭代 */ 
    c[i] = a[i] + b[i]; /* 分別計算 c[0], c[1], c[2], c[3] */

在 ICC 編譯器中可以使用 Ivec 類來表示 (需要添加頭文件 dvec.h)

Is32vec4 ivecA, ivecB, ivecC; /* 只需要一次迭代 */ 
ivecC = ivecA + ivecB;        /* 計算 ivecC0, ivecC1, ivecC2, ivecC3 */

所以可以把前面示例的程序改造如下：

[root@wuhan ~]# cat vector_class.cpp 
#include <iostream>
#include <dvec.h>
using namespace std;
#define N 100000
#define M 1024
Is32vec4 a[M/4], b[M/4], c[M/4];
int main()
{
    for (int i = 0; i < N; i ++)
        for (int j = 0; j < M/4; j ++)
            c[j] = a[j] + b[j];
    return 0;
}

把程序變成匯編碼，檢查一下是否用了 XMM 寄存器，C++ 匯編碼編譯出來太長了，就不貼全部源碼。果然 paddd 指令和 XMM 寄存器也都使用了。

[root@wuhan ~]# icpx vector_class.cpp -O1 -S
[root@wuhan ~]# cat vector_class.s | grep -i xmm
        movaps  %xmm0, c(%rbx)
        movdqa  %xmm0, (%rsp)                   # 16-byte Spill
        paddd   (%rsp), %xmm0                   # 16-byte Folded Reload
        movaps  16(%rsp), %xmm0
        movaps  (%rdi), %xmm0
        movaps  %xmm0, (%rdi)