Уже немало лет прошло, как я познакомился с инструкциями MMX, SSE, а позже и AVX на процессорах Intel. В своё время они казались какой-то магией на фоне x86 ассемблера, который уже давно стал чем-то обыденным. Они меня настолько зацепили, что пару лет назад у меня появилась идея написать свой собственный софт рендерер для одной известной игры. Сподвигло меня на это то, какую производительность обещали эти инструкции. В какой-то момент я даже думал об этом написать. Но писать текст оказалось куда сложнее кода.
В то время я хотел избежать проблем с поддержкой на разных процессорах. Хотелось иметь возможность проверить мой рендерер на максимально доступном количестве. У меня до сих пор остались знакомые со старыми AMD процессорами, и их потолок был SSE3. Поэтому на тот момент я решил ограничиться максимум SSE3. Так появилась векторная математическая библиотека, чуть менее, чем полностью реализованная на SSE, с редким включением до SSE3. Однако в какой-то момент мне стало интересно, какую максимальную производительность я смогу выжать из процессора для ряда критичных операций векторной математики. Одной из таких операций является умножение матриц float 4 на 4.
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];
}
}
}
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];
}
}
}
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(a-b) < t;
}
bool operator == (vec4 const& v) const {
return equ(x, v.x) && equ(y, v.y) && equ(z, v.z) && equ(w, v.w);
}
};
struct alignas(sizeof(__m256)) mtx4 {
// тут всё больше для наглядности хранения в памяти и регистрах
union {
struct {
float
_30, _20, _10, _00,
_31, _21, _11, _01,
_32, _22, _12, _02,
_33, _23, _13, _03;
};
__m128 r[4];
__m256 s[2];
vec4 v[4];
};
// добавляем простейшие инициализаторы
mtx4() {}
mtx4(
float i00, float i01, float i02, float i03,
float i10, float i11, float i12, float i13,
float i20, float i21, float i22, float i23,
float i30, float i31, float i32, float i33)
: _00(i00), _01(i01), _02(i02), _03(i03)
, _10(i10), _11(i11), _12(i12), _13(i13)
, _20(i20), _21(i21), _22(i22), _23(i23)
, _30(i30), _31(i31), _32(i32), _33(i33)
{}
// для передачи в реализацию умножения
operator __m128 const* () const { return r; }
operator __m128* () { return r; }
// для тестов
bool operator == (mtx4 const& m) const {
return v[0]==m.v[0] && v[1]==m.v[1] && v[2]==m.v[2] && v[3]==m.v[3];
}
// инициализаторы
static mtx4 identity() {
return mtx4(
1.f, 0.f, 0.f, 0.f,
0.f, 1.f, 0.f, 0.f,
0.f, 0.f, 1.f, 0.f,
0.f, 0.f, 0.f, 1.f);
}
static mtx4 zero() {
return mtx4(
0.f, 0.f, 0.f, 0.f,
0.f, 0.f, 0.f, 0.f,
0.f, 0.f, 0.f, 0.f,
0.f, 0.f, 0.f, 0.f);
}
};
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;
}
// представление матрицы в регистрах
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[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[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 для каждой группы строк алгоритма снова лежат в одной строке матрицы.
_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;
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)))));
}
// превращаем имена функций в обычные удобочитаемые операции (которые всё-таки лучше прятать в namespace)
__m128 operator + (__m128 const a, __m128 const b) { return _mm_add_ps(a, b); }
__m128 operator - (__m128 const a, __m128 const b) { return _mm_sub_ps(a, b); }
__m128 operator * (__m128 const a, __m128 const b) { return _mm_mul_ps(a, b); }
__m128 operator / (__m128 const a, __m128 const b) { return _mm_div_ps(a, b); }
//_mm_shuffle_ps(u, v, _MM_SHUFFLE(3,2,1,0)) превращается в shuf<3,2,1,0>(u, v)
template <int a, int b, int c, int d> __m128 shuf(__m128 const u, __m128 const v)
{ return _mm_shuffle_ps(u, v, _MM_SHUFFLE(a, b, c, d)); }
template <int a, int b, int c, int d> __m128 shuf(__m128 const v)
{ return _mm_shuffle_ps(v, v, _MM_SHUFFLE(a, b, c, d)); }
// облегчённый одноиндексный вариант
template <int i> __m128 shuf(__m128 const u, __m128 const v)
{ return _mm_shuffle_ps(u, v, _MM_SHUFFLE(i, i, i, i)); }
template <int i> __m128 shuf(__m128 const v)
{ return _mm_shuffle_ps(v, v, _MM_SHUFFLE(i, i, i, i)); }
// для float векторов можно попробовать и такой экзотический вариант,
// который иногда даёт профит, а иногда нет
template <int a, int b, int c, int d> __m128 shufd(__m128 const v)
{ return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(v), _MM_SHUFFLE(a, b, c, d))); }
template <int i> __m128 shufd(__m128 const v)
{ return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(v), _MM_SHUFFLE(i, i, i, i))); }
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]);
}
__m128 mad(__m128 const a, __m128 const b, __m128 const c) {
return _mm_add_ps(_mm_mul_ps(a, b), c);
}
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]));
}
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])));
}
//Снова представление матрицы в регистрах:
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
n0n1 = {00, 10, 20, 30} : {01, 11, 21, 31}
n2n3 = {02, 12, 22, 32} : {03, 13, 23, 33}
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);//3,3,3,3
__m256 y2 = _mm256_permute_ps(n0n1, 0xAA);//2,2,2,2
__m256 y3 = _mm256_permute_ps(n0n1, 0x55);//1,1,1,1
__m256 y4 = _mm256_permute_ps(n0n1, 0x00);//0,0,0,0
y1 = _mm256_mul_ps(y1, mm0);
y2 = _mm256_mul_ps(y2, mm1);
y3 = _mm256_mul_ps(y3, mm2);
y4 = _mm256_mul_ps(y4, mm3);
y1 = _mm256_add_ps(y1, y2);
y3 = _mm256_add_ps(y3, y4);
y1 = _mm256_add_ps(y1, y3);
__m256 n2n3 = _mm256_load_ps(&n[2].m128_f32[0]);
__m256 y5 = _mm256_permute_ps(n2n3, 0xFF);
__m256 y6 = _mm256_permute_ps(n2n3, 0xAA);
__m256 y7 = _mm256_permute_ps(n2n3, 0x55);
__m256 y8 = _mm256_permute_ps(n2n3, 0x00);
y5 = _mm256_mul_ps(y5, mm0);
y6 = _mm256_mul_ps(y6, mm1);
y7 = _mm256_mul_ps(y7, mm2);
y8 = _mm256_mul_ps(y8, mm3);
y5 = _mm256_add_ps(y5, y6);
y7 = _mm256_add_ps(y7, y8);
y5 = _mm256_add_ps(y5, y7);
_mm256_stream_ps(&r[0].m128_f32[0], y1);
_mm256_stream_ps(&r[2].m128_f32[0], y5);
}
__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]));
}
__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]));
}
__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]));
}
К сожалению, не доступен сервер mySQL