An Optimization with SIMD Instructions
本文介绍利用SIMD指令（AVX2）优化运算的一点探索，以单精度float类型的加法和乘法为例。
计算模型：
- hardware: E5 2678 v3 2.5GHz; X99; DDR4 2133 MHz 16GB x1
- software: win10 17134; mingw-W64(gcc) 8.1.0
SIMD和AVX2简介
SIMD (single instrucion multiple data)单指令多数据是用来对多个数据执行单个指令（算术或逻辑）的向量化方法，是X86-64(AMD64)指令集的拓展。由最早的SSE（SSE有多个改进版本，最后的SSE版本是SSE4.2），到后来的AVX，AVX2，最新的是AVX512，由于AVX512在较新和高端的CPU才支持，相较而言，更多的CPU支持AVX2。
AVX2指令集配合16个256位ymm寄存器（低128位为xmm寄存器），可以存储8个4字节(float, int)或4个8字节(double, long)数据，并在一个指令内同时做8个或4个运算。
Basic Solution
现在来考虑一个任务，将长度均为N = 100000的float数组a, b, c做运算(a + b) * c，存放在数组r中。将这个操作重复NTIME = 50000次，记录消耗的时间。
由于这个任务作为编程问题是如此的简单和直接，我们很快写出c程序如code 0: simd0.c，并觉得毫无优化的余地，因为仅有三行代码没有任何多余的操作，只有17,19,20三行代码是在完成任务。其中的变量sum仅仅是用来防止编译器的自动优化会直接去掉重复50000次这个步骤，虽然这个优化很厉害，但是偏离了我们的原意。
使用命令gcc simd0.c -o simd0.exe编译，运行时间为13s 803ms 846μs，取5次计算的中位数，下同。如果我们满足于这个结果，觉得已经很快了，那生活也就没有了乐趣。
1 |
|
GCC自动优化
我们可以用命令gcc -S -o simd0.s simd0.c生成上述程序simd0.exe的汇编代码，看看哪里可以优化，其实这个汇编的蠢不在于没有用SIMD指令，主要问题是有太多的访问内存指令，代码没有贴出来。
如果我们使用-O1参数进行优化，即gcc -O1 -o simd0-O1.exe simd0.c，消耗的时间下降为3s 841ms 346μs。汇编代码如code 1: simd0-O1.s, 其中略去了计时，打印和某些伪指令，下同。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21main:
movl $50000, %edx
jmp .L2
.L6:
addss %xmm1, %xmm1
subl $1, %edx
je .L4
.L2:
movl $0, %eax
.L3:
movss a(%rax), %xmm0
addss b(%rax), %xmm0
mulss c(%rax), %xmm0
movss %xmm0, r(%rax)
addq $4, %rax
cmpq $400000, %rax
jne .L3
jmp .L6
.L4:
cvttss2si %xmm1, %eax
ret
simd0-O1.s的代码和不使用优化的代码相比，主要改进是从内存中取一个数比如a[i]，只需要一次访存，因为指标i存在寄存器%rax中，这一个改进大幅降低了访问指标i和数组元素a[i]和其他数组的时间。simd0-O1.s的代码非常符合我们思维，它使用movss从内存中取一个单精度数，使用addss和mulss执行单精度浮点数的加法和乘法。执行一遍计算后再做第二遍。
使用-O2优化和-O1对这个程序没有太大差异，我们略去。
SIMD的第一尝试
如果我们使用-O3级别优化，就会放发生某些不一样的事情，汇编代码见code 2: simd0-O3.s1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17main:
movl $50000, %edx
.L2:
xorl %eax, %eax
.L3:
movaps b(%rax), %xmm0
addps a(%rax), %xmm0
addq $16, %rax
mulps c-16(%rax), %xmm0
movaps %xmm0, r-16(%rax)
cmpq $400000, %rax
jne .L3
addss %xmm1, %xmm1
subl $1, %edx
jne .L2
cvttss2si %xmm1, %eax
ret
和code 1: simd0-O2.s相比，code 2: simd0-O3.s代码启用了SSE指令集的SIMD指令，它使用完整的xmm寄存器，而不是低32位，进行浮点数的加法和乘法。其中，movaps从内存中取4个单精度浮点数放在xmm寄存器中，或者相反方向，addps和mulps分别对xmm寄存器中的4个单精度数据做加法和乘法。使用SSE指令，运行时间又降为2s 612ms 38μs。
AVX2上场
虽然我们仍然可以使用参数让编译器自动向量化，比如gcc -O3 -maxv2 -o simd1-avx2-O3.exe simd0.c，但是这样的优化结果是不可预测的，因此在编程的时候使用适当的机制或约定引导编译器进行向量化，是更加稳健的方法。
内存对齐
将变量的存储地址对齐到8/16/32/64等整数的倍数，是某些指令的要求，同时可以使寻址更加有效。在GCC环境下，使用__attribute__((aligned(__BIGGEST_ALIGNMENT__)))修饰符引导编译器进行内存对齐。比如声明长度为N的float数组，对齐到64字节1
float a[N] __attribute__((aligned(64)));
向量数据类型
在GCC环境下，可以声明基本数据类型char, int, float, double, long的数组作为向量，来引导编译器使用SIMD指令编译，比如下面的代码声明了一种8个float长的向量类型v8f，1
typedef float v8f __attribute__((vector_size(32)));
向量数据类型可以使用运算符+ - * / & ^ | %等进行逐元素运算，他们会被编译成相应的SIMD指令。
现在我们可以利用上面的技巧，来重写程序，确保我们使用了SIMD指令，见code 3: simd1.c1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
typedef float v8f __attribute__((vector_size(32)));
float a[N] __attribute__((aligned(ALIGN)));
float b[N] __attribute__((aligned(ALIGN)));
float c[N] __attribute__((aligned(ALIGN)));
float r[N] __attribute__((aligned(ALIGN)));
int main(void)
{
v8f *va, *vb, *vc, *vr;
int i, times, sum;
struct timeval starttime, endtime;
long s, ms, us, timeuse;
gettimeofday(&starttime, 0);
/* workhorse code */
for (times = 0; times < NTIMES; times++)
{
va = (v8f *)a;
vb = (v8f *)b;
vc = (v8f *)c;
vr = (v8f *)r;
for (i = 0; i < N; i += VLEN){
*vr = (*va + *vb) * (*vc);
va++, vb++, vc++, vr++;
}
sum += sum;
}
/* elapsed time */
gettimeofday(&endtime, 0);
timeuse = 1000000 * (endtime.tv_sec - starttime.tv_sec) +
endtime.tv_usec - starttime.tv_usec;
s = timeuse / 1000000;
ms = timeuse / 1000 % 1000;
us = timeuse % 1000;
printf("%ld s %ld ms %ld us\n", s, ms, us);
return sum;
}
上面的程序使用指向向量数据类型v8f的指针来进行运算，这些操作可以很好的编译成SIMD指令。并且由于向量长度是8的整数倍，所以我们不用额外的操作。使用命令gcc -mavx2 -O3 -S -o simd1-O3.s simd1.c编译，得到和我们预想一样的汇编代码1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18main:
xorl %eax, %eax
movl $50000, %ecx
.L3:
xorl %edx, %edx
.L2:
vmovaps a(%rdx), %ymm1
vaddps b(%rdx), %ymm1, %ymm0
addq $32, %rdx
vmulps c-32(%rdx), %ymm0, %ymm0
vmovaps %ymm0, r-32(%rdx)
cmpq $400000, %rdx
jne .L2
addl %eax, %eax
subl $1, %ecx
jne .L3
vzeroupper
ret
和code 2: simd0-O3.s不同的是，这里使用指令vmovaps vaddps vmulps对8个（组）float数据进行操作。可惜的是，虽然我们一顿操作猛如虎，但是效率提升并不明显，原因是因为这个例子中，最费时的是内存的读写而不是计算加法和乘法。
如果我们换一个顺序，对于某个i，我们重复NTIMES次r[i] = (a[i] + b[i]) * c[i]操作，然后在进行下一个i，由于编译器优化总是会去掉一层循环，所以我们直接修改汇编代码如下1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18main:
.L2:
movl $50000, %eax
.L3:
vmovaps a(%rdx), %ymm1
vaddps b(%rdx), %ymm1, %ymm0
vmulps c(%rdx), %ymm0, %ymm0
4,,10
addl %ebx, %ebx
vmovaps %ymm0, r(%rdx)
subl $1, %eax
jne .L3
addq $32, %rdx
cmpq $400000, %rdx
jne .L2
xorl %edx, %edx
vzeroupper
ret
用命令gcc -o simd2-0.exe simd2-0.s将汇编指令链接成二进制程序。虽然只是调换了两个循环的顺序，操作也没有少一个，但是由于L1/2/3 cache的存在，code 5显著降低了访问内存的次数，执行时间仅为0s 634ms 562μs，验证了我们上面关于效率瓶颈的论断。
如果我们显式地先读一次数据到寄存器中，做NTIMES次计算，最后再写到内存中，这样就避免了多次读写内存的操作。汇编代码如下1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20main:
.L2:
movl $50000, %eax
vmovaps a(%rdx), %ymm1
vmovaps b(%rdx), %ymm2
vmovaps c(%rdx), %ymm3
.L3:
vaddps %ymm2, %ymm1, %ymm0
vmulps %ymm3, %ymm0, %ymm0
4,,10
addl %ebx, %ebx
vmovaps %ymm0, r(%rdx)
subl $1, %eax
jne .L3
addq $32, %rdx
cmpq $400000, %rdx
jne .L2
xorl %edx, %edx
vzeroupper
ret
最后优化的程序执行时间为0s 324ms 504μs。和最初的13.8s相比，这个速度快了40多倍，虽然主要功劳并不归功于SIMD指令，但是合理的SIMD向量化有助于我们减少内存的访问，从而为程序加速。