Обзор инструкций ARM NEON для тех, кто знаком с MMX/SSE/AVX

    Мир изменился. Я чувствую это в воде, чувствую это в земле, ощущаю в воздухе.

    «Властелин колец», Джон Рональд Руэл Толкин

    Архитектура x86 долгие десятилетия была лидером по высокопроизводительным решениям. И этот факт позволял ей доминировать даже когда количество устройств на архитектуре ARM вокруг нас стало в несколько раз больше. Если вы пишете высокопроизводительный софт для серверов или рабочих станций, до недавнего времени вы могли обходиться лишь знанием x86. Но невозможно игнорировать события, которые могут поменять расклад сил уже совсем скоро. В связи с этим я решил немного поиграться с ARM и нарыл интересных фактов немного больше, чем на тред в твитере.

    У меня нет цели рассказать всё с самого начала, я буду заострять внимание лишь на отличиях и интересных моментах, которые мне встретились. Подразумевается, что вы знаете, что такое инструкции процессора, SIMD, читаете ассемблер и Си на базовом уровне.

    Мой опыт написания высокопроизводительного кода в основном связан с обработкой изображений в библиотеке Pillow-SIMD. Там я использовал интринсики в коде на Си чтобы добиться 6-8-кратного ускорения наиболее частых операций.

    Под что вообще пишем?

    Честно говоря, это самый высокий порог для вхождения в ARM архитектуру из тех, что будут. В x86 есть базовый набор команд, есть расширения (разные версии SSE, AVX, криптография или виртуализация) и есть разрядность (32 или 64 бита). В ARM же, загибайте пальцы:

    • Есть архитектуры, коих более 20. Называются они примерно так: ARMv7-M, ARMv8-R, ARMv8.3-A.

    • Есть микроархитектура. Например: Cortex-R4, Cortex-A76.

    • Есть профайл: Classic, Microcontroller, Real-time, Application.

    • Есть разные наборы команд! A32, A64, Thumb, Thumb2.

    • Наконец, расширения набора команд: SIMD, NEON, SVE.

    • Ну и никуда не делась разрядность: AArch32 и AArch64.

    Не претендуя на полноту описания, я подсвечу основные моменты и укажу, что сейчас можно опустить. Все архитектуры соответствуют одному из четырех профайлов. Classic — это прям совсем классик, такое вы вряд ли встретите. Из трёх остальных самое ходовое — это Application. Все телефоны, сервера и рабочие станции это Application. Профайл всегда отражён в названии архитектуры в виде постфикса (A, M, R). Актуальных архитектур всего две — ARMv7 и ARMv8, зато у ARMv8 вышло уже 6 минорных версий, которые тоже называются архитектурами (например, ARMv8.2-A). Причём 64-битная разрядность появилась только в ARMv8. Однако ARMv8-A не гарантирует наличие 64-битного режима у процессора, а вот ARMv8.1-A уже гарантирует.

    Если знание архитектуры чипа нужно нам, разработчикам софта, чтобы знать, какой минимальный набор функциональности возможно использовать, то микроархитектура уже нужна для разработчиков чипов, чтобы знать, сколько кеша нужно насыпать, сколько вычислительных блоков должно быть и какие опциональные технологи нужно включить в чип. Причем, бывает как микроархитектура от самой компании ARM (она обычно называется Cortex и следом снова постфикс профайла), так и кастомная, которая может называться Apple Firestorm, Neoverse N1 или никак не называться.

    Набор команд A32 используется в 32-битном режиме AArch32, а A64 в 64-битном режиме AArch64. И казалось бы, зачем выделять такие очевидные вещи. Но дело в том, что A32 — не единственный набор команд, который может быть в 32-битном режиме. A32 и A64 всегда используют 32 бита для кодирования любой инструкции, а AArch32 вышел очень давно и многим казалось, что это расточительство и тогда появились альтернативные способы кодирования Thumb и Thumb2. В них часто используемые инструкции занимали 16 бит. Для AArch64 уже ничего такого не завезли, в нём любая инструкция занимает 32 бита.

    Ну и наконец, расширения набора команд. Вообще, есть ещё расширения VFPv1-VFPv5 для работы с плавающей точкой, разницу между которыми я так и не смог понять. Как и в x86, в ARM плавающую точку завезли не сразу. В ARMv6 было добавлено расширение SIMD (так и называется), а в ARMv7 появился опциональный 128-битный advanced SIMD, он же ASIMD, он же NEON, по сути прямой аналог SSE последних версий. О нём я буду рассказывать больше всего. А вот аналога AVX в ARM нет, там пошли другим путём. Вместо того, чтобы каждые пять лет представлять новое расширение, под которые нужно будет всё переписывать, было разработано расширение Scalable Vector Extension (SVE), которое позволяет выполнять один и тот же код на чипах, реализующих разный размер векторов. Но на практике, как я понял, SVE реализован только в Fugaku supercomputer.

    Это же ужас?

    Ну, вообще, да, если вы собрались писать приложение, которое может быть выполнено на любом ARM процессоре, как это бывает с x86. Теоретически на нем может не оказаться не только NEON, но даже 64-битной арифметики с плавающей точкой. Вот только, к счастью, у ARM нет того наследия работающих систем, на которых могли бы запустить ваш код. Это в любом случае будет свежий процессор. И ещё, с максимальной вероятностью это всё же будет AArch64 система. А теперь следите за руками.

    AArch64 появился только в ARMv8. ARMv8-A уже гарантирует наличие VFPv4 (64-битный FPU), NEON и криптографии. А SVE можно даже не проверять ещё пару лет. У NEON никаких версий нет. Так же остается только один набор инструкций: A64. А микроархитектура просто ни на что не влияет.

    Получается, несмотря на огромное количество вариантов, в реальности писать код под ARM (точнее под AArch64) даже проще, чем под x86. Никакие проверки в рантайме не нужны, просто ставите `#ifdef __aarch64__` и пользуетесь всем, чем хотите.

    Знакомство с NEON

    Принципиальное устройство x86 и ARM мало чем отличаются. И там и там есть общие регистры и регистры для вычислений с плавающей точкой и SIMD. Для общей картины достаточно прочитать раздел General-purpose registers в AArch64 Instruction Set Architecture. И сразу после этого можно переходить к самому полному вводному гайду — Coding for Neon.

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

    LD3 { V0.16B, V1.16B, V2.16B }, [x0]

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

    Очень приятно удивляет кол-во вариантов всяческих сдвигов. Например, есть вариант, когда сдвиг каждого элемента вектора задается в другом векторе.

    В SSE, например, такого нет. Похожая функциональность появилась только в AVX2 с инструкциями vpsrlv[dq], vpsllv[dq], vpsravd. Но, во-первых, в NEON инструкции сдвигают в обе стороны, в зависимости от знака. Во-вторых, в AVX2 можно сдвинуть только 32-битные и 64-битные значения. На этом вкусности NEON не заканчиваются, из других вариантов сдвига есть:

    • Сдвиг и сложение с аккумулятором

    • Сдвиг вправо с округлением

    • Сдвиг с сатурацией

    • Сдвиг с уменьшением или увеличением разрядности

    Причем некоторые варианты из списка можно комбинировать, и всё это работает с любым типом данных. Поправьте, но вроде SSE ничего из этого не может предложить.

    Что касается интринсиков, в отличие от SSE/AVX, где типизированны только регистры для float и double (__m128 и __m128d), в NEON есть типы для всех целых типов и названия придерживаются конвенции stdint.h.

    uint8x16_t Sx4 = vld1q_u8(&Srgba[i]);
    uint8x16_t Dx4 = vld1q_u8(&Drgba[i]);
    uint32x4_t Sax4 = vshrq_n_u32((uint32x4_t) Sx4, 24);
    uint32x4_t Dax4 = vshrq_n_u32((uint32x4_t) Dx4, 24);

    Тут первые две переменные имеют тип «16 беззнаковых int8», вторые — «4 беззнаковых int32». Но, так как это одни и те же регистры, их можно приводить друг к другу. Интересно, что есть типы вроде uint8x16x3_t — это три регистра подряд. В основном такие типы используются для загрузки и сохранения в оперативную память.

    Справочник интринсиков

    Если вы за последнее десятилетие писали SIMD-код для x86, вы наверняка пользовались Intel Intrinsics Guide. Это прекрасный справочник с интерактивным поиском и фильтром, понятным описанием и псевдокодом для каждой инструкции. И даже есть таблицы задержек и пропускной способности по разным поколениям процессоров Intel.

    В нём можно не только искать нужное, но и изучать новое, просто выбирая какое-то поколение инструкций и читая всё подряд.

    У ARM аналогом этого гайда служит Neon Intrinsics Reference. И это просто боль и унижение.

    • Нет никаких фильтров

    • Поиск работает с перезагрузкой страницы

    • На странице выводится только 30 функций, снизу есть постраничная навигация, тоже с перезагрузкой страницы

    • Первый экран занимают бесполезные три абзаца текста, которые нужно каждый раз проматывать. А проматывать приходится очень часто, ведь любой чих — перезагрузка страницы

    • По запросу "mul" находятся 50 страниц функций! То есть 1500 штук. Знаете, почему в результатах оказалась функция, показанная на скриншоте? Потому что в описании есть слово accumulate!

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

    • Вы вообще видели этот псевдокод? Он сам по себе очень избыточен и запутан.

    Поэтому, пока я для себя нашел такой вариант: сначала я ищу что-то похожее на то, что мне нужно в алфавитном указателе инструкции с кратким описанием, а уже потом по названию инструкции ищу в справочнике интринсики.

    Подходящая задача

    Давайте попробуем NEON в деле. В качестве примера кода, на котором можно поэкспериментировать, я выбрал альфа-композитинг с premultiplied alpha. Его можно описать такой формулой:

    R_{rgba} = S_{rgba} + D_{rgba} × (1 − S_a)

    Алгоритм идеально ложится на SIMD:

    #include <stdint.h>
    #include <stddef.h>
    
    #define SHIFTFORDIV255(a)\
        ((((a) >> 8) + a) >> 8)
    
    #define DIV255(a)\
        SHIFTFORDIV255(a + 0x80)
    
    static void
    opSourceOver_premul(uint8_t* restrict Rrgba,
                        const uint8_t* restrict Srgba,
                        const uint8_t* restrict Drgba, size_t len)
    {
        size_t i = 0;
        for (; i < len*4; i += 4) {
            uint8_t Sa = Srgba[i + 3];
            Rrgba[i + 0] = DIV255(Srgba[i + 0] * 255 + Drgba[i + 0] * (255 - Sa));
            Rrgba[i + 1] = DIV255(Srgba[i + 1] * 255 + Drgba[i + 1] * (255 - Sa));
            Rrgba[i + 2] = DIV255(Srgba[i + 2] * 255 + Drgba[i + 2] * (255 - Sa));
            Rrgba[i + 3] = DIV255(Srgba[i + 3] * 255 + Drgba[i + 3] * (255 - Sa));
        }
    }

    Все операции выполняются над целочисленными значениями. 1 - Sa превращается в 255 - Sa. Так как правую часть мы умножаем на байт, то и левую нужно умножить. А после сложения нужно поделить на 255, это делает макрос DIV255 путём нескольких сдвигов.

    Запускать я буду на Raspberry Pi 4, естественно под AArch64. Чем богаты, тем и рады, как говорится. Причем в данном случае мне интересно посмотреть именно пиковую производительность, без влияния памяти. Для этого я буду тестировать на строке длиной 1000 пикселей, то есть всего будет задействовано 12 Кб данных за один вызов функции.

    Я буду пользоваться компилятором Clang-9, т.к. он в большинстве случаев выдает более быстрый код, чем GCC. Для начала интересно, как быстро работает чистый код, без векторизации.

    $ clang-9 -Wall -O2 -o run.64 main.c -fno-tree-vectorize && ./run.64
    Time elapsed: 0.189449
    Time elapsed: 0.189280
    Time elapsed: 0.189272
    Time elapsed: 0.189272

    Время указано в секундах для 20 тысяч прогонов функции с длиной строк 1000 пикселей. То есть можно сказать, что скорость работы примерно 105 МПх/с. И вообще-то это очень мало, даже для Raspberry Pi. Если включить автоматическую векторизацию, результат будет чуть лучше.

    $ clang-9 -Wall -O2 -o run.64 main.c -ftree-vectorize && ./run.64
    Time elapsed: 0.082168
    Time elapsed: 0.082341
    Time elapsed: 0.082135
    Time elapsed: 0.082147

    Ускорение в 2.3 раза существенно, но это не всё, на что можно было бы рассчитывать. Посмотрим, что можно сделать вручную.

    NEON-версия

    Первое, что нужно сделать, чтобы начать программировать на NEON — подключить заголовочный файл arm_neon.h. Согласитесь, это приятнее, чем каждый раз гуглить ничего не значащие имена заголовочных файлов, вроде smmintrin.h.

    1. Загрузка/сохранение. Хотя видно, что код хорошо векторизуется внутри цикла, можно сразу пойти чуть дальше и читать из памяти 128-битный вектор целиком и работать с четырьмя пикселями.

    #include <stdint.h>
    #include <stddef.h>
    #include <arm_neon.h>
    
    static void
    opSourceOver_premul(uint8_t* restrict Rrgba,
                        const uint8_t* restrict Srgba,
                        const uint8_t* restrict Drgba, size_t len)
    {
        size_t i = 0;
        for (; i < len*4 - 12; i += 16) {
            uint8x16_t Sx4 = vld1q_u8(&Srgba[i]);
            uint8x16_t Dx4 = vld1q_u8(&Drgba[i]);
            uint8x16_t Rx4 = vaddq_u8(Sx4, Dx4);  // Temporary stub
            vst1q_u8(&Rrgba[i], Rx4);
        }
        for (; i < len*4; i += 4) {
            uint8_t Sa = Srgba[i + 3];
            Rrgba[i + 0] = DIV255(Srgba[i + 0] * 255 + Drgba[i + 0] * (255 - Sa));
            Rrgba[i + 1] = DIV255(Srgba[i + 1] * 255 + Drgba[i + 1] * (255 - Sa));
            Rrgba[i + 2] = DIV255(Srgba[i + 2] * 255 + Drgba[i + 2] * (255 - Sa));
            Rrgba[i + 3] = DIV255(Srgba[i + 3] * 255 + Drgba[i + 3] * (255 - Sa));
        }
    }

    Тут пока что Rx4 считается намеренно неправильно. Но зато код запускается и уже можно прикинуть, сколько работает код на NEON, если он ничего не делает:

    $ clang-9 -Wall -O2 -o run.64 main.c && ./run.64
    Time elapsed: 0.008030
    Time elapsed: 0.007872
    Time elapsed: 0.007859
    Time elapsed: 0.008629

    8 мс или 2500 МПх/с! Вот мы и ускорили код с помощью NEON в 25 раз.

    2. Выделение альфа-канала. Следующая задача — нужно из вектора Sx4 отдельно вытащить все компоненты альфа-канала. Они находятся на 3, 7, 11 и 15 позиции. Причем желательно их сразу размножить на все соседние байты этого же пикселя. Несмотря на множество команд перемешивания байтов, я не нашел ничего специального и сделал через векторный поиск в таблице. Дальше нужно вычесть полученное значение из 255.

    uint8x16_t vsubq_u8 (uint8x16_t a, uint8x16_t b);
    uint8x16_t vdupq_n_u8 (uint8_t value);
    uint8x16_t vqtbl1q_u8 (uint8x16_t t, uint8x16_t idx);
    
    uint8x16_t Sax4 = vsubq_u8(
        vdupq_n_u8(255),
        vqtbl1q_u8(Sx4, (uint8x16_t){3,3,3,3, 7,7,7,7, 11,11,11,11, 15,15,15,15})
    );

    Интересно, что все компиляторы при оптимизации заменяют операцию вычитания из 255 на побитовое отрицание, что логично.

    3. Умножение. Дальше нужно все элементы Sx4 умножить на 255, а элементы Dx4 на соответствующие элементы альфы из Sax4. Все значения 8-битные.

    В NEON есть два вида умножения: либо это обычные функции vmulq_*, которые не меняют разрядность и отдают нижнюю часть результата, либо это vmull_* и vmull_high_*, которые делают операцию только над половиной вектора, но зато увеличивают разрядность и отдают результат целиком.

    uint16x8_t vmull_u8 (uint8x8_t a, uint8x8_t b);
    uint8x8_t vget_low_u8 (uint8x16_t a);
    uint8x8_t vdup_n_u8 (uint8_t value);
    uint16x8_t vmull_high_u8 (uint8x16_t a, uint8x16_t b);
    
    uint16x8_t Rx2lo = vmull_u8(vget_low_u8(Sx4), vdup_n_u8(255));
    uint16x8_t Rx2hi = vmull_high_u8(Sx4, vdupq_n_u8(255));

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

    uint16x8_t vmlal_u8 (uint16x8_t a, uint8x8_t b, uint8x8_t c);
    uint16x8_t vmlal_high_u8 (uint16x8_t a, uint8x16_t b, uint8x16_t c);
    
    Rx2lo = vmlal_u8(Rx2lo, vget_low_u8(Dx4), vget_low_u8(Sax4));
    Rx2hi = vmlal_high_u8(Rx2hi, Dx4, Sax4);

    4. Деление на 255. Ну что, пришла пора опробовать в деле крутые сдвиги. В Си-версии это происходит так:

    #define DIV255(a)\
        ((((a + 0x80) >> 8) + a + 0x80) >> 8)

    Первое, что нужно рассмотреть — сдвиги с округлением. Для этого до сдвига прибавляется константа в один бит правее сдвига. По сути это уже сильно упрощает описание деления:

    #define ROUND_SHR(a, n)\
        ((a + (1<<(n-1))) >> n)
    #define DIV255(a)\
        ROUND_SHR(ROUND_SHR(a, 8) + a, 8)

    Дальше следует обратить внимание на конструкцию ROUND_SHR(a, 8) + a. Это же сдвиг с аккумулятором vrsraq_n_u16, помните? Ну а последний сдвиг можно сделать так, чтобы он заодно уменьшал разрядность результата, ведь тот не должен превышать 255. Кроме того, можно уменьшить разрядность не только в нижнюю половину вектора, но и в верхнюю (vqrshrn_high_n_u16).

    uint8x16_t vqrshrn_high_n_u16 (uint8x8_t r, uint16x8_t a, const int n);
    uint8x8_t vqrshrn_n_u16 (uint16x8_t a, const int n);
    uint16x8_t vrsraq_n_u16 (uint16x8_t a, uint16x8_t b, const int n);
    
    uint8x16_t Rx4 = vqrshrn_high_n_u16(
        vqrshrn_n_u16(vrsraq_n_u16(Rx2lo, Rx2lo, 8), 8),
        vrsraq_n_u16(Rx2hi, Rx2hi, 8), 8);

    Всё вместе:

    static void
    opSourceOver_premul(uint8_t* restrict Rrgba,
                        const uint8_t* restrict Srgba,
                        const uint8_t* restrict Drgba, size_t len)
    {
        size_t i = 0;
    
        for (; i < len*4 - 12; i += 16) {
            uint8x16_t Sx4 = vld1q_u8(&Srgba[i]);
            uint8x16_t Dx4 = vld1q_u8(&Drgba[i]);
    
            uint8x16_t Sax4 = vsubq_u8(
                vdupq_n_u8(255),
                vqtbl1q_u8(Sx4, (uint8x16_t){3,3,3,3, 7,7,7,7, 11,11,11,11, 15,15,15,15})
            );
    
            uint16x8_t Rx2lo = vmull_u8(vget_low_u8(Sx4), vdup_n_u8(255));
            uint16x8_t Rx2hi = vmull_high_u8(Sx4, vdupq_n_u8(255));
    
            Rx2lo = vmlal_u8(Rx2lo, vget_low_u8(Dx4), vget_low_u8(Sax4));
            Rx2hi = vmlal_high_u8(Rx2hi, Dx4, Sax4);
    
            uint8x16_t Rx4 = vqrshrn_high_n_u16(
                vqrshrn_n_u16(vrsraq_n_u16(Rx2lo, Rx2lo, 8), 8),
                vrsraq_n_u16(Rx2hi, Rx2hi, 8), 8);
            vst1q_u8(&Rrgba[i], Rx4);
        }
    
        for (; i < len*4; i += 4) {
            uint8_t Sa = Srgba[i + 3];
            Rrgba[i + 0] = DIV255(Srgba[i + 0] * 255 + Drgba[i + 0] * (255 - Sa));
            Rrgba[i + 1] = DIV255(Srgba[i + 1] * 255 + Drgba[i + 1] * (255 - Sa));
            Rrgba[i + 2] = DIV255(Srgba[i + 2] * 255 + Drgba[i + 2] * (255 - Sa));
            Rrgba[i + 3] = DIV255(Srgba[i + 3] * 255 + Drgba[i + 3] * (255 - Sa));
        }
    }

    Ну и наконец, можно порадоваться результату:

    $ clang-9 -Wall -O2 -o run.64 main.c && ./run.64
    Time elapsed: 0.047613
    Time elapsed: 0.047455
    Time elapsed: 0.047452
    Time elapsed: 0.047448

    Это в 1,75 раз быстрее, чем автовекторизованная версия и в 4 раза быстрее, чем версия совсем без векторизации (кстати, для GCC ускорение получается 5,5 раза).

    Оптимизация чтения

    Полученный результат полностью доказывает жизнеспособность программирования под NEON, однако, от использованной техники (чтение и обработка 4 пикселей, 8 умножений за инструкцию) можно было ожидать больший прирост скорости. Давайте посмотрим, что можно улучшить.

    Код выполнялся на процессоре Broadcom BCM2835, который хоть и относительно современный, но супербюджетный. Можно ожидать, что на таком процессоре могут быть существенные задержки даже для доступа к кешу L1. При этом никакие инструкции предвыборки не помогут, т.к. данные уже лежат в самом близком кеше. Зато может помочь предварительное чтение. Для этого нужно на каждом шаге класть «в карман» данные, которые понадобятся на следующем шаге. А из кармана доставать то, что было выбрано на предыдущем.

        uint8x16_t Sx4_next = vld1q_u8(&Srgba[0]);
        uint8x16_t Dx4_next = vld1q_u8(&Drgba[0]);
        // for (; i < len*4 - 12; i += 16) {
        for (; i < len*4 - 12 - 16; i += 16) {
            // uint8x16_t Sx4 = vld1q_u8(&Srgba[i]);
            // uint8x16_t Dx4 = vld1q_u8(&Drgba[i]);
            uint8x16_t Sx4 = Sx4_next;
            uint8x16_t Dx4 = Dx4_next;
            Sx4_next = vld1q_u8(&Srgba[i + 16]);
            Dx4_next = vld1q_u8(&Drgba[i + 16]);
            ...

    При этом нужно не забыть, что цикл нужно уменьшить на одну итерацию, чтобы не прочитать данные за пределами указателя.

    $ clang-9 -Wall -O2 -o run.64 main.c && ./run.64
    Time elapsed: 0.038070
    Time elapsed: 0.037855
    Time elapsed: 0.037834
    Time elapsed: 0.037831

    Гипотеза оказалась верной, это дало прирост ещё 25%. Итого NEON работает ровно в 5 раз быстрее, чем код без векторизации.

    Разбор сгенерированного кода

    Ну что, пора посмотреть, что напридумывал компилятор.

    $ clang-9 -Wall -O2 -o main.s main.c -S

    Это вывод уже отформатированный, с некоторыми переименованными регистрами и с комментариями оригинального кода:

        ldr     q0, [x19]                   // Sx4 = vld1q_u8(&Srgba[0])
        ldr     q1, [x20]                   // Dx4 = vld1q_u8(&Drgba[0])
        movi    v17.2d, #0xffffffffffffffff
        mov     x9, xzr
    .LBB0_3:
        tbl     v5.16b, { v0.16b }, v16.16b // vqtbl1q_u8(Sx4, v16)
        mvn     v5.16b, v5.16b              // Sax4 = vsubq_u8(vdupq_n_u8(255), v5)
        ext     v6.16b, v0.16b, v0.16b, #8
        umull   v3.8h, v1.8b, v5.8b         // Rx2lo = vmull_u8(Dx4, Sax4);
        add     x10, x19, x9
        add     x11, x20, x9
        umull   v4.8h, v6.8b, v17.8b        // Rx2hi = vmull_high_u8(Sx4, 0xff)
        umlal   v3.8h, v0.8b, v17.8b        // vmlal_u8(Rx2lo, Sx4, 0xff)
        umlal2  v4.8h, v1.16b, v5.16b       // vmlal_high_u8(Rx2hi, Dx4, Sax4)
        ldr     q0, [x10, #16]              // Sx4 = vld1q_u8(&Srgba[i + 16])
        ldr     q1, [x11, #16]              // Dx4 = vld1q_u8(&Drgba[i + 16])
        ursra   v3.8h, v3.8h, #8            // vrsraq_n_u16(Rx2lo, Rx2lo, 8)
        ursra   v4.8h, v4.8h, #8            // vrsraq_n_u16(Rx2hi, Rx2hi, 8)
        uqrshrn v2.8b, v3.8h, #8            // Rx4 = vqrshrn_n_u16(Rx2lo, 8)
        add     x10, x9, #16            
        uqrshrn2 v2.16b, v4.8h, #8          // vqrshrn_high_n_u16(Rx4, Rx2hi, 8)
        cmp     x10, #3972                  
        str     q2, [x21, x9]               // vst1q_u8(&Rrgba[i], Rx4)
        mov     x9, x10
        b.lo    .LBB0_3

    Хочу напомнить, что q0 и v0 — это разные названия одного регистра. Первое, что стоит отметить — компилятор верно понял задумку с упреждающим чтением и, несмотря на то, что в коде мы использовали копирование из старых переменных в новые, в выводе ничего такого нет. Чтение происходит сразу после того, как данные регистров v0 и v1 больше не используются. Хоть из-за этого загрузка и сместилась с начала итерации в середину, тут стоит положиться на компилятор, т.к. скорее всего он решил так сделать основываясь на латентности операций чтения.

    Далее стоит обратить внимание на последовательности umull и umlal. Тут произошло странное, зачем-то компилятор заменил umull2 на ещё один umull. Для этого ему понадобилось в строчке №8 сделать лишнее извлечение верхней части v0 во временный регистр v6. Формально мы использовали функцию vget_low_u8, которая это и подразумевает. Однако это было сделано только для того, чтобы привести переменную к нужному типу. Если посмотреть, что генерируют компиляторы для такого кода, то видно, что они не очень понимают, что нижняя часть регистра — это и есть сам регистр. Ну а GCC вообще творит какую-то дичь: создает два разных регистра с константами, делает три копирования.

    На этом странности не заканчиваются. Для Rx2lo переставлены местами vmull_u8 и vmlal_u8. Формально это ни на что не влияет, но все равно не понятно, зачем.

    Ну и последнее, на что можно обратить внимание — это странная работа с индексами. Если для команды str q2, [x21, x9] в качестве смещения используется регистр, то для ldr смещение вычисляется заранее, причем два раза, хотя очевидно, что можно было вычислить x9 + 16 и использовать это значение в обеих загрузках.

    Пробуем всё это исправить:

        ldr     q0, [x19]                   // Sx4 = vld1q_u8(&Srgba[0])
        ldr     q1, [x20]                   // Dx4 = vld1q_u8(&Drgba[0])
        movi    v17.2d, #0xffffffffffffffff
        mov     x9, xzr                     // i = 0
    .LBB0_3:
        tbl     v5.16b, { v0.16b }, v16.16b // vqtbl1q_u8(Sx4, v16)
        mvn     v5.16b, v5.16b              // Sax4 = vsubq_u8(vdupq_n_u8(255), v5)
        add     x10, x9, #16                // x10 = i + 16
        umull   v3.8h, v0.8b, v17.8b        // Rx2lo = vmull_u8(Sx4, 0xff);
        umull2  v4.8h, v0.16b, v17.16b      // Rx2hi = vmull_high_u8(Sx4, 0xff)
        ldr     q0, [x19, x10]              // Sx4 = vld1q_u8(&Srgba[i + 16])
        umlal   v3.8h, v1.8b, v5.8b         // vmlal_u8(Rx2lo, Dx4, Sax4)
        umlal2  v4.8h, v1.16b, v5.16b       // vmlal_high_u8(Rx2hi, Dx4, Sax4)
        ldr     q1, [x20, x10]              // Dx4 = vld1q_u8(&Drgba[i + 16])
        ursra   v3.8h, v3.8h, #8            // vrsraq_n_u16(Rx2lo, Rx2lo, 8)
        ursra   v4.8h, v4.8h, #8            // vrsraq_n_u16(Rx2hi, Rx2hi, 8)
        uqrshrn v2.8b, v3.8h, #8            // Rx4 = vqrshrn_n_u16(Rx2lo, 8)
        uqrshrn2 v2.16b, v4.8h, #8          // vqrshrn_high_n_u16(Rx4, Rx2hi, 8)
        cmp     x10, #3972                  
        str     q2, [x21, x9]               // vst1q_u8(&Rrgba[i], Rx4)
        mov     x9, x10
        b.lo    .LBB0_3

    Запускаем:

    $ clang-9 -Wall -O2 -o main.o -c main.s && gcc ./main.o -o run.64 && ./run.64 
    Time elapsed: 0.033388
    Time elapsed: 0.033204
    Time elapsed: 0.033223
    Time elapsed: 0.033190

    Есть ещё 14% прироста. Итого ускорение 5,7 раз.

    Было бы интересно также посчитать, сколько тактов уходит на этот цикл. Имеем 1800МГц (тактов/с), 0.033190 с/запуск и 250 * 20 * 1000 циклов/запуск. Итого: 1800000000 * 0.033190 / (250 * 20 * 1000) ≈ 12 тактов/цикл. Учитывая, что в цикле 17 инструкций, это прекрасный результат. Я не думал, что такой простой процессор может работать так эффективно.

    Решение на SSE

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

    __m128i Sax4 = _mm_sub_epi8(
        _mm_set1_epi8((char) 255),
        _mm_shuffle_epi8(Sx4, _mm_set_epi8(
            15,15,15,15, 11,11,11,11, 7,7,7,7, 3,3,3,3)));

    Существенное отличие только в том, что порядок байтов у _mm_set_epi8 инвертирован.

    Дальше нужно 8-битное умножение. В SSE есть такое, это интринсик _mm_maddubs_epi16. И кстати, он же увеличивает разрядность результата и даже делает сложение соседних элементов. Можно было бы подумать, что дело в шляпе.

    // Это неправильный код!
    __m128i Rx2lo = _mm_maddubs_epi16(
        _mm_unpacklo_epi8(_mm_set1_epi8((char) 255), Sax4),
        _mm_unpacklo_epi8(Sx4, Dx4));
    __m128i Rx2hi = _mm_maddubs_epi16(
        _mm_unpackhi_epi8(_mm_set1_epi8((char) 255), Sax4),
        _mm_unpackhi_epi8(Sx4, Dx4));

    Количество инструкций умножения уменьшилось вдвое по сравнению с NEON-версией. Но, к сожалению, _mm_maddubs_epi16 принимает только первый 8-битный аргумент как целое без знака, а второй аргумент считается со знаком, поэтому результат будет неверным. Функции, в которой оба аргумента были бы без знака, нет. А значит, нужно использовать 16-битное умножение и распаковывать каждый аргумент с пустым регистром, что существенно увеличивает и запутывает код.

    __m128i Rx2lo = _mm_add_epi16(
        _mm_mullo_epi16(_mm_unpacklo_epi8(Sx4, _mm_setzero_si128()),
                        _mm_set1_epi16(255)),
        _mm_mullo_epi16(_mm_unpacklo_epi8(Dx4, _mm_setzero_si128()),
                        _mm_unpacklo_epi8(Sax4, _mm_setzero_si128())));
    __m128i Rx2hi = _mm_add_epi16(
        _mm_mullo_epi16(_mm_unpackhi_epi8(Sx4, _mm_setzero_si128()),
                        _mm_set1_epi16(255)),
        _mm_mullo_epi16(_mm_unpackhi_epi8(Dx4, _mm_setzero_si128()),
                        _mm_unpackhi_epi8(Sax4, _mm_setzero_si128())));

    Деление на 255 и запаковку приходится делать, описывая каждый шаг, SSE никак не помогает.

    Rx2lo = _mm_add_epi16(Rx2lo, _mm_set1_epi16(0x80));
    Rx2lo = _mm_srli_epi16(_mm_add_epi16(_mm_srli_epi16(Rx2lo, 8), Rx2lo), 8);
    Rx2hi = _mm_add_epi16(Rx2hi, _mm_set1_epi16(0x80));
    Rx2hi = _mm_srli_epi16(_mm_add_epi16(_mm_srli_epi16(Rx2hi, 8), Rx2hi), 8);
    __m128i Rx4 = _mm_packus_epi16(Rx2lo, Rx2hi);

    За вычетом загрузок/сохранений, констант и приведений типов, я насчитал 23 интринсика в SSE версии против 10 в NEON. Предложения по улучшению преветствуются.

    Моё впечатление

    NEON произвел впечатление очень продуманной и эффективной системы команд. Я нашел для себя такие плюсы:

    • В отличие от SSE, есть консистентность типов данных, с которыми работают разные инструкции.

    • Порадовало большое количество опций сдвигов, которые оказались полезными на практике.

    • Можно встраивать NEON-код в любое место приложения без проверок рантайм.

    • Использование NEON дает ощутимый прирост производительности, примерно равный такому от использования SSE.

    • Очень понятный ассемблер с типизированными аргументами.

    • Были опасения, что будет сильно не хватать инструкции _mm_madd_epi16, делающей 8 умножений и 4 сложения. Однако её функциональность покрывается парой vmull_*/vmlal_*, не требующих подготовки данных.

    Минусы я бы отметил следующие:

    • Абсолютно неюзабельный справочник интринсиков. Это не минус самого NEON, но это то, с чем придется столкнуться.

    • Производительность сильно зависит от компилятора, возможно придется залочиться на clang.

    • Имена некоторых интринсиков напоминают читы в играх: vqrshrn_n_u16, vqdmulh_s16

    Бенчмарки

    Я собрал все варианты из этой статьи в один репозиторий с make-файлом, чтобы быстро запускать на разных системах или с разным окружением. Помимо Raspberry Pi 4 я смог запустить код ещё на c6g.large инстансе в AWS, которые работают на процессорах AWS Graviton2. Вот что я намерил:

    Raspberry Pi 4

    c6g.large

    GCC 8.3.0

    Clang 9.0.1

    GCC 9.3.0

    Clang 9.0.1

    Без векторизации

    267,4 мс
    7.94x

    185,6 мс
    5.51x

    140,2 мс
    10.01x

    103,8 мс
    7.41x

    Авто векторизация

    116,8
    3.47x

    82,55
    2.45x

    140,2
    10.01x

    46,54
    3.32x

    Ручная векторизация

    46,36
    1.38x

    47,56
    1.41x

    22,85
    1.63x

    17,59
    1.26x

    Оптимизация чтения

    46
    1.37x

    36,8
    1.09x

    24,01
    1.71x

    16,86
    1.2x

    Ассемблер

    33,66 мс
    1.0x

    14 мс
    1.0x

    Коэффициентами я обозначил замедление относительно варианта на ассемблере. Таблица получилась очень интересная. Выводов можно сделать много:

    • Поведение очень сильно зависит от компилятора. Протестированные версии GCC практически везде медленнее Clang.

    • Автоматическая векторизация в целом работает, но не дает такого же эффекта, как ручная.

    • Автоматическая векторизация может внезапно не заработать, причем даже на более свежей версии компилятора. При этом функция была выбрана предельно простая для векторизации.

    • GCC также не оценил оптимизацию чтения. Я не смотрел код, но выглядит так, будто он её просто выкинул.

    • На более мощных чипах любая векторизация дает более существенный выигрыш. При этом оптимизация чтения для них уже не так критична.

    • Пока что компиляторы не умеют полностью раскрывать возможности ARM, даже при использовании интринсиков.

    Кроме этого я все же решил измерить неизмеримое и сравнить несравнимое. Запустил тесты для ARM на Apple M1, а для x86 на Intel(R) Xeon. Выбор M1 понятен — кроме него пока нет десктопных процессоров от Apple. А вот на чем запускать x86 было вопросом. На ноутбуке у меня процессор может работать в диапазоне от 2,4 до 4,1 Ггц при разных сценариях. Поэтому я решил запустить на серверном процессоре, у которого стабильная частота.

    Ещё момент, в результатах на M1 был сильнейший разброс, вплоть до двух раз. Я запускал не у себя и у меня не было возможности проверить, почему так выходит. Моё предположение, что какое-то время работали слабые ядра. В любом случае, я брал минимальное время.

    Apple M1

    c5.large

    Clang 12.0.0

    GCC 9.3.0

    Clang 9.0.1

    Без векторизации

    43,13 мс
    4.74x

    74,56 мс
    5.09x

    80,26 мс
    6.16x

    Авто векторизация

    16,71
    1.84x

    74,54
    5.09x

    80,28
    6.16x

    Ручная 128-битная векторизация

    9,09
    1.0x

    14,65
    1.0x

    13,03
    1.0x

    Ручная 256-битная векторизация

    7,76
    0.53x

    6,52
    0.5x

    Коэффициентами я обозначил замедление относительно ручной 128-битной векторизации. Хочется напомнить, что это далеко не всеобъемлющий бенчмарк и по его результатам нельзя делать выводы о производительности всей платформы. Тем не менее.

    • Скорость M1 без векторизации впечатляет. Это при том, что частота обоих чипов примерно одинаковая.

    • Упс, авто векторизация на x86 не сработала на обоих компиляторах. А случай всё ещё простейший.

    • Несмотря на огромное количество кода в цикле (напомню, 23 инструкции против 10), x86 заметно сильнее ускоряется от векторизации. Это можно объяснить большим количеством исполнительных блоков для целочисленных вычислений внутри каждого ядра или более оптимальным микрокодом.

    • Хоть 128-битная версия на M1 всё еще выполняется быстрее, чем на x86, против AVX ему нечего противопоставить.


    На этом всё. Если нашли какие-то неточности, или знаете ещё что-то интересное о NEON и архитектуре ARM, делитесь в комментариях, обсудим вместе.

    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 36

      0

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

        –1
        для замеров производительности нужно фиксировать частоту работы ядер процессора. Иногда следует отключать ядра. Почитайте про подсистема ядра ОС cpufreq.
          0
          для замеров производительности нужно

          Кому нужно, для чего?


          Если вы про это: «на ноутбуке частота процессора в диапазоне от 2,4 до 4,1 Ггц при разных сценариях», то я на нем поэтому и не тестировал. Не потому что не знал, как зафиксировать, а потому что любое измерение будет неверным.


          У всех остальных процессоров, которые я тестировал, частоты примерно одинаковые при любых сценариях.

            0
            Нужно для того, кто замеряет. Для того, чтобы измерения были валидными. Фиксация частоты обеспечивает чистоту эксперимента. Думаю, что статья будет сильнее, если избегать оговорок о разбросе результатов (как на Apple M1), и если убедиться, что частота действительно фиксированна.
            В любом случае, интересный материал.
              +2
              Для того, чтобы измерения были валидными.

              Как раз наоборот, если вы специальным образом модифицировали устройство перед измерением, ваши измерения становятся бессмысленными.


              Та же векторизация может приводить к понижению частоты, и что толку с цифр, которые вы намерите на одинаковой частоте?


              В любом случае, интересный материал.

              Спасибо

          +1
          Но на практике, как я понял, SVE реализован только в Fugaku supercomputer.

          Из существующих — да. Ещё будет SiPearl Rhea с Neoverse V1.
          Но тут давеча анонсировали ARMv9, имеющий SVE2 в базе.
          www.anandtech.com/show/16584/arm-announces-armv9-architecture
          Так что теперь заживём(с)

          А ещё будет умножение матриц (GEMM) из ARMv8.6.
          community.arm.com/developer/ip-products/processors/b/processors-ip-blog/posts/arm-architecture-developments-armv8-6-a

          так и кастомная, которая может называться Apple Firestorm, Neoverse N1 или никак не называться.

          Neoverse N1 это стоковое серверное ядро компании ARM.
          Кастомные микроархитектуры используются в Apple A6+, Samsung M1-6, ThunderX/2/3, Ampere eMAG 8180, Phytium Mars.
            +1
            Интересно почитать.
            Задам вопрос, который меня давно волнует, но я так и не находил ответа на него.
            В около компьютерных темах бывает проскакивают фразы, что x86_64 подтормаживает тот багаж, что оно должно тащить за собой. Т.е. какие-то инструкции, которые оно обязано тащить за собой. Действительно ли это так?

              0

              Я не спец по проектированию ядер процессора, но кое-что знаю.


              Само по себе наличие каких-то инструкций может быть проедает транзисторный бюджет и увеличивает площадь кристалла, но вряд ли как-то сильно тормозить процессор. Гораздо сильнее x86 тормозит плохие архитектурные решения, в частности переменная длина инструкций. В не можете параллельно выбрать 4 инструкции и сразу начать их выполнять, нужно делать по одной.

                0
                4 можно, в теории даже 6 на Tiger Lake.
                www.agner.org/optimize/microarchitecture.pdf
                12.2 Instruction fetch and decoding
                  0

                  А как это работает, если нужно знать, где начинается каждая инструкция, а это невозможно без хотя бы частичного декодирования предыдущих?


                  На ум приходит только декодировать каждый байт, а потом отбрасывать неверные варианты. Но это и есть «тормозить от плохих архитектурных решений».

                    0
                    Да, насколько я понимаю, параллельно декодируются все возможные варианты длин.
                      +2
                      Но это и есть «тормозить от плохих архитектурных решений».
                      Только не «тормозить», а «греться». Добавляется ещё одна стадия в конвеер с 16 декодерами, которые умеют только считать длину и всё. Для «больших» ядер там там даже не так сильно транзисторный бюджет вырастает.
                        –1

                        Если мы не можем сразу взять и начать декодировать несколько инструкций подряд, а нужна ещё какая-то стадия, то это тормозить.

                          +1
                          Мне страшно становится. Вы точно пытаетесь что-то там ускорить работая на уровне машинных инструкций? И не имея представления о том, что конвеер в процессоры завезли ещё полвека назад?
                          Да, дополнительная стадия увеличивает нагрузку на предсказатель ветвлений.
                          Нет, это не приводит, само по себе, к тормозам.
                          Иначе бы в соврменных процессорах не удлиняли конвеер, чтобы увеличить IPC.
                          Вот тепловыделение это повышает однозначно, так что если вам нужен не быстрый однопоток, а сотни ядер, то это плохая идея.
                +1
                В 2010(11) мне пришлось познакомиться с OpenMax и NEON. Очень понравилось. Особенно команды в которых реализовано то, что у RISC-V называют сейчас fusion. Один код делает сразу сложение, сдвиг и округление. Очень интересны команды загрузки. Правда некоторые очень неудобны в реализации, т.к. длинный операнд приходится грузить сразу в несколько разных регистров. Удобны команды с увеличением или уменьшением разрядности данных результата. Не нужно делать спец действий. Грамотно. Спасибо за информацию о последних новинках в векторной составляющей ARM, т.к. эту архитектуру я уже давно забросил :)
                  +4
                  >> Имена некоторых интринсиков напоминают читы в играх: vqrshrn_n_u16, vqdmulh_s16

                  Тут всё просто

                  — vqrshrn_n_u16:
                  v — vector
                  q — saturating (ещё с ARMv6 так)
                  r — rounded
                  shr — shift right
                  n — narrow
                  _n_ immediate
                  u16 — unsigned 16

                  — vqdmulh_s16:
                  v — vector
                  q — saturating
                  d — doubling
                  mul — multiply
                  h — high half
                  s16 — signed 16

                  >> Моё предположение, что какое-то время работали слабые ядра.
                  Из состояния сна сначала включаются мелкие ядра, потом через 30мс происходит миграция на большие. Пиковая частота достигается через 100мс.
                  www.anandtech.com/show/14892/the-apple-iphone-11-pro-and-max-review/5
                  Но это в iOS. На macOS надо бы измерить.
                    0
                    vqdmulh_s16

                    А вот кстати, чего я не понял, для чего нужно doubling?

                      +2
                      Это DSP инструкция.
                      Например, умножаем знаковый 16-битный семпл 0x7fff на fixed-point громкость 0x7fff (1.0)
                      Т.е. семпл должен оставаться на макс громкости -> в старшей части должны получить примерно те же 0x7fff.

                      integer round_const = if rounding then 1 << (esize — 1) else 0;
                      //esize == 16 -> round_const == 0x8000
                      product = (2 * element1 * element2) + round_const;
                      0x7FFF*0x7FFF
                      0x3FFF0001 * 2
                      0x7FFE0002 + 0x00008000
                      =
                      0x7FFE8002 -> 0x7FFE

                      Есть варианты без округления и с ним («if rounding»)
                      vqdmulh_s16
                      vqrdmulh_s16

                      upd: Пофиксил round_const
                        +3
                        doubling это как раз случай когда операции с одной разрядностью (например, 16), а результат уже с удвоенной (32). Далее проводим промежуточные вычисления уже с 32-разрядными данными, а в конце результат уже к меньшему формату (16). Итого более точные вычисления.
                          +3
                          Вы говорите про widening (противоположность narrowing).
                          int32x4_t vmlal_n_s16 (int32x4_t a, int16x4_t b, int16_t c)
                          Vector widening multiply accumulate with scalar

                          Инструкция, которая делает doubling (vqdmulh_s16) дополнительно умножает промежуточный результат на 2, как я расписал.
                          Да, инструкция считает с повышенной разрядностью, но слово doubling относится не к этому.
                          Signed saturating Doubling Multiply returning High half. This instruction multiplies the values of corresponding elements of the two source SIMD&FP registers, doubles the results, places the most significant half of the final results into a vector, and writes the vector to the destination SIMD&FP register.

                          Вот строка псевдокода из мануала, которая это описывает:
                          product = (2 * element1 * element2) + round_const;

                            +2
                            Спасибо за поправку. Действительно я посмотрел только vqdmul и перепутал. То, что я подразумевал реально есть widening и narrowing. Подзабыл уже :)
                      0
                      «разработано расширение Scalable Vector Extension (SVE), которое позволяет выполнять один и тот же код на чипах, реализующих разный размер векторов. Но на практике, как я понял, SVE реализован только в Fugaku supercomputer.»
                      У RISC-V векторное расширение вроде бы из этой же оперы. Правда они все никак не могут остановится на конечном варианте :)
                        0
                        у м1 есть еще AMX перевод
                        но вот напрямую доступа к нему нет, только косвенно через accelerate framework.
                        зы — 2 недели назад на досуге начал ковырять neon. интересно посмотреть, есть ли разница между ним и использование апи accelerate на моих задачах
                          0
                          Спасибо за полезный материал! Жаль, что эти инструкции всё равно останутся невостребованными, так как компиляторы не умеют полноценно векторизовать сложные манипуляции с данными (например, когда в теле цикла нам нужны значения соседних пикселей), а в сфере мобильной разработки доминируют корпорации, призывающие людей не писать нативный код.
                          Native code is primarily useful when you have an existing native codebase that you want to port to Android, not for «speeding up» parts of your Android app written with the Java language.
                          — с developer.android.com
                            0

                            Я не рассматриваю ARM как платформу для мобильного разработки, мне он интересен как серверная платформа в ближайшие несколько лет, там мне никакие корпорации не указ.

                              0
                              Корпорации могут призывать к чему угодно, но когда те же самые корпорации тупо не предоставляют никакого 3D API или NN API (neural networks, тоже очень «тяжёлая» вещь) для управляемого кода — становится ясно, что таки в тех циклах, где вы захотите брать значения соседних пикселей у вас будет явно не Java.
                              +2

                              По вашему "benchmark" кстати видно регрессию gcc. Создал баг https://gcc.gnu.org/bugzilla/show_bug.cgi?id=99856 .

                                0

                                Там на самом деле есть объективное отличие в векторизованном коде. Дело в том, что выражение Srgba[i + 0] * 255 + Drgba[i + 0] * (255 - Sa) это по-максимум 255*255 + 255*255, что переполнение для 16-битного числа: 0x1fc02. И казалось бы, переполнение будет в обоих версиях, так какая разница. Но в скаларной версии при делении на 255 самый старший бит попадёт в младший разряд 8-битного результата и увеличит его на 1 при сложении.


                                В векторизованном коде мы пользуемся знанием о том, что значения пикселей уже умножены на альфу, соответственно, если (255 - Sa) = 255, то Srgba[i + 0] = 0, а значит всё выражение всегда будет влезать в 16 бит. У компилятора таких знаний нет.


                                Кажется я пробовал искусственно ограничивать значение выражение:


                                DIV255((Srgba[i + 0] * 255 + Drgba[i + 0] * (255 - Sa)) & 0xffff)

                                Но вроде это не включало векторизацию, а скалярный код становился медленнее.

                                  0
                                  казалось бы, переполнение будет в обоих версиях, так какая разница

                                  По идее в "скалярном" никаких переполнений не будет. Там будет numeric promotions. Там сразу все в "int" будет считаться, который на данных платформах 32 битный.


                                  В векторизованном коде мы пользуемся знанием о том, что значения пикселей уже умножены на альфу, соответственно, если (255 — Sa) = 255, то Srgba[i + 0] = 0

                                  Вы хотели сказать, что Srbga[i + 3] == 0, ведь Sa это Srgba[i + 3], а не Srgba[i + 0]?

                                    0

                                    Нет, я все правильно написал. Premultiplied alpha — это формат хранения RGBA пикселей, где каждое RGB заранее умножено на A. Это упрощает расчеты в большинстве случаев. И это означает, что все Srbga[i + x] не могут быть больше Srbga[i + 3].

                                      0

                                      А assume пробовали добавлять? Например clang __builtin_assume с инвариантом, что Srbga[i + x] <= Srbga[i + 3] ?

                                        0

                                        Нет. Попробуйте, ссылка на репозиторий в топике. Там всё предельно просто: есть ридми, make-файл, нет зависимостей.

                                +1
                                Касательно оптимизации загрузки NEON векторов, я бы еще посмотрел в сторону выровненной загрузки и использования префетча:

                                #define PREFECH_SIZE 384
                                
                                        template <bool align> inline uint8x16_t Load(const uint8_t * p);
                                
                                        template <> inline uint8x16_t Load<false>(const uint8_t * p)
                                        {
                                #if defined(__GNUC__) && PREFECH_SIZE
                                            __builtin_prefetch(p + PREFECH_SIZE);
                                #endif
                                            return vld1q_u8(p);
                                        }
                                
                                        template <> inline uint8x16_t Load<true>(const uint8_t * p)
                                        {
                                #if defined(__GNUC__)
                                #if PREFECH_SIZE
                                            __builtin_prefetch(p + PREFECH_SIZE);
                                #endif
                                            uint8_t * _p = (uint8_t *)__builtin_assume_aligned(p, 16);
                                            return vld1q_u8(_p);
                                #elif defined(_MSC_VER)
                                            return vld1q_u8_ex(p, 128);
                                #else
                                            return vld1q_u8(p);
                                #endif
                                        }
                                

                                В некоторых случаях помогает (величиной PREFECH_SIZE лучше поиграться).
                                  +1
                                  В общем дошли руки посмотреть вашу программу.
                                  Это слишком маленькая нагрузка чтобы запустились быстрые ядра M1 на полной скорости.
                                  Время тоже измеряется плохо. Нужно считать минимум. У вас же считается среднее от времени работы на разной частоте.
                                  Даже если сделать 4 итерации вокруг внутреннего цикла, видно как происходит рост производительности.
                                  Ещё я увеличил на 1 итерацию внешний цикл.

                                      size_t tmin = ~0; 
                                      for (int t = 0; t < 4; t++)
                                      {   
                                          gettimeofday(&tval_before, NULL);
                                          for (size_t i = 0; i < 20 * 1000; i ++) {
                                              opSourceOver_premul(Rrgba, Srgba, Drgba, len);
                                          }
                                          gettimeofday(&tval_after, NULL);
                                          timersub(&tval_after, &tval_before, &tval_result);
                                          size_t tdelta = tval_result.tv_sec * 1000000 + tval_result.tv_usec; 
                                          if (tmin > tdelta)
                                            tmin = tdelta;
                                      }
                                      printf("Time elapsed: %ld us\n", tmin);  


                                  cc -Wall -O2 -o run.64 main.c impl.native.c && ./run.64
                                  Time elapsed: 16890 us
                                  Time elapsed: 14802 us
                                  Time elapsed: 14788 us
                                  Time elapsed: 14776 us
                                  Time elapsed: 14792 us

                                  cc -Wall -O2 -o run.64 main.c impl.native.c -fno-tree-vectorize && ./run.64
                                  Time elapsed: 43059 us
                                  Time elapsed: 43079 us
                                  Time elapsed: 43182 us
                                  Time elapsed: 43006 us
                                  Time elapsed: 42992 us

                                  cc -Wall -O2 -o run.64 main.c impl.native.c -ftree-vectorize && ./run.64
                                  Time elapsed: 16938 us
                                  Time elapsed: 14798 us
                                  Time elapsed: 14772 us
                                  Time elapsed: 14806 us
                                  Time elapsed: 14799 us

                                  cc -Wall -O2 -o run.64 main.c impl.neon.c && ./run.64
                                  Time elapsed: 9123 us
                                  Time elapsed: 6238 us
                                  Time elapsed: 5118 us
                                  Time elapsed: 4584 us
                                  Time elapsed: 4446 us

                                  cc -Wall -O2 -o run.64 main.c impl.neon_preload.c && ./run.64
                                  Time elapsed: 9077 us
                                  Time elapsed: 6194 us
                                  Time elapsed: 5113 us
                                  Time elapsed: 4580 us
                                  Time elapsed: 4429 us


                                  Видите ваши 9мс? Это примерно половина того, что делает M1 на полной частоте.

                                  Хоть 128-битная версия на M1 всё еще выполняется быстрее, чем на x86, против AVX ему нечего противопоставить.

                                  Как выясняется, вполне есть что.
                                    0

                                    Какая-то фантастика, в среднем 2,8 тактов на цикл из 17 инструкций. Как такое возможно?

                                      +1
                                      А что тут такого?
                                      Устойчивый темп 8 инструкций за такт, четыре блока NEON (512-бит).
                                      Все данные в кэше сидят.

                                      Кстати если сделать len побольше, скажем, 100000000, чтобы данные не влезали в кэши.
                                      то тогда в однопотоке M1 будет ещё более быстрым чем x86.

                                      Собрал под Windows и если ничего не напутал, получается так:
                                      Ryzen 3950X SSE Time elapsed: 1261000 us
                                      Ryzen 3950X AVX Time elapsed: 1124000 us
                                      M1: Time elapsed: 447691 us

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

                                  Only users with full accounts can post comments. Log in, please.