Векторизовать случайную инициализацию и распечатать для BigInt с десятичным массивом di git, с AVX2? - PullRequest
1 голос
/ 12 апреля 2020

Как я могу передать свой код в код AVX2 и получить тот же результат, что и раньше?

Можно ли использовать __m256i в функциях LongNumInit, LongNumPrint вместо uint8_t *L или в каком-либо подобном типе переменной?

Мои знания AVX весьма ограничены; Я немного исследовал, однако я не очень хорошо понимаю, как преобразовать мой код, любые предложения и объяснения приветствуются.

Мне действительно интересен этот код в AVX2.

void LongNumInit(uint8_t *L, size_t N )
{
  for(size_t i = 0; i < N; ++i){
      L[i] = myRandom()%10;
  }

}
void LongNumPrint( uint8_t *L, size_t N, uint8_t *Name )
{
  printf("%s:", Name);
  for ( size_t i=N; i>0;--i )
  {
    printf("%d", L[i-1]);
  }
  printf("\n");
}
int main (int argc, char **argv)
{
  int i, sum1, sum2, sum3, N=10000, Rep=50;

  seed = 12345;

  // obtain parameters at run time
  if (argc>1) { N    = atoi(argv[1]); }
  if (argc>2) { Rep  = atoi(argv[2]); }

 // Create Long Nums
  unsigned char *V1= (unsigned char*) malloc( N);
  unsigned char *V2= (unsigned char*) malloc( N);
  unsigned char *V3= (unsigned char*) malloc( N);
  unsigned char *V4= (unsigned char*) malloc( N);

  LongNumInit ( V1, N ); LongNumInit ( V2, N ); LongNumInit ( V3, N );

//Print last 32 digits of Long Numbers
  LongNumPrint( V1, 32, "V1" );
 LongNumPrint( V2, 32, "V2" );
  LongNumPrint( V3, 32, "V3" );
  LongNumPrint( V4, 32, "V4" );

  free(V1); free(V2); free(V3); free(V4);
  return 0;
}

Результат, который я получаю в своем исходном коде, таков:

V1:59348245908804493219098067811457
V2:24890422397351614779297691741341
V3:63392771324953818089038280656869
V4:00000000000000000000000000000000

1 Ответ

4 голосов
/ 13 апреля 2020

Это ужасный формат для BigInteger в целом, см. https://codereview.stackexchange.com/a/237764 для обзора кода проекта fl aws с использованием одного десятичного числа * git на байт для BigInteger, и что вы может / должен сделать вместо этого.

И посмотрите Могут ли длинные целочисленные подпрограммы извлечь выгоду из SSE? для заметок @ Mysticial о способах хранения ваших данных, которые делают SIMD для BigInteger математикой практичным, в частности частичное слово арифметика c, где ваши временные значения не могут быть "нормализованы", что позволяет вам выполнять ленивую обработку переноса.


Но, очевидно, вы просто спрашиваете о этом коде, функции random-init и print, а не как делать математику между двумя числами в этом формате.

Мы можем довольно хорошо векторизовать оба из них. Мой LongNumPrintName() является заменой для вас.

Для LongNumInit Я просто показываю строительный блок, который хранит два 32-байтовых блока и возвращает увеличенный указатель. Назовите это в al oop. (Естественно, он производит 2 вектора на вызов, поэтому для маленьких N вы можете сделать альтернативную версию.)

LongNumInit

Какой самый быстрый способ создания текстового файла объемом 1 ГБ, содержащего случайные цифры? генерирует разделенные пробелами случайные десятичные цифры ASCII со скоростью около 33 ГБ / с на 4 ГГц Skylake, включая издержки системных вызовов write() на /dev/null. (Это выше, чем пропускная способность DRAM; блокировка кэша на 128 КБ позволяет хранилищам попадать в кэш L2. Драйвер ядра для /dev/null даже не считывает буфер пространства пользователя.)

Его можно легко адаптировать в версию AVX2 void LongNumInit(uint8_t *L, size_t N ). Мой ответ там использует AVX2 xorshift128 + PRNG (векторизованный с 4 независимыми PRNG в 64-битных элементах __m256i), такой как AVX / SSE версия xorshift128 + . Это должно быть качество случайности, аналогичное вашему rand() % 10.

. Оно разбивает его на десятичные цифры с помощью мультипликативного обратного, чтобы разделить и по модулю на 10 со сдвигами и vpmulhuw, используя Почему G CC использовать умножение на странное число при реализации целочисленного деления? . (На самом деле, используя собственный синтаксис GNU C, чтобы G CC определял постоянную magi c и генерировал умножения и сдвиги для удобного синтаксиса, такого как v16u dig1 = v % ten; и v /= ten;)

. Вы можете использовать _mm256_packus_epi16 для упаковки двух векторов 16-битных цифр в 8-битные элементы вместо преобразования нечетных элементов в ASCII ' ' и четных элементов в ASCII '0'..'9'. (Поэтому измените vec_store_digit_and_space, чтобы упаковать пары векторов вместо ORing с константой.)

Скомпилируйте это с помощью g cc, clang или I CC (или, надеюсь, любого другого компилятора, который понимает GNU C диалект C99 и внутренние характеристики Intel.

См. https://gcc.gnu.org/onlinedocs/gcc/Vector-Extensions.html для части __attribute__((vector_size(32))) и https://software.intel.com/sites/landingpage/IntrinsicsGuide/ для _mm256_* вещи. Также { ссылка }.

#include <immintrin.h>

// GNU C native vectors let us get the compiler to do stuff like %10 each element
typedef unsigned short v16u __attribute__((vector_size(32)));

// returns p + size of stores.  Caller should use outpos = f(vec, outpos)
// p must be aligned
__m256i* vec_store_digit_and_space(__m256i vec, __m256i *restrict p)
{
    v16u v = (v16u)vec;
    v16u ten = (v16u)_mm256_set1_epi16(10);

    v16u divisor = (v16u)_mm256_set1_epi16(6554);  // ceil((2^16-1) / 10.0)
    v16u div6554 = v / divisor;      // Basically the entropy from the upper two decimal digits: 0..65.
    // Probably some correlation with the modulo-based values, especially dig3, but we do this instead of
    // dig4 for more ILP and fewer instructions total.

    v16u dig1 = v % ten;
    v /= ten;
    v16u dig2 = v % ten;
    v /= ten;
    v16u dig3 = v % ten;
    //  dig4 would overlap much of the randomness that div6554 gets

    // __m256i or v16u assignment is an aligned store
    v16u *vecbuf = (v16u*)p;
    vecbuf[0] = _mm256_packus_epi16(div6554, dig1);
    vecbuf[1] = _mm256_packus_epi16(dig2, dig3)
    return p + 2;  // always a constant number of full vectors
}

Лог c в random_decimal_fill_buffer, который вставляет символы новой строки, может быть полностью удален, потому что вы просто хотите получить плоский массив десятичных цифр. Просто вызывайте вышеуказанную функцию в al oop, пока вы не заполните свой буфер.

Обработка небольших размеров (меньше полного вектора):

Было бы удобно дополнить ваш mallo c до следующего кратного 32 байта, поэтому всегда безопасно выполнить 32-байтовую загрузку, не проверяя, возможно, переход на не отображенную страницу.

И используйте C11 aligned_alloc для получить 32-байтовое выровненное хранилище. Так, например, aligned_alloc(32, (size+31) & -32). Это позволяет нам делать полные 32-байтовые хранилища, даже если N нечетно. Логически только первые N байтов буфера содержат наши реальные данные, но удобно иметь отступы, которые мы можем набросать, чтобы избежать каких-либо дополнительных условных проверок, если N меньше 32 или не кратно 32.

К сожалению, ISO C и glib c отсутствуют aligned_realloc и aligned_calloc. MSV C на самом деле предоставляет такие: Почему на большинстве платформ нет 'align_realloc'? , что позволяет иногда выделять больше места в конце выровненного буфера, не копируя его. «Try_reallo c» идеально подходит для C ++, которому может потребоваться запуск конструкторов копирования, если нетривиально копируемые объекты меняют адрес. Неэкспрессивные API распределителя, которые вызывают иногда ненужное копирование, - моя любимая мозоль.


LongNumPrint

Принятие аргумента uint8_t *Name - плохой дизайн. Если вызывающий хочет сначала напечатать строку "something:", он может это сделать. Ваша функция должна просто делать то, что printf "%d" делает для int.

Так как вы сохраняете свои цифры в обратном порядке печати, вы захотите сделать обратный байт в буфер tmp и преобразовать 0,9-байтовые значения в '0'..'9' ASCII-символьные значения, используя OR '0'. Затем передайте этот буфер в fwrite.

В частности, используйте alignas(32) char tmpbuf[8192]; в качестве локальной переменной.

Вы можете работать с кусками фиксированного размера (например, 1kiB или 8kiB), вместо этого выделяя потенциально Огромный буфер. Возможно, вы захотите по-прежнему go через stdio (вместо write() напрямую и управлять своей собственной буферизацией ввода / вывода). С буфером 8 КБ эффективный fwrite может просто передать это непосредственно в write() вместо memcpy в буфер stdio. Возможно, вы захотите поэкспериментировать с настройкой этого, но если размер буфера tmp будет меньше половины кеша L1d, это будет означать, что он остается горячим в кеше, когда он перечитывается после того, как вы его записали.

Блокировка кеша делает l oop ограничивает намного более сложным, но это стоит того для очень больших N.

Байт-реверсирование 32 байта за раз :

Вы можете избежать этой работы, решив, что ваши цифры хранятся в MSD-первом порядке, но затем, если вы хотите реализовать сложение, потребуется l oop от конца назад.

Ваша функция может быть реализована с помощью SIMD _mm_shuffle_epi8 чтобы перевернуть 16-байтовые чанки, начиная с конца вашего массива di git и записи в начало вашего буфера tmp.

Или лучше загрузить vmovdqu / vinserti128 16-байтовые загрузки

для процессоров Intel, vinserti128 декодирует до загрузки + ALU uop, но может работать на любом векторном порте ALU , а не просто порт случайного воспроизведения. Таким образом, две 128-разрядные загрузки более эффективны, чем 256-разрядная загрузка -> vpshufb -> vpermq, что, вероятно, является узким местом для пропускной способности shuffle-порта, если данные были горячими в кэше. Процессоры Intel могут выполнять до 2 загрузок + 1 запоминание за такт (или в IceLake 2 загрузки + 2 запоминания). Вероятно, мы столкнемся с узким местом на внешнем интерфейсе, если нет узких мест в памяти, поэтому на практике не нужно насыщать порты загрузки + хранения и перемешивания. (https://agner.org/optimize/ и https://uops.info/)

Эта функция также упрощена, если предположить, что мы всегда можем прочитать 32 байта из L без перехода в не нанесенная на карту страница. Но после 32-байтового обратного преобразования для маленьких N первые N байтов ввода становятся последними N байтами в 32-байтовом фрагменте. Было бы наиболее удобно, если бы мы всегда могли безопасно выполнить 32-байтовую загрузку , заканчивающуюся в конце буфера, но нецелесообразно ожидать заполнения перед объектом.

#include <immintrin.h>
#include <stdalign.h>
#include <stddef.h>
#include <stdio.h>
#include <stdint.h>

// one vector of 32 bytes of digits, reversed and converted to ASCII
static inline
void ASCIIrev32B(void *dst, const void *src)
{
    __m128i hi = _mm_loadu_si128(1 + (const __m128i*)src);  // unaligned loads
    __m128i lo = _mm_loadu_si128(src);
    __m256i v = _mm256_set_m128i(lo, hi);    // reverse 128-bit hi/lo halves

    // compilers will hoist constants out of inline functions
    __m128i byterev_lane = _mm_set_epi8(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);      
    __m256i byterev = _mm256_broadcastsi128_si256(byterev_lane);  // same in each lane
    v = _mm256_shuffle_epi8(v, byterev);               // in-lane reverse

    v = _mm256_or_si256(v, _mm256_set1_epi8('0'));     // digits to ASCII
    _mm256_storeu_si256(dst, v);                       // Will usually be aligned in practice.
}

// Tested for N=32; could be bugs in the loop bounds for other N
// returns bytes written, like fwrite: N means no error, 0 means error in all fwrites
size_t LongNumPrint( uint8_t *num, size_t N)
{
    // caller can print a name if it wants

    const int revbufsize = 8192;      // 8kiB on the stack should be fine
    alignas(32) char revbuf[revbufsize];

    if (N<32) {
        // TODO: maybe use a smaller revbuf for this case to avoid touching new stack pages
        ASCIIrev32B(revbuf, num);   // the data we want is at the *end* of a 32-byte reverse
        return fwrite(revbuf+32-N, 1, N, stdout);
    }

    size_t bytes_written = 0;
    const uint8_t *inp = num+N;  // start with last 32 bytes of num[]
    do {
        size_t chunksize = (inp - num >= revbufsize) ? revbufsize : inp - num;

        const uint8_t *inp_stop = inp - chunksize + 32;   // leave one full vector for the end
        uint8_t *outp = revbuf;
        while (inp > inp_stop) {        // may run 0 times
            inp -= 32;
            ASCIIrev32B(outp, inp);
            outp += 32;
        }
        // reverse first (lowest address) 32 bytes of this chunk of num
        // into last 32 bytes of this chunk of revbuf
        // if chunksize%32 != 0 this will overlap, which is fine.
        ASCIIrev32B(revbuf + chunksize - 32, inp_stop - 32);
        bytes_written += fwrite(revbuf, 1, chunksize, stdout);
        inp = inp_stop - 32;
    } while ( inp > num );

    return bytes_written;
    // caller can putchar('\n') if it wants
}


// wrapper that prints name and newline
void LongNumPrintName(uint8_t *num, size_t N, const char *name)
{
    printf("%s:", name);
    //LongNumPrint_scalar(num, N);
    LongNumPrint(num, N);
    putchar('\n');
}

// main() included on Godbolt link that runs successfully

Это компилирует и запускает ( на Годболте ) с gcc -O3 -march=haswell и выдает идентичный вывод вашему скаляру l oop для N = 32, который проходит main. (Я использовал rand() вместо MyRandom(), поэтому мы могли бы протестировать с тем же начальным числом и получить те же числа, используя вашу функцию инициализации.)

Не тестировалось для больших N, но общая идея chunksize = min (ptrdiff, 8k) и использование этого значения до l oop вниз от конца num[] должно составлять solid.

Мы можем загружать (а не просто хранить) выровненные векторы, если преобразуем первый N%32 байт и передал это fwrite перед запуском основного l oop. Но это, вероятно, приведет либо к дополнительному системному вызову write(), либо к неуклюжему копированию внутри stdio. (Если не было уже буферизованного текста, еще не напечатанного, как Name:, в этом случае у нас уже есть такое наказание.)

Обратите внимание, что технически C UB уменьшается inp после старта num. Так что inp -= 32 вместо inp = inp_stop-32 будет иметь тот UB для итерации, который оставляет внешний l oop. Я на самом деле избегаю этого в этой версии, но в целом это работает в любом случае, потому что я думаю, что G CC предполагает плоскую модель памяти, а де-фактор определяет поведение сравнений указателя достаточно. И обычные операционные системы резервируют нулевую страницу, поэтому num определенно не может быть в пределах 32 байт от начала физической памяти (поэтому inp не может переноситься по старшему адресу.) Этот абзац по большей части оставлен после первого совершенно непроверенная попытка, как мне показалось, уменьшать указатель во внутренней l oop, чем она была на самом деле.

...