C++ 编译器优化与SIMD指令集

Reference:

1. 【公开课】编译器优化与SIMD指令集（#4）
2. 10分钟速览 C++20 新增特性
3. 编译器实验用：https://godbolt.org/

相关文章：

1. C++ 字节对齐
2. 向量化运算 和 EIGEN_MAKE_ALIGNED_OPERATOR_NEW

本次的主题是编译器优化。编译器就是从源代码生成汇编语言。所以这里要从汇编语言的角度来理解编译器是怎么优化的，以及我们如何用好它。

1. 汇编语言

首先讲解下什么是汇编语言。

下图为某 x86 的 $64$ 位架构，它有以下寄存器：

寄存器在 CPU 内，它的速度比内存快。如果读到寄存器里，然后再计算就比较高效。所以上图的 x64 就提供了这么多的寄存器，其中白色的这些是 $32$ 位模式下就有的；而 $64$ 位以后，不仅把原来 $32$ 位的这几个扩充到 $64$ 位，它还额外追加了 $8$ 个寄存器，这样我们用起来就更方便了。

图中 RIP 是它当前执行的代码的地址，也是扩充到了 $64$ 位。还有 MMX、TMM、XMM，都是用于存储浮点数的寄存器。一个 XMM 有 $128$ 位宽(YMM $256$ 位，后面还有 ZMM $512$ 位，适用于 AVX-512 寄存器)，而 float 有 $32$ 位宽，所以这里可以塞下四个 float 或者塞两个 double。这样它们在运算加法的时候就可以两个 double 一起来算，这样它效率就更高了。

1.1 通用寄存器：32位时代

开始说的 x86 架构中 $32$ 位，只有以下 $8$ 个寄存器：

• eax, ecx, edx, ebx, esi, edi, esp, ebp

其中：

• esp 是堆栈指针寄存器，和函数的调用与返回相关；
• eax 是用于保存返回值的寄存器。

1.2 通用寄存器：64位时代

到了 $64$ 位以后，又新增了 $8$ 个寄存器。$64$ 位 x86 架构中的通用寄存器有：

• rax, rcx, rdx, rbx, rsi, rdi, rsp, rbp, r8, r9, r10, r11, …, r15

其中：

• r8 到 r15 是 $64$ 位 x86 新增的寄存器，给了汇编程序员更大的空间，降低了编译器处理寄存器翻车(register spill)的压力；

寄存器多有什么好处呢？比如局部变量有 $16$ 个，那它们就都能够存进寄存器里，而不需要存到内存上，这样它读写就更快了。因此 $64$ 位比 $32$ 位机器，除了内存突破 $4$GB 的限制外，也有一定性能优势。

1.3 8位，16位，32位，64位版本

从前两小节可以看到，$32$ 位的叫 eax，而到了 $64$ 位就变成了 rax，它们之间有什么关系呢？它们是共用前 $32$ 位的关系，就是 eax 这个地方的值和 rax 的低 $32$ 位是共用的。同理还有 $16$ 位的 ax 和 eax 共用的低 $16$ 位，然后 ah 是 ax 的高 $8$ 位、al 是低 $8$ 位。

当然以 r 开头的这些数字编号寄存器也是一样的，它们通过 b、w、d、都没有 来区分的。

1.4 AT&T 汇编语言

汇编语言大致分为两种：

1. 一种是 Intel 模式的汇编，它写寄存器就直接写 eax；写操作数，比如 eax 赋值给 edx，就把 edx 写在前面，eax 写在后面，然后用 mov 指令把 eax 赋给 edx。
2. 另一种为 AT&T 汇编。它恰恰相反，它读的数字写在前面，写的数字写在后面，也就是说下图内的 movl %eax, %edx 是往右移动的，Intel 是往左移动的。

• 操作数长度：Intel 就直接 mov，而 AT&T 后面要加一个 movl，代表是 l(long)，也就是 $32$ 位、b(byte) 代表的是 $8$ 位、w(word) 代表的是 $16$ 位。
• 访问地址：Intel 是更直观的一个 []，然后后面接偏移量 lea ax, [dx + 0x10]，而 AT&T 把偏移量写在 () 前面，后面来一个 (%dx) 代表它偏移的寄存器。

总之 AT&T 非常复杂，但是 GCC 用的就是这一种。

1.5 返回值：通过 eax 传出

比如说有一个函数 fun()，右图中 ret 是函数的返回指令，movl 是赋值指令。所以可以看到函数被编译成了把 $42$ 复制给 eax 然后返回，就可以看出这个返回值是通过 eax 传出去的。

int func() {
    return 42;
}

这里是使用下面指令来编译的一个 cpp 文件，变成一个 .S 文件，-S 代表生成的不是 .o 文件，而是 .S 的汇编语言文件。

-fomit-frame-pointer 可以让生成代码更简洁。（作用为：如果函数不需要 frame pointer，就不要将 frame pointer 保留在寄存器中。当打开优化选项：-O, -O2, -O3, -Os 时或者对某些平台不打开任何优化选项时，-fomit-frame-pointer 会被默认打开）

-fverbose-asm 让上图的每一条汇编前面有一个注释，来提示这一条指令代表的是源文件当中的第几行。比如 return 42 就代表了 movl $42 %eax 这一行。

gcc -fomit-frame-pointer -fverbose-asm -S main.cpp -o /tmp/main.S

编译器会自动设置寄存器优化，不需要手动指定。

1.6 前6个参数：分别通过 edi，esi，edx，ecx，r8d，r9d 传入

下面可以看到有一个有很多参数的函数，它有 $6$ 个参数，它们通过寄存器传入的：

可以看到这里显示把所有传入的寄存器先给存到一个堆栈上面，rsp 就代表堆栈，而这个 - 就代表是堆栈上的某一个地址。它的返回值不是返回了 a 嘛，而 edi 它不是存到了 a 变量吗？它现在又取出来设为这个返回值，调用它的人一看 eax 就知道它返回了多少。

movl %edi, -4(%rsp) 相当于： *(rsp - 4) = edi。

1.7 开启优化：-O3

对先前的指令稍加改进，通过加一个 -O3 选项，加了这个 flag 以后，它就比刚才更优化：

gcc -fomit-frame-pointer -fverbose-asm -S main.cpp -o /tmp/main.S

刚才还把 b、c、d、e、f 都传到栈里了，现在直接简单了，直接把 a 传过去，就给送到 eax。它只需要一行指令就够了，这就是编译器的自动优化。它发现这几个存到栈，后来没使用过哇，所以就直接把这些指令全砍了。然后它又一想：欸？我这写入了一遍栈，又读了一遍栈，没意思啊，直接把 edi 复制给 eax 不就好了嘛，所以它就优化成了这样，movl %edi, %eax 相当于： eax = edi：

int func(int a, int b, int c, int d, int e, int f) {
    return a;
}

所以开启这个选项之后，它就能够生成更精炼更高效的代码。

1.8 32位乘法运算：imull

除了可以 movl 来复制之外，汇编语言还支持乘法、加法之类的指令。比如乘法就是 imul，和 movl 是一样的，代表它的操作数是 $32$ 位的。看看下面的这个例子做了些什么：

int func(int a, int b, int c, int d, int e, int f) {
   return a;
}

edi 我们知道代表的第一个参数 a，esi 代表的第二个参数 b。x86 汇编就是这样，它就是就地乘，没有写入到另一个寄存器，它是一种 CISC 指令的设计。这里首先把 edi 赋值给 eax，所以它就相当于，先把 eax 变成了 a，再把 b 乘到 a 里面，就变成所需的返回值了。

1.9 64位乘法运算：imulq

然后就是 $64$ 位的乘法。之前提到了 l 代表 $32$ 位，而 q 代表 $64$ 位，所以这里的 long long 它是 $64$ 位的。$64$ 位寄存器也从 e 寄存器变成 r 系列寄存器，就更宽了，能存储 $64$ 位。

int func(int a, int b, int c, int d, int e, int f) {
    return a;
}

1.10 整数加法：被优化成 leal 了

这里的整数加法我们以为它会生成一个 addl，也就是 add，没想到它却变成了一个 lea 指令。

lea 是什么呢？可以看到 (%rdi,%rsi) 是一个地址的操作数，也就是这个地址其实相当于 *(rdi+rsi)，而 lea 是加载表达式的地址，相当于 &，不是把表达式的值给取出来。它取的是这个地址，所以这个就相当于 eax = &*(rdi + rsi)，也就是这两个会抵消，也就相当于 eax = rdi + rsi。所以说这是在妙用加载地址指令，把它当作 add 来操作，就可以更简练。

1.11 整数加常数乘整数：都可以被优化成 leal

leal 除了可以直接执行加法，它还可以执行任何一次函数。比如下面的 (%rdi,%rsi,8)，就代表把第二个数乘以了 $8$。因为这种线性变换经常用于地址访问中，所以 x86 把它做进了地址访问的操作符内。同上一节，lea 却不读这个地址的值，而是将这个地址读出来。所以说这就变成了，我们可以把任何一个一次函数写在一条指令里。所以说一次函数在 x86 上都是很高效的。

1.12 指针访问对象：线性访问地址

那么线性变换为什么要做成指令呢？就是为了针对下面这种情况，比如 func(int* a, int*b)，a 是一个指针，这里要返回 a[b]，你可能会想 a 是怎么知道 a[b] 的地址，然后读取它的：

这里不能简单的 a+b，因为 a 是一个 int 类型的指针，如果直接在汇编里面 a+b，它会变成一个 char 类型的加。所以说这里要 char 类型的加，加上 sizeof(int) 乘以 b。因为 int 类型的大小是 $4$，所以这个偏移量也要乘以 $4$ 才能访问到正确的地址。所以这里就用了 x86 提供的这个很方便的线性函数来访问。

这里的 moveslq 语句又是什么意思呢？因为我们刚刚使用的 int 是 $32$ 位的，而指针和地址都是 $64$ 位的，所以说要用这个 $64$ 位的 a 去加上一个 $32$ 位的 b 之前，要把 b 转换成 $64$ 位的。可以看到这里是把 esi 这个 $32$ 位变成 rsi 这个 $64$ 位的了。这样以后才能进行线性访问。

1.13 指针的索引：尽量用 size_t

如果要避免上一节的转换呢？就可以使用 <sdtdint> 里提供的 std::size_t。这个 size_t 能够保证在 $64$ 位系统上就是 int64；而 $32$ 位系统上就是 int32，所以就不需要再进行刚才先从 $32$ 扩展到 $64$，它直接就是 $64$ 位的。所以说它有可能更搞笑，而且它还有一个好处：就比如数组的大小超过 INT_MAX，也就是数据大到超过 $2^{31}$，这时候 int 就不仅是效率问题，它是会出错的。所以我们就要用 $64$ 位的 size_t 来表示超过了那个大小的索引。

所以建议如果用于下标的话，最好用 size_t，或者说循环体的 int i=0; i<xxx，也推荐使用 size_t。

1.14 浮点作为参数和返回：xmm 系列寄存器

高性能编程中经常会用到浮点数，所以说 CPU 也对它们做了专门的指令。比如 add 指令，刚才提到 addl 代表两个 $32$ 位整数相加。这里的 $addss$ 的 s 代表 single，也就是单精度浮点数，所以 $addss$ 就代表两个数相加。它传参数也不是 edi、esi 了，它直接用 xmm 系列，就是 xmm0 传 a；xmm1 传 b，然后我们把 xmm1 加到 xmm0，因为正好 xmm0 是用于返回值的，所以我们这里直接返回了，就得到 a+b 了。

1.15 什么是 xmm 系列寄存器？

刚才提到了 xmm 这个系列的寄存器，它们都有 $128$ 位宽，可以容纳 $4$ 个 float 或者 $2$ 个 double。

刚才的例子中，因为只有一个 float 存在一个 $128$ 位的寄存器内，所以只用到了它最低的 $32$ 位。但是这样也没问题，因为我们刚才说的是，addss 它只会加最低位。这就要说到下一节将提到的 addss 了。

1.16 addss 是什么意思？

addss 可以拆成三个部分：add、s、s

1. add 表示执行加法操作。
2. 第一个 s 表示标量(scalar)，只对 xmm 的最低位进行运算；也可以是 p 表示矢量(packed)，一次对 xmm 中所有位进行运算。
3. 第二个 s 表示单精度浮点数(single)，即 float 类型；也可以是 d 表示双精度浮点数(double)，即 double 类型。

举个例子：

• addss：一个 float 加法。
• addsd：一个 double 加法。
• addps：四个 float 加法。
• addpd：两个 double 加法。

比如下图的 xmm0 和 xmm1：

addps 像左边这样同时计算四个的加；但如果是 addss 的话，它只会对最低位进行一个加法，而其他位都保留 xmm0 原来的值。

如果你看到编译器生成的汇编里，有大量 ss 结尾的指令则说明矢量化失败；如果看到大多数都是 ps 结尾则说明矢量化成功，这样生成代码就是比较高效的。

1.17 为什么需要 SIMD？单个指令处理四个数据

这种单个指令处理多个数据的技术称为 SIMD（single-instruction multiple-data），它可以大大增加计算密集型程序的吞吐量。

为什么我们需要这些指令，将四个打包到一起？这样速度更快。

因为 SIMD 把 $4$ 个 float 打包到一个 xmm 寄存器里同时运算，很像数学中矢量的逐元素加法。因此 SIMD 又被称为矢量，而原始的一次只能处理 $1$ 个 float 的方式，则称为标量。

在一定条件下，编译器能够把一个处理标量 float 的代码，转换成一个利用 SIMD 指令的，处理矢量 float 的代码，从而增强你程序的吞吐能力！

通常认为利用同时处理 $4$ 个 float 的 SIMD 指令可以加速 $4$ 倍。但是如果你的算法不适合 SIMD，则可能加速达不到 $4$ 倍；也有因为 SIMD 让访问内存更有规律，节约了指令解码和指令缓存的压力等原因，出现加速超过 $4$ 倍的情况。

2. 化简

然后来看一看编译器还有哪些能做的优化。

2.1 编译器优化：代数化简

首先一件事是，要使它能够执行一些基本代数化。比如看下图左边的例子，这里 $c$ 是 $a+b$，$d$ 是 $a-b$，那么它们相加起来是不是 $b$ 就被抵消了？所以编译器很聪明，它直接把 $a$ 作为返回值。

2.2 编译器优化：常量折叠

如果我看到 $a$ 和 $b$ 都是常量，然后它们还相加起来。我不是已经知道两个常量了，把它们拿过来，然后相加变成 $42$，它就直接优化成 return 42 了。

2.3 编译器优化：举个例子

它更疯狂的是，甚至可以弄一个 for 循环，它看到初始值是一个常数，而 i 每次的值也都是常数，那么它可以把 $1$ 加到 $100$ 给优化成直接返回 $5050$ 都没有问题。

2.4 编译器优化：我毕竟不是万能的

当然编译器也没那么聪明，就是如果把代码写的很复杂，比如用了一些 STL 容器库。当它在看到 arr.push_back() 的时候，就停止思考了。这时它就会放弃优化，干脆去一个个调用 vector。所以说像 vector 这种会在堆上分配内存的容器，编译器通常是不会去进行优化的。所以说尽量把代码写的简单一点，编译器能够看得懂，它才能帮你优化；但如果用的很复杂的话，编译器就认为优化很复杂的抽象语法树需要花很长时间，编译器就觉得不划算，于是放弃优化-----不要把答案藏的太深，它的搜索范围是有限度的，它不会无限制的搜索。所以要把代码简化，从而它的搜索时间就能变短，就能找到正确的优化。

结论：尽量避免代码复杂化，把代码写的简单一点，避免使用会造成 new/delete 的容器。简单的代码，比什么优化手段都强。

2.5 造成 new/delete 的容器：即内存分配在堆上的容器

下表可以看到哪些内存分配在堆上，哪些又在栈上：

存储在堆上（妨碍优化）	存储在栈上（利于优化）
vector, map, set, string, function, any	array, bitset, glm::vec, string_view
unique_ptr, shared_ptr, weak_ptr	pair, tuple, optional, variant

比如 vector、map、set 这种，一般来说这些老的容器都是在堆上的。

为什么要区分堆和栈呢，全部存在栈上不好嘛？事实是有区别的：比如说 array，它是固定大小的，这是它的一个缺点，但是它的优点是它可以存在栈上。那为什么 vector 不能被存在栈上？因为它可以动态的 push_back()，也就是它的大小是变化的，而栈上只能一次性开辟一定的大小，而不能动态扩充大小。所以这时候像 map、set、string 这种可以动态扩展的，得存储在堆上。像 unique_ptr 这种，它肯定是去调用了 new 和 delete 的，所以是在堆上。 (vector 这种肯定也是存在堆上的)

2.6 那改用 array 试试？

既然 array 存储在堆上，那么将之前的示例改写成 array 的，是不是就能优化成功了呢？

从图中可以看到，还是优化失败。

2.6.1 那改用手写的 reduce？

这里再改一下，改用手写的 reduce，它还是优化失败：

2.6.2 那改小到 10？成功了！

那么再修改一下，把这里的 $100$ 改成 $10$，它优化成功了：

这是因为如果代码过于复杂，涉及的语句数量过多时，编译器会放弃优化！

还是那句话，如果代码很简单，那么其实不需要什么优化手段，编译器能够理解你在干什么，那么它就能优化了。

2.7 constexpr：强迫编译器在编译期求值

当然上面的例子中，如果的确想要 $100$ 这么大的能自动在编译器求值，可以考虑使用 constexpr 语法。这时候这个函数强制编译器在运行时求值，这会让编译器花费较长的时间，即编译变慢。所以说如果想迫使编译器优化，即进行常量折叠，就使用 constexpr 关键字。

不过，constexpr 函数中无法使用非 constexpr 的容器，比如前面提到的 vector、map、set、string。因为要让编译器理解全部，它必须全部分配在栈上，而不能分配在堆上。据说 C++20 开始 vector 可以在编译期分配了，new 和 delete 也可以 constexpr 了 (10分钟速览 C++20 新增特性)。但是 C++17 还是不行的，C++17 的 constexpr 还是只能用 std::array 这些在栈上的容器。
你甚至可以把它作为一个模板，即使这个数再大，它也会去算，但是会让编译变得很慢，因为下面示例中的 $50000$ 次迭代实在编译期进行的。

3. 内联

3.1 调用外部函数：call 指令

函数分为两种：外部的函数、内部的函数。

什么是外部函数？就是像下面这种，只要一个申明，然后实际的实现在另一个文件的这种就叫外部函数。这种编译期没办法优化，它只能生成一个 call 指令 来调用这个函数：

@PLT 是 Procedure Linkage Table 的缩写，即函数链接表。链接器会查找其他 .o 文件中是否定义了 _Z5otheri 这个符号，如果定义了则把这个 @PLT 替换为它的地址。可以看到 PLT 在 call _Z5otheri@PLT 出现，代表这个函数是外部函数。它会在链接的时候，如果另一个文件定义了 _Z5otheri 这个符号，它就会把那个符号的地址填充到这里。所以如果别的函数里没有定义 _Z5otheri，它这里就没办法填充，从而链接器就会报错。

gcc -fomit-frame-pointer -fverbose-asm -S main.cpp -o /tmp/main.S

3.1.1 编译器优化：call 变 jmp

上面一小节是没开优化的结果，开优化以后就没有 call 指令了。它直接跳转到 jmp 的 other 这个地址。在前面一小节，调用完 func 一样是要返回的，不如让 other 代为返回，这样也是可以的

gcc -O3 -fomit-frame-pointer -fverbose-asm -S main.cpp -o /tmp/main.S

3.2 多个函数定义在同一个文件中

刚才提到外部函数会让编译器无法优化，但如果是内部函数呢？内部函数是声明和定义在同一个文件，就是它定义在 func 调用它的相同文件。编译器在编译 func 的同时，它是看得到 other 的定义的。也就是说它知道 other 里面什么都没做，直接返回了 $a$，从而它就可以优化了。下图内上面一个红框为 other 的函数，下面的是 func 的函数。
也就是说，这样的话，它本来不优化的话是会直接去调用的，而开了优化之后，定义在同一个文件里是没有 @PLT 的，这样链接器也不用到时候把它替换了。

gcc -fomit-frame-pointer -fverbose-asm -S main.cpp -o /tmp/main.S

3.2.1 编译器优化：内联化

如果在同一个文件里，而且开启了 -O3，这时候我们发现 func() 根本没有调用 other()，它直接返回了 233，这是为什么呢？

因为编译器看到 other 就是直接返回 $a$，所以这里就直接替换成了 return 233 了。

这就是内联：当编译器看得到被调用函数（other）实现的时候，会直接把函数实现贴到调用它的函数（func）里。

没使用关键字 inline 而开启了 -O3 也会内联，编译器只要看得到这个函数体就行了，它自动帮你 inline，不需要去指定。

gcc -O3 -fomit-frame-pointer -fverbose-asm -S main.cpp -o /tmp/main.S

注意：

1. 只有定义在同一个文件的函数可以被内联！否则编译器看不见函数体里的内容是没法内联的。
2. 为了效率我们可以尽量把常用函数定义在头文件里，然后声明为 static(这里的重点是要声明为 static)。这样调用它们的时候编译器看得到它们的函数体，从而有机会内联。

3.2.2 局部可见函数：static

如果 other 声明为 static，就代表只有这文件可见，编译器就干脆不生成 other 函数了，这样的话还是可以内联。如果不是 static，就像上一小节那样，它也是可以内联的，不过会将额外生成的 other 暴露出来。

与上一小节相比，可以看见编译器内少了 _Z5otheri 的定义。这是因为因为 static 声明表示不会暴露 other 给其他文件，而且 func 也已经内联了 other，所以编译器干脆不定义 other 了。

gcc -O3 -fomit-frame-pointer -fverbose-asm -S main.cpp -o /tmp/main.S

3.2.3 inline 关键字？不需要！

在现代编译器中，它们都足够只智能，知道哪些函数要内联，不需要去提醒它，下图为加与不加 inline 的代码：
、
可以看到提醒前和提醒后的代码是完全一样的，编译器都帮着内联了。

结论：
在现代编译器的高强度优化下，加不加 inline 无所谓。编译器不是傻子，只要它看得见 other 的函数体定义，就会自动内联。
内联与否和 inline 没关系，内联与否只取决于是否在同文件，且函数体够小。要性能的，定义在头文件声明为 static 即可，没必要加 inline。static 纯粹是为了避免多个 .cpp 引用同一个头文件造成冲突，并不是必须 static 才内联。如果你不确定某修改是否能提升性能，那你最好实际测一下，不要脑内模拟。
inline 在现代 C++ 中有其他含义，但和内联没有关系，它是一个迷惑性的名字。

在线做编译器实验推荐这个网站：https://godbolt.org/，比如在下面写了一个平方函数，可以在右侧看到它自动变成了图上所示的这些指令。
当然也可以在 Compiler options 处加上优化开关，可以看到它就优化成下面这个样子了：

4. 指针

4.1 指针别名现象（pointer aliasing）

我们知道 C 语言的历史遗产就是指针。下面弄一个 func 函数，可以看到：它将 $a$ 指向的值取出来，然后写入到 $c$，再把 $b$ 指向的值写入到 $c$。那你肯定想：先把 $a$ 复制到 $c$，再把 $b$ 复制到 $c$，那第一个语句不就浪费了嘛。我直接将 $b$ 复制给 $c$ 不就行了，因为 $a$ 的值被覆盖了。那么为啥编译器明明是开了优化它还是复制了两遍呢？
考虑这种调用方式，$b$ 和 $c$ 指向同一个变量：
优化前相当于：$b=a;b=b;$ 最后 $b$ 变成了 $a$；
如果优化了：$b=b;$ 这样 $b$ 就没有发生改变，这就是编译器放弃优化的原因。

编译器并不知道 $a$ 和 $b$ 是否指向了同一个地址，所以如果它优化掉的话，结果就不对了。

编译器是保守的，它宁愿变慢也不要出错。

4.1.1 告诉编译器别怕指针别名：__restrict 关键字

如果你想让编译器放心大胆，告诉它 $b$ 永远不会等于 $c$ 的，可以这样写：

这时 gcc 特有的 __restrict 关键字，它在 Ｃ语言里是标准，而在 C++ 里却没有标准，所以需要加上 __ 前缀代表它是当前编译器特定的。__restrict 是一个提示性的关键字，是程序员向编译器保证：这些指针之间不会发生重叠。所以 __restrict 提示以后它不会再重复两遍了， 它只会重复一遍。也就是说，向编译器保证 $b$ 和 $c$ 是不同的值，从而它就可以安全的把 *c=*a 给优化掉

4.1.2 __restrict 关键字：只需加在非 const 的即可

那是不是在程序里每一个出现指针的地方都得加个 __restrict 呢？不一定，只需要非 const 的加一个就行了。因为只有写入才会造成指针别名的问题，而如果只读的话就没有任何影响。即：所有非 const 的指针都声明 __restrict。

4.1.3 禁止优化：volatile

除此之外，还有一个让编译器不要优化的关键字，就是 volatile。

可以看到上面示例中，先是把 $a$ 赋值为 $42$，然后有返回了 *a。也就是说，这时候先把 $42$ 写入，然后这里直接被编译器优化成 return 42 了。因为它想刚刚我写入的 $42$ 嘛，但是万一这个 $a$ 指针是指向一个多线程同时在写入的地方，也就是比如 *a=42 我赋值 $42$ 了以后，又有某个线程进来，把 $a$ 修改了，这时候这个返回值就可能是不一样的。所以这时候可以告诉编译器不要优化这个变量，就是 int volatile：

这时候它就会硬生生去写一遍再读一遍，生怕函数内再插入了什么东西。如果想优化的话，就不要加 volatile 了。如果在使用多线程后想要测试这个内存读写有多快，就可以使用关键字 volatile，它不会优化读写，从而可以测试它的性能。

4.1.4 注意一下区别

• volatile：
  volatile int *a 或 int volatile *a

• __restrict：
  int *__restrict a

• 语法上区别：volatile 在 * 前面而 __restrict 在 * 后面。

• 功能上区别：volatile 是禁用优化，__restrict 是帮助优化。

• 是否属于标准上区别：volatile 和 const 一样是 C++ 标准的一部分。restrict 是 C99 标准关键字，但不是 C++ 标准的关键字。__restrict 其实是编译器的“私货”，好在大多数主流编译器都支持。

• volatile 是可以针对不是指针的，而 __restrict 必定是针对指针。

4.2 编译器优化：合并写入

说道指针，还有一点就是，下面示例中有两个写入，第一个是针对 a[0]，然后是 a[1]，这两个分别是 int，也就是 $32$ 位的写入。在编译出来之后就变成了一个 q，也就是变成了一个 $64$ 位的写入。因为这里的两个 $32$ 位是连续的地址，所以它能够自动优化成一个 $64$ 位的写入：

4.2.1 合并写入：不能跳跃

刚才说的写入能够变成单独一个 $64$ 位的，但如果地址中间跳了，变成了 a[0] 和 a[2]，中间 a[1] 空了，那它没办法优化了，只能变成写入两个 $32$ 位的了。
所以说我们在设计数据结构的时候，尽可能把它设计的紧凑点，不要把很多无关的变量交错排列，这样它就很难优化了。

5. 矢量化

5.1 更宽的合并写入：矢量化指令（SIMD）

前面说道可以将 $2$ 个 int 变成一个int64，其实还可以把 $4$ 个 int 优化成一个 __m128，也就是 xmm 寄存器。比如下图中有 $4$ 个 int，因为 $4$ 个 $32$ 位就是 $128$ 位，所以就可以用 xmm 来实现一次写入 $4$ 个 int：
xmm 是 SSE 引入的，它是 $128$ 位的。它本来是用于存储 $4$ 个 float 的，但它也可以存 $4$ 个 int。

这里还有一个 movups：move unaligned packed single，之前都是没有这个 u 的，这个 u 是特有的，它代表上面 %rdi 这个地址不一定对齐到 $16$ 字节。也可以是 movaps：move aligned packed single，这时候就是向它保证 %rdi 一定是对齐到 $16$ 字节的，这样可能有微乎其微的提升。

5.1.1 SIMD 指令：敢不敢再宽一点？

刚才说的是 $4$ 个可以变成 $128$，然后 $8$ 个就可以变成 $256$，这个是 ymm 寄存器。但我为什么明明像下面这样写，它却用了两个 xmm 的寄存器，而没有优化成一个 ymm 呢？

这是因为编译器不敢保证你的电脑支持 AVX。也就是说编译器它只能生成、只能保证因为所有 $64$ 位的电脑都支持 xmm，所以它只能保证 cmm 是能用的，而 $256$ 位的 ymm 只有一些比较新的电脑上才有，它就不敢用了。

5.1.2 让编译器自动检测当前硬件支持的指令集

所以说你要让它敢用，可以加以下 flag，跟它说，你来检测一下我的电脑是否支持 AVX，如果我的电脑支持，那你就可以用 ymm 优化了：

gcc -march=native -O3

比如上面的示例中写了 $8$ 个，编译器也的确用 ymm 读取了一次，然后写入。这样就更加高效了。

-march=native：让编译器自动判断当前硬件支持的指令。不过注意这样编译出的程序，可能放到别人不支持 AVX 的电脑上没法运行。

所以如果用于发布到特定平台的话，最好还是不要用这个 flag。

5.2 数组清零：自动调用标准库的 memset

下面为数组清零的代码，可以看到代码内除了清零什么都没做，也没有调任何标准库函数，而在右边汇编内，它缺变成了一个对 memset 的调用，这是为什么呢？

这是因为编译器是和标准库绑定的。所以说在这里生成一个像是把整个数组清零操作，它会识别到，然后就能够一个对标准库的调用了。这样就不必为了高效，手动改写成对 memcpy/memset 的调用，影响可读性。编译器会自动分析是在做拷贝还是清零，并优化成对标准库 memcpy/memset 的调用，标准库里面的实现通常是更高效的。

5.3 从 0 到 1024 填充：SIMD 加速

现在有一个从 $0$ 到 $1024$，然后填入到数组 a[] 中:

void func(int *a) {
    for (int i = 0; i < 1024; i++) {
        a[i] = i;
    }
}

它的汇编如下图所示：

其中 paddd 表示 $4$ 个 int 的加法；movdqa 表示加载 $4$ 个 int。

是不是挺难看懂的。实际上它就是把代码优化成下面这个方便理解的伪代码：

void func(int *a) {
    __m128i curr = {0, 1, 2, 3};
    __m128i delta = {4, 4, 4, 4};
    for (int i = 0; i < 1024; i += 4) {
        a[i : i + 4] = curr;
        curr += delta;
    }
}

这里 i++ 变成 i+=4 了，就是一次可以写入 $4$ 个，然后我们这里在 $a[i]$ 写入 $0,1,2,3$ 这四个数。在下一步的时候把 curr 给全部加上 $4$，curr 就变成 $4,5,6,7$ 了；这时候再写入就是 $4,5,6,7$，再加一下，就是 $8,9,10,11$。这样和原来的运行结果完全一致，但却能更加高效，因为它是用 SIMD 指令并行去加的。

5.4 如果不是 4 的倍数？边界特判法

上面的示例中有一个缺点，因为这里是 i+=4，如果上面不是 $1024$ 而是 $1023$ 的话，那它这里就会多写入一个 int，而这个地址可能是别人的数据就把它覆盖掉，这该怎么办呢？

void func(int *a, int n) {
    for (int i = 0; i < n; i++) {
        a[i] = i;
    }
}

可以用边界特判法。刚才是一个固定的 $1024$，编译器看得到，它是一个 $4$ 的倍数，可以直接像上一小节那样写。但在边界不是 $4$ 的倍数的时候，它就生成了一个超复杂的代码：

它做的事是，假如 n 是 $1023$ 的话，它会先把 n 编程一个刚才优化过的 $1020$ 个元素的填入，这样每次可以处理四个。然后还剩下三个是没办法直接一次写入了，它就变成了传统的标量模式，来一个处理一个。所以这种思想就是，如果 n 不是 $4$ 的倍数，就先取其中和 $4$ 整除的部分，然后剩下不整除的部分，再用低效的标量代码去处理。这样大部分 $1020$ 个数据都是矢量化的加，而剩下的 $3$ 个数据才是标量加。所以它既高效又不会犯错，这就是边界特判。

编译器做优化时，它会自己去判断。但如果要自己手动去写 mm 什么的 SIMD 指令的话，需要手动去判断这个边界。

5.4.1 n 总是 4 的倍数？避免边界特判

如果能保证 $n$ 总归是 $4$ 的倍数，不想让它这么复杂的生成两个版本，可以这样写：

void func(int *a, int n) {
    n = n / 4 * 4;
    for (int i = 0; i < n; i++) {
        a[i] = i;
    }
}

用 n 先除以 $4$，因为除是整除，就让 n 会落入 $4$ 的倍数中，这样的话编译器能够发现 n 始终是 $4$ 的倍数，即 $n\%4=0$，从而它就不会生成边界特判的分支了，相当智能：

5.5 假定指针是 16 字节对齐的：assume_aligned

GCC 内置的函数，都是以两个下划线 __，然后 builtin 开头的。(int *)__builtin_assume_aligned(a, 16) 是在告诉编译器，保证 $a$ 是 $16$ 组合对齐的。

void func(int *a, int n) {
    n = n / 4 * 4;
    a = (int *)__builtin_assume_aligned(a, 16);
    for (int i = 0; i < n; i++) {
        a[i] = i;
    }
}

来看看这与上面的示例有什么区别：

先前这里是 movups，这里变成了 movaps。这就让 CPU 知道，保证是 $16$ 字节对齐的，传 CPU 可能更快一点，也可能不快。但是这样的话就是调用 GCC 才有的函数，如果在微软编译器上，就不能通过。

所以为了通用，C++20 在 <memory> 头文件内引入了这个模板函数 std::assume_aligned，然后通过模板参数来制定对齐到 $16$ 字节：

#include <memory>
void func(int *a, int n) {
    n = n / 4 * 4;
    a = std::assume_aligned<16>(a);
    for (int i = 0; i < n; i++) {
        a[i] = i;
    }
}

5.9 数组求和：reduction 的优化

float func(float *a) {
    float ret = 0;
    for (int i = 0; i < 1024; i++) {
        ret += a[i];
    }
    return ret;
}

优化成了下面这个样子：

伪代码如下：

#include <x86intrin.h>

float func(float *a) {
    __m128 ret = _mm_setzero_ps();
    for (int i = 0; i < 1024; i += 4) {
        __m128 a_i = _mm_loadu_ps(&a[i]);
        ret = _mm_add_ps(ret, a_i);
    }
    float r[4];
    _mm_storeu_ps(r, ret);
    return r[0] + r[1] + r[2] + r[3];
}

总之就是非常高效。

6. 循环

我们的程序中有 热 和 冷 这两个术语，这里的热是指经常被访问的代码。比如上面一节中 float ret=0，它只会调用一次，而这个 for 循环里面被调用了 $1024$ 次。这个调用次数多的地方称为热，而调用数少的地方称为冷。一般都是在热的地方优化，因为热的地方通常它执行很多遍，更容易成为瓶颈，而往往循环的地方都能执行很多遍。所以说优化通常都是针对循环的优化，而循环也往往是优化的重点。

6.1 循环中的矢量化：还在伺候指针别名

这里提供了一个用循环的 func 函数，它的工作就是把 b[] 里面的数 $+1$，然后复制到 a[]。

void func(float *a, float *b) {
    for (int i = 0; i < 1024; i++) {
        a[i] = b[i] + 1;
    }
}

来看看这会发生什么：

按照我们的理解，这个代码应该是可以矢量化的，但是编译器却生成了两个版本，这是为什么？编译器还在担心前面章节提到的，$a$ 和 $b$ 只想的数组是否有重合。考虑 func(a,a+1) 的情况，那样会产生数据依赖链，没法 SIMD 化。

为了优化而不失正确性，编译器索性生成两份代码。 一份是 SIMD 的，一份是传统标量的：
它在运行时检测 $a$, $b$ 指针的差是否超过 $1024$ 来判断是否有重叠现象。

1. 如果没有重叠，则跳转到 SIMD 版本高效运行。
2. 如果重叠，则跳转到标量版本低效运行，但至少不会错。

6.1.1 循环中的矢量化：解决指针别名

所以我们为了让编译器放心，会加上之前提到的关键字 __restrict，告诉编译器这两个数组保证不会重合，打消编译器的顾虑：

void func(float *__restrict a, float *__restrict b) {
    for (int i = 0; i < 1024; i++) {
        a[i] = b[i] + 1;
    }
}

这样它就只生成了一个 SIMD 版，而不需要在运行时判断是不是有重合了：

这样的话，因为不需要分支它效率就能提升，代码更为合理。

6.2 循环中的矢量化：OpenMP 强制矢量化

#pragma 的作用是设定编译器的状态或者是指示编译器完成一些特定的动作 。#pragma 指令对每个编译器给出了一个方法，在保持与 C 和 C++ 语言完全兼容的情况下，给出主机或操作系统专有的特征。依据定义，编译指示是机器或操作系统专有的，且对于每个编译器都是不同的。

除了可以用 GCC 特有的 __restrict 让编译器放心做 SIMD，也可以使用比较跨平台的 OpenMP：

void func(float *a, float *b) {
#pragma omp simd
    for (int i = 0; i < 1024; i++) {
        a[i] = b[i] + 1;
    }
}

OpenMP 是高性能计算的一个框架，它可以通过 -fopenmp 来打开：

gcc -fopenmp -O3 -fomit-frame-pointer -fverbose-asm -S main.cpp -o /tmp/main.S

使用起来只需要用 #pragma omp simd 语句就行了。它会告诉你：下面代码有可能出现数据依赖，但是不要担心。这样编译器就知道这里都显式的说 SIMD 了，那应该能默认这两个指针不会打架，从而它也是只会生成一份代码：

6.3 循环中的矢量化：编译器提示忽略指针别名

除了可以用上面的两种方法外，如果只用 GCC，还可以用 ivdep 这个 pragma，它是 ignore vector dependency 的简称，即忽略矢量依赖关系，表示忽略下方 for 循环内可能的指针别名现象。就是说在下面的代码中，它可以忽略 $b$ 可能依赖于 $a$ 的情况，这个是 GCC 特有的：

void func(float *a, float *b) {
#pragma GCC ivdep
    for (int i = 0; i < 1024; i++) {
        a[i] = b[i] + 1;
    }
}

其他编译器这个 pragma 可能不一样，需要自己去查，这里只是拿 GCC 举例。

6.4 循环中的 if 语句：挪到外面来

还有下面这种循环的优化。这里作者的用意很明显，在 is_mul 为真时执行 a*=b，否则执行 a+=b：

void func(float *__restrict a, float *__restrict b, bool is_mul) {
    for (int i = 0; i < 1024; i++) {
        if (is_mul) {
            a[i] = a[i] * b[i];
        } else {
            a[i] = a[i] + b[i];
        }
    }
}

然而有 if 分支的循环体是难以 SIMD 矢量化的。所以说编译器会把这 if 语句挪到外面，具体是过程如下：

可以看到 corl %eax, %eax 内的 eax，是参数 is_mul 的寄存器，它会先判断这个是 true 还是 false。如果为 true 就跳到 L2，可以看到生成 addps 的一个循环，反之跳到 L3，可以看到 mulps，也就是乘法。

它会去提前判断，也就是它相当于进行了这么一个优化：

void func(float *__restrict a, float *__restrict b, bool is_mul) {
 	if (is_mul) {
    	for (int i = 0; i < 1024; i++) {
            a[i] = a[i] * b[i];
       }
    } else {
  		for (int i = 0; i < 1024; i++) {
            a[i] = a[i] + b[i];
        }
    }
}

也就是编译器看到 is_mul 是一个常量，把 if 挪到了 for 外面，相当于生成了两个版本：一个乘法、一个加法。这样就可以里面没有 if 语句，从而可以矢量化了。这样写起来比较麻烦，所以有时候也像前面那样写。只要编译器知道 is_mul 是循环中不会改变的量，它就能够这样优化。

6.5 循环中的不变量：挪到外面来

除了上面可以挪到外面的 if，还有可以挪到外面的表达式。比如优化之前是长下面这样的：

void func(float *__restrict a, float *__restrict b, float dt) {
    for (int i = 0; i < 1024; i++) {
        a[i]=a[i]+b[i]*(dt*dt);
    }
}

优化之后它发现，dt*dt 在循环中是不会改变的。那如果先把 dt*dt 计算好，之后直接乘就可以了，就不要再重新计算 $1024$ 遍了，也就是这样写：

void func(float *__restrict a, float *__restrict b, float dt) {
    float dt2 = dt * dt;
    for (int i = 0; i < 1024; i++) {
        a[i]=a[i]+b[i]*dt2;
    }
}

所以再汇编中可以看到，它先是进行了一个 dt*dt 的运算 mulss %xmm0, %xmm0，后面就直接应用这个乘好的 %mm 了，就可以节省 $1024$ 次乘法运算：

这种思想就是，把热的代码尽量移到冷的代码里。

6.5.1 挪到外面来：优化失败

如果把上面代码的括号去掉会怎样？

void func(float *__restrict a, float *__restrict b, float dt) {
    for (int i = 0; i < 1024; i++) {
        a[i]=a[i]+b[i]*dt*dt;
    }
}

只要去掉 (dt*dt) 的括号就会优化失败：因为乘法是左结合的，就相当于 (b[i]*dt)*dt。编译器识别不到不变量，从而优化失败。由于浮点运算不满足结合律，因为浮点不能精确表示实数。

因此，要么帮编译器打上括号帮助它识别，要么手动提取不变量到循环体外。

6.6 调用不在另一个文件的函数：SIMD 优化失败

还有一个会导致 SIMD 优化失败的例子，就是调用外部函数。前面提到过，外部函数是优化失败的罪魁祸首，比如下面代码有一个 for 循环，本来这有这个是可以优化成功的，而在代码内调用了一个在其他文件的函数声明，编译器没办法矢量化，这里并不知道 other 是不是会改变 $a$，因为编译器并不知道 other 里干了啥哪怕 other 在定义它的文件里是个空函数，它也不敢优化。

void other();

float func(float *a) {
    float ret = 0;
    for (int i = 0; i < 1024; i++) {
        ret += a[i];
        other();
    }
    return ret;
}

可以看到汇编内全是 ss，说明它优化失败了，都是标量的。如果全是 ps 的话，才说明优化成功。

6.6.1 解决方案：放在同一个文件里

所以说，尽可能把 other 放在同一个文件里定义，这样编译器就能够看到。比如下面的代码：

void other() {
}

float func(float *a) {
    float ret = 0;
    for (int i = 0; i < 1024; i++) {
        ret += a[i];
        other();
    }
    return ret;
}

这里讲 other 放到了和 func 同一个 .cpp 文件里，这样编译器看得到 other 的函数体，就可以内联化该代码。(再次提醒，是不需要 inline 的)

所以说在热的 for 循环里，尽量不要调用外部函数，把它们移到同一个文件里，或者放在同文件声明为 static 函数。

6.7 循环中的下标

6.7.1 随机访问

像下面这种随机访问也会影响 SIMD，使优化失败：

void func(float *a, int *b) {
    for (int i = 0; i < 1024; i++) {
        a[b[i]] += 1;
    }
}

这是因为 b[i] 每次的值是不确定的，它可能是跳跃的在访问，这时候编译器就没有相应的指令能够生成优化。(理论上这种没什么好的优化方式)

6.7.2 循环中的下标：跳跃访问

下面这种情况也比较恶劣，就是 i*2，也就是说先访问 $0$、再访问 $2$、$4$…

void func(float *a) {
    for (int i = 0; i < 1024; i++) {
        a[i * 2] += 1;
    }
}

就如同刚刚说的，这种跳跃的访问是不利于 SIMD 的，但是编译器就非常自然了，它用了一大堆 shufps 指令，还是成功的对 $a$ 进行了一定的 SIMD 优化。(矢量化部分成功，但是非常艰难)

6.7.3 连续访问

所以最好的情况是什么呢？也就是直接 a[i]，也就是顺序访问。这种情况下，编译器能够只用一个指令进行读写。

void func(float *a) {
    for (int i = 0; i < 1024; i++) {
        a[i] += 1;
    }
}

也就是说，不管是编译器还是 CPU，都喜欢顺序的连续访问。

6.12 编译器指令：循环展开

还有一个优化方式是循环展开。像下面这种循环体非常简单：

void func(float *a) {
    for (int i = 0; i < 1024; i++) {
        a[i] = 1;
    }
}

生成的代码一行是 movups，接着的一行是 addq，然后是 cmpq。就是说实际执行计算的指令只有一个，而进行比较和 i++ 的指令却用了两个，还有一个跳转指令。这样就导致大部分时间都花在跳转和 + 上了，而不是在实际的计算上。

对于 GCC 编译器，可以用 #pragma GCC unroll 4，表示把循环体展开为四个。添加后的代码为：

void func(float *a) {
#pragma GCC unroll 4
    for (int i = 0; i < 1024; i++) {
        a[i] = 1;
    }
}

比如下面是把四个循环体放在一个循环内，然后 i++ 改成 i+=4。注意，这里不建议手动这样写，会妨碍编译器的 SIMD 矢量化：

void func(float *a) {
    for (int i = 0; i < 1024; i += 4) {
        a[i + 0] = 1;
        a[i + 1] = 1;
        a[i + 2] = 1;
        a[i + 3] = 1;
    }
}

这样的话它就有 $4$ 个 movups 和 $3$ 个其他指令，从而实际执行计算的时间就更多了，从而就可以避免反复比较的 overhead。(也就是本来是 $1$ 个 movups + $3$ 个其他指令，优化后变成 $4$ 个 movups + $3$ 个其他指令，其他指令使用的次数变少了)

注意，指令 #pragma GCC unroll 4 是 GCC 特定的，笔者目前没有找到什么通用的替代。

$4$ 这个数字可以改，这里代表的是把 $4$ 个循环体变成了一个，可以根据实际的情况决定要用多少。可以根据实际的情况决定要用多少，就是小的循环器一般 unroll 一下是划算的，但最好不要 unroll 大的循环体，否则会造成指令缓存的压力反而变慢。

7. 结构体

这章是重点。

7.1 字节对齐

7.1.1 两个 float：对齐到 8 字节

比如我们有一个 MyVec，它里面有 $x,y$ 两个值。

struct MyVec {
    float x;
    float y;
};

MyVec a[1024];

void func() {
    for (int i = 0; i < 1024; i++) {
        a[i].x *= a[i].y;
    }
}

func 内会使用到 $x$ 和 $y$ 两个成员，可以看到非常复杂的代码：

里面有 shufps 什么的，但最终还是用到了 mulps，也就是它的确矢量化成功了。

7.1.2 三个 float：对齐到 12 字节

但是如果我们突然想把这个改成三维的了，就凭空加一个 $z$，然后我们的 func 甚至没有任何变化：

struct MyVec {
    float x;
    float y;
    float z;
};

MyVec a[1024];

void func() {
    for (int i = 0; i < 1024; i++) {
        a[i].x *= a[i].y;
    }
}

汇编出来的就全部变成了 ss，也就是矢量化失败了，这是什么原因呢？

这是因为刚才之所以能够矢量化，是因为它把两个 MyVec 给读出来，然后读到一个 SIMD 的 $128$ 位寄存器里了。而在这里有 $3$ 个 float 的话，它如果要读到一个 SIMD 里，$128$ 位就是 xyzx，剩余的一个 $y$ 是没有办法读出来的。所以说这时候编译器就尴尬了，它没办法生成一个 SIMD 的读，所以说要怎么办呢？

7.1.3 添加一个辅助对齐的变量：对齐到 16 字节

如果想让它再变回高效，甚至可以像下面这样直接加一个 char padding[4]，这是一个没有任何用处的空变量：

struct MyVec {
    float x;
    float y;
    float z;
    char padding[4];
};

MyVec a[1024];

void func() {
    for (int i = 0; i < 1024; i++) {
        a[i].x *= a[i].y;
    }
}

可以看到这时它又能够优化成功了。可以看到这里生成的 ps 的代码，是进行矢量化的：

可以发现，如果你的 struct 是 $2$ 的整数幂大小，这时候是有利于 SIMD 优化的；而如果像刚才那样是三个一组，这时候它就没有办法对齐，从而导致 SIMD 优化失败。

结论：计算机喜欢 $2$ 的整数幂，2, 4, 8, 16, 32, 64, 128…结构体大小若不是 2 的整数幂，往往会导致 SIMD 优化失败。

7.2 C++11 新语法：alignas

除了上面那种直接加一个从来用不到的变量的方法外，C++11 还新增了 alignas，也就是对齐到 $16$ 字节。这样 MyVec 的起始地址，它会对齐到 $16$，而且它的大小也会变成 $16$ 的整数倍。(这里的 $16$ 代表的是字节而不是位)

struct alignas(16) MyVec {
    float x;
    float y;
    float z;
};

MyVec a[1024];

void func() {
    for (int i = 0; i < 1024; i++) {
        a[i].x *= a[i].y;
    }
}

可以看到这时候生成的汇编是一样的，但是又不需要手动加一个 char padding：

那是不是所有结构体打上 alignas(16) 我的程序就会变快？
错了，有可能不仅不变快，反而还变慢！SIMD 和缓存行对齐只是性能优化的一个点，又不是全部。还要考虑结构体变大会导致内存带宽的占用，对缓存的占用等一系列连锁反应，比如导致它更容易变成 memory bound，从而有可能还会变慢，总之，要根据实际情况选择优化方案。实际测了以后如果的确变快了那才去加。

7.5 结构体的内存布局：AOS 与 SOA

• AOS(Array of Struct)：单个对象的属性紧挨着存-----xyzxyzxyzxyz
• SOA(Struct of Array)：属性分离存储在多个数组-----xxxxyyyyzzzz
• AOS 必须对齐到 $2$ 的幂才高效，SOA 就不需要；
• AOS 符合直觉，不一定要存储在数组这种线性结构，而 SOA 可能无法保证多个数组大小一致；
• SOA 不符合直觉，但通常是更高效的。

刚才说的结构是不是 x、y、z、padding？如下图所示，就是先 xyz，然后图上的黑框处空一个 padding，再 xyz，然后空一个 padding。这种结构类型叫 Array of Struct，它指的是 x、y、z 单个对象的属性紧紧凑在一起。

与之相对的就是 SOA，它会把 $1$ 个对象拆成 $3$ 个数组，就是一个对象的 $x$ 在一个数组里，一个对象的 $y$ 又在另一个数组里，这样的好处就是它不需要去考虑 padding 了。因为这样，如果只访问了 $x$ 和 $y$ 而 $z$ 根本就没有用到，那这时候就相当于只用到了 $x$、$y$；而如果用 AOS 的话，$z$ 和 padding 都相当于没有用的那个属性。所以如果刚才只用 $x$ 和 $y$，那 CPU 看到就是在不停地跳着两格读，从而效率是最低的；而如果用 SOA，只访问 $x$、$y$，没用 $z$，那 CPU 看到的是在访问两个数组，从而它能够用 SIMD 优化。所以说 SOA 虽然它好像不符合我们面向对象的思想，但它其实对 CPU 更友好。

7.5.1 AOS：紧凑存储多个属性

来看看怎么做的。首先是刚才说的比较烂的 AOS，这符合一般面向对象编程(OOP)的习惯，但常常不利于性能：

struct MyVec {
    float x;
    float y;
    float z;
};

MyVec a[1024];

void func() {
    for (int i = 0; i < 1024; i++) {
        a[i].x *= a[i].y;
    }
}

我们用上面代码生成的就是一个不能 SIMD 的代码：

7.5.2 SOA：分离存储多个属性

那要怎么变成 SOA 呢？可以将上面代码写成这样：

struct MyVec {
    float x[1024];
    float y[1024];
    float z[1024];
};

MyVec a;

void func() {
    for (int i = 0; i < 1024; i++) {
        a.x[i] *= a.y[i];
    }
}

之前是 a[1024]，现在 a[1024] 不要了，把它移动到 MyVec 里，就是 $x$ 有 $1024$ 个，$y$ 有 $1024$。访问的时候本来是 a[i].x，现在变成 a.x[i]，就是括号换了下位置。它们分离开来的存储，可以看到下面就用到了 ps，从而矢量化是成功的：

这样不符合 面向对象编程(OOP) 的习惯，但常常有利于性能。又称为 面向数据编程(DOP, date-oriented programming)，这种在高性能计算中比较常见。

7.5.3 AOSOA：中间方案

但是可以发现，SOA 有时不利于优化，所以说推出了一种中间解决方案，就 $4$ 个对象打包乘以 SOA，再用一个 $n/4$ 大小的数组打包成一个 AOS，这样既有 AOS 的直观，又有 SOA 的高效：

struct MyVec {
    float x[4];
    float y[4];
    float z[4];
};

MyVec a[1024 / 4];

void func() {
    for (int i = 0; i < 1024 / 4; i++) {
        for (int j = 0; j < 4; j++) {
            a[i].x[j] *= a[i].y[j];
        }
    }
}

这里其实就相当于把 $4$ 个变成一个 m128 了。但是这种其实也不常用，最常用的还是上面的 SOA。

• AOSOA 优点：
  SOA 便于 SIMD 优化；AOS 便于存储在传统容器；AOSOA 两者得兼。
• AOSOA 缺点：
  需要两层 for 循环，不利于随机访问；需要数组大小是 $4$ 的整数倍，不过可以用边界特判法解决。在代码中可以看到，前面都只需要一个 for 循环，而这里却需要两个，就很麻烦了，而且这对随机访问也非常不友好。

7.9 测试一下加速了多少倍？

这里做个总结测试，比如刚才的程序，有 $x$、$y$、$z$ $3$ 个属性：

struct Point {
    float x;
    float y;
    float z;
};

Point ps[N];

void compute() {
    for (int i = 0; i < N; i++) {
        ps[i].x = ps[i].x + ps[i].y + ps[i].z;
    }
}

第一步将它优化成 SOA，也就是把 N 移上去。再加上每个编译器特有的 unroll 语句，如 GCC 是 #progma GCC unroll 32，MSC 也就是微软的比较霸气，直接就是 #progma unroll 32：

struct Point {
    float x[N];
    float y[N];
    float z[N];
};

Point ps;

void compute() {
    for (int i = 0; i < N; i++) {
        ps.x[i] = ps.x[i] + ps.y[i] + ps.z[i];
    }
}

下面是测试结果：

• $811620ns$：未经优化的耗时 ；
• $989485ns$：用了 $16$ 字节对齐以后的效率，它耗时反而更长了；
• $176316ns$：如果把原来代码一点不变，只是加上一个 parallel for 的效果；
• $218063ns$：使用 SOA 优化的耗时；
• $187422ns$：SOA+grpgrma omp simd 忽略指针别名后的耗时；
• $186657ns$：使用 size_t 作为下标以后的耗时；
• $158938ns$：用 size_t 以后再去 unroll，可以看到它甚至比前面直接加一个 parallel 还要快；
• $59488ns$：如果优化以后再上一个并行，这样就比谁都快了；
• $541613ns$：使用 AOSOA 的优化，它虽然快一点，但还是不如 SOA；

下面是我电脑上的测试效果，差的有点多，首先是对齐后的效率是会更快的，再者 AOSOA 甚至比 SOA 还快。

图中可以看到，一条绿线是无脑并行的，而粉线是我们深度优化，可以看到再深度优化之后，再低数据量的时候甚至超过了无脑并行，因为并行需要创建很多线程，这也需要时间，所以对于小数据甚至超过了它。

这里的结论是：SOA 是针对这个案例最高效的数据排布格式。当然不一定所有的案例都可以用 SOA 就能优化，需要看情况。

8. STL 容器

前面都是基于 C 语言的指针，有些老了。C++ 喜欢 RAII 思想，那就得用 STL 的 vector 容器。

8.1 std::vector：也有指针别名问题

首先是一个经典案例，就是第一章说到的指针别名问题，它不能保证下面代码内的 $c$ 和 $b$ 是同一个 vector，从而它得生成两遍：

#include <vector>

void func(std::vector<int> &a,
          std::vector<int> &b,
          std::vector<int> &c) {
    c[0] = a[0];
    c[0] = b[0];
}

汇编结果如下：

8.1.1 __restrict：能否用于 std::vector？

那我想把这个引用声明 __restrict 行不行：

#include <vector>

void func(std::vector<int> &__restrict a,
          std::vector<int> &__restrict b,
          std::vector<int> &__restrict c) {
    c[0] = a[0];
    c[0] = b[0];
}

可以从下面汇编看到，没有任何效果。因为这里只是把 vector 本身做了一个 __restrict ，而这个 vector 里的指针，它没有上 __restrict，这时候编译器还是不知道 c 和 a 它是不是同一地址的：

8.1.2 解决方案：pragma omp simd 或 pragma GCC ivdep

所以结论就是：要么用 #pragma omp simd 让它忽视可能存在的指针别名，或者用 #pragma gcc ivdep 也可以：
#pragma omp simd：

void func(std::vector<int> &a,
          std::vector<int> &b) {
#pragma omp simd
    for (int i = 0; i < 1024; i++) {
        a[i] = b[i] + 1;
    }
}

#pragma gcc ivdep：

void func(std::vector<int> &a,
          std::vector<int> &b) {
#pragma GCC ivdep
    for (int i = 0; i < 1024; i++) {
        a[i] = b[i] + 1;
    }
}

C/C++ 缺点：指针自由度过高，允许多个 immutable reference 指向同一个对象。

所以我们在适当的时候要加上一些提示来帮助编译器优化。

8.4 std::vector：也能实现 SOA！

先前不是说 SOA 比较好嘛，但是如果使用了 vector 该怎么办？可以像下面这样把 vector 移上去：

• 优化前(AOS)：

  struct MyVec {
      float x;
      float y;
      float z;
  };
  std::vector<MyVec> a;
  void func() {
      for (std::size_t i = 0; i < a.size(); i++) {
          a[i].x *= a[i].y;
      }
  }
  

• 优化后(SOA)：

  struct MyVec {
  	std::vector<float> x;
  	std::vector<float> y;
  	std::vector<float> z;
  };
  MyVec a;
  void func() {
  	for (std::size_t i = 0; i < a.x.size(); i++) {
      	a.x[i] *= a.y[i];
  	}
  }
  

  不过需要保证，$x$、$y$、$z$ 当中一个被 push 的时候，$3$ 个都被 push。否则的话下面的 for 循环会出错。

9. 数学运算

9.1 数学优化：除法变乘法

下面示例接受 $a$ 参数，然后变成 a/2：

float func(float a) {
    return a / 2;
}

返回 a/2，生成的汇编如下所示：

这里 mulss 是一个乘，.long 1056964608 是浮点数 $0.5$ 的意思，也就是说它把一个除法优化成了更快的乘法。

9.2 编译器放弃的优化：分离公共除数

虽然乘法更快，但这个示例中它却没有把除法提取成一个倒数出来：

void func(float *a, float b) {
    for (int i = 0; i < 1024; i++) {
        a[i] /= b;
    }
}

为什么它不能优化呢？因为它害怕 $b$ 是 $0$ 的情况，b 是 $0$ 的话公式 a[i]/=b 会变成 inf 或者说它会报一个错。

9.2.1 解决方案1：手动优化

这时候编译器就不敢优化，所以我们可以手动进行优化，就把 1/b 预先计算出来， 然后在前面公式内除法改成乘法：

void func(float *a, float b) {
    float inv_b = 1 / b;
    for (int i = 0; i < 1024; i++) {
        a[i] *= inv_b;
    }
}

这样就能够把 $1024$ 次除法变成 $1024$ 次乘法加一次除法。这样其实是更快的。

也就是说，乘法比除法更快。提前计算好 b 的倒数避免重复求除法。

9.2.2 解决方案2：-ffast-math

既然前面放弃优化了，这里还有一个解决方案，就是加选项 gcc -ffast-math -O3：

void func(float *a, float b) {
    for (int i = 0; i < 1024; i++) {
        a[i] /= b;
    }
}

-ffast-math 选项让 GCC 更大胆地尝试浮点运算的优化，有时能带来 $2$ 倍左右的提升。作为代价，它对 NaN 和无穷大的处理，可能会和 IEEE 标准规定的不一致。如果能保证程序中永远不会出现 NaN 和无穷大，那么可以放心打开 -ffast-math。（测试发现有时候加了 -ffast-math 是快，但是算出来的值差很多）

9.3 数学函数请加 std:: 前缀！

还有一个问题在于，数学函数 sqrt 只接受 double 的，如果使用 float 去调用 sqrt，它也会返回一个 double，而且计算精度也是基于 double 的，计算速度会慢一些。而如果想要调用 float 的开根号，其实应该用 sqrtf。所以说最好的解决方案应该是 C++ 的 sqrt 函数，它重载了double 和 float 两个版本，就能够自动根据拿什么调用它就自动返回什么类型：

#include <cmath>
#include <type_traits>

void func(float *a) {
    for (int i = 0; i < 1024; i++) {
        auto x = sqrt(a[i]);      // wrong and bad
        static_assert(std::is_same_v<decltype(x), double>);
        auto y = sqrtf(a[i]);     // correct but bad
        static_assert(std::is_same_v<decltype(y), float>);
        auto z = std::sqrt(a[i]); // correct and good
        static_assert(std::is_same_v<decltype(z), float>);
    }
}

也就是推荐不要再用 C 语言全局的 sqrt，多用用 std 里面的。

更坑的是 abs，如果使用全局命名空间的 abs，会发现 abs(1.4f) = 1，这是因为 abs 只接受 int，所以 $1.4$ 被隐式的转换成了 int。要用 fabs 才能调用 double 的；fabsf 调用 float 的。而 std::abs 它重载了 int、 float、double ，就不会出现上面奇怪的现象了。

9.4 -ffast-math 的另一优点

-ffast-math 还有一个优点，比如下面是对 a 的每一个数求根号：

void func(float *a) {
    for (int i = 0; i < 1024; i++) {
        a[i] = std::sqrt(a[i]);
    }
}

9.4.1 开启前：sqrt 矢量化失败

开启前可以看到，这里是没有矢量化的：

为什么这里没有矢量化呢？因为它生怕 a[i] 是一个负数，那么计算到 a[i] 就会出错，然后停止。但是如果矢量化了，那就变成 a[i] 出错时它已经写入了，这里不能保证写入的顺序 ，所以它就不敢优化。

9.4.2 开启后：sqrt 矢量化成功

在开启 -ffast-math 以后，就告诉它不要怕，我能保证 a 绝对不会出错，它就能优化成下面这个：

9.5 嵌套循环：直接累加，有指针别名问题

下面的嵌套循环在计算卷积或者矩阵乘法都会用到：

void func(float *a, float *b, float *c) {
    for (int i = 0; i < 1024; i++) {
        for (int j = 0; j < 1024; j++) {
            c[i] += a[i] * b[j];
        }
    }
}

但直接这样写是有问题的，可以看到这里 $a$ 和 $c$ 是有可能重合的，甚至 $b$ 和 $c$ 是重合的。编译器担心 $c$ 和 $a$ 可能会指向同一个地址，而连续判断三个指针是否有重合又过于复杂，所以放弃了矢量化。

从这里的汇编可以看到，它与前面的章节不同，前面章节还生成了两个版本，而这里只生成了一个版本，即没有生成 SIMD 版本。这就是因为当编译器碰到问题复杂化时，它就直接放弃优化了。所以要么像前面的方式添加 __restrict，要么使用下面的方式：

9.5.1 解决方案1：先读到局部变量，累加完毕后，再写入

可以先把 $c$ 读到 tmp，然后对 tmp 进行追加，读完之后再写回到 $c$。

void func(float *a, float *b, float *c) {
    for (int i = 0; i < 1024; i++) {
        float tmp = c[i];
        for (int j = 0; j < 1024; j++) {
            tmp += a[i] * b[j];
        }
        c[i] = tmp;
    }
}

这样对寄存器更加友好，这时候编译器认为不会存在指针别名的问题，矢量化成功：

9.5.2 解决方案2：先累加到初始为 0 的局部变量，再累加到 c

另一种更好的方法是将 tmp 初始化为 $0$。这与上面方法是有区别的：上面一种加法是在 c[i] 已经有值的情况下去加，这与如果 c[i] 是一个比较大的值，比如 $1000$，而 $a$ 是一个比较小的值，比如 $0.01$。浮点数的构造有一个指数有一个底数，而指数相差太大的话，它们两个之间相加的精度就会受损。就比如 $1000+0.000001$，可能还是 $1000$，就导致加出来的结果可能不对，所以我们最好是把 tmp 初始化为 $0$，然后这样再去 c[i] += tmp 的话就会更精确。因为这样的话 tmp 就足够大，从而不会出现指数位差太远，导致精度误差的问题。

void func(float *a, float *b, float *c) {
    for (int i = 0; i < 1024; i++) {
        float tmp = 0;
        for (int j = 0; j < 1024; j++) {
            tmp += a[i] * b[j];
        }
        c[i] += tmp;
    }
}

这里是精度上的提升，而不是性能上的提升。

10. 优化手法总结

1. 函数尽量写在同一个文件内；
2. 避免在 for 循环内调用外部函数；
3. 非 const 指针加上 __restrict 修饰；
4. 试着用 SOA 取代 AOS；
5. 对齐到 16 或 64 字节($16$ 是 SIMD 宽度，$64$ 是缓存行宽度)；
6. 简单的代码，不要复杂化(把代码写的简单点让编译器看懂；太复杂的看不懂就没法优化)；
7. 试试看 #pragma omp simd；
8. 循环中不变的常量挪到外面来；
9. 对小循环体用 #pragma unroll；
10. -ffast-math 和 -march=native(对于 GCC 编译器可以加上这两个，让它快 $2$ 倍甚至 $4$ 倍都有可能)。

在 CMake 中开启的几个开关

• CMake 中开启 -O3
  默认的 set(CMAKE_BUILD_TYPE Release) 是 Debug，它会生成一些调试，但是它的优化就会被关掉。而 Release 模式，它就相当于 -O3 了

• CMake 中开启 -fopenmp
  因为 fopenmp 是只针对 GCC 的开启方式，微软的可能不一样，所以可以使用这种通用的方式来开启：

  find_package(OpenMP REQUIRED)
  target_link_libraries(testbench PUBLIC OpenMP::OpenMP_CXX)
  

  它是做成了一个库的形式。

• CMake 中开启 -ffast-math 和 -march=native

  target_compile_options(testbench PUBLIC -ffast-math -march=native)
  

  这个好像没办法变成一个跨平台的，只能针对 GCC 开启这两个。