Самый быстрый способ умножить массив int64_t?

Ответ 1

Кажется, вы считаете, что long - 64 бит в вашем коде, но затем используйте __uint64_t. В 32 бит x32 ABI, а в Windows long - 32-битный тип. В вашем названии упоминается long long, но тогда ваш код игнорирует его. Некоторое время я интересовался, если ваш код предполагал, что long был 32 бит.

Вы полностью стреляете в ногу, используя нагрузки AVX256, но затем накладывая указатель на __m256i для выполнения скалярных операций. gcc просто сдается и дает вам ужасный код, который вы просили: векторная нагрузка, а затем куча инструкций extract и insert. Ваш способ записи означает, что оба вектора должны быть распакованы для выполнения sub в скаляре, а не для использования vpsubq.

Современные процессоры x86 имеют очень быстрый кеш L1, который может обрабатывать две операции за такт. (Хасуэлл и позже: две нагрузки и один магазин за часы). Выполнение нескольких скалярных нагрузок из одной и той же строки кэша лучше, чем векторная загрузка и распаковка. (Несовершенное вычисление uop уменьшает пропускную способность примерно до 84% от этого: см. Ниже)

gcc 5.3 -O3 -march = haswell (Godbolt compiler explorer) автоматически сгенерирует простую скалярную реализацию. Когда AVX2 недоступен, gcc по-дурацки по-прежнему авто-векторизует с векторами 128b: на Haswell это будет примерно на 1/2 скорости идеального скалярного 64-битного кода. (см. анализ perf ниже, но замените 2 элемента на вектор вместо 4).

#include <stdint.h>    // why not use this like a normal person?
#define BASE_VEX_STOP 1024
#define BASE_DIMENSION 1028

// restrict lets the compiler know the arrays don't overlap,
// so it doesn't have to generate a scalar fallback case
void multiply_simple(uint64_t *restrict Gi_vec, uint64_t q, const uint64_t *restrict Gj_vec){
    for (intptr_t i=0; i<BASE_DIMENSION; i++)   // gcc doesn't manage to optimize away the sign-extension from 32bit to pointer-size in the scalar epilogue to handle the last less-than-a-vector elements
        Gi_vec[i] -= Gj_vec[i] * q;
}

внутренний цикл:

.L4:
    vmovdqu ymm1, YMMWORD PTR [r9+rax]        # MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B]
    add     rcx, 1    # ivtmp.30,
    vpsrlq  ymm0, ymm1, 32      # tmp174, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B],
    vpmuludq        ymm2, ymm1, ymm3        # tmp173, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B], vect_cst_.25
    vpmuludq        ymm0, ymm0, ymm3        # tmp176, tmp174, vect_cst_.25
    vpmuludq        ymm1, ymm4, ymm1        # tmp177, tmp185, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B]
    vpaddq  ymm0, ymm0, ymm1    # tmp176, tmp176, tmp177
    vmovdqa ymm1, YMMWORD PTR [r8+rax]        # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B]
    vpsllq  ymm0, ymm0, 32      # tmp176, tmp176,
    vpaddq  ymm0, ymm2, ymm0    # vect__13.24, tmp173, tmp176
    vpsubq  ymm0, ymm1, ymm0    # vect__14.26, MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__13.24
    vmovdqa YMMWORD PTR [r8+rax], ymm0        # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__14.26
    add     rax, 32   # ivtmp.32,
    cmp     rcx, r10  # ivtmp.30, bnd.14
    jb      .L4 #,

Переведите это обратно в intrinsics, если хотите, но это будет намного проще просто дать автогенератору компилятор. Я не пытался анализировать его, чтобы убедиться, что он оптимален.

Если вы обычно не компилируете с помощью -O3, вы можете использовать #pragma omp simd перед циклом (и -fopenmp).

Конечно, вместо скалярного эпилога это было бы проблематично. быстрее выполнить невыровненную нагрузку последних 32B Gj_vec и сохранить в последние 32B Gi_vec, потенциально перекрывающиеся с последним хранилищем из цикла. (Скалярный резерв по-прежнему необходим, если массивы меньше 32B.)

Улучшенная внутренняя версия вектора для Haswell

Из моих комментариев по поводу ответа Z Boson. На основе код библиотеки векторного класса Agner Fog.

Версия Agner Fog сохраняет инструкцию, но узкие места в порту тасования, используя phadd + pshufd, где я использую psrlq/paddq/pand.

Так как один из ваших операндов является постоянным, убедитесь, что он прошел set1(q) как b, а не a, поэтому можно поднять тасовку "bswap".

// replace hadd -> shuffle (4 uops) with shift/add/and (3 uops)
// The constant takes 2 insns to generate outside a loop.
__m256i mul64_haswell (__m256i a, __m256i b) {
    // instruction does not exist. Split into 32-bit multiplies
    __m256i bswap   = _mm256_shuffle_epi32(b,0xB1);           // swap H<->L
    __m256i prodlh  = _mm256_mullo_epi32(a,bswap);            // 32 bit L*H products

    // or use pshufb instead of psrlq to reduce port0 pressure on Haswell
    __m256i prodlh2 = _mm256_srli_epi64(prodlh, 32);          // 0  , a0Hb0L,          0, a1Hb1L
    __m256i prodlh3 = _mm256_add_epi32(prodlh2, prodlh);      // xxx, a0Lb0H+a0Hb0L, xxx, a1Lb1H+a1Hb1L
    __m256i prodlh4 = _mm256_and_si256(prodlh3, _mm256_set1_epi64x(0x00000000FFFFFFFF)); // zero high halves

    __m256i prodll  = _mm256_mul_epu32(a,b);                  // a0Lb0L,a1Lb1L, 64 bit unsigned products
    __m256i prod    = _mm256_add_epi64(prodll,prodlh4);       // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32
    return  prod;
}

Посмотрите на Godbolt.

Обратите внимание, что это не включает окончательный вычитание, только умножение.

Эта версия должна работать лучше на Haswell, чем на gcc autovectorized. (например, один вектор за 4 цикла вместо одного вектора за 5 циклов, узкий по пропускной способности порта 0. Я не рассматривал другие узкие места для полной проблемы, поскольку это было последним дополнением к ответу.)

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

Перф-анализ автосекторизованного кода GCC (не встроенная версия)

Фон: см. Agner Fog insn tables и руководство микроархива и другие ссылки в x86 теги wiki.

До AVX512 (см. ниже) это, вероятно, всего лишь быстрее, чем скалярный 64-битный код: imul r64, m64 имеет пропускную способность одного за такт на процессорах Intel (но один на 4 такта на семействе AMD Bulldozer). load/imul/sub-with-memory-dest - это 4 процессора с интегрированными доменами на процессорах Intel (с режимом адресации, который может быть с микро-предохранителем, который gcc не может использовать). Ширина конвейера - это 4 сдобренных домена на каждом часе, поэтому даже большой разворот не может заставить его выдавать один раз в часы. При достаточном разворачивании мы будем сталкиваться с пропускной способностью загрузки/хранения. 2 аккумулятора и один магазин за часы возможны на Haswell, но магазины-адреса uops, загружающие порты загрузки, будут снизить пропускную способность примерно до 81/96 = 84% от этого, согласно руководству Intel.

Поэтому, возможно, лучший способ для Haswell будет загружать и умножать на скаляр, (2 uops), затем vmovq/pinsrq/vinserti128, чтобы вы могли сделать вычитание с помощью vpsubq. Это 8 uops для загрузки и умножения всех 4 скаляров, 7 shuffle uops, чтобы получить данные в __m256i (2 (movq) + 4 (pinsrq is 2 uops) + 1 vinserti128) и еще 3 раза для загрузки векторной нагрузки/vpsubq/векторный магазин. Таким образом, 18 совпадающих доменов fused-domain на 4 умножения (4.5 цикла для выдачи), но 7 shuffle uops (7 циклов для выполнения). Итак, nvm, это нехорошо по сравнению с чистым скаляром.

Автовекторный код использует 8 векторных инструкций ALU для каждого вектора из четырех значений. На Haswell, 5 из этих uops (умножения и сдвиги) могут работать только на порте 0, поэтому независимо от того, как вы разворачиваете этот алгоритм, он достигнет в лучшем случае одного вектора за 5 циклов (т.е. Умножается на 5/4 циклов).

Сдвиги могут быть заменены на pshufb (порт 5), чтобы переместить данные и сдвинуть нули. (Другие тасования не поддерживают обнуление вместо копирования байта с входа, и на входе, который мы могли бы скопировать, не было каких-либо известных нулей).

paddq/psubq может работать на портах 1/5 на Haswell или p015 на Skylake.

Skylake запускает pmuludq, а вектор быстрого счета сдвигается на p01, поэтому он может теоретически управлять пропускной способностью одного вектора на максимальный (5/2, 8/3, 11/4) = 11/4 = 2,75 циклы. Таким образом, это узкие места на общей пропускной способности пропущенного домена (включая 2 векторных нагрузки и 1 векторный магазин). Поэтому немного разворачивается цикл. Вероятно, конфликты ресурсов из-за несовершенного планирования будут сглаживать его до немногим меньше, чем 4 сдобренных домена на каждые часы. Накладные расходы на цикл могут, надеюсь, выполняться на порте 6, который может обрабатывать только некоторые скалярные операционные системы, включая add и compare-and-branch, оставляя порты 0/1/5 для векторных ALU ops, поскольку они близки к насыщению (8/3 = 2,666 часов). Однако порты загрузки/хранения нигде не близки к насыщению.

Таким образом, Skylake может теоретически управлять одним вектором за 2,75 цикла (плюс накладные расходы на цикл) или одним умножить на 0,7 цикла, по сравнению с наилучшим вариантом Haswell (по одному на 1,2 цикла в теории со скаляром, или один на 1,25 цикла в теории с векторами). Скалярный один за 1,2 цикла, вероятно, потребует ручной цикл asm, потому что компиляторы не знают, как использовать режим однорежимной адресации для хранилищ, и режим регистрации двух регистров для нагрузок (dst + (src-dst) и приращение dst).

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

AVX-512DQ вводит 64bx64b- > 64b векторное умножение

gcc может автоматически использовать его, если вы добавите -mavx512dq.

.L4:
    vmovdqu64       zmm0, ZMMWORD PTR [r8+rax]    # vect__11.23, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B]
    add     rcx, 1    # ivtmp.30,
    vpmullq zmm1, zmm0, zmm2  # vect__13.24, vect__11.23, vect_cst_.25
    vmovdqa64       zmm0, ZMMWORD PTR [r9+rax]    # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B]
    vpsubq  zmm0, zmm0, zmm1    # vect__14.26, MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__13.24
    vmovdqa64       ZMMWORD PTR [r9+rax], zmm0    # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__14.26
    add     rax, 64   # ivtmp.32,
    cmp     rcx, r10  # ivtmp.30, bnd.14
    jb      .L4 #,

Таким образом, AVX512DQ (ожидаемый в составе многоядерного Xeon (Purley) Skylake в ~ 2017 году) даст намного больше, чем 2x ускорение (из более широких векторов), если эти инструкции конвейерны на один за такт.

Обновление: Skylake-AVX512 (aka SKL-X или SKL-SP) запускает VPMULLQ за один раз в 1,5 цикла для векторов xmm, ymm или zmm. Это 3 часа с задержкой 15 с. (Возможно, задержка 1 с для версии zmm, если это не ошибка измерения в результатах AIDA.)

vpmullq намного быстрее, чем все, что вы можете построить из 32-битных кусков, поэтому очень важно иметь инструкцию для этого, даже если текущие процессоры не имеют аппаратное ускорение с использованием векторов с 64-разрядным элементом. (Предположительно, они используют множители мантиссы в единицах FMA.)

Ответ 2

Если вам интересны операции с SIMD 64bx64b до 64b (ниже), вы можете найти решения AVX и AVX2 от Agner Fog Vector Class Library, Я бы проверил их с массивами и посмотрел, как он сравнивается с тем, что делает GCC с общим циклом, например, в ответе .

AVX (используйте SSE - вы все еще можете скомпилировать с -mavx, чтобы получить vex-кодировку).

// vector operator * : multiply element by element
static inline Vec2q operator * (Vec2q const & a, Vec2q const & b) {
#if INSTRSET >= 5   // SSE4.1 supported
    // instruction does not exist. Split into 32-bit multiplies
    __m128i bswap   = _mm_shuffle_epi32(b,0xB1);           // b0H,b0L,b1H,b1L (swap H<->L)
    __m128i prodlh  = _mm_mullo_epi32(a,bswap);            // a0Lb0H,a0Hb0L,a1Lb1H,a1Hb1L, 32 bit L*H products
    __m128i zero    = _mm_setzero_si128();                 // 0
    __m128i prodlh2 = _mm_hadd_epi32(prodlh,zero);         // a0Lb0H+a0Hb0L,a1Lb1H+a1Hb1L,0,0
    __m128i prodlh3 = _mm_shuffle_epi32(prodlh2,0x73);     // 0, a0Lb0H+a0Hb0L, 0, a1Lb1H+a1Hb1L
    __m128i prodll  = _mm_mul_epu32(a,b);                  // a0Lb0L,a1Lb1L, 64 bit unsigned products
    __m128i prod    = _mm_add_epi64(prodll,prodlh3);       // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32
    return  prod;
#else               // SSE2
    int64_t aa[2], bb[2];
    a.store(aa);                                           // split into elements
    b.store(bb);
    return Vec2q(aa[0]*bb[0], aa[1]*bb[1]);                // multiply elements separetely
#endif
}

AVX2

// vector operator * : multiply element by element
static inline Vec4q operator * (Vec4q const & a, Vec4q const & b) {
    // instruction does not exist. Split into 32-bit multiplies
    __m256i bswap   = _mm256_shuffle_epi32(b,0xB1);           // swap H<->L
    __m256i prodlh  = _mm256_mullo_epi32(a,bswap);            // 32 bit L*H products
    __m256i zero    = _mm256_setzero_si256();                 // 0
    __m256i prodlh2 = _mm256_hadd_epi32(prodlh,zero);         // a0Lb0H+a0Hb0L,a1Lb1H+a1Hb1L,0,0
    __m256i prodlh3 = _mm256_shuffle_epi32(prodlh2,0x73);     // 0, a0Lb0H+a0Hb0L, 0, a1Lb1H+a1Hb1L
    __m256i prodll  = _mm256_mul_epu32(a,b);                  // a0Lb0L,a1Lb1L, 64 bit unsigned products
    __m256i prod    = _mm256_add_epi64(prodll,prodlh3);       // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32
    return  prod;
}

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