C++中使用SIMD的方法

写在前面

        最近在看《游戏引擎架构》这本书，第四章介绍了利用simd来加速向量运算，感到十分有趣，在此记录下simd的使用方法。

1.SIMD是什么

        SIMD全称为Single Instruction Multiple Data，即单指令多数据流。x86指令集中提供了一些SIMD指令，使用SIMD指令可以将多个操作数打包到一个专用的寄存器中进行运算。这样一来可以方便地进行小到向量内积，大到矩阵运算音视频解码等操作，同时获得一定的性能提升。
        另外值得一提的是x86指令集中的SIMD指令集自发布以来并不是一成不变的，从最早的MMX（1996年）到SSE，再到后来的AVX、AVX2等，在使用时要加以区分尽量避免混用。

2.示例代码

2.1 主要数据结构

        前文中也提到了，simd指令是对多个打包的操作数进行运算，所以C++中提供了相应的数据结构来存储这些打包的数据，用浮点数来举例说，一个单精度浮点数占4个字节也就是32位，我们要同时运算四个浮点数的话就需要将其封装为长为16字节（128位）的数据包。这样的数据结构在程序中以__m128的形式来表示。除此之外还有很多封装例如

__m64				// 64位紧缩整数（MMX）
__m128d				// 128位紧缩双精度（SSE2）
__m128i				// 128位紧缩整数（SSE2）
__m256				// 256位紧缩单精度（AVX） 
__m256i				// 256位紧缩整数（AVX）
__m256d				// 256位紧缩双精度（AVX）

等等。

2.2 字节对齐

        在封装的时候需要保证原始数据是字节对齐的，否则在程序运行过程中会造成segfault的结果，需要注意的是如果是在windows环境下编程，由于MSVC有其自己的方言存在，所以可以通过一下方式进行对齐(以__m128为例)：
__declspec(align(16)) float Arr[] = {1.0f, 2.0f, 3.0f, 4.0f};
在Unix/Linux环境下则需要通过以下的方式：
__attribute__((aligned(16))) float Arr[] = {1.0f, 2.0f, 3.0f, 4.0f};
当然不对齐也是有解决方法的，具体后面会提到。

2.3 示例：向量内积运算

        想要使用各种simd指令可以通过内联汇编的形式使用，也可以通过编译器提供的内部函数来使用，后者需要引入xmmintrin.h头文件，SSE、SSE2、AVX等不同的内部函数分布在不同的头文件中，具体请移步参考资料[3]查阅。以下是参考代码。

#include <iostream>
#include <xmmintrin.h>

__attribute__((aligned(16))) float A[] = {1.0f, 2.0f, 3.0f, 4.0f};
__attribute__((aligned(16))) float B[] = {3.0f, 4.0f, -1.0f, 2.0f};
__attribute__((aligned(16))) float C[] = {0.0f, 0.0f, 0.0f, 0.0f};

int main()
{   
    //load to reg
    __m128 a = _mm_load_ps(&A[0]);
    __m128 b = _mm_load_ps(&B[0]);
    //calculate
    __m128 c = _mm_mul_ps(a, b);
    //store result to C array
    _mm_store_ps(&C[0], c);

    std::cout << "result:" << C[0] << " " << C[1] << " " << C[2] << " " << C[3] << std::endl;
    return 0;
}

结果如下：

这里就要提到不用对齐的方法了，观察代码我们可以发现在进行运算之前原始数据的封装是利用_mm_load_ps函数实现的，这里编译器提供了另一个名叫_mm_loadu_ps的方法，对于没有对齐的数组，直接将原方法替换掉就可以了。

2.4 示例：矢量与矩阵相乘

        有线性代数基础的朋友都知道，矩阵乘法运算口诀是是左行乘右列，但是按照这种传统方法来做，会产生许多列向量形式的的中间结果，只有将一个列向量内的所有元素相加才能得到最终结果的其中一维，这显然是十分麻烦的。所以这里的技巧是使用M的行向量相乘，即让向量V中的元素如Vx与矩阵M中的第一行乘，得到的结果再与下一次运算(Vy与M的第二行乘)的结果相加，如此循环得到最终结果。所以代码中利用_mm_shuffle_ps定义了一些方法，来进行辅助运算，该方法可以调整SSE寄存器中分量顺序。

#include <iostream>
#include <xmmintrin.h>

// 以下代码可以简略的理解为将其中一个分量暂时分配到另外三个分量的位置
// 例如将[x, y, z, w] 通过_mm_replicate_x_ps(v)变为 [x, x, x, x]
#define SHUFFLE_PARAM(x, y, z, w) \
    ((x) | ((y) << 2) | ((z) << 4) | ((w) << 6))
#define _mm_replicate_x_ps(v) \
    _mm_shuffle_ps((v), (v), SHUFFLE_PARAM(0, 0, 0, 0))
#define _mm_replicate_y_ps(v) \
    _mm_shuffle_ps((v), (v), SHUFFLE_PARAM(1, 1, 1, 1))
#define _mm_replicate_z_ps(v) \
    _mm_shuffle_ps((v), (v), SHUFFLE_PARAM(2, 2, 2, 2))
#define _mm_replicate_w_ps(v) \
    _mm_shuffle_ps((v), (v), SHUFFLE_PARAM(3, 3, 3, 3))

// 将乘法和加法整合到一起
#define _mm_madd_ps(a, b, c) \
    _mm_add_ps(_mm_mul_ps((a), (b)), (c))

// 运算向量和矩阵相乘的子方法
__m128 mulVecMatrix(__m128 v, 
    __m128 Mrow1, __m128 Mrow2, __m128 Mrow3, __m128 Mrow4)
{
    __m128 result;
    result = _mm_mul_ps(_mm_replicate_x_ps(v), Mrow1);
    result = _mm_madd_ps(_mm_replicate_y_ps(v), Mrow2, result);
    result = _mm_madd_ps(_mm_replicate_z_ps(v), Mrow3, result);
    result = _mm_madd_ps(_mm_replicate_w_ps(v), Mrow4, result);

    return result;
}

// 定义测试数据
__attribute__((aligned(16))) float row1[] = {1.0f, 2.0f, 3.0f, -4.0f};
__attribute__((aligned(16))) float row2[] = {3.0f, 4.0f, -1.0f, 2.0f};
__attribute__((aligned(16))) float row3[] = {-3.0f, 2.0f, 4.0f, 1.0f};
__attribute__((aligned(16))) float row4[] = {4.0f, -1.0f, 2.0f, 3.0f};
__attribute__((aligned(16))) float vec[] = {7.0f, 8.0f, 6.0f, 5.0f};
__attribute__((aligned(16))) float C[] = {0.0f, 0.0f, 0.0f, 0.0f};

int main()
{   
    //load to reg
    __m128 Mrow1 = _mm_load_ps(&row1[0]);
    __m128 Mrow2 = _mm_load_ps(&row2[0]);
    __m128 Mrow3 = _mm_load_ps(&row3[0]);
    __m128 Mrow4 = _mm_load_ps(&row4[0]);
    __m128 v = _mm_load_ps(&vec[0]);
    //calculate
    __m128 c = mulVecMatrix(v, Mrow1, Mrow2, Mrow3, Mrow4);
    //store result to C array
    _mm_store_ps(&C[0], c);

    std::cout << "result:" << C[0] << " " << C[1] << " " << C[2] << " " << C[3] << std::endl;
    return 0;
}

运行结果：

3. 总结

        一次很有意思的体验，对之后的高性能计算程序的开发会起到一定的指导作用。

参考资料

1. SIMD(MMX/SSE/AVX)变量命名规范心得 http://www.cnblogs.com/zyl910/archive/2012/04/23/simd_var_name.html
2. SIMD指令集发展历程表(MMX、SSE、AVX等) https://www.cnblogs.com/zyl910/archive/2012/02/26/x86_simd_table.html
3. Intrinsics头文件与SIMD指令集、Visual Studio版本对应表 https://www.cnblogs.com/zyl910/archive/2012/04/26/md00.html
4. https://stackoverflow.com/questions/10366670/how-to-compile-simd-code-with-gcc