本文是一篇向量化编程入门文章。
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).
与通常的标量程序不同，向量化的程序可以一个指令处理一个向量（若干个相邻值），是指令级并行的技术。向量化可以分为由编译器自动完成的自动向量化和由开发者显式调用向量化指令完成的显式向量化。
自动向量化
自动向量化由编译器尝试自动使用 SIMD 指令。用户也可以提供一些附加信息（如 hint 和 pragma）来提示编译器。
编程指引
- 适用简单的 for 循环，避免复杂的循环终止条件
- straight-line code，避免分支语句
- 避免数据 dependency
- 尽量使用数组类型，不要通过指针来访问数组
- 通过 loop index 访问数组，不要通过单独的自增变量作为下标访问数组
- 16-byte 对齐（针对 SSE 指令集）
- 使用 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;
}
例子是一个简单的求和函数，我们通过这个例子来简单验证自动向量化的性能表现。通过设置 -fopt-info-vec-optimized，在屏幕上打印编译器的优化信息。实际上，-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
通过不同的编译选项，得到最终的执行结果如下，可以看到向量化对于性能的影响还是非常显著的。
| 编译选项 | 执行时间 | 优化方式 |
|---|---|---|
| -O3 -fno-tree-vectorize | 10.31s | 无 |
| -O3 | 2.62s | loop vectorized using 16 byte vectors |
| -msse4.2 -O3 | 2.53s | loop vectorized using 16 byte vectors |
| -mavx2 -O3 | 1.44s | loop vectorized using 32 byte vectors |
| -mavx512f -O3 | 0.79s | loop vectorized using 64 byte vectors |
无法向量化的场景
- 不连续的内存访问
4个相邻的整形或者浮点型数值，或2个相邻的双精度浮点型数值，可以通过一条 SSE 指令加载。在内存中的不相邻的数据，若通过 SIMD 指令访问可能会加载额外的数据。甚至，对于循环步长超过指令集位宽的场景，无法通过向量化加速。
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
- 无依赖，每次循环写的变量都不一样
- read-after-write dependency，变量在一次循环迭代中被写入，在随后的迭代中被读取，无法被安全地向量化
for (j=1; j<MAX; j++) A[j]=A[j-1]+1;
- write-after-read dependency，变量在一次循环迭代中被读取，在随后地循环中被写入，可以安全地被向量化
for (j=1; j<MAX; j++) A[j-1]=A[j]+1;
- Read-after-read dependency，实际上都不构成数据依赖
- Write-after-write dependency，同一个变量在多次循环迭代中被写入，不能安全地被向量化
指示编译器向量化
- pragma 预处理命令
- #pragma ivdep
指示编译器忽略可能的 data dependency - #pragma loop count (n)
指示编译器循环次数的典型值，帮助其评估向量化的收益 - #pragma vector always
指示编译器向量化循环 - #pragma vector align
断言循环中的变量都是内存对齐的 - #pragma novector
指示编译器不要向量化循环 - #pragma vector nontemporal
提示编译器数据不会被再次使用，可以使用 stream store 来绕过 cache
注：这里所有的 pragma 指令都是针对 Intel 编译器 ICC 的，对于 GCC 可能会有不同。如 #pragma ivdep 在 GCC 中对应的使用方式是 #pragma GCC ivdep。
- restrict 关键词
restrict 关键词是C99标准引入的，用于指示编译器通过只能通过该指针访问指针指向的内存。类似于 #pragma ivdep，编译器可以假定通过指针访问的内存是互斥的，不存在交叠。
显式向量化
有多种方式可以实现向量化，需要在易用性和程序的自定义程度上进行取舍。
相比于 Auto Vectorization，SIMD vectorization、SIMD intrinsic、Vector intrinsics 都是更为底层的向量化编程方式。SIMD vectorization 使用 #pragma omp simd 预处理指令命令编译器进行向量化。SIMD vectorization 在代码形式上与 Auto vectorization 非常类似。不同的是，前者是程序员显示命令编译器进行向量化，当无法向量化代码时编译器会返回错误导致编译失败；后者只作为一个hint，即使是 #pragma vector always 下，无法向量化代码时依然可以编译成功。
SIMD-enabled 函数
SIMD-enabled 函数通过对一个元素的标量操作，描述向量化的算法实现。在程序调用时，与普通函数一样操作单个元素。编译器实际上会编译生成针对多个元素的变体，采用 SIMD 指令并行处理。
可以通过 #pragma omp declare simd clauses
的预处理指令，声明一个 SIMD-enabled 函数。
// 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
