如何使用SSE来优化算子核心
一：准备
1.背景资料
SSE 全称Intel® Streaming SIMD Extensions (Intel® SSE),本质是Intel公司封装汇编语句提供的底层操作指令函数集。同样属于底层操作指令集的还有著名的Intel® AVX(Advanced Vector Extensions), 及 Intel® AVX2(Intel® Advanced Vector Extensions 2)。基于同样原理封装的还有Arm 对应Arm Intrinsic，MIPS中对应MIPS Intrinsic。
SSE的版本包含：SSE/SSE2/SSE3/SSE4.1/SSE4.2。下文中在描述CPU特性上统称为SSE系列指令集。在描述具体使用指令函数中的CPUID Flags，才会具体区分SSE不同版本。
自从MSVC不再支持x64的汇编指令后（虽然可以强制使用，但不推荐不安全）。SSE，AVX等成为MSVC 支持的最佳底层优化方法。
本文将从SSE的使用出发，以ncnn实现为例，展示如何使用SSE优化深度学习中算子。
优化算子工作需要三方面的准备事项：
- 测试正确的原生代码
- 快速测试验证环境
- 基准统计程序
2.确认硬件是否支持SSE
在开始SSE优化之前，首先请确保您硬件支持SSE指令集，对于大多数Intel CPU都支持SSE指令集。但在各种系统环境下，查看方式不同。我们有：
1.windows环境
windows环境下推荐简单使用GPU-Z来检测当前处理器是否支持SSE扩展。在GPU-Z官网下载后，运行，在“处理器”-“CPU支持的特性”项目下，若包含SSE系列指令集，即当前CPU支持SSE。
2.Linux环境和类Unix环境
Linux环境和类Unix环境下，使用查看cpuinfo文件来确认CPU特性；
cat /proc/cpuinfo
...
flags: *** sse sse2 *** #在cpu flags中即可检查是否支持sse扩展
2.macOS环境：
macOS本质是像Unix环境，所以同样使用sysctl 来查看CPU特性.(注意Mac的 M1 M2系列芯片是arm架构，不支持SSE)
sysctl machdep.cpu # 结果同Linux环境
二：编写原生代码
使用SSE来优化算法的过程本质就是代码重构的一种情况。代码重构的首要条件是完成完备的代码行为测试集合。所以，这部分将从测试代码的编写开始。
其次优化过程的目标是调优某些性能指标的过程。所以第二部分将讨论性能指标的选定和优先级；
1.编写测试代码
在大多数情况下，看到这篇文章的人肯定是比笔者更会写算法，所以我在这里只谈一些编写测试的注意事项（这里的测试指验证你算法满足你的要求所编写的代码行为，跟其他人无关）。
编写测试代码主要注意事项：
- 思考如何构造基础数据结构才能满足算法行为的输入要求。举例来说，如果你准备为ncnn贡献算法，请阅读ncnn中关于Mat结构的函数。最好编写相关测试来验证该数据结构满足你的需要。（笔者的建议是可以先从简单结构来验证，比如需要做一个支持f32任意大小的矩阵加法算子，可以先从支持固定矩阵int8类型的加法开始编写测试代码）。
- 保持结果的正确性。首先考虑，你所编写的原生代码行为上，是否满足你所需要的结果（不论这个结果是手算的，numpy算的或者pytorch算的）。其次要考虑，结果在内存上结构如何排布。以ncnn为例，思考你的结果该如何放入到一个Mat中，Mat的size该如何设定。（在后续SSE优化中，我们将多次以原生代码结果作为target结果，验证每次优化后的正确性，原生代码能够稳定输出正确结果非常重要）
- 不用过早考虑算法的完备性，应该随着每次测试结果的正确来迭代重构算法和测试代码。二者同样重要。如果能够自动化测试，请尽量让一个简单的脚本执行来完成所有你慢吞吞的命令行。
2.考虑性能指标
性能指标的主要作用是随着每次优化的迭代，告诉我们所采取的措施在什么方面取得效果，是正面优化还是负面优化。
性能指标很多，包括吞吐量，还有类似计算稳定度，时间延迟，视频方面还有fps 等等。无法确认有效的性能指标也是大多数优化算法的困难点之一。
随着简单粗暴地叠晶体管数量来解决电脑运行问题，性能指标似乎变得越来越不重要。这是一种错误观念，如果在单核上编写非常烂的代码，增加N个核心只是把烂代码重复N次而已。另一方面，性能指标有着客观性，在开发板上和集群设备上运行同样的算法，性能指标的优先级也不一样。但是，我认为应当满足最基础性能指标有这两个：
- 吞吐量：算法在单位时间内执行的次数，用Gflops表示，该值越大越好（也可以认为执行同样算法所平均占用的时间，时间越短越好）；
- 性能衰退：即随着数据规模的增加，Gflops在不同数据规模下的波动情况。更低的离散程度意味着吞吐量保持在一定范围内不发生变化。
其余有效的性能指标应当由业务环境和任务需求决定。负责技术基础设施建设的算法工程师，一方面应该理解业务所需求中的最高优先级，另一方面也应该追求做到更好。
以SSE优化算法，本质上是重构的迭代过程。不用在初期就考虑如何达到最大性能指标，而是应该考虑每次迭代中带来一定量的性能优化。
三：理解SSE
SSE主要由SSE基础数据类型 及 针对性的SSE操作函数构成。前文提到，SSE是针对汇编语句的封装，所以本身不具备错误检查和错误处理（错误检查和错误处理一般由编译器完成）。使用不当的话，诸如segmentation fault之内指针指向不存在的内存错误非常常见。我在此处建议：使用SSE优化之前，确保理解代码指针位置和移动原理，原生代码已经完成测试，输出结果正确。
1.SSE数据类型
SSE数据类型形如：
__m<bit><type> //__m适用代表申请mm寄存器
// bit 代表数据类型的字节长度，在SSE中为128 或 64
// 默认type为单精度浮点（f32），其余为int 或double
// 另外要注意所有SSE的类型除__m128和__m64外，随着版本更新有不同的类型，建议根据需要且确定硬件性能后选择合适的类型
// 举例如下：
__m128 //4xf32 含有4个单精度浮点数；SSE
__m64 //4xf32 含有2个单精度浮点数；SSE
__m128i //8个int类型（8x16) ；SSE3
__m128d //2个double类型(2x64)
2.SSE内联函数结构
SSE内联函数在线查询：Intel® Intrinsics Guide 在此 单个指令的结构如下：
- Synopsis ：摘要。描述指令的接口定义，需要引入的头文件，对应的指令，CPU必须支持的标志；
- Description：描述该指令的行为；
- Operation：逻辑层面描述指令行为；
- Performance：在不同架构中所需要的延迟和执行所需要的时钟周期数（CPI）。
值得指出的的是此处默认使用小端存储，即左边为高位，右边为低位。
相似的内联函数有很多，在使用时候一定要注意Operation中的逻辑满足您的要求。
另外，在ncnn中，ncnn已经将部分SSE内联函数以NCNN内联的方式封装。在为NCNN添加SSE优化的算法的过程中，请首先考虑搜索“NCNNINLINE”宏封装的SSE函数。
四：样例
1.一个简单的样例：4x4矩阵乘法
矩阵乘法方面，已经有很多出色的成果。值得一读的比如how-to-optimize-gemm，及 以Arm Intrinsic优化矩阵乘法。我建议感兴趣同学参考和学习这两份项目，来探究如何从0到1优化一份算法；
矩阵乘法原理很简单：
假设有A，B两个矩阵，如下：
$$
A_{[4][4]} =
\begin{bmatrix}
a_0 & a_4 & a_8 & a_{12} \
a_1 & a_5 & a_9 & a_{13} \
a_2 & a_6 & a_{10} & a_{14} \
a_3 & a_7 & a_{11} & a_{15}
\end{bmatrix}
~~
B_{[4][4]} =
\begin{bmatrix}
b_0 & b_4 & b_8 & b_{12} \
b_1 & b_5 & b_9 & b_{13} \
b_2 & b_6 & b_{10} & b_{14} \
b_3 & b_7 & b_{11} & b_{15}
\end{bmatrix}
~~
C_{[4][4]} =
\begin{bmatrix}
c_0 & c_4 & c_8 & c_{12} \
c_1 & c_5 & c_9 & c_{13} \
c_2 & c_6 & c_{10} & c_{14} \
c_3 & c_7 & c_{11} & c_{15}
\end{bmatrix}
$$
对于C 矩阵的第一列，我们有：
$$
c_0 = a_0b_0 + a_4b_1 + a_8b_2 + a_{12}b_3 \
c_1 = a_1b_0 + a_5b_1 + a_9b_2 + a_{13}b_3 \
c_2 = a_2b_0 + a_6b_1 + a_{10}b_2 + a_{14}b_3 \
c_3 = a_3b_0 + a_7b_1 + a_{11}b_2 + a_{15}b_3
$$
1.编写测试代码和基准测量程序
在该样例中，测试代码很容易编写出来，我们只需要初始化4x4的二维数组，并返回指针即可。此时，可以不考虑泛用性，初始化为固定值即可。
// <代码片段>
...
float A[16] = {0.0f}; // 此处已经将输入和输出的矩阵默认展开成im2col 后的单行（inch = 1） 宽度为h*w = 16的矩阵
float B[16] = {0.0f};
float C[16] = {0.0f};
matrix_init_rand(A, 4, 4); // 随机初始化A数组
matrix_init_rand(B, 4, 4); // 随机初始化B数组
编写验证正确性的测试代码。
// <代码片段>
...
float T[16] = {...}; // Target即为预测的C的结果数组，可用numpy或者纸笔计算
...
float error = 0.0001;
bool CheckAuc(T, C, error);
// 注意：float在计算机中不能完全表示，只能使用绝对误差的判别方法。gtest等测试框架的EXCEPT宏无法处理1.234e5这样结构的float数的对比。
同样，编写计算耗时的基准测量代码，此处使用1000次操作所占的平均时间来作为基准。
// <代码片段>
...
const int loop = 1000;
clock_gettime_(CLOCK_REALTIME, &time_start);
for(init i = 0; i < loop; i++)
{
matirx_mult_native(C, A, B);
}
clock_gettime_(CLOCK_REALTIME, &time_end);
clocks_c = (time_end.tv_sec - time_start.tv_sec) * 1000000 + (time_end.tv_sec - time_start.tv_sec) /1000;
2.编写原生代码
编写原生代码，使得正确性测试能够通过。
// <代码片段>
static void matirx_mult_native(float *C, float *A, float *B)
{
for(int i_idx = 0; i_idx < 4; i_idx++)
{
for(int j_idx = 0; j_idx < 4; j_idx++)
{
for(int k_idx = 0; k_idx < 4; k_idx++)
{
C[4*j_idx + i_idx] += A[4*k_idx + i_idx] * B[4*j_idx + k_idx];
}
}
}
}
3.优化原生代码
注意到上述代码中，先取c0 - c3 的计算作为样例考虑： $$ c_0 = a_0b_0 + a_4b_1 + a_8b_2 + a_{12}b_3 \ c_1 = a_1b_0 + a_5b_1 + a_9b_2 + a_{13}b_3 \ c_2 = a_2b_0 + a_6b_1 + a_{10}b_2 + a_{14}b_3 \ c_3 = a_3b_0 + a_7b_1 + a_{11}b_2 + a_{15}b_3 $$
1.装载寄存器
-
考虑竖排a0-a1-a2-a3 为4个f32 数据，又因为SSE可以申请mm寄存器，单次保存128bit，那么不妨把a0-a4保存在寄存器中，
-
对于b0-b3 则是，单次读取一个值，能够重复用4次，不妨考虑b0 重复4次，排满单个128bit的mm寄存器；
-
同理把c0-c3也放入寄存器，从列方向上考虑，取名为_c0
_m128 _a0 = _mm_load_ps(a_ptr); //a0 -a1 -a2 -a3 _m128 _a1 = _mm_load_ps(a_ptr + 4); //a4 -a5 -a6 -a7 _m128 _a2 = _mm_load_ps(a_ptr + 8); //a8 -a9 -a10-a11 _m128 _a3 = _mm_load_ps(a_ptr + 12); //a12-a13-a14-a15 _m128 _b0 = _mm_load_ps1(b_ptr); // b0 - b0 - b0 - b0 _m128 _b1 = _mm_load_ps1(b_ptr + 4); // b1 - b1 - b1 - b1 _m128 _b2 = _mm_load_ps1(b_ptr + 8); // b2 - b2 - b2 - b2 _m128 _b3 = _mm_load_ps1(b_ptr + 12); // b3 - b3 - b3 - b3
2.编写第一列的计算结果
对于_a0 -_a3 数据与_b0 数据相乘 ，有：
// 保存结果新建一个_c0 作为临时变量
_m128 _c0 = _mm_set_ps1(0.0f);
_c0 = _mm_mul_ps(_a0, _b0);
_c0 = _mm_add_ps(_mm_mul_ps(_a1, _b1),_c0);
_c0 = _mm_add_ps(_mm_mul_ps(_a2, _b2),_c0);
_c0 = _mm_add_ps(_mm_mul_ps(_a3, _b3),_c0);
// 把 _sum0存会以c指针开头的内存中，完美！
_mm_store_ps(c_ptr, _c0);
3.将单列输出扩展到所有列：
我们针对剩下的c中的c1 列也做相同的操作： 对于C1 列 有： $$ c_4 = a_0b_4 + a_4b_5 + a_8b_6 + a_{12}b_7 \ c_5 = a_1b_4 + a_5b_5 + a_9b_6 + a_{13}b_7 \ c_6 = a_2b_4 + a_6b_5 + a_{10}b_6 + a_{14}b_7 \ c_7 = a_3b_4 + a_7b_5 + a_{11}b_6 + a_{15}b_7 $$
// a 系列不变 b系列指针+1
_m128 _b4 = _mm_load_ps1(b_ptr + 1); // b4 - b4 - b4 - b4
_m128 _b5 = _mm_load_ps1(b_ptr + 4 + 1); // b5 - b5 - b5 - b5
_m128 _b6 = _mm_load_ps1(b_ptr + 8 + 1); // b6 - b6 - b6 - b6
_m128 _b7 = _mm_load_ps1(b_ptr + 12+ 1); // b7 - b7 - b7 - b7
// 保存结果新建一个_c0 作为临时变量
_m128 _c1 = _mm_set_ps1(0.0f);
_c1 = _mm_mul_ps(_a0, _b4);
_c1 = _mm_add_ps(_mm_mul_ps(_a1, _b5),_c1);
_c1 = _mm_add_ps(_mm_mul_ps(_a2, _b6),_c1);
_c1 = _mm_add_ps(_mm_mul_ps(_a3, _b7),_c1);
// 把 _sum0存会以c指针开头的内存中，完美！
_mm_store_ps(c_ptr, _c1);
此时我们发现，对于C1列的操作与C0列及其相似，只不过是b_ptr的指针发生移动，不妨将其放到同一个循环中，有：
// a 系列不变
_m128 _a0 = _mm_load_ps(a_ptr); //a0 -a1 -a2 -a3
_m128 _a1 = _mm_load_ps(a_ptr + 4); //a4 -a5 -a6 -a7
_m128 _a2 = _mm_load_ps(a_ptr + 8); //a8 -a9 -a10-a11
_m128 _a3 = _mm_load_ps(a_ptr + 12); //a12-a13-a14-a15
for(int i = 0; i < 4; i++)
{
_m128 _b0 = _mm_load_ps1(b_ptr); // b0 - b0 - b0 - b0
_m128 _b1 = _mm_load_ps1(b_ptr + 4); // b1 - b1 - b1 - b1
_m128 _b2 = _mm_load_ps1(b_ptr + 8); // b2 - b2 - b2 - b2
_m128 _b3 = _mm_load_ps1(b_ptr + 12); // b3 - b3 - b3 - b3
_m128 _ci = _mm_set_ps1(0.0f);
_ci = _mm_mul_ps(_a0, _b0);
_ci = _mm_add_ps(_mm_mul_ps(_a1, _b1),_ci);
_ci = _mm_add_ps(_mm_mul_ps(_a2, _b2),_ci);
_ci = _mm_add_ps(_mm_mul_ps(_a3, _b3),_ci);
_mm_store_ps(c_ptr, _ci);
b_ptr += 1; // 移动b_ptr
c_ptr += 4; // 移动保存内存的c_ptr
}
2.NCNN中以SSE优化算子的注意事项
1.线程与openmp
以上计算Benchmark 和 SSE优化的方法大多集中在单个核心中，但是在实际使用ncnn中，ncnn使用Option opt 中提供的num_threads 给openmp赋值，以实现多线程并行化，同时运行在多个核心上。
#pragma omp parallel for num_threads(opt.num_threads)
在优化成SSE代码的初期，可以考虑锁定为单线程，或者直接不用考虑线程的影响，仅对单核以SSE优化，保证单核的结果正确后，再加上opt的多线程进行结果测试。
2.展开循环
在实际ncnn实现的原生代码的算法中，循环是非常常见的。针对以SSE优化这类循环，遵循非常简单的原则：循环中，迭代器等于零时刻，整个输出的结果也是正确的。
那么，在我们使用SSE优化过程中，不妨以迭代器等于零的时刻，函数计算结果作为此时目标结果。在此基础上再利用SSE优化代码。与目标结果核对正确以后，再进一步去考虑迭代器等于1的情况（重复这个过程直到迭代器达到最大值）。在迭代器的每个元素下，SSE优化出的代码都与结果相等，那么我们可以说，该次优化是正确性，且完全覆盖了需执行代码。（一般来说不用考虑到最大值，根据数学归纳法，n有效，n+1有效，那么n的序列都是有效的）
五：总结
本文描述SSE的使用及以4x4矩阵乘法的样例来优化SSE代码。
值得注意的是，SSE只是128bit数据宽度的指令集，但是也可以用来模拟256bit 和 512bit数据宽度，来实现以pack4拼接成pack8，甚至pack16的做法，只不过在输出结果管理上更加繁琐而已。感兴趣的同学可以尝试一下。
六：引用
- SSE指令扩展快查；
- 浮点性能基准计算-浮点峰值那些事儿
- 硬件性能基准测试计算样例：M1芯片搞数据科学好使吗？5种基准测试给你答案
- 讨论矩阵乘法如何优化的系列论文：how-to-optimized-gemm
- 讨论以Arm Intrinsic 优化gemm的系列文章：OpenBLAS gemm从零入门
Source: docs/developer-guide/how-to-write-a-sse-optimized-op-kernel.zh.md