本篇文章包含的内容有SIMD指令集简介以及简短的practice环节。
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支持的数据类型
不同处理器对于SIMD指令集的支持如下图：
不同处理器对于SIMD指令集的支持
下图给出了使用SIMD的不同方法：
使用SIMD的不同方法
首先是最简单的方法是使用Intel开发的跨平台函数库（IPP，Intel Integrated Performance Primitives ），里面的函数实现都使用了SIMD指令进行优化。
其次是借助于Auto-vectorization(自动矢量化)，借助编译器将标量操作转化为矢量操作。
第三种方法是使用编译器指示符(compiler directive)，如Cilk里的#pragma simd和OpenMP里的#pragma omp simd。如下所示，使用#pragma simd强制循环矢量化：
https://zhuanlan.zhihu.com/p/55327037void add_floats（float * a，float * b，float * c，float * d，float * e，int n）
{
int i;
#pragma simd
for（i = 0; i <n; i ++）
{
a [i] = a [i] + b [i] + c [i] + d [i] + e [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指令。
a.头文件
SSE/AVX指令主要定义于以下一些头文件中：
每一个头文件都包含了之前的所有头文件，所以如果你想要使用SSE4.2以及之前SSE3, SSE2, SSE中的所有函数就只需要包含<nmmintrin.h>头文件。
b.命名规则
SSE/AVX提供的数据类型和函数的命名规则如下：
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提供的算术运算操作包括：
此外还有向下取整，向上取整等运算，这里只列出了浮点数支持的算术运算类型，还有一些整形的算术运算类型未列出。
3).比较运算
SSE/AVX提供的比较运算操作包括：
此外还有一些执行单分量对比的比较运算
4).逻辑运算
SSE/AVX提供的逻辑运算操作包括：
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 <immintrin.h>
#elif defined(_M_IX86_FP) && (_M_IX86_FP == 2)
#define SIMD_LEVEL 1 //Support SSE2
#include <emmintrin.h>
#else
#define SIMD_LEVEL 0
#endif
#endif // !SIMD_LEVEL
#ifndef AMVector
#if SIMD_LEVEL == 0
typedef Vector4<float> AMVector;
typedef const Vector4<float>& CRAMVector;
#else
typedef __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 == 0
float dot = lhs.x() * rhs.x() + lhs.y() * rhs.y() + lhs.z() * rhs.z();
return Vector4<float>{dot};
#elif SIMD_LEVEL == 1
AMVector 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);
#else
return _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 == 0
return 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 == 0
float dot = v.x() * v.x() + v.y() * v.y() + v.z() * v.z();
return Vector4<float>{ std::sqrt(dot) };
#elif SIMD_LEVEL == 1
AMVector 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);
#else
return _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 == 0
float dot = v.x() * v.x() + v.y() * v.y() + v.z() * v.z();
return vector4_divide(v, std::sqrt(dot));
#elif SIMD_LEVEL == 1
AMVector 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));
#else
return _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, 0
AM_MASK mask = _mm_movemask_ps(v); // 7
bool 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 == 1
return _mm_or_ps(_mm_andnot_ps(condition, a), _mm_and_ps(b, condition));
#else
return _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<float> AMMatrix;
typedef const Matrix4x4<float>& 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;
//Row0
AMVector 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);
//Row1
vec = 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);
//Row2
vec = 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);
//Row3
vec = 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<float*>(_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();
//Scalar
for (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<double> time_span = std::chrono::duration_cast<std::chrono::duration<double>>(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();
//SSE
for (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<std::chrono::duration<double>>(t1 - t0);
std::cout << "Time Span: " << time_span.count() << std::endl;
std::cout << sum << std::endl;
_aligned_free(vars);
}
上面的例子还可以扩展出更多内容，不过本篇文章的内容已经够多了，所以留着以后有机会再讲吧...
如果你觉得本篇文章的内容对你有帮助，请点赞！！！
引用：
如果你想查看所有的SSE/AVX指令集可以查看Intrinsics Guide
https://db.in.tum.de/~finis/x86%20intrinsics%20cheat%20sheet%20v1.0.pdf 这个pdf也罗列出了大部分的SSE/AVX指令集，以及对应的功能