本文是一篇向量化編程入門文章。

Vectorization is the process of converting an algorithm from a scalar implementation, which does an operation one pair of operands at a time, to a vector process where a single instruction can refer to a vector (a series of adjacent values).

與通常的標量程序不同，向量化的程序可以一個指令處理一個向量（若干個相鄰值），是指令級并行的技術(shù)。向量化可以分為由編譯器自動完成的自動向量化和由開發(fā)者顯式調(diào)用向量化指令完成的顯式向量化。

自動向量化

自動向量化由編譯器嘗試自動使用 SIMD 指令。用戶也可以提供一些附加信息（如 hint 和 pragma）來提示編譯器。

編程指引

1. 適用簡單的 for 循環(huán)，避免復(fù)雜的循環(huán)終止條件
2. straight-line code，避免分支語句
3. 避免數(shù)據(jù) dependency
4. 盡量使用數(shù)組類型，不要通過指針來訪問數(shù)組
5. 通過 loop index 訪問數(shù)組，不要通過單獨的自增變量作為下標訪問數(shù)組
6. 16-byte 對齊（針對 SSE 指令集）
7. 使用 structure of arrays(SoA)，而不是 array of structures(AoS)

性能測試

const int repeat_times = 1000000;
const int sz = 8192;

int a[sz], b[sz], c[sz];
int foo () {
  for (int i=0; i<sz; i++){ // loop vectorized
    a[i] = b[i] + c[i];
  }

  int sum = 0;
  for (int i=0; i<sz; i++){ // loop vectorized
    sum += a[i];
  }
  return sum;
}

int main() {
  int sum = 0;
  for(int i = 0; i < repeat_times; ++i) {
    sum += foo();
  }
  return sum;
}

例子是一個簡單的求和函數(shù)，我們通過這個例子來簡單驗證自動向量化的性能表現(xiàn)。通過設(shè)置 -fopt-info-vec-optimized，在屏幕上打印編譯器的優(yōu)化信息。實際上，-fopt-info-vec- 可以支持 all、note、optimized 多個級別。

> g++ main.cc -fopt-info-vec-optimized -mavx512f  -O3 -o a5.out
main.cc:11:18: optimized: loop vectorized using 64 byte vectors
main.cc:6:18: optimized: loop vectorized using 64 byte vectors
main.cc:11:18: optimized: loop vectorized using 64 byte vectors
main.cc:6:18: optimized: loop vectorized using 64 byte vectors
> time ./a5.out
./a5.out  0.79s user 0.00s system 99% cpu 0.800 total

通過不同的編譯選項，得到最終的執(zhí)行結(jié)果如下，可以看到向量化對于性能的影響還是非常顯著的。

無法向量化的場景

• 不連續(xù)的內(nèi)存訪問
  4個相鄰的整形或者浮點型數(shù)值，或2個相鄰的雙精度浮點型數(shù)值，可以通過一條 SSE 指令加載。在內(nèi)存中的不相鄰的數(shù)據(jù)，若通過 SIMD 指令訪問可能會加載額外的數(shù)據(jù)。甚至，對于循環(huán)步長超過指令集位寬的場景，無法通過向量化加速。

obstacle.cc
const int SIZE = 1024;
int fun(double a[1024], double b[1024], double x[1024], int index[1024]) {
  // arrays accessed with stride 2
  for (int i = 0; i < SIZE; i += 3) b[i] += a[i] * x[i];
  return b[10];
}

gcc -c obstacle.cc -fopt-info-vec-all -O3 -msse
obstacle.cc:5:21: missed: couldn't vectorize loop
obstacle.cc:5:21: missed: Loop costings may not be worthwhile.
obstacle.cc:3:5: note: vectorized 0 loops in function.

• Data Dependency

1. 無依賴，每次循環(huán)寫的變量都不一樣
2. read-after-write dependency，變量在一次循環(huán)迭代中被寫入，在隨后的迭代中被讀取，無法被安全地向量化

for (j=1; j<MAX; j++)  A[j]=A[j-1]+1; 

3. write-after-read dependency，變量在一次循環(huán)迭代中被讀取，在隨后地循環(huán)中被寫入，可以安全地被向量化

for (j=1; j<MAX; j++)  A[j-1]=A[j]+1; 

4. Read-after-read dependency，實際上都不構(gòu)成數(shù)據(jù)依賴
5. Write-after-write dependency，同一個變量在多次循環(huán)迭代中被寫入，不能安全地被向量化

指示編譯器向量化

• pragma 預(yù)處理命令

1. #pragma ivdep
   指示編譯器忽略可能的 data dependency
2. #pragma loop count (n)
   指示編譯器循環(huán)次數(shù)的典型值，幫助其評估向量化的收益
3. #pragma vector always
   指示編譯器向量化循環(huán)
4. #pragma vector align
   斷言循環(huán)中的變量都是內(nèi)存對齊的
5. #pragma novector
   指示編譯器不要向量化循環(huán)
6. #pragma vector nontemporal
   提示編譯器數(shù)據(jù)不會被再次使用，可以使用 stream store 來繞過 cache

注：這里所有的 pragma 指令都是針對 Intel 編譯器 ICC 的，對于 GCC 可能會有不同。如 #pragma ivdep 在 GCC 中對應(yīng)的使用方式是 #pragma GCC ivdep。

• restrict 關(guān)鍵詞
  restrict 關(guān)鍵詞是C99標準引入的，用于指示編譯器通過只能通過該指針訪問指針指向的內(nèi)存。類似于 #pragma ivdep，編譯器可以假定通過指針訪問的內(nèi)存是互斥的，不存在交疊。

顯式向量化

各種向量化的方式

有多種方式可以實現(xiàn)向量化，需要在易用性和程序的自定義程度上進行取舍。
相比于 Auto Vectorization，SIMD vectorization、SIMD intrinsic、Vector intrinsics 都是更為底層的向量化編程方式。SIMD vectorization 使用 #pragma omp simd 預(yù)處理指令命令編譯器進行向量化。SIMD vectorization 在代碼形式上與 Auto vectorization 非常類似。不同的是，前者是程序員顯示命令編譯器進行向量化，當無法向量化代碼時編譯器會返回錯誤導(dǎo)致編譯失??；后者只作為一個hint，即使是 #pragma vector always 下，無法向量化代碼時依然可以編譯成功。

SIMD-enabled 函數(shù)

SIMD-enabled 函數(shù)通過對一個元素的標量操作，描述向量化的算法實現(xiàn)。在程序調(diào)用時，與普通函數(shù)一樣操作單個元素。編譯器實際上會編譯生成針對多個元素的變體，采用 SIMD 指令并行處理。
可以通過 #pragma omp declare simd clauses
的預(yù)處理指令，聲明一個 SIMD-enabled 函數(shù)。

// routine prototype
#pragma omp declare simd                           // universal but slowest definition matches the use in all three loops
#pragma omp declare simd linear(in1) linear(ref(in2)) uniform(mul) // matches the use in the first loop
#pragma omp declare simd linear(ref(in2))                            // matches the use in the second and the third loops
#pragma omp declare simd linear(ref(in2)) linear(mul)              // matches the use in the second loop
#pragma omp declare simd linear(val(in2:2))                          // matches the use in the third loop
extern int func(int* in1, int& in2, int mul);

參考文獻

Vectorization (intel.com)
gcc Multiply.c Driver.c -lm -O3 -fopt-info-vec-all -o a2.out