OpenMP4.0: #pragma openmp simd實現SIMD指令優化(ARM,X86,MIPS)
考慮一下，CPU一般都是32或64位的暫存器，一次處理的資料長度達到32或64位，對於影象處理來說，一般是每個畫素以8位為單位，那麼我們在對一幅影象每個畫素做處理時，用32位或64位的暫存器來處理8位的資料，其實就是一效能上的浪費。有沒有辦法充分利用CPU 32/64位的處理能能力,讓CPU一次處理多個8位資料呢？這就是本文要說的SIMD.
向量化( Vectorization)
向量化( Vectorization)是一種單指令多資料( Single Instruction Mutiple Data,簡稱SIMD)的並行執行方式。具體而言,向量化是指相同指令在硬體向量處理單元( Vector Processing Unit簡稱VPU)上對多個數據流進行操作。這些硬體向量處理單元也被稱為SIMD單元。
例如,兩個向量的加法形成的第三個向量就是一個典型的SMD操作。許多處理器具有可同時執行2、4、8或更多的SIMD(向量)單元執行相同的操作。
它通過迴圈展開、資料依賴分析、指令重排等方式充分挖掘程式中的並行性,將程式中可以並行化的部分合成處理器支援的向量指令,通過複製多個運算元並把它們直接打包在暫存器中,從而完成在同一時間內採用同步方式對多個數據執行同一條指令,有效地提高程式效能。
還以前面影象處理的應用場景為例，向量化( Vectorization)可以允許一條SIMD指令一次實現多個8位畫素的運算處理。以intel CPU的SSE指令為例，SSE的暫存器達到128bit寬度，一次可以實現16個byte的算術運算。(SSE是Intelr SIMD指令集，進一步，還有升級版的AVX 256bit,和AVX512)。可想而知，在不增加硬體裝置投入的前提下，SIMD對於密集運算程式的效能會帶來數倍乃至數十倍的提升。所以向量化可以充分挖掘處理器並行處理能力,非常適合於處理並行程度高的程式程式碼.
不同的CPU體系的有不同的SIMD指令集標準，比如:
Intel有的x86體系有SSE以及後續的升級版的AVX,AVX2,AVX512 等(參見
arm 平臺也有自己的SIMD指令集,叫NEON(參見《NEON》).
mips體系的SIMD指令集叫MSA(參見《MIPS SIMD》).
看到這裡估計你該頭痛了，SIMD好是好，但這麼多互不相容SIMD指令標準。實際開發中該怎麼用呢？
向量化的實現通常可採用兩種方式:自動向量化和手動向量化.
手動向量化
通過內嵌手工編寫的彙編程式碼或目標處理器的內部函式來新增SIMD指令從而實現程式碼的向量化。
說白了，就是開發者要手工編寫彙編程式使用CPU的SIMD指令來實現向量化( Vectorization)。這要求開發者具備很高的底層彙編開發能力，這個過程對於開發者而言痛苦而低效。而且只能針對特定平臺編寫程式，程式碼不能跨平臺使用，總之代價很高，吃力不討好。
自動向量化
編譯器通過分析程式中控制流和資料流的特徵,識別並選出可以向量化執行的程式碼,並將標量指令自動轉換為相應的SMD指令的過程。
也就是說，向量化的過程由編譯器自動完成，開發者只要編寫正常的C程式碼就好，編譯器會自動分析程式碼結構，將適合向量化的C程式碼部分自動生成SIMD指令的向量化程式碼。而且這些C程式碼可以跨平臺編譯，針對不同的平臺生成不同的SIMD指令。開發者不需要詳細瞭解SIMD指令的用法。也不需要具備彙編程式的編寫能力。
2013年, OpenMP4.0提供了預處理指令simd對函式和迴圈進行向量化。現在主流編譯器都支援了OpenMP4.0(比如gnu,intel Compiler，參見 https://www.openmp.org/resources/openmp-compilers-tools/)。感謝OpenMP4.0，為SIMD指令的跨平臺應用提供了可能。
OpenMP又是啥?
按照Wiki的解釋，OpenMP（Open Multi-Processing）是一套支援跨平臺共享記憶體方式的多執行緒併發的程式設計API，使用C,C++和Fortran語言，可以在大多數的處理器體系和作業系統中執行，包括Solaris, AIX, HP-UX, GNU/Linux, Mac OS X, 和Microsoft Windows。包括一套編譯器指令、庫和一些能夠影響執行行為的環境變數。參見（https://zh.wikipedia.org/wiki/OpenMP）
OpenMP早期是用來實現跨平臺的多執行緒併發程式設計的一套標準。到了OpenMP4.0加入了對SIMD指令的支援，以實現跨平臺的向量化支援。
那麼如何使用OpenMP來實現SIMD指令優化呢（向量化）呢？簡單說只要在程式碼的迴圈邏輯前加入#pragma omp simd預處理指令就可以，不需要任何依賴庫。簡單吧？
#pragma omp simd指令應用於程式碼中的迴圈邏輯，可以讓多個迭代的迴圈利用simd指令實現併發執行。
示例
多說無益，還是舉個栗子吧！
下面就是一個簡單BGRA轉RGB影象的程式，沒有什麼複雜的邏輯，就是把4位元組BGRA格式畫素轉為3位元組的RGB格式畫素。與普通的C程式沒有任何不同，只是在for迴圈前面多了一個#pragma omp simd預處理指令。
這個預處理令告訴編譯器下面這個迴圈要使用SIMD指令來實現向量化。
test_simd.c
/*
* test_simd.c
*
* Created on: Nov 27, 2018
* Author: gyd
*/
#if 1
void bgra2rgb(const char *src,char*dst,int w,int h)
{
#pragma omp simd
for(int y=0;y<h;++y)
{
for(int x=0;x<w;++x)
{
dst[(y*w+x)*3 ] = src[(y*w+x)*4 + 2];
dst[(y*w+x)*3+1] = src[(y*w+x)*4 + 1];
dst[(y*w+x)*3+2] = src[(y*w+x)*4 + 0];
}
}
}
int main()
{
char bgra_mat[480*640*4];
char rgb_mat[480*640*3];
bgra2rgb(bgra_mat,rgb_mat,480,640);
}
#endif
程式部分就這樣了，只是多了一行預處理指令而已，夠簡單吧。重要的是程式碼的編譯方式，以gcc編譯器為例，下面是命令列編譯test_simd.c的過程：
$ gcc -O3 -fopt-info -fopenmp -mavx2 -o test_simd test_simd.c
test_simd.c:13:3: note: loop vectorized
test_simd.c:13:3: note: loop versioned for vectorization because of possible aliasing
上面編譯命令執行時輸出test_simd.c:13:3: note: loop vectorized,就顯示line 13的迴圈程式碼已經實現了迴圈向量化.下面詳細解釋幾個特別的編譯選項的意義：
-
-fopenmp開啟OpenMP預處理指令支援開關，使用此選項，程式碼中的#pragma omp simd預處理指令才有效。
參見 https://gcc.gnu.org/onlinedocs/gcc/C-Dialect-Options.html#C-Dialect-Options -
-mavx2指定使用intel AVX2指令集。如果目標CPU不支援AVX，也可以根據目標CPU的型別改為低版本的-msse4.1 -msse4.2 -msse4 -mavx
參見 https://gcc.gnu.org/onlinedocs/gcc/Option-Summary.html#Option-Summary -
-fopt-info顯示優化過程的輸出，該選項只是用於輸出顯示，指示哪些程式碼已經被優化了，可以不用，就沒有上面的輸出顯示。
參見 https://gcc.gnu.org/onlinedocs/gcc/Developer-Options.html#Developer-Options -
-O33級優化選項，參見 https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#Optimize-Options
對於mips平臺，編譯方式是這樣的,與x86平臺唯一的不同就是-mavx2改為-mmsa(參見 《Option-Summary》)：
$ mips-linux-gnu-gcc -O3 -fopt-info -fopenmp -mmsa -o test_simd_msa test_simd.c
test_simd.c:13:3: note: loop vectorized
test_simd.c:13:3: note: loop versioned for vectorization because of possible aliasing
如果是arm平臺，編譯方式應該是這樣的（我還沒有試過）,參見參考資料5,6：
arm-none-linux-gnueabi-gcc -mfpu=neon -ftree-vectorize -ftree-vectorizer-verbose=1 -c test_simd.c
驗證
如何驗證程式碼是SIMD指令實現的呢？
最直接的辦法 就是檢視生成的可執行檔案的反彙編程式碼。
可以用gdb開啟生成的可執行檔案test_simd，通過檢視生成的指令來驗證是否對迴圈實現了向量化優化。
執行gdb test_simd開啟gdb,再執行disassemble /m bgra2rgb顯示bgra2rgb函式的彙編程式碼，翻幾頁就可以看到類似vmovdqa 0x52f(%rip),%ymm11這樣的指令，像vmovdqa這種v開頭的指令就是AVX2的SIMD指令。代表SIMD指令已經被用於程式中
$ gdb test_simd
GNU gdb (Ubuntu 7.11.1-0ubuntu1~16.5) 7.11.1
Copyright (C) 2016 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law. Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
Type "show configuration" for configuration details.
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<http://www.gnu.org/software/gdb/documentation/>.
For help, type "help".
Type "apropos word" to search for commands related to "word".
(gdb) disassemble /m bgra2rgb
Dump of assembler code for function bgra2rgb:
0x0000000000400660 <+0>: test %ecx,%ecx
0x0000000000400662 <+2>: jle 0x400a1a <bgra2rgb+954>
0x0000000000400668 <+8>: lea 0x8(%rsp),%r10
0x000000000040066d <+13>: and $0xffffffffffffffe0,%rsp
0x0000000000400671 <+17>: lea 0x0(,%rdx,4),%eax
0x0000000000400678 <+24>: xor %r11d,%r11d
0x000000000040067b <+27>: pushq -0x8(%r10)
0x000000000040067f <+31>: push %rbp
0x0000000000400680 <+32>: mov %rsp,%rbp
0x0000000000400683 <+35>: push %r15
0x0000000000400685 <+37>: push %r14
0x0000000000400687 <+39>: push %r13
0x0000000000400689 <+41>: push %r12
0x000000000040068b <+43>: xor %r13d,%r13d
0x000000000040068e <+46>: push %r10
0x0000000000400690 <+48>: push %rbx
0x0000000000400691 <+49>: xor %r10d,%r10d
0x0000000000400694 <+52>: xor %ebx,%ebx
0x0000000000400696 <+54>: mov %eax,-0x34(%rbp)
0x0000000000400699 <+57>: lea (%rdx,%rdx,2),%eax
0x000000000040069c <+60>: vmovdqa 0x41c(%rip),%ymm8 # 0x400ac0
0x00000000004006a4 <+68>: mov %eax,-0x38(%rbp)
---Type <return> to continue, or q <return> to quit---
0x00000000004006a7 <+71>: mov %edx,%eax
0x00000000004006a9 <+73>: lea (%rax,%rax,2),%r15
0x00000000004006ad <+77>: shl $0x2,%rax
0x00000000004006b1 <+81>: mov %rax,-0x40(%rbp)
0x00000000004006b5 <+85>: lea -0x21(%rdx),%eax
0x00000000004006b8 <+88>: shr $0x5,%eax
0x00000000004006bb <+91>: add $0x1,%eax
0x00000000004006be <+94>: mov %eax,-0x54(%rbp)
0x00000000004006c1 <+97>: shl $0x5,%eax
0x00000000004006c4 <+100>: mov %eax,-0x48(%rbp)
0x00000000004006c7 <+103>: lea -0x1(%rdx),%eax
0x00000000004006ca <+106>: mov %eax,-0x44(%rbp)
0x00000000004006cd <+109>: lea (%rax,%rax,2),%rax
0x00000000004006d1 <+113>: mov %rax,-0x50(%rbp)
0x00000000004006d5 <+117>: nopl (%rax)
0x00000000004006d8 <+120>: test %edx,%edx
0x00000000004006da <+122>: jle 0x4009ac <bgra2rgb+844>
0x00000000004006e0 <+128>: movslq %r11d,%r9
0x00000000004006e3 <+131>: movslq %ebx,%r12
0x00000000004006e6 <+134>: lea (%rdi,%r9,1),%r8
0x00000000004006ea <+138>: add -0x40(%rbp),%r9
0x00000000004006ee <+142>: lea (%rsi,%r12,1),%rax
0x00000000004006f2 <+146>: add %rdi,%r9
---Type <return> to continue, or q <return> to quit---
0x00000000004006f5 <+149>: cmp %r9,%rax
0x00000000004006f8 <+152>: lea (%r15,%r12,1),%r9
0x00000000004006fc <+156>: setae %r14b
0x0000000000400700 <+160>: add %rsi,%r9
0x0000000000400703 <+163>: cmp %r9,%r8
0x0000000000400706 <+166>: setae %r9b
0x000000000040070a <+170>: or %r9b,%r14b
0x000000000040070d <+173>: je 0x4009e0 <bgra2rgb+896>
0x0000000000400713 <+179>: cmp $0x1f,%edx
0x0000000000400716 <+182>: jbe 0x4009e0 <bgra2rgb+896>
0x000000000040071c <+188>: xor %r9d,%r9d
0x000000000040071f <+191>: cmpl $0x1f,-0x44(%rbp)
0x0000000000400723 <+195>: jbe 0x40095c <bgra2rgb+764>
0x0000000000400729 <+201>: vmovdqa 0x52f(%rip),%ymm11 # 0x400c60
0x0000000000400731 <+209>: vmovdqa 0x547(%rip),%ymm10 # 0x400c80
0x0000000000400739 <+217>: vmovdqa 0x55f(%rip),%ymm9 # 0x400ca0
0x0000000000400741 <+225>: vmovdqa 0x577(%rip),%ymm7 # 0x400cc0
0x0000000000400749 <+233>: vmovdqa 0x58f(%rip),%ymm6 # 0x400ce0
0x0000000000400751 <+241>: vmovdqa 0x5a7(%rip),%ymm5 # 0x400d00
0x0000000000400759 <+249>: vmovdqa 0x5bf(%rip),%ymm4 # 0x400d20
0x0000000000400761 <+257>: vmovdqu (%r8),%xmm1
0x0000000000400766 <+262>: add $0x1,%r9d
0x000000000040076a <+266>: sub $0xffffffffffffff80,%r8
---Type <return> to continue, or q <return> to quit---
0x000000000040076e <+270>: add $0x60,%rax
0x0000000000400772 <+274>: vmovdqu -0x60(%r8),%xmm13
0x0000000000400778 <+280>: vinserti128 $0x1,-0x70(%r8),%ymm1,%ymm1
0x000000000040077f <+287>: vmovdqu -0x40(%r8),%xmm3
0x0000000000400785 <+293>: vinserti128 $0x1,-0x50(%r8),%ymm13,%ymm13
0x000000000040078c <+300>: vmovdqu -0x20(%r8),%xmm12
0x0000000000400792 <+306>: vinserti128 $0x1,-0x30(%r8),%ymm3,%ymm3
0x0000000000400799 <+313>: vinserti128 $0x1,-0x10(%r8),%ymm12,%ymm12
0x00000000004007a0 <+320>: vpand %ymm13,%ymm8,%ymm2
0x00000000004007a5 <+325>: vpsrlw $0x8,%ymm13,%ymm13
0x00000000004007ab <+331>: vpand %ymm1,%ymm8,%ymm0
0x00000000004007af <+335>: vpsrlw $0x8,%ymm1,%ymm1
0x00000000004007b4 <+340>: vpackuswb %ymm13,%ymm1,%ymm13
0x00000000004007b9 <+345>: vpand %ymm12,%ymm8,%ymm14
0x00000000004007be <+350>: vpsrlw $0x8,%ymm12,%ymm1
0x00000000004007c4 <+356>: vpackuswb %ymm2,%ymm0,%ymm0
0x00000000004007c8 <+360>: vpand %ymm3,%ymm8,%ymm2
0x00000000004007cc <+364>: vpsrlw $0x8,%ymm3,%ymm3
0x00000000004007d1 <+369>: vpackuswb %ymm1,%ymm3,%ymm1
0x00000000004007d5 <+373>: vpermq $0xd8,%ymm13,%ymm13
0x00000000004007db <+379>: vpackuswb %ymm14,%ymm2,%ymm14
0x00000000004007e0 <+384>: vpermq $0xd8,%ymm1,%ymm1
0x00000000004007e6 <+390>: vpand %ymm13,%ymm8,%ymm3
---Type <return> to continue, or q <return> to quit---
0x00000000004007eb <+395>: vpermq $0xd8,%ymm0,%ymm0
0x00000000004007f1 <+401>: vpermq $0xd8,%ymm14,%ymm14
0x00000000004007f7 <+407>: vpand %ymm1,%ymm8,%ymm1
0x00000000004007fb <+411>: vpand %ymm0,%ymm8,%ymm2
0x00000000004007ff <+415>: vpsrlw $0x8,%ymm0,%ymm0
0x0000000000400804 <+420>: vpand %ymm14,%ymm8,%ymm15
0x0000000000400809 <+425>: vpsrlw $0x8,%ymm14,%ymm14
0x000000000040080f <+431>: vpackuswb %ymm1,%ymm3,%ymm1
0x0000000000400813 <+435>: vpackuswb %ymm14,%ymm0,%ymm0
0x0000000000400818 <+440>: vpackuswb %ymm15,%ymm2,%ymm2
0x000000000040081d <+445>: vmovdqa 0x41b(%rip),%ymm15 # 0x400c40
0x0000000000400825 <+453>: vpermq $0xd8,%ymm1,%ymm1
0x000000000040082b <+459>: vpermq $0xd8,%ymm0,%ymm0
0x0000000000400831 <+465>: vpermq $0xd8,%ymm2,%ymm2
0x0000000000400837 <+471>: vpshufb 0x2c0(%rip),%ymm1,%ymm12 # 0x400b00
0x0000000000400840 <+480>: vpshufb 0x297(%rip),%ymm0,%ymm3 # 0x400ae0
0x0000000000400849 <+489>: vpermq $0x4e,%ymm12,%ymm13
0x000000000040084f <+495>: vpermq $0x4e,%ymm3,%ymm14
0x0000000000400855 <+501>: vpshufb 0x2e2(%rip),%ymm1,%ymm12 # 0x400b40
0x000000000040085e <+510>: vpshufb 0x2b9(%rip),%ymm0,%ymm3 # 0x400b2---Type <return> to continue, or q <return> to quit---
如果你不習慣用命令列的gdb工具，也可以用eclipse來檢視反彙編程式碼，如下，在程式中加個斷點，除錯執行到指定的斷點，在Disassembly視窗就可以檢視到對應的彙編程式碼
總結
上面的例子非常簡單，說明#pragma omp simd預處理指令的強大，但這並不是全部,也並不是表面看的那麼簡單，#pragma omp simd不是萬能的，一段迴圈程式碼是不是能被向量化，有不少的限制條件。並不是所有的迴圈都可以直接用#pragma omp simd來向量化優化。關於#pragma omp simd更詳細的說明請參見參考資料2,3。如果你覺得英文看得吃力，建議找本書翻翻，系統化的資料比網上零散的文章看起來更有效率，比如這本《多核異構平行計算（OpenMP4.5C\C++篇）》 ,我也是前幾天從京東買的，寫得一般，不夠通俗，但這樣的系統化中文書籍本身就不多，也只有它了，看看就成。
參考資料：
1.《#pragma omp simd - IBM》
2.《PDF:SIMD Vectorization with OpenMP》
3.《Options Controlling C Dialect.》
4. 《GCC Developer Options》
5. 《ARM NEON Development》
6. 《1.4.3. Automatic vectorization》
7. 《OpenMP in Visual C++》