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

Я хочу, чтобы векторизовать умножение двух массивов с выравниванием по памяти. Я не нашел способ размножить 64 * 64 бит в AVX / AVX2, поэтому я просто сделал цикл-разворот и загрузил / хранилища AVX2. Есть ли более быстрый способ сделать это?

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

void multiply_vex(long *Gi_vec, long q, long *Gj_vec){ int i; __m256i data_j, data_i; __uint64_t *ptr_J = (__uint64_t*)&data_j; __uint64_t *ptr_I = (__uint64_t*)&data_i; for (i=0; i<BASE_VEX_STOP; i+=4) { data_i = _mm256_load_si256((__m256i*)&Gi_vec[i]); data_j = _mm256_load_si256((__m256i*)&Gj_vec[i]); ptr_I[0] -= ptr_J[0] * q; ptr_I[1] -= ptr_J[1] * q; ptr_I[2] -= ptr_J[2] * q; ptr_I[3] -= ptr_J[3] * q; _mm256_store_si256((__m256i*)&Gi_vec[i], data_i); } for (; i<BASE_DIMENSION; i++) Gi_vec[i] -= Gj_vec[i] * q; } 

UPDATE: Я использую микроархитектуру Haswell с компиляторами ICC / GCC. Таким образом, оба AVX и AVX2 в порядке. Я заменяю -= на C intrisic _mm256_sub_epi64 после цикла размножения умножения, где он получает некоторое ускорение. В настоящее время это ptr_J[0] *= q; ... ptr_J[0] *= q; ...

Я использую __uint64_t но это ошибка . Правильный тип данных – __int64_t .

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

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

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


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

 #include  // 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 

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

 .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 Босона. На основе кода библиотеки векторного classа Agner Fog .

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

Поскольку один из ваших операндов является постоянным, обязательно передайте set1(q) как b , а не a , чтобы можно было поднять «shswle» «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. (например, один вектор за 4 цикла вместо одного вектора за 5 циклов, узкий по пропускной способности порта 0. Я не рассматривал другие узкие места для полной проблемы, поскольку это было последним дополнением к ответу.)

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


Анализ Perf автосекторизуемого кода GCC (а не внутренней версии)

Справочная информация: см . Таблицы Insn Agner Fog и руководство по микроархитектуре и другие ссылки в вики-теге x86 .

До AVX512 (см. Ниже) это, вероятно, всего лишь быстрее, чем скалярный 64- imul r64, m64 код: imul r64, m64 имеет пропускную способность по одному за такт на процессорах Intel (но один на 4 такта на семействе AMD Bulldozer). load / imul / sub-with-memory-dest - это 4 процессора с интегрированными доменами на процессорах Intel (с режимом адресации, который может быть с микро-предохранителем, который gcc не может использовать). Ширина конвейера - это 4 сдобренных домена на каждом часе, поэтому даже большой разворот не может заставить его выдавать один раз в часы. При достаточном разворачивании мы будем сталкиваться с пропускной способностью загрузки / хранения. В Haswell есть 2 нагрузки и один магазин за такт, но магазины-адреса 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 / vector хранить. Итак, это 18 совпадающих доменов на 4 умножения (4.5 цикла для выпуска), но 7 перетасовки (7 циклов для выполнения). Итак, nvm, это нехорошо по сравнению с чистым скаляром.


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

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

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 такта , по сравнению с лучшим вариантом Хасуэлла (по одному на 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.)

Если вас интересуют операции SIMD 64bx64b до 64b (ниже), здесь представлены решения AVX и AVX2 из библиотеки векторных classов Agner Fog. Я бы проверил их с массивами и посмотрел, как он сравнивается с тем, что 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 является постоянным в цикле, вам не нужно пересчитывать некоторые вещи на каждую итерацию, но ваш компилятор, вероятно, все равно это увидит.