自从我熟悉Intel处理器上的MMX，SSE和后来的AVX指令以来，已经过去了很多年。 一次，在x86汇编程序的背景下，它们似乎像是一种魔术，而这早已是世俗的事物。 他们深深吸引了我，几年前，我想到了为一个著名游戏编写自己的软件渲染器的想法。 答应这些指示的事情向我保证了这一点。 在某些时候，我什至考虑过编写它。 但是事实证明，编写文本比编写代码要复杂得多。
当时，我想避免在不同处理器上支持的问题。 我希望能够检查我的渲染器的最大可用数量。 我仍然和旧的AMD处理器有朋友，他们的上限是SSE3。 因此，当时我决定将自己限制为SSE3的最大值。 因此，有一个向量数学库，比在SSE上完全实现的要少一些，在SSE3之前很少包含。 但是，在某些时候，我想知道对于一些关键的矢量数学运算，我可以从处理器中挤出什么最大性能。 一种这样的操作是将4个矩阵乘以4。
实际上，为了娱乐起见，我决定更多地从事这项业务。 我已经写过并且一直在我在SSE上的软件渲染中使用矩阵乘法，这对我来说似乎足够了。 但是后来我决定看看将2个float4x4矩阵相乘可以从原则上挤出多少个度量。 在我当前的SSE上，这是16个时钟周期。 的确，自从我开始为某些指令写1 *而不是0 *以来，最近到
IACA 3的过渡开始显示19。 显然，这只是分析仪中的一个缺陷。
简要介绍所使用的实用程序
对于代码分析，我使用了著名的
英特尔架构代码分析器实用程序。 为了进行分析，我使用Haswell体系结构（HSW）作为支持AVX2的最低要求。 编写代码也非常方便使用：《
英特尔技术指南》和《
英特尔优化手册》 。
对于组装，我从控制台使用MSVS 2017社区。 我在带有内在函数的版本中编写代码。 您只需编写一次，通常它可以立即在不同平台上运行。 另外，x64 VC ++编译器不支持嵌入式汇编程序，但我希望它也能在x64下工作。
由于本文已经超出了SIMD编程的初学者水平，因此我将不介绍寄存器，指令，绘制（或刮削）精美的图片，而是尝试使用SIMD指令学习编程。 英特尔网站上有许多出色，清晰和详细的文档。
我想让一切变得更轻松。。。
这是时刻的开始，这使实现和本文都变得很复杂。 因此，我将对此稍作介绍。 对我来说，用元素的标准行布局编写矩阵乘法并不有趣。 谁需要它，所以他们就读大学或自己学习。 我们的目标是提高生产力。 首先，很久以前我切换到列布局。 我的软件渲染器基于OpenGL API，因此，为了避免不必要的转置，我开始将元素存储在列中。 这也很重要，因为矩阵乘法不是那么关键。 乘以2-5-10矩阵。 就是这样。 然后，我们将完成的矩阵乘以成千上万个顶点。 而且此操作更为关键。 当然，您每次都可以转置。 但是为什么，如果可以避免的话。
但是只回到矩阵。 我们确定了列的存储量。 但是，这可能会更加复杂。 对于我来说，将向量和矩阵行的高级元素存储在SIMD寄存器中比较方便，以便
x处于最高浮点数（索引3），而
w处于次要浮点（索引0）。 显然，在这里，我们将不得不退缩。
事实是，在矢量的软件渲染器中，您必须更频繁地操作
w组件（
1 / z存储在其中），并且通过
_ss版本的操作（仅在
xmm寄存器的下浮点中的组件进行操作）非常方便，而无需触摸
x，y， 。 因此，在SSE寄存器中，向量以易于理解的顺序
x，y，z，w存储，并以相反的
w，z，y，x顺序存储在存储器中。
此外，所有乘法选项也由单个功能实现。 这样做是因为我根据支持的指令类型使用它们来替代所需的选项。
好在这里描述。我们实现基本功能
带循环的乘法，行有序
元素的行布局选项for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { r[i][j] = 0.f; for (int k = 0; k < 4; ++k) { r[i][j] += m[i][k] * n[k][j]; } } }
这里的一切都很简单明了。 对于每个元素，我们进行4次乘法和3次加法。 总共有64个乘法和48个加法。 而且这没有考虑到记录元素的读取。
简而言之，一切都很难过。 对于此选项，对于内部周期，IACA发出：
x86组装为3.65个时钟周期，x64组装为2.97个时钟 。 不要问为什么小数。 我不知道 IACA 2.1并未遭受此困扰。 无论如何，这些数字应乘以约4 * 4 * 4 =64。即使您使用x64，结果也约为192小节。 显然，这是一个粗略的估计。 对于该选项，我看不出更准确地评估性能的意义。
循环实现，列有序
转置矩阵，重新排列行和列索引 for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { r[j][i] = 0.f; for (int k = 0; k < 4; ++k) { r[j][i] += m[k][i] * n[j][k]; } } }
周期乘法，面向SIMD的存储
以相反的顺序将行存储在内存中 for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { r[j][3-i] = 0.f; for (int k = 0; k < 4; ++k) { r[j][3-i] += m[k][3-i] * n[j][3-k]; } } }
这种实现在某种程度上简化了对内部发生的事情的理解，但显然还不够。
助手类
为了方便理解和编写参考和调试代码，实现两个辅助类很方便。 仅此而已，一切都只是为了理解。 我注意到，完整的向量和矩阵类的实现是一个单独的难题，并且未包含在本文的主题中。
向量和矩阵类 struct alignas(sizeof(__m128)) vec4 { union { struct { float w, z, y, x; }; __m128 fmm; float arr[4]; }; vec4() {} vec4(float a, float b, float c, float d) : w(d), z(c), y(b), x(a) {} static bool equ(float const a, float const b, float t = .00001f) { return fabs(ab) < t; } bool operator == (vec4 const& v) const { return equ(x, vx) && equ(y, vy) && equ(z, vz) && equ(w, vw); } }; struct alignas(sizeof(__m256)) mtx4 {
测试参考功能
由于矩阵中元素的可接受顺序使理解变得很复杂，因此我们也不会被引用
清除功能所打扰，该功能将在以后的实现中显示一切正常。 我们将与之比较后续结果。
要创建它，只需采取并扩大周期 void mul_mtx4_mtx4_unroll(__m128* const _r, __m128 const* const _m, __m128 const* const _n) { mtx4 const& m = **reinterpret_cast<mtx4 const* const*>(&_m); mtx4 const& n = **reinterpret_cast<mtx4 const* const*>(&_n); mtx4& r = **reinterpret_cast<mtx4* const*>(&_r); r._00 = m._00*n._00 + m._01*n._10 + m._02*n._20 + m._03*n._30; r._01 = m._00*n._01 + m._01*n._11 + m._02*n._21 + m._03*n._31; r._02 = m._00*n._02 + m._01*n._12 + m._02*n._22 + m._03*n._32; r._03 = m._00*n._03 + m._01*n._13 + m._02*n._23 + m._03*n._33; r._10 = m._10*n._00 + m._11*n._10 + m._12*n._20 + m._13*n._30; r._11 = m._10*n._01 + m._11*n._11 + m._12*n._21 + m._13*n._31; r._12 = m._10*n._02 + m._11*n._12 + m._12*n._22 + m._13*n._32; r._13 = m._10*n._03 + m._11*n._13 + m._12*n._23 + m._13*n._33; r._20 = m._20*n._00 + m._21*n._10 + m._22*n._20 + m._23*n._30; r._21 = m._20*n._01 + m._21*n._11 + m._22*n._21 + m._23*n._31; r._22 = m._20*n._02 + m._21*n._12 + m._22*n._22 + m._23*n._32; r._23 = m._20*n._03 + m._21*n._13 + m._22*n._23 + m._23*n._33; r._30 = m._30*n._00 + m._31*n._10 + m._32*n._20 + m._33*n._30; r._31 = m._30*n._01 + m._31*n._11 + m._32*n._21 + m._33*n._31; r._32 = m._30*n._02 + m._31*n._12 + m._32*n._22 + m._33*n._32; r._33 = m._30*n._03 + m._31*n._13 + m._32*n._23 + m._33*n._33; }
此处清楚地画出了经典算法，很难犯错（但您可以:-）。 IACA在其上发布了：
x86-69.95措施，x64-64措施 。 这里大约有64个周期，我们将在将来查看此操作的加速情况。
SSE实施
经典SSE算法
为什么要经典？ 因为它早已在
FVec的实现中作为MSVS的一部分。 首先，我们将编写如何在SSE寄存器中显示矩阵元素。 在这里看起来已经很简单了。 只是一个转置矩阵。
// 00, 10, 20, 30 // m[0] - SIMD / 01, 11, 21, 31 // m[1] 02, 12, 22, 32 // m[2] 03, 13, 23, 33 // m[3]
我们采用上述变体的
展开代码。 不知何故他对上交所并不友好。 第一组行由结果矩阵的列的结果组成：
r._00，r._01，r._02，r._03 。 我们有此列，但需要一行。 是的
，M，N出现计算不便。 因此，我们重新排列算法的行，以使结果
r是逐行的。
// , r[0] r00 = m00*n00 + m01*n10 + m02*n20 + m03*n30; r10 = m10*n00 + m11*n10 + m12*n20 + m13*n30; r20 = m20*n00 + m21*n10 + m22*n20 + m23*n30; r30 = m30*n00 + m31*n10 + m32*n20 + m33*n30; // , r[1] r01 = m00*n01 + m01*n11 + m02*n21 + m03*n31; r11 = m10*n01 + m11*n11 + m12*n21 + m13*n31; r21 = m20*n01 + m21*n11 + m22*n21 + m23*n31; r31 = m30*n01 + m31*n11 + m32*n21 + m33*n31; // , r[2] r02 = m00*n02 + m01*n12 + m02*n22 + m03*n32; r12 = m10*n02 + m11*n12 + m12*n22 + m13*n32; r22 = m20*n02 + m21*n12 + m22*n22 + m23*n32; r32 = m30*n02 + m31*n12 + m32*n22 + m33*n32; // , r[3] r03 = m00*n03 + m01*n13 + m02*n23 + m03*n33; r13 = m10*n03 + m11*n13 + m12*n23 + m13*n33; r23 = m20*n03 + m21*n13 + m22*n23 + m23*n33; r33 = m30*n03 + m31*n13 + m32*n23 + m33*n33;
但这已经好多了。 实际上，我们看到了什么？ 根据每组中算法的列，我们具有矩阵
m的行：
m [0] = {00,10,20,30}，m [1] = {01,11,21,31}，m [2] = {02,12,22,32}，m [3] = {03,13,23,33}，
它们乘以矩阵
n的相同元素。 例如，对于第一个组，它是：
n._00，n._10，n._20，n._30 。 对于算法的每一行行，矩阵
n的元素再次位于矩阵的一行中。
然后一切都变得很简单：我们只需按索引取矩阵
m的行，但对于元素
n取其行，并通过
shuffle指令将其分散到寄存器的所有4个元素，以便乘以寄存器中矩阵
m的行。 例如，对于元素
n._00 （记住，它在寄存器中的移位具有索引3），它将是：
_mm_shuffle_ps（n [0]，n [0]，_ MM_SHUFFLE（3,3,3,3））
以简化形式，该算法如下所示：
// n[0]={00,10,20,30} r[0] = m[0] * n00 + m[1] * n10 + m[2] * n20 + m[3] * n30; // n[1]={01,11,21,31} r[1] = m[0] * n01 + m[1] * n11 + m[2] * n21 + m[3] * n31; // n[2]={02,12,22,32} r[2] = m[0] * n02 + m[1] * n12 + m[2] * n22 + m[3] * n32; // n[3]={03,13,23,33} r[3] = m[0] * n03 + m[1] * n13 + m[2] * n23 + m[3] * n33;
基本的SSE实施 void mul_mtx4_mtx4_sse_v1(__m128* const r, __m128 const* const m, __m128 const* const n) { r[0] = _mm_add_ps( _mm_add_ps( _mm_mul_ps(m[0], _mm_shuffle_ps(n[0], n[0], _MM_SHUFFLE(3,3,3,3))), _mm_mul_ps(m[1], _mm_shuffle_ps(n[0], n[0], _MM_SHUFFLE(2,2,2,2)))), _mm_add_ps( _mm_mul_ps(m[2], _mm_shuffle_ps(n[0], n[0], _MM_SHUFFLE(1,1,1,1))), _mm_mul_ps(m[3], _mm_shuffle_ps(n[0], n[0], _MM_SHUFFLE(0,0,0,0))))); r[1] = _mm_add_ps( _mm_add_ps( _mm_mul_ps(m[0], _mm_shuffle_ps(n[1], n[1], _MM_SHUFFLE(3,3,3,3))), _mm_mul_ps(m[1], _mm_shuffle_ps(n[1], n[1], _MM_SHUFFLE(2,2,2,2)))), _mm_add_ps( _mm_mul_ps(m[2], _mm_shuffle_ps(n[1], n[1], _MM_SHUFFLE(1,1,1,1))), _mm_mul_ps(m[3], _mm_shuffle_ps(n[1], n[1], _MM_SHUFFLE(0,0,0,0))))); r[2] = _mm_add_ps( _mm_add_ps( _mm_mul_ps(m[0], _mm_shuffle_ps(n[2], n[2], _MM_SHUFFLE(3,3,3,3))), _mm_mul_ps(m[1], _mm_shuffle_ps(n[2], n[2], _MM_SHUFFLE(2,2,2,2)))), _mm_add_ps( _mm_mul_ps(m[2], _mm_shuffle_ps(n[2], n[2], _MM_SHUFFLE(1,1,1,1))), _mm_mul_ps(m[3], _mm_shuffle_ps(n[2], n[2], _MM_SHUFFLE(0,0,0,0))))); r[3] = _mm_add_ps( _mm_add_ps( _mm_mul_ps(m[0], _mm_shuffle_ps(n[3], n[3], _MM_SHUFFLE(3,3,3,3))), _mm_mul_ps(m[1], _mm_shuffle_ps(n[3], n[3], _MM_SHUFFLE(2,2,2,2)))), _mm_add_ps( _mm_mul_ps(m[2], _mm_shuffle_ps(n[3], n[3], _MM_SHUFFLE(1,1,1,1))), _mm_mul_ps(m[3], _mm_shuffle_ps(n[3], n[3], _MM_SHUFFLE(0,0,0,0))))); }
现在我们将算法中的元素
n更改为对应的
shuffle ，乘以
_mm_mul_ps ，乘以
_mm_add_ps ，就可以了。 可以用 但是，代码看起来比算法本身看起来差很多。 IACA对此代码发布：
x86-18.89，x64-16个周期 。 这比上一个快4倍。 在SSE中注册第四个浮点数 几乎线性关系。
装饰SSE实施
尽管如此，在代码中看起来还是很糟糕。 我们将尝试通过编写一些语法糖来改善这一点。
编译器可以完美地内联这些函数（尽管有时不带__forceinline）。
因此代码正在转向... void mul_mtx4_mtx4_sse_v2(__m128* const r, __m128 const* const m, __m128 const* const n) { r[0] = m[0]*shuf<3>(n[0]) + m[1]*shuf<2>(n[0]) + m[2]*shuf<1>(n[0]) + m[3]*shuf<0>(n[0]); r[1] = m[0]*shuf<3>(n[1]) + m[1]*shuf<2>(n[1]) + m[2]*shuf<1>(n[1]) + m[3]*shuf<0>(n[1]); r[2] = m[0]*shuf<3>(n[2]) + m[1]*shuf<2>(n[2]) + m[2]*shuf<1>(n[2]) + m[3]*shuf<0>(n[2]); r[3] = m[0]*shuf<3>(n[3]) + m[1]*shuf<2>(n[3]) + m[2]*shuf<1>(n[3]) + m[3]*shuf<0>(n[3]); }
因此，它已经更好，更易读。 IACA它发布了关于预期的结果
：86 - 19，64 - 16（什么是不是小数？）。 实际上，性能没有改变，但是代码更加漂亮和易于理解。
对未来优化的贡献很小
让我们在铁版本中最近出现的功能级别上进行另一项改进。
多次加法（fma）操作 。
fma（a，b，c）= a * b + c 。
多次添加实施 __m128 mad(__m128 const a, __m128 const b, __m128 const c) { return _mm_add_ps(_mm_mul_ps(a, b), c); }
为什么这是必要的？ 首先，为了将来的优化。 例如，您可以根据需要通过相同的宏在完成的代码中简单地用
mad替换
mad 。 但是，我们现在将为优化奠定基础：
有多次添加的变体 void mul_mtx4_mtx4_sse_v3(__m128* const r, __m128 const* const m, __m128 const* const n) { r[0] = mad(m[0], shuf<3>(n[0]), m[1]*shuf<2>(n[0])) + mad(m[2], shuf<1>(n[0]), m[3]*shuf<0>(n[0])); r[1] = mad(m[0], shuf<3>(n[1]), m[1]*shuf<2>(n[1]) + mad(m[2], shuf<1>(n[1]), m[3]*shuf<0>(n[1])); r[2] = mad(m[0], shuf<3>(n[2]), m[1]*shuf<2>(n[2])) + mad(m[2], shuf<1>(n[2]), m[3]*shuf<0>(n[2])); r[3] = mad(m[0], shuf<3>(n[3]), m[1]*shuf<2>(n[3])) + mad(m[2], shuf<1>(n[3]), m[3]*shuf<0>(n[3])); }
IACA：
x86-18.89，x64-16 。 再次分数。 不过，IACA有时仍会产生奇怪的结果。 代码变化不大。 可能更糟。 但是优化有时需要做出这样的牺牲。
我们通过_mm_stream进行保存
各种优化指南再次建议不要为大容量保存操作拉入缓存。 当您处理成千上万个顶点时，通常这是合理的。 但是对于矩阵而言，这可能并不那么重要。 不过，我还是会添加它。
流保存选项 void mul_mtx4_mtx4_sse_v4(__m128* const r, __m128 const* const m, __m128 const* const n) { _mm_stream_ps(&r[0].m128_f32[0], mad(m[0], shuf<3>(n[0]), m[1]*shuf<2>(n[0])) + mad(m[2], shuf<1>(n[0]), m[3]*shuf<0>(n[0]))); _mm_stream_ps(&r[1].m128_f32[0], mad(m[0], shuf<3>(n[1]), m[1]*shuf<2>(n[1])) + mad(m[2], shuf<1>(n[1]), m[3]*shuf<0>(n[1]))); _mm_stream_ps(&r[2].m128_f32[0], mad(m[0], shuf<3>(n[2]), m[1]*shuf<2>(n[2])) + mad(m[2], shuf<1>(n[2]), m[3]*shuf<0>(n[2]))); _mm_stream_ps(&r[3].m128_f32[0], mad(m[0], shuf<3>(n[3]), m[1]*shuf<2>(n[3])) + mad(m[2], shuf<1>(n[3]), m[3]*shuf<0>(n[3]))); }
从字面上看，这里的时间没有任何变化。 但是，根据建议，我们现在不再再次触摸缓存。
AVX实施
基本AVX选项
接下来，我们进入下一步的优化。 SSE寄存器中包含第4个浮点数，而AVX中已经包含第4个浮点数。也就是说，理论上有机会减少执行的操作数量，如果不是一半，则提高生产率，至少是1.5倍。 但是有一点告诉我，过渡到AVX并非一切都会如此简单。 我们可以从双寄存器中获得必要的数据吗？
让我们尝试找出答案。 同样，我们写出上面使用的乘法算法。 您无法执行此操作，但是当所有内容都在附近并且无需向上滚动半页时，使用该代码更方便。
// : 00, 10, 20, 30, 01, 11, 21, 31, 02, 12, 22, 32, 03, 13, 23, 33 // SSE: r0 = m0*n00 + m1*n10 + m2*n20 + m3*n30 r1 = m0*n01 + m1*n11 + m2*n21 + m3*n31 r2 = m0*n02 + m1*n12 + m2*n22 + m3*n32 r3 = m0*n03 + m1*n13 + m2*n23 + m3*n33
在输出处，我们期望得到
ymm = {r0：r1}和
ymm = {r2：r3}的结果 。 如果在SSE版本中我们的算法被泛化为列，那么现在我们需要将其泛化为行。 因此，在SSE选项的情况下无法正常工作。
如果考虑寄存器
ymm中的矩阵
m ，则分别得到
ymm = {m0：m1}和
ymm = {m2：m3} 。 以前，我们在寄存器中只有矩阵列，现在只有列和行。
如果您尝试像以前那样操作，则必须将
ymm = {m0：m1}乘以寄存器
ymm = {n00，n00，n00，n00}：{n10，n10，n10，n10} 。 由于
n00和
n01在矩阵
n的同一行中，根据可用的AVX指令集判断，将它们分散
ymm会很昂贵。
ymm寄存器内的float的两个四分之一（高
xmm和低
xmm ）中的每一个都分别进行
混洗和
置换 。
如果我们假设
n个 YMM矩阵
，我们得到
N00和
N10在较旧的2内部
YMM xmm寄存器中两个元件。
{n00，n10，n20，n30}：{n01，n11，n21，n31} 。 通常，现有指令的索引为0到3。它的地址仅在两个内部
ymm寄存器中的一个
xmm寄存器中浮动。 无法
廉价地将
n10从较旧的
xmm转移到较年轻的
xmm 。 然后必须多次重复此焦点。 我们不能忍受这种措施的损失。 有必要提出其他建议。
我们曾经对列进行泛化，但现在对行进行泛化。 因此，我们将尝试以一些不同的方式进行。 我们需要在
{r0：r1}中得到结果。 这意味着必须对算法进行改进，而不是在算法的不同行中进行改进，而应同时进行两次改进。 在这里，
洗牌和
置换工作的
弊端对我们来说将是一个加分。 当我们考虑矩阵
n时，我们将看看
ymm寄存器中的
内容 。
n0n1 = {00, 10, 20, 30} : {01, 11, 21, 31} n2n3 = {02, 12, 22, 32} : {03, 13, 23, 33}
是的，我们注意到在
ymm寄存器的不同
xmm部分中，我们有元素
00和
01 。 可以通过
{_00，_00，_00，_00}中的permute命令通过寄存器乘以大小写
：{_ 01，_01，_01，_01} ，这两个
xmm部件仅指示一个索引3。 这正是我们所需要的。 实际上，系数也用在不同的行中。 仅现在才在对应的
ymm寄存器中进行乘法运算，有必要保留
{m0：m0} ，即矩阵
m的第一行重复。
因此，我们将更详细地描述算法。 我们在
ymm寄存器中读取矩阵
m的双行：
mm[0] = {m0:m0} mm[1] = {m1:m1} mm[2] = {m2:m2} mm[3] = {m3:m3}
然后我们将乘法计算为：
r0r1 = mm[0] * {n00,n00,n00,n00:n01,n01,n01,n01} + // permute<3,3,3,3>(n0n1) mm[1] * {n10,n10,n10,n10:n11,n11,n11,n11} + // permute<2,2,2,2>(n0n1) mm[2] * {n20,n20,n20,n20:n21,n21,n21,n21} + // permute<1,1,1,1>(n0n1) mm[3] * {n30,n30,n30,n30:n31,n31,n31,n31} // permute<0,0,0,0>(n0n1) r2r3 = mm[0] * {n02,n02,n02,n02:n03,n03,n03,n03} + // permute<3,3,3,3>(n2n3) mm[1] * {n12,n12,n12,n12:n13,n13,n13,n13} + // permute<2,2,2,2>(n2n3) mm[2] * {n22,n22,n22,n22:n23,n23,n23,n23} + // permute<1,1,1,1>(n2n3) mm[3] * {n32,n32,n32,n32:n33,n33,n33,n33} // permute<0,0,0,0>(n2n3)
我们重写得更清楚：
r0r1 = mm[0]*n0n1<3,3,3,3>+mm[1]*n0n1<2,2,2,2>+mm[2]*n0n1<1,1,1,1>+mm[3]*n0n1<0,0,0,0> r2r3 = mm[0]*n2n3<3,3,3,3>+mm[1]*n2n3<2,2,2,2>+mm[2]*n2n3<1,1,1,1>+mm[3]*n2n3<0,0,0,0>
或以简化形式：
r0r1 = mm[0]*n0n1<3> + mm[1]*n0n1<2> + mm[2]*n0n1<1> + mm[3]*n0n1<0> r2r3 = mm[0]*n2n3<3> + mm[1]*n2n3<2> + mm[2]*n2n3<1> + mm[3]*n2n3<0>
一切似乎都很清楚。
剩下的只是写一个实现 void mul_mtx4_mtx4_avx_v1(__m128* const r, __m128 const* const m, __m128 const* const n) { __m256 mm0 = _mm256_set_m128(m[0], m[0]); __m256 mm1 = _mm256_set_m128(m[1], m[1]); __m256 mm2 = _mm256_set_m128(m[2], m[2]); __m256 mm3 = _mm256_set_m128(m[3], m[3]); __m256 n0n1 = _mm256_load_ps(&n[0].m128_f32[0]); __m256 y1 = _mm256_permute_ps(n0n1, 0xFF);
它一直有趣的数字从
IACA：86 - 12.53，64 - 12。 虽然，当然，我想要更好。 错过了什么。
AVX优化加语法糖
似乎在上面的代码中，AVX并未充分发挥其潜力。 我们发现不用在
ymm寄存器中设置两条相同的行，而是可以使用
broadcast ，它可以用两个相同的
xmm值填充
ymm寄存器。 同样，在此过程中，为AVX功能添加一些“语法糖”。
增强的AVX实施 __m256 operator + (__m256 const a, __m256 const b) { return _mm256_add_ps(a, b); } __m256 operator - (__m256 const a, __m256 const b) { return _mm256_sub_ps(a, b); } __m256 operator * (__m256 const a, __m256 const b) { return _mm256_mul_ps(a, b); } __m256 operator / (__m256 const a, __m256 const b) { return _mm256_div_ps(a, b); } template <int i> __m256 perm(__m256 const v) { return _mm256_permute_ps(v, _MM_SHUFFLE(i, i, i, i)); } template <int a, int b, int c, int d> __m256 perm(__m256 const v) { return _mm256_permute_ps(v, _MM_SHUFFLE(a, b, c, d)); } template <int i, int j> __m256 perm(__m256 const v) { return _mm256_permutevar_ps(v, _mm256_set_epi32(i, i, i, i, j, j, j, j)); } template <int a, int b, int c, int d, int e, int f, int g, int h> __m256 perm(__m256 const v) { return _mm256_permutevar_ps(v, _mm256_set_epi32(a, b, c, d, e, f, g, h)); } __m256 mad(__m256 const a, __m256 const b, __m256 const c) { return _mm256_add_ps(_mm256_mul_ps(a, b), c); } void mul_mtx4_mtx4_avx_v2(__m128* const r, __m128 const* const m, __m128 const* const n) { __m256 const mm[] { _mm256_broadcast_ps(m+0), _mm256_broadcast_ps(m+1), _mm256_broadcast_ps(m+2), _mm256_broadcast_ps(m+3) }; __m256 const n0n1 = _mm256_load_ps(&n[0].m128_f32[0]); _mm256_stream_ps(&r[0].m128_f32[0], mad(perm<3>(n0n1), mm[0], perm<2>(n0n1)*mm[1])+ mad(perm<1>(n0n1), mm[2], perm<0>(n0n1)*mm[3])); __m256 const n2n3 = _mm256_load_ps(&n[2].m128_f32[0]); _mm256_stream_ps(&r[2].m128_f32[0], mad(perm<3>(n2n3), mm[0], perm<2>(n2n3)*mm[1])+ mad(perm<1>(n2n3), mm[2], perm<0>(n2n3)*mm[3])); }
结果在这里已经更加有趣了。 IACA产生数字：x86-10
，x64-8.58 ，看起来好多了，但仍然不是2倍。
AVX + FMA选项（最终）
让我们再尝试一次。 现在，再次调用FMA指令集是合乎逻辑的，因为它是在AVX之后添加到处理器中的。 只需更改单个
mul +添加一次操作即可。 尽管我们仍然使用乘法指令为编译器提供更多的优化机会，并为处理器提供并行执行乘法的机会。 通常，我会在汇编器中查看生成的代码，以确保哪个选项更好。
在这种情况下，我们需要计算
a * b + c * d + e * f + g * h 。 您可以这样做：
fma（a，b，fma（c，d，fma（e，f，g * h））） 。 但是，正如我们看到的，如果不完成前一个操作就不可能在这里执行一个操作。 这意味着我们将无法使用成对乘法的功能，因为SIMD流水线允许我们这样做。 如果我们转换位计算
FMA（A，B，C * d）+ FMA（E，F，G * H），我们看到，计算可以被并行化。 首先执行两个独立的乘法，然后执行两个独立的
fma运算。
AVX + FMA实施 __m256 fma(__m256 const a, __m256 const b, __m256 const c) { return _mm256_fmadd_ps(a, b, c); } void mul_mtx4_mtx4_avx_fma(__m128* const r, __m128 const* const m, __m128 const* const n) { __m256 const mm[]{ _mm256_broadcast_ps(m + 0), _mm256_broadcast_ps(m + 1), _mm256_broadcast_ps(m + 2), _mm256_broadcast_ps(m + 3) }; __m256 const n0n1 = _mm256_load_ps(&n[0].m128_f32[0]); _mm256_stream_ps(&r[0].m128_f32[0], fma(perm<3>(n0n1), mm[0], perm<2>(n0n1)*mm[1])+ fma(perm<1>(n0n1), mm[2], perm<0>(n0n1)*mm[3])); __m256 const n2n3 = _mm256_load_ps(&n[2].m128_f32[0]); _mm256_stream_ps(&r[2].m128_f32[0], fma(perm<3>(n2n3), mm[0], perm<2>(n2n3)*mm[1])+ fma(perm<1>(n2n3), mm[2], perm<0>(n2n3)*mm[3])); }
IACA：
x86-9.21，x64-8 。 现在很好。 可能有人会说可以做得更好的方法，但我不知道如何做。
基准测试
我立即注意到，这些数字不应被视为最终真理。 即使进行了固定测试，它们也会在一定范围内游泳。 而且，它们在不同平台上的行为也有所不同。 .
- Function: . s — , mov ( ). , .
- IACA cycles: IACA
- Measured cycles: ( , )
- IACA speedup: /
- Measured speedup: / ( , )
loop_m 64. . .
i3-3770:
i7-8700K:
. , .
BONUS
, , AVX512, . , . , AVX+FMA. , .
, __m512 operator + (__m512 const a, __m512 const b) { return _mm512_add_ps(a, b); } __m512 operator - (__m512 const a, __m512 const b) { return _mm512_sub_ps(a, b); } __m512 operator * (__m512 const a, __m512 const b) { return _mm512_mul_ps(a, b); } __m512 operator / (__m512 const a, __m512 const b) { return _mm512_div_ps(a, b); } template <int i> __m512 perm(__m512 const v) { return _mm512_permute_ps(v, _MM_SHUFFLE(i, i, i, i)); } template <int a, int b, int c, int d> __m512 perm(__m512 const v) { return _mm512_permute_ps(v, _MM_SHUFFLE(a, b, c, d)); } __m512 fma(__m512 const a, __m512 const b, __m512 const c) { return _mm512_fmadd_ps(a, b, c); } void mul_mtx4_mtx4_avx512(__m128* const r, __m128 const* const m, __m128 const* const _n) { __m512 const mm[]{ _mm512_broadcast_f32x4(m[0]), _mm512_broadcast_f32x4(m[1]), _mm512_broadcast_f32x4(m[2]), _mm512_broadcast_f32x4(m[3]) }; __m512 const n = _mm512_load_ps(&_n[0].m128_f32[0]); _mm512_stream_ps(&r[0].m128_f32[0], fma(perm<3>(n), mm[0], perm<2>(n)*mm[1])+ fma(perm<1>(n), mm[2], perm<0>(n)*mm[3])); }
:
x86 — 4.79, x64 — 5.42 (IACA SKX). , 64 48 .
PS
. . .