SSE学习-一个小栗子

指令集优化在高性能计算中至关重要，所以用 C/C++ 到后面感觉难免要用上指令集。虽然指令集学习和使用都不太容易，但想一想能够在不增加 CPU 占用的情况下提升数倍计算速度，确实挺诱人的。

具体指令集是啥就不多介绍了，几个名词：

SIMD: 单指令多数据，即指令集加速技术
SSE: Streaming SIMD Extensions, 使用128位寄存器的指令集（Intel）
AVX: Advanced Vector Extensions, 使用256位寄存器的指令集（Intel）
Neon: ARM 上使用 128 位寄存器的指令集

基本上 00 年后的机器都是支持指令集的，但256位指令集（AVX）出现的较晚，而且 ARM 上貌似还只有128位的 Neon，考虑兼容性尤其是移植问题的话，还是用128位指令集为主吧。博主还是个初学者，本文只是介绍 SSE 入门的一个小例子。

配置 C++ 工程

想在 C++ 代码中使用 SSE，只需要两个简单的步骤：
1. 包含相关头文件
具体头文件如下(stackoverflow)：

<mmintrin.h>  MMX
<xmmintrin.h> SSE
<emmintrin.h> SSE2
<pmmintrin.h> SSE3
<tmmintrin.h> SSSE3
<smmintrin.h> SSE4.1
<nmmintrin.h> SSE4.2
<ammintrin.h> SSE4A
<wmmintrin.h> AES
<immintrin.h> AVX, AVX2, FMA

目前只需要包含<immintrin.h>就相当于包含了上述全部头文件。
2. 添加编译选项
加入-msse4.2(或者-msse2等，如果不支持高版本 SSE)。比如用 g++ 编译：

g++ -msse4.2 -std=c++17 -O3 main.cpp -o main

或者使用 CMake 时，在 CMakeLists.txt 中加入：

target_compile_options(MyProject PRIVATE -msse4.2)

C++ 代码

SSE 开方计算测试：

#include <cmath>
#include <immintrin.h>
#include <iostream>

#include <TestFuncs/TicToc.hpp> // 时间测试功能

void sqrtNormal(float* a, int N)
{
    for (int i = 0; i < N; ++i)
        a[i] = sqrt(a[i]);
}

void sqrtSSE(float* a, int N)
{
    for (__m128 *ptr = reinterpret_cast<__m128*>(a),
                *end = reinterpret_cast<__m128*>(a + N);
         ptr < end; ptr++) {
        __m128 val = _mm_sqrt_ps(*ptr);
        *ptr = val;
    }
}

int main(int argc, char* argv[])
{
    // prepare data
    constexpr int N = 4000000;
    float* buf1 = static_cast<float*>(std::aligned_alloc(16, N * 4));
    float* buf2 = static_cast<float*>(std::aligned_alloc(16, N * 4));
    for (int i = 0; i < N; ++i) {
        float val = rand() % 10000;
        buf1[i] = buf2[i] = val;
    }

    bxg::TicToc::start("normal");
    sqrtNormal(buf1, N);
    auto t1 = bxg::TicToc::showTime("normal");

    bxg::TicToc::start("sse");
    sqrtSSE(buf2, N);
    auto t2 = bxg::TicToc::showTime("sse");
    std::cout << "speed up: " << t1 / t2 << std::endl;

    // deallocate memory
    std::free(buf1);
    std::free(buf2);
    return 0;
}

首先要了解，内存数据的对齐（数据地址为 16 Byte 的倍数）对于 SIMD 操作很重要。大部分 SIMD 指令都要求数据必须对齐，否则运行时报错。即使有些指令不强制要求对齐，但不对齐的数据很可能会导致性能下降。

内存对齐可以像上面例子那样使用std::aligned_alloc，也可以用alignas说明符，或者编译器相关的属性标识符，如 MSVC 中__declspec(align(16))，或者 gcc 中__attribute__ ((aligned (16))。需注意的是，即使没有使用上述方法显式对齐，new 出来的内存空间也很有可能是按16字节对齐的，所以就算你忘记对齐，程序可能也不报错，但这有很多不可控因素，使用 SIMD 时还是应当手动对齐数据！

上面例子中，使用 SSE 的就是函数sqrtSSE()，具体过程为：

1. 将内存数据加载到128位寄存器
2. 在128位寄存器中作运算：_mm_sqrt_ps同时计算四个单精度浮点数的开方
3. 将计算结果取出到内存

相应的代码：

// float* a 存放了已对齐的内存数据
__m128 *ptr = reinterpret_cast<__m128*>(a); // 从内存加载4个 float 到128位寄存器
__m128 val = _mm_sqrt_ps(*ptr); // 4个数同时计算，结果保存在128位寄存器中
*ptr = val; // 将128位寄存器的数据存入内存

__m128 *ptr = reinterpret_cast<__m128*>(a)等同于__m128 *ptr = (__m128*)(a)，通过指针的转化自动进行内存与寄存器的数据交换，更容易理解的使用方法是：

__m128 v1 = _mm_load_ps(a); // 从内存加载4个 float 到128位寄存器
__m128 v2 = _mm_sqrt_ps(v1); // 4个数同时计算，结果保存在128位寄存器中
_mm_store_ps(a, v2); // 将128位寄存器的数据存入内存

这两种方法效果一样。

通过反汇编理解 SSE 指令

C++ 查看反汇编指令的方法（stackoverflow）：

• gcc/g++ 加入编译选项 -g -Wa,-alh，可将汇编指令输出到终端（与 C++ 代码一起）
• gcc/g++ 加入编译选项 -S，可将编译指令输出到xxx.s文件（不含 C++ 代码）
• objdump -d xxx.out（UNIX 或 cygwin）
• dumpbin /DISASM xxx.exe（Windows）

更方便的方法是借助 IDE和调试器(gdb/cdb) 查看反汇编指令：

• QtCreator: 调试命中断点后，Debug菜单 -> Operate by Instruction
• VisualStudio: 调试命中断点后，右键菜单 -> 转到反汇编

这里需要注意的是，我们想分析的是 release 模式生成的代码，但 release 编译优化会带来一些调试上的困难（主要是因为变量优化、函数自动内联等）。如果一个函数在 release 编译过程中被自动内联，那么我们打在里面的断点就不会被命中。为此，我们需要禁止编译器内联，具体做法是给函数加入 noinline 属性标识符：

__declspec(noinline) void foo(){} // MSVC
void __attribute__ ((noinline)) foo(){} // gnu
[[gnu::noinline]] void foo(){} // gnu, with c++11 style

至此，VisualStudio 已经可以正常调试 release 模式下的函数了，但 QtCreator 还必须将编译模式切换成 Release with Debug Information 才行。

例子中 SSE 部分代码的反编译指令：

        14 [1]	{
0x5555555568e0                  f3 0f 1e fa  endbr64
        15 [1]	    for (__m128 *ptr = reinterpret_cast<__m128*>(a),
0x5555555568e4  <+    4>        48 63 f6     movslq %esi,%rsi
0x5555555568e7  <+    7>        48 8d 04 b7  lea    (%rdi,%rsi,4),%rax
        17 [1]	         ptr < end; ptr++) {
0x5555555568eb  <+   11>        48 39 c7     cmp    %rax,%rdi
0x5555555568ee  <+   14>        73 10        jae    0x555555556900 <sqrtSSE(float*, int)+32>
        210 [1]	  return (__m128) __builtin_ia32_sqrtps ((__v4sf)__A);
0x5555555568f0  <+   16>        0f 51 07     sqrtps (%rdi),%xmm0
        19 [1]	        *ptr = val;
0x5555555568f3  <+   19>        48 83 c7 10  add    $0x10,%rdi
0x5555555568f7  <+   23>        0f 29 47 f0  movaps %xmm0,-0x10(%rdi)
        17 [1]	         ptr < end; ptr++) {
0x5555555568fb  <+   27>        48 39 f8     cmp    %rdi,%rax
0x5555555568fe  <+   30>        77 f0        ja     0x5555555568f0 <sqrtSSE(float*, int)+16>
0x555555556900  <+   32>        c3           retq

其中，sqrtps和movaps都是 SIMD 指令，一次操作4个 float 数。但这儿我还不太理解的是，为什么开头是movslq而非movaps（从内存加载数据到寄存器），希望有大佬能解答一下。

加速比

上述例子的测试结果，debug 模式下：

timer [normal]: 0.0095135 s
timer [sse]: 0.00239473 s
speed up: 3.97269

release 模式下：

timer [normal]: 0.00336943 s
timer [sse]: 0.00147091 s
speed up: 2.29071

debug 模式确实能够达到接近4的加速比，但为什么 release 模式下加速比只有2呢？我们分别看一下 debug 和 release 模式下 sqrtNormal() 函数的反汇编指令：

// debug
        10 [1]	        a[i] = sqrt(a[i]);
0x55555555666b  <+   34>        8b 45 fc                 mov    -0x4(%rbp),%eax
0x55555555666e  <+   37>        48 98                    cltq
0x555555556670  <+   39>        48 8d 14 85 00 00 00 00  lea    0x0(,%rax,4),%rdx
0x555555556678  <+   47>        48 8b 45 e8              mov    -0x18(%rbp),%rax
0x55555555667c  <+   51>        48 01 d0                 add    %rdx,%rax
0x55555555667f  <+   54>        f3 0f 10 00              movss  (%rax),%xmm0
0x555555556683  <+   58>        f3 0f 5a c0              cvtss2sd %xmm0,%xmm0
0x555555556687  <+   62>        e8 24 fd ff ff           callq  0x5555555563b0 <sqrt@plt>
0x55555555668c  <+   67>        8b 45 fc                 mov    -0x4(%rbp),%eax
0x55555555668f  <+   70>        48 98                    cltq
0x555555556691  <+   72>        48 8d 14 85 00 00 00 00  lea    0x0(,%rax,4),%rdx
0x555555556699  <+   80>        48 8b 45 e8              mov    -0x18(%rbp),%rax
0x55555555669d  <+   84>        48 01 d0                 add    %rdx,%rax
0x5555555566a0  <+   87>        f2 0f 5a c0              cvtsd2ss %xmm0,%xmm0
0x5555555566a4  <+   91>        f3 0f 11 00              movss  %xmm0,(%rax)

// release
        10 [1]	        a[i] = sqrt(a[i]);
0x555555556890  <+   32>        f3 0f 10 07           movss  (%rdi),%xmm0
0x555555556894  <+   36>        0f 2e d0              ucomiss %xmm0,%xmm2
0x555555556897  <+   39>        0f 28 c8              movaps %xmm0,%xmm1
0x55555555689a  <+   42>        f3 0f 51 c9           sqrtss %xmm1,%xmm1
0x55555555689e  <+   46>        77 19                 ja     0x5555555568b9 <sqrtNormal(float*, int)+73>
0x5555555568a0  <+   48>        f3 0f 11 0f           movss  %xmm1,(%rdi)

可以看到，release 模式下，编译器自动进行了大量的优化，甚至还自动使用了指令集(movaps)！

参考资料

Intel Intrinsics Guide: SSE 指令查询，非常有用！
C++ SSE Optimization - Lesson: 油管视频教程
SIMD vector instructions: How to Write Fast Numerical Code 课程 ppt