SIMD简介
SIMD简介 - 知乎本篇文章包含的内容有SIMD指令集简介以及简短的practice环节。 1.SIMD的历史与分类SIMD( Single Instruction Multiple Data)即单指令流多数据流，是一种采用一个控制器来控制多个处理器，同时对一组数据（又称“数…
https://zhuanlan.zhihu.com/p/55327037
本篇文章包含的内容有SIMD指令集简介以及简短的practice环节。
1.SIMD的历史与分类
SIMD(Single Instruction Multiple Data)即单指令流多数据流，是一种采用一个控制器来控制多个处理器，同时对一组数据（又称“数据向量”）中的每一个分别执行相同的操作从而实现空间上的并行性的技术。简单来说就是一个指令能够同时处理多个数据。
标量运算与SIMD运算对比
如上图所示，使用标量运算一次只能对一对数据执行乘法操作，而采用SIMD乘法指令，则一次可以对四对数据同时执行乘法操作。
SIMD于20世纪70年代首次引用于ILLIAC IV大规模并行计算机上。而大规模应用到消费级计算机则是在20实际90年代末。
1996年Intel推出了X86的MMX(MultiMedia eXtension)指令集扩展，MMX定义了8个寄存器，称为MM0到MM7，以及对这些寄存器进行操作的指令。每个寄存器为64位宽，可用于以“压缩”格式保存64位整数或多个较小整数，然后可以将单个指令一次应用于两个32位整数，四个16位整数或8个8位整数。
intel在1999年又推出了全面覆盖MMX的SSE(Streaming SIMD Extensions， 流式SIMD扩展)指令集，并将其应用到Pentium III系列处理器上，SSE添加了八个新的128位寄存器(XMM0至XMM7)，而后来的X86-64扩展又在原来的基础上添加了8个寄存器(XMM8至XMM15)。SSE支持单个寄存器存储4个32位单精度浮点数，之后的SSE2则支持单个寄存器存储2个64位双精度浮点数，2个64位整数或4个32位整数或8个16位短整形。SSE2之后还有SSE3，SSE4以及AVX，AVX2等扩展指令集。
AVX引入了16个256位寄存器(YMM0至YMM15)，AVX的256位寄存器和SSE的128位寄存器存在着相互重叠的关系(XMM寄存器为YMM寄存器的低位)，所以最好不要混用AVX与SSE指令集，否在会导致transition penalty(过渡处罚)，两种寄存器的关系如下图：
AVX256位寄存器与SSE128位寄存器的关系
AVX与SSE支持的数据类型如下：
AVX与SSE支持的数据类型
不同处理器对于SIMD指令集的支持如下图：
不同处理器对于SIMD指令集的支持
如果想知道CPU的SIMD支持等级可以使用cpuid指令 ，或者直接使用cpuz软件查看。
2.如何使用SIMD
下图给出了使用SIMD的不同方法：
使用SIMD的不同方法
首先是最简单的方法是使用Intel开发的跨平台函数库（IPP，Intel Integrated Performance Primitives ），里面的函数实现都使用了SIMD指令进行优化。
其次是借助于Auto-vectorization(自动矢量化)，借助编译器将标量操作转化为矢量操作。
第三种方法是使用编译器指示符(compiler directive)，如Cilk里的#pragma simd和OpenMP里的#pragma omp simd。如下所示，使用#pragma simd强制循环矢量化：
void add_floats（float * a，float * b，float * c，float * d，float * e，int n）
{int i;
#pragma simdfor（i = 0; i 第四种方法则是使用内置函数(intrinsics)的方式，如下所示，使用SSE _mm_add_ps 内置函数，一次执行8个单精度浮点数的加法：
int main()
{__m128 v0 = _mm_set_ps(1.0f, 2.0f, 3.0f, 4.0f);__m128 v1 = _mm_set_ps(1.0f, 2.0f, 3.0f, 4.0f);__m128 result = _mm_add_ps(v0, v1);
}
最后一种方法则是使用汇编直接操作寄存器，当然直接使用汇编有点太麻烦了，所以本篇文章主要介绍的方法是使用intrinsics的方式使用SIMD指令。
3.SSE/AVX Intrinsics介绍
a.头文件
SSE/AVX指令主要定义于以下一些头文件中：
: SSE, 支持同时对4个32位单精度浮点数的操作。 : SSE 2, 支持同时对2个64位双精度浮点数的操作。 : SSE 3, 支持对SIMD寄存器的水平操作(horizontal operation)，如hadd, hsub等...。 : SSSE 3, 增加了额外的instructions。 : SSE 4.1, 支持点乘以及更多的整形操作。 : SSE 4.2, 增加了额外的instructions。 : AVX, 支持同时操作8个单精度浮点数或4个双精度浮点数。
每一个头文件都包含了之前的所有头文件，所以如果你想要使用SSE4.2以及之前SSE3, SSE2, SSE中的所有函数就只需要包含
b.命名规则
SSE/AVX提供的数据类型和函数的命名规则如下：
- 数据类型通常以_mxxx(T)的方式进行命名，其中xxx代表数据的位数，如SSE提供的__m128为128位，AVX提供的__m256为256位。T为类型，若为单精度浮点型则省略，若为整形则为i，如__m128i，若为双精度浮点型则为d，如__m256d。
- 操作浮点数的内置函数命名方式为：_mm(xxx)_name_PT。 xxx为SIMD寄存器的位数，若为128m则省略，如_mm_addsub_ps，若为_256m则为256，如_mm256_add_ps。 name为函数执行的操作的名字，如加法为_mm_add_ps，减法为_mm_sub_ps。 P代表的是对矢量(packed data vector)还是对标量(scalar)进行操作，如_mm_add_ss是只对最低位的32位浮点数执行加法，而_mm_add_ps则是对4个32位浮点数执行加法操作。 T代表浮点数的类型，若为s则为单精度浮点型，若为d则为双精度浮点，如_mm_add_pd和_mm_add_ps。
- 操作整形的内置函数命名方式为：_mm(xxx)_name_epUY。xxx为SIMD寄存器的位数，若为128位则省略。 name为函数的名字。U为整数的类型，若为无符号类型则为u，否在为i，如_mm_adds_epu16和_mm_adds_epi16。Y为操作的数据类型的位数，如_mm_cvtpd_pi32。
c.instructions的分类
instructions按执行操作类别的不同主要分为以下几类：
1).存取操作(load/store/set)
load系列可以用来从内存中载入数据到SSE/AVX提供的类型中，如：
void test()
{__declspec(align(16)) float p[] = { 1.0f, 2.0f, 3.0f, 4.0f };__m128 v = _mm_load_ps(p);
}
_mm_load_ps可以从16字节对齐的连续内存中加载4个32位单精度浮点数到__m128数据类型中（若不对齐则加载会出错）。
_mm_loadu_ps同_mm_load_ps的作用相同，但不要求提供的内存地址对齐。
_mm_load_ps1是从内存中载入一个32位浮点数，并重复存储到__m128中的4个浮点数中，即：m[0] = p[0], m[1] = p[0], m[2] = p[0], m[3] = p[0]。
_mm_load_ss则是从内存中载入一个32位浮点数，并将其赋值给__m128中的最低位的浮点数，并将高位的3个浮点数设置为0，即：m[0] = p[0], m[1] = 0, m[2] = 0, m[3] = 0。
_mm_loadr_ps位载入4个32位浮点数并将其反向赋值给__m128中的4个浮点数，即：m[0] = p[3], m[1] = p[2], m[2] = p[1], m[3] = p[0]。
除此之外还有_mm_loadh_pd，_mm_loadl_pi等...
store系列可以将SSE/AVX提供的类型中的数据存储到内存中，如：
void test()
{__declspec(align(16)) float p[] = { 1.0f, 2.0f, 3.0f, 4.0f };__m128 v = _mm_load_ps(p);__declspec(align(16)) float a[] = { 1.0f, 2.0f, 3.0f, 4.0f };_mm_store_ps(a, v);
}
_mm_store_ps可以__m128中的数据存储到16字节对齐的内存。
_mm_storeu_ps不要求存储的内存对齐。
_mm_store_ps1则是把__m128中最低位的浮点数存储为4个相同的连续的浮点数，即：p[0] = m[0], p[1] = m[0], p[2] = m[0], p[3] = m[0]。
_mm_store_ss是存储__m128中最低位的位浮点数到内存中。
_mm_storer_ps是按相反顺序存储__m128中的4个浮点数。
set系列可以直接设置SSE/AVX提供的类型中的数据，如：
__m128 v = _mm_set_ps(0.5f, 0.2f, 0.3f, 0.4f);
_mm_set_ps可以将4个32位浮点数按相反顺序赋值给__m128中的4个浮点数，即：_mm_set_ps(a, b, c, d) : m[0] = d, m[1] = c, m[2] = b, m[3] = a。
_mm_set_ps1则是将一个浮点数赋值给__m128中的四个浮点数。
_mm_set_ss是将给定的浮点数设置到__m128中的最低位浮点数中，并将高三位的浮点数设置为0.
_mm_setzero_ps是将__m128中的四个浮点数全部设置为0.
2). 算术运算
SSE/AVX提供的算术运算操作包括：
- _mm_add_ps，_mm_add_ss等加法系列
- _mm_sub_ps，_mm_sub_pd等减法系列
- _mm_mul_ps，_mm_mul_epi32等乘法系列
- _mm_div_ps，_mm_div_ss等除法系列
- _mm_sqrt_pd，_mm_rsqrt_ps等开平方系列
- _mm_rcp_ps，_mm_rcp_ss等求倒数系列
- _mm_dp_pd，_mm_dp_ps计算点乘
此外还有向下取整，向上取整等运算，这里只列出了浮点数支持的算术运算类型，还有一些整形的算术运算类型未列出。
3).比较运算
SSE/AVX提供的比较运算操作包括：
- _mm_max_ps逐分量对比两个数据，并将较大的分量存储到返回类型的对应位置中。
- _mm_min_ps逐分量对比两个数据，并将较小的分量存储到返回类型的对应位置中。
- _mm_cmpeq_ps逐分量对比两个数据是否相等。
- _mm_cmpge_ps逐分量对比一个数据是否大于等于另一个是否相等。
- _mm_cmpgt_ps逐分量对比一个数据是否大于另一个是否相等。
- _mm_cmple_ps逐分量对比一个数据是否小于等于另一个是否相等。
- _mm_cmplt_ps逐分量对比一个数据是否小于另一个是否相等。
- _mm_cmpneq_ps逐分量对比一个数据是否不等于另一个是否相等。
- _mm_cmpnge_ps逐分量对比一个数据是否不大于等于另一个是否相等。
- _mm_cmpngt_ps逐分量对比一个数据是否不大于另一个是否相等。
- _mm_cmpnle_ps逐分量对比一个数据是否不小于等于另一个是否相等。
- _mm_cmpnlt_ps逐分量对比一个数据是否不小于另一个是否相等。
此外还有一些执行单分量对比的比较运算
4).逻辑运算
SSE/AVX提供的逻辑运算操作包括：
- _mm_and_pd对两个数据逐分量and
- _mm_andnot_ps先对第一个数进行not，然后再对两个数据进行逐分量and
- _mm_or_pd对两个数据逐分量or
- _mm_xor_ps对两个数据逐分量xor
5).Swizzle运算
包含_mm_shuffle_ps，_mm_blend_ps, _mm_movelh_ps等。
这里主要介绍以下_mm_shuffle_ps:
void test()
{__m128 a = _mm_set_ps(1, 2, 3, 4);__m128 b = _mm_set_ps(5, 6, 7, 8);__m128 v = _mm_shuffle_ps(a, b, _MM_SHUFFLE(1, 0, 3, 2)); // 2, 1, 8, 7
}
_mm_shuffle_ps读取两个__m128类型的数据a和b，并按照_MM_SHUFFLE提供的索引将返回的__m128类型数据的低两位设置为a中按索引值取得到的对应值，将高两位设置为按索引值从b中取得到的对应值。索引值在0到3之间，分别以相反的顺序对应__m128中的四个浮点数。
SSE/AVX还提供了类型转换等操作，这里就不做介绍了。
PS:补充一个操作指令：_mm_cvtss_f32，可以获取__m128中的最低位浮点数并返回。
4. Practice环节
现在可以使用之前所介绍的instructions来实现一个简单的数学运算库。
首先需要定义一些宏：
#ifndef SIMD_LEVEL#if defined(__AVX__) || defined(__AVX2__)#define SIMD_LEVEL 2 //Support SSE4, AVX, AVX2#include #elif defined(_M_IX86_FP) && (_M_IX86_FP == 2)#define SIMD_LEVEL 1 //Support SSE2#include #else#define SIMD_LEVEL 0 #endif
#endif // !SIMD_LEVEL#ifndef AMVector#if SIMD_LEVEL == 0typedef Vector4 AMVector;typedef const Vector4& CRAMVector;#elsetypedef __m128 AMVector;typedef const AMVector CRAMVector;#endif
#endif // !AMVector#define SHUFFLE4(V, X,Y,Z,W) (_mm_shuffle_ps(V, V, _MM_SHUFFLE(W,Z,Y,X)))
#define SHUFFLE3(V, X,Y,Z) (_mm_shuffle_ps(V, V, _MM_SHUFFLE(3,Z,Y,X)))
#define SHUFFLE2(V, X,Y) (_mm_shuffle_ps(V, V, _MM_SHUFFLE(3,2,Y,X)))#define AM_INLINEF __forceinline
#define AM_CALLCONV __vectorcall
visual studio上可以通过__AVX__，__AVX2__等宏检测是否支持AVX指令集，_M_IX86_FP为2则表示支持SSE2，为1则表示支持SSE，否则不支持SSE。
这里将__m128类型定义为AMVector，并且定义了SHUFFLE4等宏方便对__m128类型进行Swizzle运算。
因为我们定义的函数都比较短，所以有必要把所有函数都定义为_forceinline来减少函数调用开销。x64上的默认函数调用约定为__fastcall，前两个参数通过寄存器传递，其他参数通过堆栈传递，而__vectorcall能使用比__fastcall更多的寄存器传递参数，并且支持__m128矢量类型，在可能的情况下可以通过寄存器返回函数返回值。（关于__vectorcall的更多信息可以查看Introducing ‘Vector Calling Convention’）
首先需要定义加，减，乘，除等算术运算与赋值操作：
AM_INLINEF AMVector AM_CALLCONV operator + (CRAMVector lhs, CRAMVector rhs)
{return _mm_add_ps(lhs, rhs);
}AM_INLINEF AMVector& AM_CALLCONV operator += (AMVector &lhs, CRAMVector rhs)
{lhs = _mm_add_ps(lhs, rhs);return lhs;
}AM_INLINEF AMVector AM_CALLCONV operator - (CRAMVector lhs, CRAMVector rhs)
{return _mm_sub_ps(lhs, rhs);
}AM_INLINEF AMVector& AM_CALLCONV operator -= (AMVector &lhs, CRAMVector rhs)
{lhs = mm_sub_ps(lhs, rhs);return lhs;
}AM_INLINEF AMVector AM_CALLCONV operator * (CRAMVector lhs, CRAMVector rhs)
{return _mm_mul_ps(lhs, rhs);
}AM_INLINEF AMVector AM_CALLCONV operator * (CRAMVector lhs, float rhs)
{return _mm_mul_ps(lhs, _mm_set1_ps(rhs));
}AM_INLINEF AMVector AM_CALLCONV operator * (float lhs, CRAMVector rhs)
{return _mm_mul_ps(_mm_set1_ps(lhs), rhs);
}AM_INLINEF AMVector& AM_CALLCONV operator *= (AMVector& lhs, float rhs)
{lhs = _mm_mul_ps(lhs, _mm_set1_ps(rhs));return lhs;
}AM_INLINEF AMVector& AM_CALLCONV operator *= (AMVector& lhs, CRAMVector rhs)
{lhs = _mm_mul_ps(lhs, rhs);return lhs
}AM_INLINEF AMVector AM_CALLCONV operator / (CRAMVector lhs, CRAMVector rhs)
{return _mm_div_ps(lhs, rhs);
}AM_INLINEF AMVector AM_CALLCONV operator / (CRAMVector lhs, float rhs)
{return _mm_div_ps(lhs, _mm_set1_ps(rhs
}AM_INLINEF AMVector AM_CALLCONV operator / (float lhs, CRAMVector rhs)
{return _mm_div_ps(_mm_set1_ps(lhs), rhs);
}AM_INLINEF AMVector& AM_CALLCONV operator /= (AMVector& lhs, float rhs)
{lhs = _mm_div_ps(lhs, _mm_set1_ps(rhs));return lhs;
}AM_INLINEF AMVector& AM_CALLCONV operator /= (AMVector& lhs, CRAMVector rhs)
{lhs = _mm_div_ps(lhs, rhs);return lh
}
然后可以定义一些Set，get 操作方便我们存储和读取__m128中的值：
AM_INLINEF AMVector AM_CALLCONV am_vector_set(float x, float y, float z, float w)
{return _mm_set_ps(w, z, y, x);
}AM_INLINEF float AM_CALLCONV am_vector_get_y(CRAMVector v)
{return _mm_cvtss_f32(_mm_shuffle_ps(v, v, _MM_SHUFFLE(1, 1, 1, 1)));
}
之后就是必不可少的点乘和叉乘操作了：
AM_INLINEF AMVector AM_CALLCONV am_vector3_dot(CRAMVector lhs, CRAMVector rhs)
{
#if SIMD_LEVEL == 0float dot = lhs.x() * rhs.x() + lhs.y() * rhs.y() + lhs.z() * rhs.z();return Vector4{dot};
#elif SIMD_LEVEL == 1AMVector dot = _mm_mul_ps(lhs, rhs);AMVector temp = SHUFFLE4(dot, 1, 1, 2, 2);dot = _mm_add_ps(dot, temp);temp = SHUFFLE4(dot, 3, 3, 3, 3);dot = _mm_add_ps(dot, temp);return SHUFFLE4(dot, 0, 0, 0, 0);
#elsereturn _mm_dp_ps(lhs, rhs, 0x7f);
#endif // SIMD_LEVEL == 0
}
如果支持SSE4则可以直接采用_mm_dp_ps实现点乘操作，否则若支持SSE则可以将lhs和rhs相乘，然后通过shuffle操作打乱乘法操作的结果，并将对应的的三个分量加起来。
AM_INLINEF AMVector AM_CALLCONV am_vector3_cross(CRAMVector lhs, CRAMVector rhs)
{
#if SIMD_LEVEL == 0return vector4_cross3(lhs, rhs);
#else//http://threadlocalmutex.com/?p=8 Investigating SSE Cross Product Performance//return _mm_sub_ps(_mm_mul_ps(SHUFFLE3(lhs, 1, 2, 0), SHUFFLE3(rhs, 2, 0, 1)), _mm_mul_ps(SHUFFLE3(lhs, 2, 0, 1), SHUFFLE3(rhs, 1, 2, 0)));AMVector result = _mm_sub_ps(_mm_mul_ps(lhs, SHUFFLE3(rhs, 1, 2, 0)), _mm_mul_ps(SHUFFLE3(lhs, 1, 2, 0), rhs));return SHUFFLE3(result, 1, 2, 0);
#endif // SIMD_LEVEL == 0
}
三维向量的叉乘操作定义为：cross(a, b) = a.yzx * b.zxy - a.zxy * b.yzx;
可以表示为：
cross(a, b).x = a.y * b.z - a.z * b.y;cross(a, b).y = a.z * b.x - a.x * b.z;cross(a, b).z = a.x * b.y - a.y * b.x;
可以将上面的式子重新排序如下：
cross(a, b).z = a.x * b.y - a.y * b.x;cross(a, b).x = a.y * b.z - a.z * b.y;cross(a, b).y = a.z * b.x - a.x * b.z;
然后我们就可以得到：
cross(a, b).zxy = a * b.yzx - a.yzx * b;
即：
cross(a, b) = (a * b.yzx - a.yzx * b).yzx;
这样以来叉乘的实现就只需要3个shuffle操作即可完成，比起代码中我原来实现的方法少了一个shuffle指令，效率更高(上面介绍的叉乘操作方法来自于Investigating SSE Cross Product Performance)。
随后可以定义向量的取模和归一化操作：
AM_INLINEF AMVector AM_CALLCONV am_vector3_length_sq(CRAMVector v)
{return am_vector3_dot(v, v);
}AM_INLINEF AMVector AM_CALLCONV am_vector3_length(CRAMVector v)
{
#if SIMD_LEVEL == 0float dot = v.x() * v.x() + v.y() * v.y() + v.z() * v.z();return Vector4{ std::sqrt(dot) };
#elif SIMD_LEVEL == 1AMVector dot = _mm_mul_ps(v, v);AMVector temp = SHUFFLE4(dot, 1, 1, 2, 2);dot = _mm_add_ps(dot, temp);temp = SHUFFLE4(dot, 3, 3, 3, 3);dot = _mm_sqrt_ss(_mm_add_ps(dot, temp));return SHUFFLE4(dot, 0, 0, 0, 0);
#elsereturn _mm_sqrt_ps(_mm_dp_ps(v, v, 0x7f));
#endif // SIMD_LEVEL == 0
}
计算向量的模可以先通过对向量和它本身做点击计算出模的平方，然后再对其开平方求得。
AM_INLINEF AMVector AM_CALLCONV am_vector3_normalize_est(CRAMVector v)
{
#if SIMD_LEVEL == 0float dot = v.x() * v.x() + v.y() * v.y() + v.z() * v.z();return vector4_divide(v, std::sqrt(dot));
#elif SIMD_LEVEL == 1AMVector dot = _mm_mul_ps(v, v);AMVector temp = SHUFFLE4(dot, 1, 1, 2, 2);dot = _mm_add_ps(dot, temp);temp = SHUFFLE4(dot, 3, 3, 3, 3);dot = _mm_rsqrt_ss(_mm_add_ps(dot, temp));return _mm_mul_ps(v, SHUFFLE4(dot, 0, 0, 0, 0));
#elsereturn _mm_mul_ps(v, _mm_rsqrt_ps(_mm_dp_ps(v, v, 0x7f)));
#endif // SIMD_LEVEL == 0
}
对向量进行归一化即 v/v.length，需要求向量模的倒数，可以通过_mm_rsqrt_ps取倒数，但是会有一些误差，也可以通过_mm_div_ps(1, length)取得精确的倒数值。
此外还定义了sqrt, lerp, clamp, saturate等函数。
AM_INLINEF AMVector AM_CALLCONV am_vector_sqrt_est(CRAMVector v)
{return _mm_rsqrt_ps(v);
}AM_INLINEF AMVector AM_CALLCONV am_vector_lerp(float t, CRAMVector lhs, CRAMVector rhs)
{return _mm_add_ps(lhs, am_vector_scale(_mm_sub_ps(rhs, lhs), t));
}AM_INLINEF AMVector AM_CALLCONV am_vector_clamp(CRAMVector v, CRAMVector min, CRAMVector max)
{return _mm_max_ps(_mm_min_ps(v, max), min);
}AM_INLINEF AMVector AM_CALLCONV am_vector_saturate(CRAMVector v)
{return _mm_max_ps(_mm_min_ps(v, am_const_one), am_const_zero);
}
还有一些算术运算如log，exp，pow等函数可以借助标准库std::log, std::exp, std::pow实现
(pow实现整数指数计算比较简单，但是要支持浮点数就有些麻烦了...，exp虽然可以通过泰勒展开式计算，但是要计算到比较精确的值需要计算100项以上，速度比较慢，而且误差也比加大...)
三角函数可以通过泰勒展开式计算，只需要计算5项左右就可以得到比较好的精度了：
AM_INLINEF AMVector AM_CALLCONV am_vector_acos(CRAMVector v)
{AMVector v2 = _mm_mul_ps(v, v);AMVector v3 = _mm_mul_ps(v2, v);AMVector v5 = _mm_mul_ps(v3, v2);AMVector v7 = _mm_mul_ps(v5, v2);AMVector v9 = _mm_mul_ps(v7, v2);AMVector result = _mm_sub_ps(_mm_sub_ps(_mm_sub_ps(_mm_sub_ps(_mm_sub_ps(am_const_half_pi, v), _mm_mul_ps(v3, _mm_set1_ps(0.16666666666667f))),_mm_mul_ps(v5, _mm_set1_ps(0.075f))), _mm_mul_ps(v7, _mm_set1_ps(0.04464285714f))), _mm_mul_ps(v9, _mm_set1_ps(0.030381944444444f)));return result;
}
常见的三角函数泰勒展开式都可以在wiki上找到：
常见三角函数的泰勒展开式
除了这些算术运算外还定义了floor, ceil等函数直接借助于SSE的_mm_floor_ps和_mm_ceil_ps函数进行计算。
然后就是比较运算：
typedef int AM_MASK;AM_INLINEF AM_MASK AM_CALLCONV am_vector_less(CRAMVector lhs, CRAMVector rhs)
{AMVector temp = _mm_cmplt_ps(lhs, rhs);return _mm_movemask_ps(temp);
}AM_INLINEF bool am_mask_all3(AM_MASK m)
{return (m & 7) == 7;
}AM_INLINEF bool am_mask_any3(AM_MASK m)
{return (m & 7) != 0;
}
_mm_cmplt_ps函数返回的是__m128类型的数据，其中存储了逐分量比较的结果，如果条件成立则为-nan，否则为0：
void test()
{__m128 a = am_vector_set(5, 6, 7, 4);__m128 b = am_vector_set(5, 6, 7, 8);__m128 v = _mm_cmpeq_ps(a, b); // -nan, -nan, -nan, 0AM_MASK mask = _mm_movemask_ps(v); // 7bool all = am_mask_all3(mask); //true
}
可以使用_mm_movemask_ps将比较结果转化为int类型的值，此int值从低到高的4位分别对应两个__m128从低到高逐分量比较的结果，所以这个int类型的值的范围在0-15之间。可以通过检查这个int类型的值判断条件是否都成立，或其中一个是否成立。
AM_INLINEF AMVector AM_CALLCONV am_vector_select(CRAMVector a, CRAMVector b, CRAMVector condition)
{
//.....
#elif SIMD_LEVEL == 1return _mm_or_ps(_mm_andnot_ps(condition, a), _mm_and_ps(b, condition));
#elsereturn _mm_blendv_ps(a, b, condition);
#endif // SIMD_LEVEL == 0
}
可以定义一个select函数，根据比较结果，选择a或b中的对应分量作为返回值的对应分量(如果为true则选择b，否则选择a)，如果支持sse则select函数通过位操作实现，若支持sse4则可以通过_mm_blendv_ps指令实现。select函数的效果如下：
void test()
{__m128 a = am_vector_set(1, 9, 3, 10);__m128 b = am_vector_set(5, 6, 7, 8);__m128 v = _mm_cmplt_ps(a, b); // -nan, 0, -nan, 0__m128 result = am_vector_select(a, b, v); // 5, 9, 7, 10
}
以上差不多介绍了所有向量支持的操作。接下来可以实现下简单的矩阵操作。
矩阵的定义如下（这里矩阵采用的是行主序(Row Major)）：
struct alignas(16) SIMDMatrix
{AMVector m_rows[4];SIMDMatrix()= default;SIMDMatrix(const SIMDMatrix &m){m_rows[0] = m.m_rows[0];m_rows[1] = m.m_rows[1];m_rows[2] = m.m_rows[2];m_rows[3] = m.m_rows[3];}SIMDMatrix(float a00, float a01, float a02, float a03,float a10, float a11, float a12, float a13,float a20, float a21, float a22, float a23,float a30, float a31, float a32, float a33){m_rows[0] = am_vector_set(a00, a01, a02, a03);m_rows[1] = am_vector_set(a10, a11, a12, a13);m_rows[2] = am_vector_set(a20, a21, a22, a23);m_rows[3] = am_vector_set(a30, a31, a32, a33);}
//......
}#if SIMD_LEVEL == 0
typedef Matrix4x4 AMMatrix;
typedef const Matrix4x4& CRAMMatrix;
#else
typedef SIMDMatrix AMMatrix;
typedef const SIMDMatrix CRAMMatrix;
#endif
矩阵也定义了一些set函数，和AMVector定义的set函数差不多，这里就不再给出了。
接下来定义矩阵与矩阵的相乘操作，考虑两个矩阵a和b相乘得到c：
矩阵相乘
根据矩阵相乘的定义(a的行乘以b的列)，可以得到c的第一行为：
对上面的式子重新排列后可以得到：
即c的第0行为a的第零行第零列的元素乘以b的第0行元素加上a的第零行第一列的元素乘以b的第2行元素加上a的第零行第二列的元素乘以b的第3行元素最后再加上a的第零行第三列的元素乘以b的第3行元素。c的其他行与第一行类似，代码实现如下：
AM_INLINEF AMMatrix AM_CALLCONV am_matrix_multiply(CRAMMatrix lhs, CRAMMatrix rhs)
{AMMatrix result;//Row0AMVector vec = lhs.m_rows[0];AMVector vec_x = SHUFFLE4(vec, 0, 0, 0, 0);AMVector vec_y = SHUFFLE4(vec, 1, 1, 1, 1);AMVector vec_z = SHUFFLE4(vec, 2, 2, 2, 2);AMVector vec_w = SHUFFLE4(vec, 3, 3, 3, 3);vec_x = _mm_mul_ps(vec_x, rhs.m_rows[0]);vec_y = _mm_mul_ps(vec_y, rhs.m_rows[1]);vec_z = _mm_mul_ps(vec_z, rhs.m_rows[2]);vec_w = _mm_mul_ps(vec_w, rhs.m_rows[3]);vec_x = _mm_add_ps(vec_x, vec_y);vec_z = _mm_add_ps(vec_z, vec_w);result.m_rows[0] = _mm_add_ps(vec_x, vec_z);//Row1vec = lhs.m_rows[1];vec_x = SHUFFLE4(vec, 0, 0, 0, 0);vec_y = SHUFFLE4(vec, 1, 1, 1, 1);vec_z = SHUFFLE4(vec, 2, 2, 2, 2);vec_w = SHUFFLE4(vec, 3, 3, 3, 3);vec_x = _mm_mul_ps(vec_x, rhs.m_rows[0]);vec_y = _mm_mul_ps(vec_y, rhs.m_rows[1]);vec_z = _mm_mul_ps(vec_z, rhs.m_rows[2]);vec_w = _mm_mul_ps(vec_w, rhs.m_rows[3]);vec_x = _mm_add_ps(vec_x, vec_y);vec_z = _mm_add_ps(vec_z, vec_w);result.m_rows[1] = _mm_add_ps(vec_x, vec_z);//Row2vec = lhs.m_rows[2];vec_x = SHUFFLE4(vec, 0, 0, 0, 0);vec_y = SHUFFLE4(vec, 1, 1, 1, 1);vec_z = SHUFFLE4(vec, 2, 2, 2, 2);vec_w = SHUFFLE4(vec, 3, 3, 3, 3);vec_x = _mm_mul_ps(vec_x, rhs.m_rows[0]);vec_y = _mm_mul_ps(vec_y, rhs.m_rows[1]);vec_z = _mm_mul_ps(vec_z, rhs.m_rows[2]);vec_w = _mm_mul_ps(vec_w, rhs.m_rows[3]);vec_x = _mm_add_ps(vec_x, vec_y);vec_z = _mm_add_ps(vec_z, vec_w);result.m_rows[2] = _mm_add_ps(vec_x, vec_z);//Row3vec = lhs.m_rows[3];vec_x = SHUFFLE4(vec, 0, 0, 0, 0);vec_y = SHUFFLE4(vec, 1, 1, 1, 1);vec_z = SHUFFLE4(vec, 2, 2, 2, 2);vec_w = SHUFFLE4(vec, 3, 3, 3, 3);vec_x = _mm_mul_ps(vec_x, rhs.m_rows[0]);vec_y = _mm_mul_ps(vec_y, rhs.m_rows[1]);vec_z = _mm_mul_ps(vec_z, rhs.m_rows[2]);vec_w = _mm_mul_ps(vec_w, rhs.m_rows[3]);vec_x = _mm_add_ps(vec_x, vec_y);vec_z = _mm_add_ps(vec_z, vec_w);result.m_rows[3] = _mm_add_ps(vec_x, vec_z);return result;
}
向量与矩阵相乘采用的是左乘，计算方式与上面计算c的第一行的方式完全相同。
接下来是定义矩阵转置操作，转置操作用到了_mm_unpacklo_ps，_mm_unpackhi_ps，_mm_movelh_ps，_mm_movehl_ps等操作。这里简单介绍下_mm_unpacklo_ps和_mm_movelh_ps。
void test()
{__m128 a = am_vector_set(1, 2, 3, 4);__m128 b = am_vector_set(5, 6, 7, 8);__m128 v = _mm_unpacklo_ps(a, b); // 1, 5, 2, 6//v[0] = a[0]//v[1] = b[0]//v[2] = a[1]//v[3] = b[1]v = _mm_movelh_ps(a, b); // 1, 2, 5, 6//v[0] = a[0]//v[1] = a[1]//v[2] = b[0]//v[3] = b[1]
}
对矩阵A做转置操作可以先对第0行和第1行做_mm_unpacklo_ps操作得到temp1 : {a00, a10, a01, a11}， 再对第2行和第3行做_mm_unpacklo_ps操作得到temp2 : {a20, a30, a21, a31}，再对temp1和temp2做_mm_movelh_ps操作就可以得到转置后矩阵的第一行元素{a00, a10, a20, a30}了。其他行的操作与第一行类似。
转置操作的完整实现如下：
AM_INLINEF AMMatrix AM_CALLCONV am_matrix_transpose(CRAMMatrix m)
{AMVector t0 = _mm_unpacklo_ps(m.m_rows[0], m.m_rows[1]);AMVector t1 = _mm_unpacklo_ps(m.m_rows[2], m.m_rows[3]);AMVector t2 = _mm_unpackhi_ps(m.m_rows[0], m.m_rows[1]);AMVector t3 = _mm_unpackhi_ps(m.m_rows[2], m.m_rows[3]);return AMMatrix{ _mm_movelh_ps(t0, t1), _mm_movehl_ps(t1, t0), _mm_movelh_ps(t2, t3), _mm_movehl_ps(t3, t2) };
}
至于矩阵的求逆操作，采用的是分块求逆：
矩阵分块求逆，A代表左上角2x2的submatrix，B为右上角2x2submatrix，C,D同理
分块求逆的原理和实现可以参考：Fast 4x4 Matrix Inverse with SSE SIMD, Explained
到这里整个数学运算库差不多介绍完毕了。
当然使用SSE/AVX的方法也不止使用数学运算库这一种方式，下面给出了一个10000个数字求和的例子，使用SSE指令的求和运算，其速度为标量求和运算的3倍多：
void test()
{constexpr size_t num = 10000;float *vars = static_cast(_aligned_malloc(sizeof(float) * num, 16));for (size_t i = 0; i < num; i++){vars[i] = 1;}float sum = 0.0f;auto t0 = std::chrono::steady_clock::now();//Scalarfor (size_t i = 0; i < 10000; i++){sum += vars[i];}auto t1 = std::chrono::steady_clock::now();std::cout << "Scalar sum！" << std::endl;std::chrono::duration time_span = std::chrono::duration_cast>(t1 - t0);std::cout << "Time Span: " << time_span.count() << std::endl;std::cout << sum << std::endl;AMVector vector_sum = am_vector_set(0, 0, 0, 0);t0 = std::chrono::steady_clock::now();//SSEfor (size_t i = 0; i < 10000; i += 4){vector_sum += am_vector_load_aligned(vars + i);}sum = am_vector4_hadd(vector_sum);t1 = std::chrono::steady_clock::now();std::cout << "SSE sum！" << std::endl;time_span = std::chrono::duration_cast>(t1 - t0);std::cout << "Time Span: " << time_span.count() << std::endl;std::cout << sum << std::endl;_aligned_free(vars);
}
上面的例子还可以扩展出更多内容，不过本篇文章的内容已经够多了，所以留着以后有机会再讲吧...
如果你觉得本篇文章的内容对你有帮助，请点赞！！！
引用：
- https://software.intel.com/zh-cn/articles/ticker-tape-part-2
- Writing C++ Wrappers for SIMD Intrinsics (1)
- How To Write A Maths Library In 2016
- Fast 4x4 Matrix Inverse with SSE SIMD, Explained
- 在C/C++代码中使用SSE等指令集的指令(3)SSE指令集基础 - 。。。。 - CSDN博客
- Easy SIMD through Wrappers
- C++中使用SIMD的几种方法 - 道道道人间道 - CSDN博客
如果你想查看所有的SSE/AVX指令集可以查看Intrinsics Guide
https://db.in.tum.de/~finis/x86%20intrinsics%20cheat%20sheet%20v1.0.pdf 这个pdf也罗列出了大部分的SSE/AVX指令集，以及对应的功能
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！