Считайте каждую битовую позицию отдельно для множества 64-битных битовых масок, используя AVX, но не AVX2 - PullRequest
12 голосов
/ 09 марта 2019

Я работаю над проектом на C, где мне нужно пройти через десятки миллионов масок (типа ulong (64-бит)) и обновить массив (называемый target) из 64 коротких целых чисел (uint16) на основе по простому правилу:

// for any given mask, do the following loop
for (i = 0; i < 64; i++) {
    if (mask & (1ull << i)) {
        target[i]++
    }
}

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

В настоящее время я использую gcc 4.8.4 в Ubuntu 14.04 (i7-2670QM, поддерживающий AVX, а не AVX2) для компиляции и запуска следующего кода, и это заняло около 2 секунд. Хотелось бы, чтобы он работал под 200 мс.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <sys/stat.h>

double getTS() {
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return tv.tv_sec + tv.tv_usec / 1000000.0;
}
unsigned int target[64];

int main(int argc, char *argv[]) {
    int i, j;
    unsigned long x = 123;
    unsigned long m = 1;
    char *p = malloc(8 * 10000000);
    if (!p) {
        printf("failed to allocate\n");
        exit(0);
    }
    memset(p, 0xff, 80000000);
    printf("p=%p\n", p);
    unsigned long *pLong = (unsigned long*)p;
    double start = getTS();
    for (j = 0; j < 10000000; j++) {
        m = 1;
        for (i = 0; i < 64; i++) {
            if ((pLong[j] & m) == m) {
                target[i]++;
            }
            m = (m << 1);
        }
    }
    printf("took %f secs\n", getTS() - start);
    return 0;
}

Заранее спасибо!

Ответы [ 4 ]

12 голосов
/ 10 марта 2019

В моей системе, 4-летнем MacBook (Intel Core i5 с частотой 2,7 ГГц) с clang-900.0.39.2 -O3, ваш код работает за 500 мс.

Простое изменение внутреннего теста на if ((pLong[j] & m) != 0) экономит 30%, работаяза 350 мс.

Дальнейшее упрощение внутренней части до target[i] += (pLong[j] >> i) & 1; без теста приводит к уменьшению ее до 280 мс.

Для дальнейших улучшений требуются более совершенные методы, такие как распаковка битов в блоки по 8ulongs и добавляя их параллельно, обрабатывая 255 ulongs за раз.

Вот улучшенная версия, использующая этот метод.он работает в моей системе за 45 мс.

#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <sys/stat.h>

double getTS() {
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return tv.tv_sec + tv.tv_usec / 1000000.0;
}

int main(int argc, char *argv[]) {
    unsigned int target[64] = { 0 };
    unsigned long *pLong = malloc(sizeof(*pLong) * 10000000);
    int i, j;

    if (!pLong) {
        printf("failed to allocate\n");
        exit(1);
    }
    memset(pLong, 0xff, sizeof(*pLong) * 10000000);
    printf("p=%p\n", (void*)pLong);
    double start = getTS();
    uint64_t inflate[256];
    for (i = 0; i < 256; i++) {
        uint64_t x = i;
        x = (x | (x << 28));
        x = (x | (x << 14));
        inflate[i] = (x | (x <<  7)) & 0x0101010101010101ULL;
    }
    for (j = 0; j < 10000000 / 255 * 255; j += 255) {
        uint64_t b[8] = { 0 };
        for (int k = 0; k < 255; k++) {
            uint64_t u = pLong[j + k];
            for (int kk = 0; kk < 8; kk++, u >>= 8)
                b[kk] += inflate[u & 255];
        }
        for (i = 0; i < 64; i++)
            target[i] += (b[i / 8] >> ((i % 8) * 8)) & 255;
    }
    for (; j < 10000000; j++) {
        uint64_t m = 1;
        for (i = 0; i < 64; i++) {
            target[i] += (pLong[j] >> i) & 1;
            m <<= 1;
        }
    }
    printf("target = {");
    for (i = 0; i < 64; i++)
        printf(" %d", target[i]);
    printf(" }\n");
    printf("took %f secs\n", getTS() - start);
    return 0;
}

Методика раздувания байта до 64-битной длины исследуется и объясняется в ответе: https://stackoverflow.com/a/55059914/4593267.Я сделал массив target локальной переменной, а также массив inflate и распечатал результаты, чтобы компилятор не оптимизировал вычисления.В рабочей версии вы будете вычислять массив inflate отдельно.

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

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

11 голосов
/ 10 марта 2019

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

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


См. , есть ли обратная инструкция к инструкции movemask в intel avx2? для краткого изложения методов растрового изображения -> векторных для различных размеров. Предложение Ext3h в другом ответе хорошо: распаковка битов в нечто более узкое, чем массив окончательного подсчета, дает вам больше элементов на инструкцию. Байты эффективны с SIMD, и затем вы можете сделать до 255 вертикальных paddb без переполнения перед распаковкой для накопления в массив 32-битных счетчиков.

Для хранения всех 64 uint8_t элементов требуется только 4x 16-байтовых __m128i векторов, поэтому эти аккумуляторы могут оставаться в регистрах, добавляя их в память только при расширении до 32-битных счетчиков во внешнем цикле.

Распаковка не обязательно должна быть в порядке : вы всегда можете перетасовать target[] один раз в самом конце, собрав все результаты.

Внутренний цикл можно развернуть для запуска с 64- или 128-битной векторной загрузкой и распаковать 4 или 8 различными способами, используя pshufb (_mm_shuffle_epi8).


Еще лучшая стратегия - постепенно расширяться

Начиная с 2-разрядных аккумуляторов, затем маскируйте / сдвигайте, чтобы расширить их до 4-разрядных. Таким образом, в самом внутреннем цикле большинство операций работают с «плотными» данными, а не «разводят» их слишком много сразу. Более высокая плотность информации / энтропии означает, что каждая инструкция выполняет больше полезной работы.

Использование SWAR методов для 32-битного 2-битного добавления внутри скалярных или SIMD-регистров легко / дешево, потому что мы все равно должны избегать возможности выполнения вершины элемента. При правильном SIMD мы потеряли бы это количество, а при использовании SWAR мы испортили бы следующий элемент.

uint64_t x = *(input++);        // load a new bitmask
const uint64_t even_1bits = 0x5555555555555555;  // 0b...01010101;

uint64_t lo = x & even_1bits;
uint64_t hi = (x>>1) & even_1bits;            // or use ANDN before shifting to avoid a MOV copy

accum2_lo += lo;   // can do up to 3 iterations of this without overflow
accum2_hi += hi;   // because a 2-bit integer overflows at 4

Затем вы повторяете до 4 векторов 4-битных элементов, затем 8 векторов 8-битных элементов, затем вам нужно расширить до 32 и накапливать в массиве в памяти, потому что вы все равно исчерпаете регистры и эта работа с внешним внешним циклом достаточно редка, поэтому нам не нужно беспокоиться о переходе на 16-битный режим. (Особенно если мы вручную векторизируем).

Самый большой недостаток: этот не автоматически векторизован, в отличие от версии @njuffa. Но с gcc -O3 -march=sandybridge для AVX1 (затем запускается код на Skylake), этот запущенный скаляр 64-разрядная на самом деле все еще немного быстрее , чем 128-разрядная AVX, автоматически векторизованная asm из кода @ njuffa.

Но это время для Skylake, у которого есть 4 скалярных порта ALU (и mov-elmination), в то время как Sandybridge не имеет mov-el сокращений и имеет только 3 порта ALU, поэтому скалярный код, вероятно, достигнет узких мест внутреннего порта выполнения. (Но SIMD-код может быть почти таким же быстрым, потому что есть много AND / ADD, смешанных со сдвигами, и у SnB есть исполнительные блоки SIMD на всех 3 его портах, на которых есть какие-либо ALU. Haswell только что добавил порт 6 для скалярного -только включая смены и филиалы.)

При хорошей ручной векторизации это должно быть почти в 2 или 4 раза быстрее.

Но если вам нужно выбрать между этим скаляром или @ njuffa с автовекторизацией AVX2, @ njuffa быстрее на Skylake с -march=native

Если построение на 32-битной цели возможно / требуется, это сильно сказывается (без векторизации из-за использования uint64_t в 32-битных регистрах), в то время как векторизованный код почти не страдает (потому что вся работа происходит в векторных регистрах такой же ширины).

// TODO: put the target[] re-ordering somewhere
// TODO: cleanup for N not a multiple of 3*4*21 = 252
// TODO: manual vectorize with __m128i, __m256i, and/or __m512i

void sum_gradual_widen (const uint64_t *restrict input, unsigned int *restrict target, size_t length)
{
    const uint64_t *endp = input + length - 3*4*21;     // 252 masks per outer iteration
    while(input <= endp) {
        uint64_t accum8[8] = {0};     // 8-bit accumulators
        for (int k=0 ; k<21 ; k++) {
            uint64_t accum4[4] = {0};  // 4-bit accumulators can hold counts up to 15.  We use 4*3=12
            for(int j=0 ; j<4 ; j++){
                uint64_t accum2_lo=0, accum2_hi=0;
                for(int i=0 ; i<3 ; i++) {  // the compiler should fully unroll this
                    uint64_t x = *input++;    // load a new bitmask
                    const uint64_t even_1bits = 0x5555555555555555;
                    uint64_t lo = x & even_1bits; // 0b...01010101;
                    uint64_t hi = (x>>1) & even_1bits;  // or use ANDN before shifting to avoid a MOV copy
                    accum2_lo += lo;
                    accum2_hi += hi;   // can do up to 3 iterations of this without overflow
                }

                const uint64_t even_2bits = 0x3333333333333333;
                accum4[0] +=  accum2_lo       & even_2bits;  // 0b...001100110011;   // same constant 4 times, because we shift *first*
                accum4[1] += (accum2_lo >> 2) & even_2bits;
                accum4[2] +=  accum2_hi       & even_2bits;
                accum4[3] += (accum2_hi >> 2) & even_2bits;
            }
            for (int i = 0 ; i<4 ; i++) {
                accum8[i*2 + 0] +=   accum4[i] & 0x0f0f0f0f0f0f0f0f;
                accum8[i*2 + 1] +=  (accum4[i] >> 4) & 0x0f0f0f0f0f0f0f0f;
            }
        }

        // char* can safely alias anything.
        unsigned char *narrow = (uint8_t*) accum8;
        for (int i=0 ; i<64 ; i++){
            target[i] += narrow[i];
        }
    }
    /* target[0] = bit 0
     * target[1] = bit 8
     * ...
     * target[8] = bit 1
     * target[9] = bit 9
     * ...
     */
    // TODO: 8x8 transpose
}

Нас не волнует порядок, поэтому accum4[0] имеет 4-битные аккумуляторы для каждого 4-го бита, например. Последним исправлением, необходимым (но еще не реализованным) в самом конце, является транспонирование 8x8 массива uint32_t target[64], , которое может быть эффективно выполнено с использованием unpck и vshufps только с AVX1.( Транспонировать поплавок 8x8 с использованием AVX / AVX2 ).А также цикл очистки для последних до 251 маски.

Мы можем использовать любую ширину элемента SIMD для реализации этих сдвигов;мы все равно должны маскироваться для ширины ниже 16-битной (SSE / AVX не имеет сдвигов гранулярности, только 16-битный минимум.)

Результаты тестов в Arch Linux i7-6700k из тестового ремня @ njuffa, с этим добавлено( Godbolt ) N = (10000000 / (3*4*21) * 3*4*21) = 9999864 (т. Е. 10000000, округленное до кратного 252 итерационного коэффициента "развертывания", поэтому моя упрощенная реализация выполняет тот же объем работы, не считая повтор-заказ target[], что не делает, поэтому выводит результаты о несоответствии. Но количество отпечатков совпадает с другой позицией ссылочного массива.)

Я запускал программу 4 раза подряд (чтобы убедиться,ЦП прогрелся до максимальной скорости турбокомпрессора) и взял один из прогонов, который выглядел хорошо (ни один из 3-х кратно аномально высоких).

ref: лучший битовый цикл (следующий раздел)
быстро:код @ нюффа.(автоматическая векторизация с использованием 128-битных целочисленных инструкций AVX).
постепенный: моя версия (не векторизована gcc или clang, по крайней мере, не во внутреннем цикле.) gcc и clang полностью развертывают внутренние 12 итераций.

  • gcc8.2 -O3 -march=sandybridge -fpie -no-pie
    ref: 0,331373 секунды, быстрый: 0,011387 секунд, постепенный: 0,009966 секунд
  • gcc8.2 -O3 -march=sandybridge -fno-pie -no-pie
    ref: 0,397175 секунд, быстрый: 0,011255 секунд,постепенное: 0,010018 с
  • clang7.0 -O3 -march=sandybridge -fpie -no-pie
    ref: 0,352381 с, быстрое: 0,011926 с, постепенное: 0,009269 с (очень низкое значение для порта 7 мопов, clang использует индексированную адресацию для магазинов)
  • clang7.0 -O3 -march=sandybridge -fno-pie -no-pie
    ref: 0,293014 с , быстрый: 0,011777 с, постепенный: 0,009235 с

-марш = skylake (допуская AVX2 для 256-битных целочисленных векторов) помогает обоим, но больше всего @ njuffa, потому что большая часть его векторизована (включая его самый внутренний цикл):

  • gcc8.2 -O3 -march=skylake -fpie -no-pie
    ref: 0,328725 с, быстрое: 0,007621 с, постепенное: 0,010054 с (gcc не показывает усиления для «постепенного», только «быстрое»)
  • gcc8.2 -O3 -march=skylake -fno-pie -no-pie
    ref: 0,333922 с, быстрый: 0,007620 с, постепенный: 0,009866 с

  • clang7.0 -O3 -march=skylake -fpie -no-pie
    ref: 0,260616 с,быстрый: 0,007521 с, постепенный: 0,008535 с (ИДК, почему постепенный быстрее, чем -march = Sandybridge;он не использует BMI1 andn.Я думаю, потому что он использует 256-битный AVX2 для внешнего цикла k = 0..20 с vpaddq)

  • clang7.0 -O3 -march=skylake -fno-pie -no-pie
    ref: 0.259159 с , быстро: 0,007496 с , постепенное: 0,008671 с

Без AVX только SSE4.2: (-march=nehalem), причудливое ускорение ускоренияс AVX / мелодия = песчаный мост."fast" только чуть медленнее, чем с AVX.

  • gcc8.2 -O3 -march=skylake -fno-pie -no-pie
    ref: 0,337178 с, быстрый: 0,011983 с, постепенный: 0,010587 с
  • clang7.0 -O3 -march=skylake -fno-pie -no-pie
    ref: 0,293555 с , быстрый: 0,012549 с, постепенный: 0,008697 с

-fprofile-generate / -fprofile-use помогают некоторым для GCC, особеннодля версии «ref», где по умолчанию она вообще не разворачивается.

Я выделил лучшее, но часто они находятся в пределах допустимых помех друг для друга.Неудивительно, что -fno-pie -no-pie иногда был быстрее: индексирование статических массивов с помощью [disp32 + reg] является , а не режимом индексированной адресации, просто base + disp32, поэтому он никогда не запускается на процессорах семейства Sandybridge.

Но с gcc иногда -fpie был быстрее;Я не проверял, но я предполагаю, что gcc просто как-то выстрелил себе в ногу, когда возможна абсолютная 32-битная адресация.Или просто невинно выглядящие различия в code-gen вызвали проблемы с выравниванием или uop-кешем;Я не проверял подробно.


Для SIMD мы можем просто делать 2 или 4x uint64_t параллельно, накапливая только горизонтально на последнем шаге, где мы расширяем байты до 32-битных элементов. (Возможно, путем перетасовки в очереди изатем используя pmaddubsw с множителем _mm256_set1_epi8(1) для добавления горизонтальных пар байтов в 16-битные элементы.)

TODO: векторизация вручную __m128i и __m256i__m512i) версии этого.Должно быть примерно в 2, 4, или даже 8 раз быстрее, чем «постепенное» время, указанное выше. Вероятно, предварительная выборка HW все еще может идти в ногу с этим, за исключением, может быть, версии AVX512 с данными, поступающими из DRAM, особенно если есть конфликт от другихпотоки.Мы выполняем значительный объем работы с каждым прочитанным словом.


Устаревший код: улучшения в битовом цикле

Ваша портативная скалярная версия тоже может быть улучшена ускорение с ~ 1,92 секунды ( с общей ошибкой прогнозирования ветвления 34% , с быстрыми зацикленными комментариями!) До ~ 0,35 с (clang7.0 -O3 -march=sandybridge) с надлежащим случайным входом на 3,9 ГГцSkylake.Или 1,83 с для версии с ветвлением с != 0 вместо == m, поскольку компиляторы не могут доказать, что m всегда имеет ровно 1 установленный бит, и / или оптимизировать соответственно.

(против 0,01 с для@ njuffa или моя быстрая версия выше, так что это совершенно бесполезно в абсолютном смысле, но стоит упомянуть в качестве общего примера оптимизации использования кода без ответвлений.)

Если вы ожидаете случайное сочетание нулей ите, которые вы хотите что-то без разветвлений, которые не будут неправильно прогнозировать.Выполнение += 0 для элементов, которые были равны нулю, избегает этого, а также означает, что абстрактная машина C определенно касается этой памяти независимо от данных.

Компиляторам не разрешается изобретать записи, поэтому, если они хотят автоматически-векторизируйте вашу if() target[i]++ версию, им придется использовать хранилище в маске, например x86 vmaskmovps, чтобы избежать неатомарного чтения / перезаписи неизмененных элементов target.Поэтому некоторым гипотетическим будущим компиляторам, которые могут автоматически векторизовать простой скалярный код, было бы легче с этим.

В любом случае, один из способов написать это - target[i] += (pLong[j] & m != 0);, используя преобразование bool-> int, чтобы получить 0/ 1 целое число.

Но мы получим лучший ассемблер для x86 (и, вероятно, для большинства других архитектур), если просто сдвинуть данные и изолировать младший бит с помощью &1.Компиляторы довольно глупы и, похоже, не замечают этой оптимизации.Они хорошо оптимизируют счетчик дополнительных циклов и превращают m <<= 1 в add same,same для эффективного смещения влево, но они все еще используют xor-zero / test / setne для создания целого числа 0/1.

Внутренний цикл, подобный этому, компилируется немного эффективнее (но все же намного хуже, чем мы можем сделать с SSE2 или AVX, или даже скалярный, используя таблицу поиска @ chrqlie, которая будет оставаться горячей в L1d при многократном использованиивот так, разрешив SWAR в uint64_t):

    for (int j = 0; j < 10000000; j++) {
#if 1  // extract low bit directly
        unsigned long long tmp = pLong[j];
        for (int i=0 ; i<64 ; i++) {   // while(tmp) could mispredict, but good for sparse data
            target[i] += tmp&1;
            tmp >>= 1;
        }
#else // bool -> int shifting a mask
        unsigned long m = 1;
        for (i = 0; i < 64; i++) {
            target[i]+= (pLong[j] & m) != 0;
            m = (m << 1);
        }
#endif

Обратите внимание, что unsigned long не гарантированно является 64-битным типом и не поддерживается в x86-64 System V x32 (ILP32 в64-битный режим) и Windows x64.Или в 32-разрядных интерфейсах ABI, таких как i386 System V.

Скомпилировано в проводнике компилятора Godbolt с помощью gcc, clang и ICC , это меньше на 1 моп в цикле с gcc.Но все они просто скаляры, clang и ICC развернуты на 2.

# clang7.0 -O3 -march=sandybridge
.LBB1_2:                            # =>This Loop Header: Depth=1
   # outer loop loads a uint64 from the src
    mov     rdx, qword ptr [r14 + 8*rbx]
    mov     rsi, -256
.LBB1_3:                            #   Parent Loop BB1_2 Depth=1
                                    # do {
    mov     edi, edx
    and     edi, 1                              # isolate the low bit
    add     dword ptr [rsi + target+256], edi   # and += into target

    mov     edi, edx
    shr     edi
    and     edi, 1                              # isolate the 2nd bit
    add     dword ptr [rsi + target+260], edi

    shr     rdx, 2                              # tmp >>= 2;

    add     rsi, 8
    jne     .LBB1_3                       # } while(offset += 8 != 0);

Это немного лучше, чем мы получаем от test / setnz.Без развертывания bt / setc могли бы быть равными, но компиляторы плохо используют bt для реализации bool (x & (1ULL << n)) или bts для реализации x |= 1ULL << n.

Еслиу многих слов самый высокий установленный бит намного ниже бита 63, зацикливание на while(tmp) может быть выигрышем .Неправильные прогнозы ветвей делают его не стоящим, если он экономит от ~ 0 до 4 итераций большую часть времени, но если он часто экономит 32 итерации, это может стоить того.Возможно, разверните исходный код, чтобы цикл проверял только tmp каждые 2 итерации (потому что компиляторы не будут выполнять это преобразование за вас), но тогда ветвь цикла может быть shr rdx, 2 / jnz.

В семействе Sandybridge это 11 мопов с плавким доменом для внешнего интерфейса на 2 бита ввода.(add [mem], reg с неиндексированным режимом адресации микросопрягает нагрузку + ALU и адрес хранилища + данные хранилища, все остальное - одиночные uop. Add / jcc macro-fuses. См. Руководство Agner Fog и https://stackoverflow.com/tags/x86/info). Таким образом, он должен работать примерно с 3 циклами на 2 бита = один uint64_t на 96 циклов (Sandybridge не «разворачивается» внутри своего буфера цикла, так что число мопов, не кратное 4, в основном округляетсявверх, в отличие от Haswell и более поздних версий).

против не развернутой версии gcc: 7 мопов на 1 бит = 2 цикла на бит. Если вы скомпилировали с gcc -O3 -march=native -fprofile-generate / test-run / gcc -O3 -march=native -fprofile-use,Оптимизация по профилю позволила бы развернуть цикл.

Это, вероятно, медленнее, чем ветвящаяся версия для полностью предсказуемых данных, как вы получаете из memset с любым повторяющимся байтовым шаблоном . Я бы предложилзаполнение массива случайно сгенерированными данными из быстрого PRNG, например, xorshift + SSE2, или, если вы просто синхронизируете цикл подсчета, тогда используйте все, что захотите, например rand().

8 голосов
/ 10 марта 2019

Один из способов значительно ускорить это, даже без AVX, состоит в том, чтобы разбить данные на блоки до 255 элементов и накапливать счетчики битов в байтах в обычных uint64_t переменных. Поскольку исходные данные имеют 64 бита, нам нужен массив из 8 байтовых аккумуляторов. Первый аккумулятор считает биты в позициях 0, 8, 16, ... 56, второй аккумулятор считает биты в позициях 1, 9, 17, ... 57; и так далее. После того, как мы закончили обработку блока данных, мы переносим счетчики из байтового аккумулятора в счетчики target. Функция обновления счетчиков target для блока до 255 номеров может быть закодирована простым способом в соответствии с описанием выше, где BITS - это количество битов в исходных данных:

/* update the counts of 1-bits in each bit position for up to 255 numbers */
void sum_block (const uint64_t *pLong, unsigned int *target, int lo, int hi)
{
    int jj, k, kk;
    uint64_t byte_wise_sum [BITS/8] = {0};
    for (jj = lo; jj < hi; jj++) {
        uint64_t t = pLong[jj];
        for (k = 0; k < BITS/8; k++) {
            byte_wise_sum[k] += t & 0x0101010101010101;
            t >>= 1;
        }
    }
    /* accumulate byte sums into target */
    for (k = 0; k < BITS/8; k++) {
        for (kk = 0; kk < BITS; kk += 8) {
            target[kk + k] += (byte_wise_sum[k] >> kk) & 0xff;
        }
    }
}

Вся программа ISO-C99, которая должна работать как минимум на платформах Windows и Linux, показана ниже. Он инициализирует исходные данные с помощью PRNG, выполняет проверку на соответствие эталонной реализации запрашивающего и сравнивает эталонный код и ускоренную версию. На моей машине (Intel Xeon E3-1270 v2 @ 3,50 ГГц) при компиляции с MSVS 2010 при полной оптимизации (/Ox) вывод программы составляет:

p=0000000000550040
ref took 2.020282 secs, fast took 0.027099 secs

, где ref относится к исходному решению аскера. Ускорение здесь примерно в 74 раза. Различные ускорения будут наблюдаться с другими (и особенно более новыми) компиляторами.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>

#if defined(_WIN32)
#if !defined(WIN32_LEAN_AND_MEAN)
#define WIN32_LEAN_AND_MEAN
#endif
#include <windows.h>
double second (void)
{
    LARGE_INTEGER t;
    static double oofreq;
    static int checkedForHighResTimer;
    static BOOL hasHighResTimer;

    if (!checkedForHighResTimer) {
        hasHighResTimer = QueryPerformanceFrequency (&t);
        oofreq = 1.0 / (double)t.QuadPart;
        checkedForHighResTimer = 1;
    }
    if (hasHighResTimer) {
        QueryPerformanceCounter (&t);
        return (double)t.QuadPart * oofreq;
    } else {
        return (double)GetTickCount() * 1.0e-3;
    }
}
#elif defined(__linux__) || defined(__APPLE__)
#include <stddef.h>
#include <sys/time.h>
double second (void)
{
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return (double)tv.tv_sec + (double)tv.tv_usec * 1.0e-6;
}
#else
#error unsupported platform
#endif

/*
  From: geo <gmars...@gmail.com>
  Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
  Subject: 64-bit KISS RNGs
  Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)

  This 64-bit KISS RNG has three components, each nearly
  good enough to serve alone.    The components are:
  Multiply-With-Carry (MWC), period (2^121+2^63-1)
  Xorshift (XSH), period 2^64-1
  Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64  (kiss64_t = (kiss64_x << 58) + kiss64_c, \
                kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
                kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64  (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
                kiss64_y ^= (kiss64_y << 43))
#define CNG64  (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

#define N          (10000000)
#define BITS       (64)
#define BLOCK_SIZE (255)

/* cupdate the count of 1-bits in each bit position for up to 255 numbers */
void sum_block (const uint64_t *pLong, unsigned int *target, int lo, int hi)
{
    int jj, k, kk;
    uint64_t byte_wise_sum [BITS/8] = {0};
    for (jj = lo; jj < hi; jj++) {
        uint64_t t = pLong[jj];
        for (k = 0; k < BITS/8; k++) {
            byte_wise_sum[k] += t & 0x0101010101010101;
            t >>= 1;
        }
    }
    /* accumulate byte sums into target */
    for (k = 0; k < BITS/8; k++) {
        for (kk = 0; kk < BITS; kk += 8) {
            target[kk + k] += (byte_wise_sum[k] >> kk) & 0xff;
        }
    }
}

int main (void) 
{
    double start_ref, stop_ref, start, stop;
    uint64_t *pLong;
    unsigned int target_ref [BITS] = {0};
    unsigned int target [BITS] = {0};
    int i, j;

    pLong = malloc (sizeof(pLong[0]) * N);
    if (!pLong) {
        printf("failed to allocate\n");
        return EXIT_FAILURE;
    }
    printf("p=%p\n", pLong);

    /* init data */
    for (j = 0; j < N; j++) {
        pLong[j] = KISS64;
    }

    /* count bits slowly */
    start_ref = second();
    for (j = 0; j < N; j++) {
        uint64_t m = 1;
        for (i = 0; i < BITS; i++) {
            if ((pLong[j] & m) == m) {
                target_ref[i]++;
            }
            m = (m << 1);
        }
    }
    stop_ref = second();

    /* count bits fast */
    start = second();
    for (j = 0; j < N / BLOCK_SIZE; j++) {
        sum_block (pLong, target, j * BLOCK_SIZE, (j+1) * BLOCK_SIZE);
    }
    sum_block (pLong, target, j * BLOCK_SIZE, N);
    stop = second();

    /* check whether result is correct */
    for (i = 0; i < BITS; i++) {
        if (target[i] != target_ref[i]) {
            printf ("error @ %d: res=%u ref=%u\n", i, target[i], target_ref[i]);
        }
    }

    /* print benchmark results */
    printf("ref took %f secs, fast took %f secs\n", stop_ref - start_ref, stop - start);
    return EXIT_SUCCESS;
}
7 голосов
/ 09 марта 2019

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

Так что просто следуйте следующей стратегии для распаковки битов в байты вектора: https://stackoverflow.com/a/24242696/2879325

Теперь, когда вы добавили каждый бит к 8 битам, вы можете просто сделать это для блоков до 255 битовых масок за раз и собрать их все в один векторный регистр. После этого вам следует ожидать потенциального переполнения, поэтому вам необходимо выполнить перевод.

После каждого блока из 255, распакуйте снова в 32bit и добавьте в массив. (Вам не нужно делать ровно 255, просто какое-то удобное число меньше 256, чтобы избежать переполнения байтовых аккумуляторов).

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


typedef uint64_t T;
const size_t bytes = 8;
const size_t bits = bytes * 8;
const size_t block_size = 128;

static inline __m256i expand_bits_to_bytes(uint32_t x)
{
    __m256i xbcast = _mm256_set1_epi32(x);    // we only use the low 32bits of each lane, but this is fine with AVX2

    // Each byte gets the source byte containing the corresponding bit
    const __m256i shufmask = _mm256_set_epi64x(
        0x0303030303030303, 0x0202020202020202,
        0x0101010101010101, 0x0000000000000000);
    __m256i shuf = _mm256_shuffle_epi8(xbcast, shufmask);

    const __m256i andmask = _mm256_set1_epi64x(0x8040201008040201);  // every 8 bits -> 8 bytes, pattern repeats.
    __m256i isolated_inverted = _mm256_andnot_si256(shuf, andmask);

    // this is the extra step: byte == 0 ? 0 : -1
    return _mm256_cmpeq_epi8(isolated_inverted, _mm256_setzero_si256());
}

void bitcount_vectorized(const T *data, uint32_t accumulator[bits], const size_t count)
{
    for (size_t outer = 0; outer < count - (count % block_size); outer += block_size)
    {
        __m256i temp_accumulator[bits / 32] = { _mm256_setzero_si256() };
        for (size_t inner = 0; inner < block_size; ++inner) {
            for (size_t j = 0; j < bits / 32; j++)
            {
                const auto unpacked = expand_bits_to_bytes(static_cast<uint32_t>(data[outer + inner] >> (j * 32)));
                temp_accumulator[j] = _mm256_sub_epi8(temp_accumulator[j], unpacked);
            }
        }
        for (size_t j = 0; j < bits; j++)
        {
            accumulator[j] += ((uint8_t*)(&temp_accumulator))[j];
        }
    }
    for (size_t outer = count - (count % block_size); outer < count; outer++)
    {
        for (size_t j = 0; j < bits; j++)
        {
            if (data[outer] & (T(1) << j))
            {
                accumulator[j]++;
            }
        }
    }
}

void bitcount_naive(const T *data, uint32_t accumulator[bits], const size_t count)
{
    for (size_t outer = 0; outer < count; outer++)
    {
        for (size_t j = 0; j < bits; j++)
        {
            if (data[outer] & (T(1) << j))
            {
                accumulator[j]++;
            }
        }
    }
}

В зависимости от выбранного компилятора, векторизованная форма достигла ускорения примерно в 25 раз по сравнению с наивным.

В Ryzen 5 1600X векторизованная форма примерно достигла прогнозируемой пропускной способности ~ 600 000 000 элементов в секунду.

Удивительно, но на самом деле это все еще на 50% медленнее, чем решение, предложенное @ njuffa.

...