Ускоряем умножение матриц float 4x4 с помощью SIMD

Уже немало лет прошло, как я познакомился с инструкциями MMX, SSE, а позже и AVX на процессорах Intel. В своё время они казались какой-то магией на фоне x86 ассемблера, который уже давно стал чем-то обыденным. Они меня настолько зацепили, что пару лет назад у меня появилась идея написать свой собственный софт рендерер для одной известной игры. Сподвигло меня на это то, какую производительность обещали эти инструкции. В какой-то момент я даже думал об этом написать. Но писать текст оказалось куда сложнее кода.

В то время я хотел избежать проблем с поддержкой на разных процессорах. Хотелось иметь возможность проверить мой рендерер на максимально доступном количестве. У меня до сих пор остались знакомые со старыми AMD процессорами, и их потолок был SSE3. Поэтому на тот момент я решил ограничиться максимум SSE3. Так появилась векторная математическая библиотека, чуть менее, чем полностью реализованная на SSE, с редким включением до SSE3. Однако в какой-то момент мне стало интересно, какую максимальную производительность я смогу выжать из процессора для ряда критичных операций векторной математики. Одной из таких операций является умножение матриц float 4 на 4.


Собственно, этим делом решил заняться больше ради развлечения. Давно уже написал и использую умножение матриц для моего софт рендера на SSE и вроде мне хватает. Но тут решил посмотреть, сколько тактов я смогу выжать в принципе из умножения 2-х матриц float4x4. На моём текущем SSE варианте это 16 тактов. Правда недавний переход на IACA 3 стал показывать 19, так как на некоторые инструкции вместо 0* стал писать 1*. Видимо раньше это было просто недоработкой анализатора.

Коротко об используемых утилитах


Для анализа кода использовал известную утилиту Intel Architecture Code Analyzer. Для анализа использую архитектуру Haswell (HSW), как минимальную с поддержкой AVX2. Для написания кода также очень удобно пользоваться: Intel Intrinsics Guide и Intel optimization manual.

Для сборки использую MSVS 2017 Community с консоли. Код пишу в варианте с интринсиками. Пишешь один раз, и обычно он сразу работает на разных платформах. К тому же x64 компилятор на VC++ не поддерживает инлайн ассемблер, а хочется чтобы и под x64 работало.

Поскольку данная статья уже несколько выходит за рамки уровня начинающего в SIMD программировании, я не буду описывать регистры, инструкции, рисовать (или тырить) красивые картинки и пытаться учить программировать с использованием SIMD инструкций. На сайте Intel полно отличной, понятной и подробной документации.

Хотел сделать всё проще.… А получилось как всегда


Вот тут начинается момент, который немало усложняет как реализацию, так и статью. Поэтому немного на нём остановлюсь. Писать умножение матриц со стандартным построчным расположением элементов мне не интересно. Кому было надо, и так изучили в ВУЗах или самостоятельно. Наша же цель — производительность. Во-первых, я давно перешёл на расположение по столбцам. Мой софт рендерер базируется на OpenGL API и поэтому, дабы избежать лишних транспонирований, я начал хранить элементы по столбцам. Также это важно потому, что умножение матриц не так критично. Ну умножили 2-5-10 матриц. И всё. А потом умножаем готовую матрицу на тысячи-миллионы вершин. И вот эта операция уже куда критичнее. Можно, конечно, каждый раз транспонировать. Но зачем, если этого можно избежать.

Но вернёмся исключительно к матрицам. С хранением по столбцам определились. Однако можно ещё усложнить. Мне удобнее хранить старшие элементы векторов и матричных строк в SIMD регистрах так, что x в старшем float (индекс 3), а w в младшем (индекс 0). Тут, видимо, придётся снова сделать отступление на тему почему так.

Дело в том, что в софт рендерере в векторе компонентой w приходится манипулировать чаще (там хранится 1/z), и очень удобно делать это через _ss вариант операции (операции исключительно с компонентой в младшем float регистра xmm), не трогая x, y, z. Поэтому в SSE регистре вектор хранится в понятном порядке x, y, z, w, а в памяти в обратном w, z, y, x.

Далее, все варианты умножения также реализованы отдельными функциями. Так сделано потому, что их я использую для подстановки нужного варианта в зависимости от поддерживаемого типа инструкций. Хорошо описано тут.

Реализуем базовый функционал


Умножение с циклами, row ordered


Вариант для построчного расположения элементов
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 выдала: 3.65 тактов под x86 сборку и 2.97 под x64 сборку. Не спрашивайте, почему дробные цифры. Не знаю. IACA 2.1 подобным не страдала. В любом случае, эти цифры надо умножить примерно на 4*4*4 = 64. Даже если взять x64, в результате получается около 192 тактов. Понятно, что это примерная оценка. Оценивать производительность точнее для данного варианта не вижу смысла.

Реализация с циклами, column ordered


транспонированная матрица, переставляем индексы строк и столбцов
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(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;
}


Здесь наглядно расписан классический алгоритм, ошибиться сложно (но можно :-) ). На него 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]

Берём код unroll варианта выше. Какой-то он недружелюбный для SSE. Первая группа строк состоит из результатов по столбцу результирующей матрицы: 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 регистре 4-е float. Почти линейная зависимость.

Украшаем SSE реализацию


И всё-таки в коде это выглядит ужасно. Попытаемся это улучшить, написав немного синтаксического сахара.

Операторы и улучшаторы
// превращаем имена функций в обычные удобочитаемые операции (которые всё-таки лучше прятать в 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))); }


Данные функции компилятор умеет отлично инлайнить (хотя иногда без __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 выдала примерно ожидаемый результат: x86 — 19 (а чего не дробный?), x64 — 16. По сути производительность не изменилась, но код намного красивее и понятнее.

Небольшой вклад в будущую оптимизацию


Введём ещё одно улучшение на уровне функции, не так уж давно появившейся в железном варианте. Операцию multiple-add (fma). fma(a, b, c) = a * b + c.

Реализация multiple-add
__m128 mad(__m128 const a, __m128 const b, __m128 const c) {
    return _mm_add_ps(_mm_mul_ps(a, b), c);
}


Зачем это надо? Прежде всего, для оптимизации в будущем. Например можно просто в готовом коде заменить mad на fma через те же макросы, кому как удобно. Но основу под оптимизацию мы заложим уже сейчас:

Вариант с multiple-add
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-е float, а в AVX уже 8. То есть, есть теоретический шанс уменьшить число выполняемых операций, и увеличить производительность если не вдвое, то хотя бы раза в 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 вариантом уже не выйдет.

Если считать матрицу m в регистры ymm, мы получим ymm = {m0:m1} и ymm = {m2:m3} соответственно. Раньше у нас были только столбцы матрицы в регистре, а теперь и столбцы и строки.

Если попробовать действовать как раньше, то надо ymm={m0:m1} умножать на регистр ymm={n00,n00,n00,n00}:{n10,n10,n10,n10}. Так как n00 и n01 в одной строке матрицы n, то, судя по имеющемуся набору AVX инструкций, раскидывать их по ymm будет дорого. И shuffle, и permute работают отдельно для каждой из двух четверок float (старший и младший xmm) внутри ymm регистров.

Если считать ymm из матрицы n, то получим оба элемента n00 и n10 в старшем из 2-х xmm внутри ymm регистра. {n00,n10,n20,n30}:{n01,n11,n21,n31}. Обычно индекс у имеющихся инструкций от 0 до 3. И адресует float лишь внутри одного xmm регистра из двух внутри ymm регистра. Перекинуть n10 из старшего xmm в младший дёшево не получится. А тут ещё этот фокус надо повторить несколько раз. С такой потерей тактов мы смириться не можем. Надо придумать что-то другое.

Раньше мы обобщали столбцы, а теперь строки. Поэтому попробуем зайти немного с другой стороны. Нам нужно получить результат в {r0:r1}. Значит и улучшать алгоритм надо не по отдельным строчкам алгоритма, а сразу по две. И тут то, что было минусом в работе shuffle и permute, станет для нас плюсом. Смотрим, что у нас будет в регистрах ymm, когда мы считаем матрицу n.

n0n1 = {00, 10, 20, 30} : {01, 11, 21, 31}
n2n3 = {02, 12, 22, 32} : {03, 13, 23, 33}

Ага, замечаем, что в разных xmm частях ymm регистра у нас есть элементы 00 и 01. Их можно размножить по регистру через команду permute в {_00,_00,_00,_00}:{_01,_01,_01,_01}, указав лишь один индекс 3 для обоих xmm частей. Это именно то, что нам надо. Ведь коэффициенты тоже используются в разных строчках. Только теперь в соответствующем 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);//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);
}


Вот уже интересные цифры от IACA: x86 — 12.53, x64 — 12. Хотя, конечно, хотелось получше. Что-то упустили.

AVX оптимизация плюс «синтаксический сахар»


Похоже в коде выше AVX использовался не на полную мощь. Находим, что вместо установки двух одинаковых строк в регистр ymm можно задействовать broadcast, который умеет заполнять регистр ymm двумя одинаковыми значениями xmm. Также попутно добавим немного «синтаксического сахара» для 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+add на одну операцию. Хотя инструкцию умножения мы всё равно задействуем, чтобы дать компилятору больше возможностей для оптимизации, а процессору для параллельного выполнения умножений. Обычно я смотрю генерируемый код на ассемблере, чтобы убедиться какой вариант лучше.

В данном случае нам требуется вычислить 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. Вот теперь совсем хорошо. Наверное, кто-то скажет, что можно сделать ещё лучше, но я уже не знаю, как.

Benchmarks


Сразу отмечу, что эти цифры не стоит воспринимать как истину в последней инстанции. Даже при зафиксированном тесте они плавают в некоторых пределах. И уж тем более ведут себя по разному на разных платформах. При любой оптимизации делайте замеры конкретно для вашего случая.

Содержание таблицы


  • Function: название функции. Окончание на s — функции со стримингом, иначе обычный mov (без стриминга). Добавил для наглядности, так как это достаточно важно.
  • IACA cycles: количество тактов на функцию вычисленное IACA
  • Measured cycles: замеренное количество тактов (чем меньше, тем лучше)
  • IACA speedup: количество тактов в нулевой строке / количество тактов в строке
  • Measured speedup: количество тактов в нулевой строке / количество тактов в строке (чем больше, тем лучше)


Для loop_m такты из статьи были умножены на 64. То есть это сильно примерное значение. По факту так и оказалось.

i3-3770:


x86
Function IACA cycles Measured cycles IACA speedup Measured speedup
unroll_m 70.00 50.75 1.00 1.00
loop_m 233.60 119.21 0.30 0.43
sse_v1m 18.89 27.51 3.70 1.84
sse_v2m 19.00 27.61 3.68 1.84
sse_v3m 18.89 27.22 3.70 1.86
sse_v4s 18.89 27.18 3.70 1.87
avx_v1m 13.00 19.21 5.38 2.64
avx_v1s 13.00 20.03 5.38 2.53
avx_v2m 10.00 12.91 6.99 3.93
avx_v2s 10.00 17.34 6.99 2.93

x64
Function IACA cycles Measured cycles IACA speedup Measured speedup
unroll_m 70 68.60 1.00 1.00
loop_m 233.60 119.37 0.30 0.57
sse_v1m 18.89 21.98 3.70 3.12
sse_v2m 19.00 21.09 3.68 3.25
sse_v3m 18.89 22.19 3.70 3.09
sse_v4s 18.89 22.39 3.70 3.06
avx_v1m 13.00 9.61 5.38 7.13
avx_v1s 13.00 16.90 5.38 4.06
avx_v2m 10.00 9.20 6.99 7.45
avx_v2s 10.00 14.64 6.99 4.68

i7-8700K:


x86
Function IACA cycles Measured cycles IACA speedup Measured speedup
unroll_m 69.95 40.25 1.00 1.00
loop_m 233.60 79.49 0.30 0.51
sse_v1m 18.89 19.31 3.70 2.09
sse_v2m 19.00 19.98 3.68 2.01
sse_v3m 18.89 19.69 3.70 2.04
sse_v4s 18.89 19.67 3.70 2.05
avx_v1m 13.00 14.22 5.38 2.83
avx_v1s 13.00 14.13 5.38 2.85
avx_v2m 10.00 11.73 6.99 3.43
avx_v2s 10.00 11.81 6.99 3.41
AVX+FMAm 9.21 10.38 7.60 3.88
AVX+FMAs 9.21 10.32 7.60 3.90

x64
Function IACA cycles Measured cycles IACA speedup Measured speedup
unroll_m 69.95 57.11 1.00 1.00
loop_m 233.60 75.73 0.30 0.75
sse_v1m 18.89 15.83 3.70 3.61
sse_v2m 19.00 17.22 3.68 3.32
sse_v3m 18.89 15.92 3.70 3.59
sse_v4s 18.89 16.18 3.70 3.53
avx_v1m 13.00 7.03 5.38 8.12
avx_v1s 13.00 12.98 5.38 4.40
avx_v2m 10.00 5.40 6.99 10.57
avx_v2s 10.00 11.39 6.99 5.01
AVX+FMAm 9.21 9.73 7.60 5.87
AVX+FMAs 9.21 9.81 7.60 5.82

Код тестов в исходниках. Если есть обоснованные предложения как их улучшить, пишите в комментариях.

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 сложений.

P.S. Код из статьи




Это мой первый опыт в написании статьей. Всем большое спасибо за комментарии. Они помогают сделать код и статью лучше.
Поделиться публикацией

Комментарии 72

    +2
    Не помешала бы сводная таблица в конце статьи.
      0
      Согласен. Пока набираюсь опыта. Собираю пожелания.
      +6
      А дон знает толк в ветряных мельницах! И результат весьма впечатляет, почти как у пресловутого дона — о нем мир помнит, а мельницы сгинули.

      Можно поинтересоваться судьбой софтового рендера?
        +4
        Конечно. Пока посмотреть можно тут: www.youtube.com/watch?v=gBp1PoWzeO8
        Сначала писал статью. Потом понял, что одна не получится. Потом понял, что сразу всё писать сложно. Эта статья как раз вышла из неё. Решил посмотреть как пойдёт и будет ли интерес.
        Если кратко, то: Quake2 1600x1200, MT (multythreaded), SSE3, i7-3770, 30 fps.
        Переделываю с нуля под AVX. Цикл статьей планируется по SSE варианту.
        +1
        Минутку. А где unit test из 100500 повторений вызова функции, замером времени, и расчётом реального количества инструкций на функцию?
        Вы настолько верите IACA? Тогда вас может ждать много удивительных открытий при погружении в реальный мир.
          +5
          1. Да, верю. По опыту. Результаты как минимум коррелируют. Ну будут у меня цифры немного другие, сути это не меняет. А IACA разбирается в архитектуре Intel лучше меня уж точно.
          2. Про реальный мир вы зря: я давно уже использую SIMD и тестирую и так, и так. Тут решил обойтись.
          3. Не могу протестировать всё. У меня максимум AVX (i7-3770).
          4. Это будет синтетика. Я уже отмечал, что тысячами матрицы я не умножаю.
          5. Однако вы правы: пункт 4 имеет место быть, но IACA не может учитывать реальную работу с данными, памятью. Поэтому стоит написать. Интересно будет не только вам.
            0
            3. Не могу протестировать всё. У меня максимум AVX (i7-3770).
            Однако как раз AVX512 это особенно важно. Дело в том, что увидев AVX512-инструкцию процессор «пугается» и снижает частоту. Сразу, превентивно. Потому скорость в тактах становится куда менее важной чем для предыдущих SIMD-команд.
          +2
          > На каждый элемент мы делаем 4 умножения и 3 сложения. В сумме это 256 умножений и 192 сложения.

          В матрице 16 элементов, в сумме это 64 умножения и 48 сложений.
            0
            Спасибо!!! Как проглядел непонятно. Исправлю.
              0
              Кстати я понял почему так вышло. Я взял те самые 4*4*4 из текста, предназначенные для одного умножения во внутреннем цикле, и перенёс их на mulps, в котором 4 умножения.
              –1
              Хорошее велосипедостроительство. Но реализация libxsmm мне нравиться больше:
              Small matrix multiplication kernels (SMMs) are generated for Intel SSE, Intel AVX, Intel AVX2, IMCI (KNCni) for Intel Xeon Phi coprocessors (KNC), and Intel AVX‑512 as found in the Intel Xeon Phi processor family (KNL, KNM) and Intel Xeon processors (SKX). Highly optimized code for small convolutions is targeting Intel AVX2 and Intel AVX‑512, whereas other targets can automatically leverage specialized SMMs to perform convolutions.
                +5
                Универсальные решения как правило проигрывают специализированным а libxsmm даже не самая быстрая реализация из существующих. Если нужно постоянно работать с матрицами 4x4 а не зоопарком где фигурируют перемножения матриц 5x12 на 12x7 то использование кастомного кода выглядит разумнее. Но было бы конечно любопытно взглянуть на код который libxsmm сгенерирует для 4x4 матрицы из примера автора
                  0
                  А можно утверждение про «не самая быстрая» чем то подкрепить?
                  +2
                  Незачем тащить с собой монструозную некроссплатформенную библиотеку с нетривиальной установкой и подключением для решения мелкой проблемы, если нет в этом особой необходимости. К тому же, насколько я понял, libxsmm принуждает к использованию собственных типов для хранения матриц, что может потребовать либо переписывание кода, либо ненужное копирование данных.
                  +2

                  Интересно, удалось ли обогнать оптимизатор GCC с флагами "-m64 -O3 -march=native -funroll-loops -ffast-math"? Перемножение двух целочисленных uint32_t матриц 4096x4096 у меня заняло ~6.8 сек на одном ядре Intel Xeon E3-1220 v5, задача здесь. Алгоритм — слегка оптимизированный O(n^3), отданный на растерзание компилятору.

                    +3
                    GCC плохо/либо же ужасно умеет использовать такие инструкции. Хотя по времени выигрыш может оказаться и не очень большим.

                    Статью стоит рассматривать как мини введение в векторные инструкции, а не как конкретную реализацию умножения матриц либо еще какой конкретной задачи.
                      0
                      Согласен. Если бы компиляторы так умели, я бы точно не стал заморачиваться.
                        0
                        А у Вас нет случаем бенчмарков для сравнения Вашего кода и того, что нагенерировал GCC/Clang/ICC для этого случая?

                        Потому что мне кажется, что тот же ICC очень неплохо занимается векторизацией. Было бы интересно посмотреть на результаты бенчмарка.
                          0
                          Бенчмарки делаю, но пока не очень получается. Пока погоду показывают. Ни линукса, ни ICC (он вроде платный) у меня под руками нет. Потом когда будут бенчмарки, затестю на том, что подвернётся.
                          А пока разве только разные сайты использовать, где можно онлайн скомпилить код любым компилятором, и посмотреть результат. Но кому интересно, уже наверняка это сделали. Код то я выложил.
                            0
                            Какие именно проблемы с бенчмарками? Может что подсказать стоит?
                            ICC компилятор имеет уже давно бесплатную(возобновляемую каждые 90 дней) даже коммерческую лицензию. Просто заходите и качаете. К тому же Intel System Studio работает и на Windows, так что можете провести измерения msvc vs icc.

                            Как уже отметили, иногда сложно предсказать, какие способы будут работать быстрее, глядя просто на то, что там нагенерировал компилятор — а вот по бенчмаркам уже проще смотреть.
                              0
                              i7-3770:
                              x86: тормозят все реализации раза в 2 от ожидаемых
                              x64: все SSE быстрые, все AVX медленнее в 20-30 от ожидаемых
                              На 8700K сказали не проявляется.
                              AVX_v1 быстрее без стриминга. AVX_v2 быстрее со стримингом.
                        0
                        Справедливости ради стоит отметить, что в GCC 8.2, который вышел прям совсем недавно, улучшили генерацию кода с AVX-512
                        0
                        4096x4096 это по сути 4k экран. Тут уже SIMD + MT (multithreading). А у меня речь про чистый SIMD пока. А сама по себе задача интересная.
                        +1
                        AVX+FMA вариант (финальный)… В данном случае нам требуется вычислить a*b + c*d + e*f + g*h.

                        Выполняется одной командой _mm_dp_ps из SSE4 — она гораздо быстрее
                        0

                        У вас покрасивее, чем у Intel-а код выглядит https://software.intel.com/en-us/articles/benefits-of-intel-avx-for-small-matrices.
                        А почему вы не используете __restrict__ и __builtin_assume_aligned/__assume_aligned? Может быть тогда и инлайниться лучше будет.

                          0
                          В моём рабочем коде классов классе пишу вроде
                          __forceinline __m128 __vectorcall operator + (__m128 const) const noexcept;

                          Однако здесь всё итак работает хорошо. Всегда проверяю сгенерированный ассемблерный код. Насчёт *assume_aligned я посмотрю что это. Скорее всего то, чего мне сильно и давно не хватает. Но если оно linux/intel_compiler/… specific, то мне, к сожалению, не подойдёт.
                            0
                            Посмотрел. Да, это что мне нужно. Увы, ICC specific.
                              0

                              __assume_aligned — это Intel, __builtin_assume_aligned — это GCC.

                            0
                            а почему бы не использовать DPPS _mm_dp_ps (SSE4.1/AVX)? готовое скалярное произведение пары векторов
                              +3
                              Очень разумный вопрос. Статья по факту представляет чистый положительный опыт. Если собрать весь отрицательный, то она будет минимум в 3 раза больше. Конечно dpps я пробовал когда то. Не оправдала от слова совсем.

                              Когда SSE4.1 только рекламировали, я уже надеялся, что такая классная инструкция dpps всё упростит. И ждал. Но увы…

                              1. Она дорогая: 2 такта (и судя по всему не параллелится как mulps), то есть даже без учёта транспонирования, их надо 16 по 2 такта, получается 32!!! (IACA согласна)
                              1.1 dpps — 4 умножения, 3 сложения
                              1.2 mulps — 4 умножения, параллелится, это 8 умножений за такт, 16 за два такта, или, грубо, 8 умножений + 3 сложения за 2.5 такта
                              2. Одну матрицу надо предварительно трансформировать: dpps(строка M, столбец N)
                              3.1 дополнительная инструкция movss для извлечения результата и работа с памятью на уровне 1-го float
                              ИЛИ
                              3.2 дополнительный код для компоновки результата из 4 float в одном xmm

                              Не стоит думать об инструкции, как вещи в себе. Инструкции нужны данные, они должны быть в нужном представлении как на входе, так и на выходе.

                              По совокупности факторов её даже Intel не использует в своих примерах матричного умножения.
                                0
                                … Очень разумный вопрос. Статья по факту представляет чистый положительный опыт. Если собрать весь отрицательный, то она будет минимум в 3 раза больше…
                                А может будет это вторая статья?

                                По моему опыту, dpps имеет приличный выигрыш, но случаи бывают всякие. И подготовка данных, полностью согласен, стоит времени ЦПУ
                                  0
                                  Как вариант приложить в коде. Но на статью, пожалуй, это не тянет. Простой набросок без транспонирования дал уже 32 такта.
                                0
                                Я уверен, что всё равно найдутся несогласные. Но для меня это пройденный опыт, и я не хочу на нём останавливаться. Если вы считаете, что можете лучше с dpps, киньте сюда свою реализацию вычисления умножения с ней, которая будет лучше хотя бы 20 тактов (с учётом чтения матриц m/n и записью в результирующую матрицу r), и я публично признаю свою ошибку.
                                +1
                                сравнивать число тактов в AVX-реализации с числом тактов в SSE-реализации некорректно.
                                Процессоры Intel снижают частоту ниже базовой при большом числе AVX инструкций.
                                https://en.wikichip.org/wiki/intel/frequency_behavior
                                  0
                                  Спасибо за информацию. Тем более прав был old_bear, когда настаивал на замерах времени: habr.com/post/418247/#comment_18921601.
                                    +1
                                    Спасибо за информацию. Погуглил, пишут что снижение частоты зависит от степени нагрева процессора. При равной тепловой нагрузке AVX-код будет выполняться на частоте до ~20% ниже (т.е. 8 реальных тактов автора превратятся в 10 «эквивалентных»), это все равно намного быстрее чем 16 в SSE-варианте.
                                      +1
                                      Подстава заключается в том, что частота падает не только во время выполнения этих инструкций, но и некоторое время после. В реальности отдельная функция может начать работать быстрее, но программа в целом замедлится. Этот эффект наблюдается с новым OpenSSL — по всем бенчмаркам самой библиотеки стало быстрее, однако веб-сервер начинает обслуживать https запросы медленнее.
                                        +1
                                        Это вы ещё AVX512 не видели. Там вообще нет смысла в использовании, если приложение ускоряется меньше чем на 20% по тактам, т.к. частота на эти же 20% может сбрасываться. Относительно версии с AVX/AVX2.
                                      0
                                      Дополнительно ещё стоит учесть работу с памятью.
                                      Например, разницы при сложении двух очень длинных векторов с использованием SSE и AVX при максимальном распараллеливании не будет, т.к. всё упрётся в пропускную способность памяти.
                                        0
                                        А AMD не снижают. Стало быть корректно?
                                        Зачем быть заложником одной компании?
                                          +1
                                          У современных процессоров AMD в два раза меньше вычислительных блоков AVX, чем у процессоров Intel. Так что AMD пока в проигрыше.
                                            0
                                            Смотря каких.
                                            www.agner.org/optimize/blog/read.php?i=838

                                            AMD has four 128-bit units for floating point and vector operations. Two of these can do addition and two can do multiplication. Intel has two 256-bit units, both of which can do addition as well as multiplication. This means that floating point code with scalars or vectors of up to 128 bits will execute on the AMD processor at a maximum rate of four instructions per clock (two additions and two multiplications), while the Intel processor can do only two. With 256-bit vectors, AMD and Intel can both do two instructions per clock.

                                            Т.е. тоже самое.
                                            Но FMA быстрее на Интел.

                                            Intel beats AMD on 256-bit fused multiply-and-add instructions, where AMD can do one while Intel can do two per clock. Intel is also better than AMD on 256-bit memory writes, where Intel has one 256-bit write port while the AMD processor has one 128-bit write port.
                                        0

                                        Очень интересная статья! Кстати, кроме как тем, что нужен свой формат хранения, матрички из GLM не такие же шустрые? Кажется, там была какая-то поддержка SIMD.

                                        0
                                        А можно такой вопрос — от неспециалиста в SIMD инструкциях.
                                        По поводу разрядности. Я так понимаю что float в данном случае 32 бита?
                                        А если точность нужна выше? 64 бита на Double. Этот же код сможет работать с double — данные влезут в регистры? Или надо будет модифицировать?
                                        И возможен ли вариант на 80 бит (extended)? Или SIMD не поддерживает 80-битный float сопроцессора x86?
                                          +1
                                          1 регистр SSE — это 4 float или 2 double
                                          1 регистр AVX — это 8 float или 4 double
                                          Код просто будет аналогичным.

                                          80-битные числа SSE/AVX не поддерживаются
                                            0
                                            Благодарю!
                                            А вот насчет 80 бит жаль :-( Странно что Intel не дала возможности использовать его в своих же расширениях. Пусть бы и не так быстро как float32.
                                              0
                                              Странно что Intel не дала возможности использовать его в своих же расширениях.
                                              Нет в этом ничего странного. Операции с 80битами сильно дороже, чем с 64мя. А используются редко. Тратить драгоценные транзисторы на то, что не используется никто не будет.

                                              Собственно изначально MMX работал с целыми числами, SSE — с 32-битными float'ами, SSE2 — с 64битными double'ами. И каждый раз собиралась статистика: кто будет использовать эти команды и зачем. Очевидно, что когда «выпекали» SSE3/4/etc толп людей, желающих воспользоваться 90битными флоатами у маркетиногово отдела Intel — не стояло…
                                              0
                                              Код просто будет аналогичным.

                                              Если бы! AVX — это не в два раза больший SSE модуль, это два SSE модуля, работающих в паре. Если для сложения и умножения это ничего не меняет, то например _mm256_shuffle_ps ведёт себя совсем иначе.

                                                0
                                                На AVX2 есть _mm256_permute4x64_pd, на AVX действительно придётся допиливать напильником. Но в целом, признаюсь, я далёк от double тематики.
                                                  0
                                                  это не в два раза больший SSE модуль

                                                  В 2 раза больший.
                                                  это два SSE модуля, работающих в паре.

                                                  Амд-проблемы.

                                                    0
                                                    Я бы не назвал это прям уж проблемой. Иногда ymm нужен именно как спаренный sse (для float индекс будет в диапазоне 0-3), а иногда как один большой регистр с независимыми числами (в этом случае для float уже нужен индекс в диапазоне 0-7). Проблема, что они вторую часть реализуют в следующем наборе, или не реализуют вообще. Как пример, выше упомянутая _mm256_permute4x64_pd.
                                                      +1

                                                      Вы вообще не поняли о чем я, видимо ни разу не использовали AVX. Я не про какие-то реализации каких-то процессоров, я про весь набор команд и их архитектуру. Посмотрите как работает например _mm256_shuffle_ps — в dst[127:0] могут попасть только биты из a[127:0] и b[127:0], но не из a[255:128] и b[255:128]. То есть нижние 128 бит работают только с нижними 128 битами, старшие со старшими — это больше похоже на два модуля SSE, чем на один модуль SSE, но в два раза больше. И это касается большинства команд. То же самое можно сказать про AVX512 — это скорее 4 модуля SSE, работающих параллельно, чем один очень длинный модуль SSE.

                                                        –2
                                                        Посмотрите как работает например _mm256_shuffle_ps — в dst[127:0] могут попасть только биты из a[127:0] и b[127:0], но не из a[255:128] и b[255:128]. То есть нижние 128 бит работают только с нижними 128 битами, старшие со старшими — это больше похоже на два модуля SSE, чем на один модуль SSE, но в два раза больше. И это касается большинства команд. То же самое можно сказать про AVX512 — это скорее 4 модуля SSE, работающих параллельно, чем один очень длинный модуль SSE.

                                                        Похоже вы ничего не знаете об микроархитектуре, странно как-то вы используете avx и ничего о ней не знаете.

                                                        Во-первых, что значит «например» — если это единственный пример(вернее там много похожих инструкций, но суть одна). Зачем вы делаете вид, будто бы у вас много примеров?

                                                        Во-вторых ваша аналитика ничего не имеет общего с реальность. И если бы вы действительно использовали тот самых avx, то вы бы знали в чём именно заключается проблема, а заключается проблема в 8битном селекторе.

                                                        Опять же, ваши рассуждения имели бы смысл, если данные операции работали в рамках 128 бит. Т.е. если бы из А можно было перемещать в рамках первой части dst, а из Б в рамках второй. Это было бы так, если бы они попросту склеили два 128битных шафла, но это не так. Они почему-то шафлят aa bb aa bb. Неужели там mmx?

                                                        И да, никаких «модулей» не существует. Это обывательская терминология, ничего общего с реальностью не имеющая.

                                                          0

                                                          Нет. Видимо, вы путаете AVX и AVX2. В AVX все команды работают с регистрами как с двумя 128-битными независимыми половинками (за исключением пары команд, которые обращаются с 128-битными частями как единым целым).


                                                          И если бы вы действительно использовали тот самых avx, то вы бы знали в чём именно заключается проблема, а заключается проблема в 8битном селекторе.

                                                          Когда 8-битного селектора не хватает, то в качестве аргумента команды может использоваться другой регистр, например, в команде AVX2 _mm256_permutevar8x32_ps.


                                                          Они почему-то шафлят aa bb aa bb

                                                          Потому что lower(result) = op(lower(aa), lower(bb)), higher(result) = op(higher(aa), higher(bb)).


                                                          SELECT4(src, control)
                                                          {
                                                            CASE(control[1:0])
                                                              0: tmp[31:0] := src[31:0]
                                                              1: tmp[31:0] := src[63:32]
                                                              2: tmp[31:0] := src[95:64]
                                                              3: tmp[31:0] := src[127:96]
                                                            ESAC RETURN tmp[31:0]
                                                          }
                                                          
                                                          dst[31:0] := SELECT4(a[127:0], imm8[1:0])
                                                          dst[63:32] := SELECT4(a[127:0], imm8[3:2])
                                                          dst[95:64] := SELECT4(b[127:0], imm8[5:4])
                                                          dst[127:96] := SELECT4(b[127:0], imm8[7:6])
                                                          dst[159:128] := SELECT4(a[255:128], imm8[1:0])
                                                          dst[191:160] := SELECT4(a[255:128], imm8[3:2])
                                                          dst[223:192] := SELECT4(b[255:128], imm8[5:4])
                                                          dst[255:224] := SELECT4(b[255:128], imm8[7:6])
                                                          dst[MAX:256] := 0
                                                            0
                                                            Во-первых, что значит «например» — если это единственный пример(вернее там много похожих инструкций, но суть одна). Зачем вы делаете вид, будто бы у вас много примеров?

                                                            Ооох, да пожалуйста:


                                                            __m256 _mm256_unpacklo_ps (__m256 a, __m256 b)
                                                            
                                                            INTERLEAVE_DWORDS(src1[127:0], src2[127:0]){
                                                                dst[31:0] := src1[31:0] 
                                                                dst[63:32] := src2[31:0] 
                                                                dst[95:64] := src1[63:32] 
                                                                dst[127:96] := src2[63:32] 
                                                                RETURN dst[127:0]
                                                            }   
                                                            
                                                            dst[127:0] := INTERLEAVE_DWORDS(a[127:0], b[127:0])
                                                            dst[255:128] := INTERLEAVE_DWORDS(a[255:128], b[255:128])
                                                            dst[MAX:256] := 0

                                                            __m256i _mm256_packs_epi32 (__m256i a, __m256i b)
                                                            
                                                            dst[15:0] := Saturate_Int32_To_Int16 (a[31:0])
                                                            dst[31:16] := Saturate_Int32_To_Int16 (a[63:32])
                                                            dst[47:32] := Saturate_Int32_To_Int16 (a[95:64])
                                                            dst[63:48] := Saturate_Int32_To_Int16 (a[127:96])
                                                            dst[79:64] := Saturate_Int32_To_Int16 (b[31:0])
                                                            dst[95:80] := Saturate_Int32_To_Int16 (b[63:32])
                                                            dst[111:96] := Saturate_Int32_To_Int16 (b[95:64])
                                                            dst[127:112] := Saturate_Int32_To_Int16 (b[127:96])
                                                            dst[143:128] := Saturate_Int32_To_Int16 (a[159:128])
                                                            dst[159:144] := Saturate_Int32_To_Int16 (a[191:160])
                                                            dst[175:160] := Saturate_Int32_To_Int16 (a[223:192])
                                                            dst[191:176] := Saturate_Int32_To_Int16 (a[255:224])
                                                            dst[207:192] := Saturate_Int32_To_Int16 (b[159:128])
                                                            dst[223:208] := Saturate_Int32_To_Int16 (b[191:160])
                                                            dst[239:224] := Saturate_Int32_To_Int16 (b[223:192])
                                                            dst[255:240] := Saturate_Int32_To_Int16 (b[255:224])
                                                            dst[MAX:256] := 0

                                                            __m256 _mm256_dp_ps (__m256 a, __m256 b, const int imm8)
                                                            
                                                            DP(a[127:0], b[127:0], imm8[7:0]) {
                                                                FOR j := 0 to 3
                                                                    i := j*32
                                                                    IF imm8[(4+j)%8]
                                                                        temp[i+31:i] := a[i+31:i] * b[i+31:i]
                                                                    ELSE
                                                                        temp[i+31:i] := 0
                                                                    FI
                                                                ENDFOR
                                                            
                                                                sum[31:0] := (temp[127:96] + temp[95:64]) + (temp[63:32] + temp[31:0])
                                                            
                                                                FOR j := 0 to 3
                                                                    i := j*32
                                                                    IF imm8[j%8]
                                                                        tmpdst[i+31:i] := sum[31:0]
                                                                    ELSE
                                                                        tmpdst[i+31:i] := 0
                                                                    FI
                                                                ENDFOR
                                                                RETURN tmpdst[127:0]
                                                            }
                                                            
                                                            dst[127:0] := DP(a[127:0], b[127:0], imm8[7:0])
                                                            dst[255:128] := DP(a[255:128], b[255:128], imm8[7:0])
                                                            dst[MAX:256] := 0

                                                            И в все остальные команды арифметических и логических действий соблюдают это правило: результат работы с нижними 128 битами записывается в нижние 128 битов. А теперь ваша очередь показать все остальные инструкции спокойно перекидываются битами из нижних и верхних частей? Такие команды конечно есть, но работают они со всеми 128 битами одновременно, то есть являются инструкциями копирования, а никак не SIMD-операциями.

                                                              +1
                                                              В основном да. Но исключения всё-таки есть. Та же _mm256_permute4x64_pd.

                                                              SELECT4(src, control){
                                                              	CASE(control[1:0])
                                                              	0:	tmp[63:0] := src[63:0]
                                                              	1:	tmp[63:0] := src[127:64]
                                                              	2:	tmp[63:0] := src[191:128]
                                                              	3:	tmp[63:0] := src[255:192]
                                                              	ESAC
                                                              	RETURN tmp[63:0]
                                                              }
                                                              
                                                              dst[63:0] := SELECT4(a[255:0], imm8[1:0])
                                                              dst[127:64] := SELECT4(a[255:0], imm8[3:2])
                                                              dst[191:128] := SELECT4(a[255:0], imm8[5:4])
                                                              dst[255:192] := SELECT4(a[255:0], imm8[7:6])
                                                              dst[MAX:256] := 0
                                                                0
                                                                Вот только это команда AVX2, а не AVX.
                                                                  0
                                                                  Увы, да. Однако это проблема тех, кому нужен double.
                                                                    0
                                                                    Почему же? Есть аналогичная команда для float _mm256_permutevar8x32_ps.

                                                                    Это проблема тех, кто до сих пор пользуется Sandy Bridge, в котором нет AVX2 (например, я).
                                                    +1
                                                    При миграции с SSE на AVX (повышение разрядности инструкций), сможет с минимальными изменениями. В SSE регистр вмещалось 4 float, в AVX 4 double. Если на той же разрядности регистров нужен double, тогда придётся переписывать.

                                                    80 только FPU. SSE: 32*4=64*2=128bit, AVX: 64*4=32*8=256bit.
                                                    По сути точность несколько меньше, плата за производительность.
                                                      0
                                                      Благодарю, а такой вопрос — если код написан под AVX — а процессор поддерживает только SSE, то в этом случае надо писать два варианта кода для поддержки обоих архитектур или есть возможность рассчитывать на эмуляцию с какой-нибудь стороны?
                                                        0
                                                        Для эмуляции есть способ. Где то находил в интернетах SDK для эмуляции. С практической точки зрения SDK интересно лишь для отладки.

                                                        Здесь вскользь упомянут вариант, как делать для разных платформ. Пишем разные реализации, и потом после анализа процессора через cpuid подставляем нужный указателю на функцию. И потом работаем с данным функционалом через указатель. Данный подход полезен для больших, ёмких операций. Если использовать его для одной инструкции — точно будет только хуже.

                                                        Как пример: функция для умножения матрицы на вектор — станет медленнее, функция для умножения матрицы на массив векторов — имеет смысл.

                                                        Гибкий вариант есть здесь optimizing_cpp.pdf. Искать по: CriticalFunction_Dispatch.

                                                          0
                                                          Хорошо, буду думать, а то тоже хотелось ускорить умножение матриц 4*4, но немного — не более 2-5 десятков. Но в идеале extended на 80 бит. Значит или многопоточный вариант с FPU или переходить на double для float64.
                                                    +3
                                                    Разные руководства по оптимизации рекомендуют лишний раз не дёргать кэш для массовых операций сохранения.

                                                    Это явно не ваш случай. Эти инструкции полезны, когда вам нужно обработать огромный массив данных (по крайней мере больше кеша) и вы уверены, что к нему не будет обращений до конца вашей работы. В вашем же случае вы оперируете 64 байтами, к которым с большой вероятностью в ближайшее время будет дальнейшая работа.


                                                    По опыту использования в обработке изображений стриминг давал незначительный прирост на больших изображениях и проигрыш в несколько раз на изображениях поменьше.

                                                      0
                                                      Про это было отмечено в статье, что большие массивы данных нужны скорее на умножении матрицы на массив вершин.
                                                        0

                                                        Не понял к чему вы это. Финальный код содержит стриминговые инструкции, которые в данном случае сильно вредны. Ваше нежелание или невозможность тестировать на реальных кейсах подводит вас.

                                                          0
                                                          В чём то вы оказались правы. Но как то странно получилось. Stream никак не сказалось на SSE версии, зато сильно портила жизнь AVX1 (ускорение 6-7 против 8-9). Архитектура IvyBridge. В целом корректно было бы IVB такты измерить. А это лишь IACA 2 может. Но они как то по разному мерят вторая с третьей.
                                                            0
                                                            А теперь нет. Видимо над тестом придётся ещё поработать.
                                                              0

                                                              У вас данные выровнены по линиям кеша? 4x4x4 — это как раз одна линия. Возможно, если данные пересекают линии, результат будет отличаться.

                                                    Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                                                    Самое читаемое