Как транспонировать матрицу 16x16 с помощью инструкций SIMD?

В настоящее время я пишу код, предназначенный для будущих инструкций Intel AVX-512 SIMD, которые поддерживают 512-битные операции.

Теперь, предполагая, что есть матрица, представленная 16 регистрами SIMD, каждый из которых содержит 16 32-битных целых чисел (соответствует строке), как я могу транспонировать матрицу с помощью чисто инструкций SIMD?

Уже есть решения для транспонирования матриц 4x4 или 8x8 с SSE и AVX2 соответственно. Но я не мог понять, как расширить его до 16x16 с помощью AVX-512.

Любые идеи?


person lei_z    schedule 08.04.2015    source источник
comment
Часто самый быстрый способ что-то сделать - вместо этого ничего не делать - по сути, дать каждой матрице транспонированный флаг и просто инвертировать этот флаг. Конечно, это означает, что вам нужно проверить транспонированный флаг и поменять местами индекс столбца и индекс строки в любом другом коде, который может иметь дело с транспонированными матрицами. Например. если у вас есть функция для добавления 2 матриц, вы можете получить 3 случая (ни один транспонированный, один транспонированный, оба транспонированные), где результатом сложения всегда является матрица, которая не транспонируется.   -  person Brendan    schedule 12.04.2015
comment
Не могли бы вы из любопытства объяснить, почему вас интересует транспонирование 16x16? Например. Это для ядра для большего транспонирования? Имеет ли значение чтение / запись для вас или это сгенерированные данные?   -  person Z boson    schedule 12.04.2015
comment
@Zboson Это часть алгоритма шифрования, который мы пытаемся оптимизировать с помощью AVX512. Фактически, мы можем использовать инструкцию сбора для транспонирования матрицы при загрузке из памяти. Но нам удалось сделать это с помощью SSE / AVX2, когда нет инструкций по сбору / разбросу, поэтому мне просто любопытно, как мы можем сделать то же самое с AVX512, то есть транспонирование в регистре.   -  person lei_z    schedule 13.04.2015
comment
@ lei.april, хорошо, в этом случае мое решение выполняет транспонирование в регистре, поэтому, когда выйдет AVX512, вы можете сравнить его с производительностью сбора / разброса. Было бы здорово, если бы производительность сбора / разброса была лучше, но я бы не стал на это рассчитывать.   -  person Z boson    schedule 13.04.2015
comment
@Zboson Некоторые приблизительные значения задержки / пропускной способности для KNL отсутствуют. Как и ожидалось, сбор / разброс все еще медленный. 2 элемента / цикл нагрузки, 1 / цикл запоминания. Итак, 8 циклов / поплавок-сборка и 16 циклов / разброс поплавка. IOW, инструкции сбора / разброса все еще разбиваются на отдельные мопы для каждого элемента и переходят в соответствующие порты. Это просто более эффективно, чем в предыдущих поколениях, когда у них была масса других дополнительных мопов.   -  person Mysticial    schedule 17.08.2016
comment
Я не видел способа добиться большего, чем это год назад. И похоже, что у Intel пока нет ничего лучше. Поскольку по сути, сбор и разброс представляют собой наборы более мелких независимых операций. Я не совсем уверен, как кеши графического процессора могут обрабатывать массивные параллельные нагрузки от сотен потоков одновременно, но у них, похоже, есть проблемы с разбросом.   -  person Mysticial    schedule 17.08.2016
comment
@Mysticial, группа HPC на работе дала мне аккаунт на своей карте Knights Landing с AVX512. Я попробовал свой код, и он сработал с первой попытки. Это хорошо знать. Я еще не проводил тестов производительности. Я получил аккаунт около 30 минут назад.   -  person Z boson    schedule 19.12.2016
comment
@Zboson Тебе повезло. Я вряд ли получу какое-либо оборудование AVX512 до Skylake Purley.   -  person Mysticial    schedule 19.12.2016
comment
@ lei.april Я только что добавил новый ответ на этот вопрос, основанный на тестах с реальным оборудованием AVX512. У вас есть оборудование AVX512, чтобы это проверить. Если да, то какие результаты вы получили? Я получаю, что сбор все еще медленнее, чем мой метод.   -  person Z boson    schedule 21.12.2016


Ответы (2)


Для двух инструкций операндов с использованием SIMD вы можете показать, что количество операций, необходимых для транспонирования матрицы nxn, равно n*log_2(n), тогда как при использовании скалярных операций это O(n^2). Фактически, позже я покажу, что количество операций чтения и записи с использованием скалярных регистров равно 2*n*(n-1). Ниже приведена таблица, показывающая количество операций транспонирования матриц 4x4, 8x8, 16x16 и 32x32 с использованием SSE, AVX, AVX512 и AVX1024 по сравнению со скалярными операциями.

n            4(SSE)          8(AVX)    16(AVX512)    32(AVX1024)  
SIMD ops          8              24           64            160
SIMD +r/w ops    16              40           96            224     
Scalar r/w ops   24             112          480           1984

где SIMD + r / w ops включает операции чтения и записи (n*log_2(n) + 2*n).

Причина, по которой транспонирование SIMD может быть выполнено в n*log_2(n) операциях, заключается в том, что алгоритм:

permute n 32-bit rows
permute n 64-bit rows
...
permute n simd_width/2-bit rows

Например, для 4x4 есть 4 строки, и поэтому вам нужно переставить 32-битные полосы 4 раза, а затем 64-битные полосы 4 раза. Для 16x16 вам нужно переставить 32-битные полосы, 64-битные полосы, 128-битные полосы и, наконец, 256 полос, по 16 раз для каждой.

Я уже показал, что 8x8 можно выполнить с помощью 24 операций с AVX < / а>. Итак, вопрос в том, как это сделать для 16x16 с использованием AVX512 в 64 операциях? Общий алгоритм:

interleave 32-bit lanes using 
    8x _mm512_unpacklo_epi32
    8x _mm512_unpackhi_epi32
interleave 64-bit lanes using
    8x _mm512_unpacklo_epi64 
    8x _mm512_unpackhi_epi64 
permute 128-bit lanes using
   16x _mm512_shuffle_i32x4
permute 256-bit lanes using again
   16x _mm512_shuffle_i32x4

Вот непроверенный код, делающий это

    //given __m512i r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, ra, rb, rc, rd, re, rf;
    __m512i t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, ta, tb, tc, td, te, tf;

    t0 = _mm512_unpacklo_epi32(r0,r1); //   0  16   1  17   4  20   5  21   8  24   9  25  12  28  13  29 
    t1 = _mm512_unpackhi_epi32(r0,r1); //   2  18   3  19   6  22   7  23  10  26  11  27  14  30  15  31
    t2 = _mm512_unpacklo_epi32(r2,r3); //  32  48  33  49 ...
    t3 = _mm512_unpackhi_epi32(r2,r3); //  34  50  35  51 ...
    t4 = _mm512_unpacklo_epi32(r4,r5); //  64  80  65  81 ...  
    t5 = _mm512_unpackhi_epi32(r4,r5); //  66  82  67  83 ...
    t6 = _mm512_unpacklo_epi32(r6,r7); //  96 112  97 113 ...
    t7 = _mm512_unpackhi_epi32(r6,r7); //  98 114  99 115 ...
    t8 = _mm512_unpacklo_epi32(r8,r9); // 128 ...
    t9 = _mm512_unpackhi_epi32(r8,r9); // 130 ...
    ta = _mm512_unpacklo_epi32(ra,rb); // 160 ...
    tb = _mm512_unpackhi_epi32(ra,rb); // 162 ...
    tc = _mm512_unpacklo_epi32(rc,rd); // 196 ...
    td = _mm512_unpackhi_epi32(rc,rd); // 198 ...
    te = _mm512_unpacklo_epi32(re,rf); // 228 ...
    tf = _mm512_unpackhi_epi32(re,rf); // 230 ...

    r0 = _mm512_unpacklo_epi64(t0,t2); //   0  16  32  48 ...
    r1 = _mm512_unpackhi_epi64(t0,t2); //   1  17  33  49 ...
    r2 = _mm512_unpacklo_epi64(t1,t3); //   2  18  34  49 ...
    r3 = _mm512_unpackhi_epi64(t1,t3); //   3  19  35  51 ...
    r4 = _mm512_unpacklo_epi64(t4,t6); //  64  80  96 112 ...  
    r5 = _mm512_unpackhi_epi64(t4,t6); //  65  81  97 114 ...
    r6 = _mm512_unpacklo_epi64(t5,t7); //  66  82  98 113 ...
    r7 = _mm512_unpackhi_epi64(t5,t7); //  67  83  99 115 ...
    r8 = _mm512_unpacklo_epi64(t8,ta); // 128 144 160 176 ...  
    r9 = _mm512_unpackhi_epi64(t8,ta); // 129 145 161 178 ...
    ra = _mm512_unpacklo_epi64(t9,tb); // 130 146 162 177 ... 
    rb = _mm512_unpackhi_epi64(t9,tb); // 131 147 163 179 ...
    rc = _mm512_unpacklo_epi64(tc,te); // 192 208 228 240 ... 
    rd = _mm512_unpackhi_epi64(tc,te); // 193 209 229 241 ...
    re = _mm512_unpacklo_epi64(td,tf); // 194 210 230 242 ...
    rf = _mm512_unpackhi_epi64(td,tf); // 195 211 231 243 ...

    t0 = _mm512_shuffle_i32x4(r0, r4, 0x88); //   0  16  32  48   8  24  40  56  64  80  96  112 ...
    t1 = _mm512_shuffle_i32x4(r1, r5, 0x88); //   1  17  33  49 ...
    t2 = _mm512_shuffle_i32x4(r2, r6, 0x88); //   2  18  34  50 ...
    t3 = _mm512_shuffle_i32x4(r3, r7, 0x88); //   3  19  35  51 ...
    t4 = _mm512_shuffle_i32x4(r0, r4, 0xdd); //   4  20  36  52 ...
    t5 = _mm512_shuffle_i32x4(r1, r5, 0xdd); //   5  21  37  53 ...
    t6 = _mm512_shuffle_i32x4(r2, r6, 0xdd); //   6  22  38  54 ...
    t7 = _mm512_shuffle_i32x4(r3, r7, 0xdd); //   7  23  39  55 ...
    t8 = _mm512_shuffle_i32x4(r8, rc, 0x88); // 128 144 160 176 ...
    t9 = _mm512_shuffle_i32x4(r9, rd, 0x88); // 129 145 161 177 ...
    ta = _mm512_shuffle_i32x4(ra, re, 0x88); // 130 146 162 178 ...
    tb = _mm512_shuffle_i32x4(rb, rf, 0x88); // 131 147 163 179 ...
    tc = _mm512_shuffle_i32x4(r8, rc, 0xdd); // 132 148 164 180 ...
    td = _mm512_shuffle_i32x4(r9, rd, 0xdd); // 133 149 165 181 ...
    te = _mm512_shuffle_i32x4(ra, re, 0xdd); // 134 150 166 182 ...
    tf = _mm512_shuffle_i32x4(rb, rf, 0xdd); // 135 151 167 183 ...

    r0 = _mm512_shuffle_i32x4(t0, t8, 0x88); //   0  16  32  48  64  80  96 112 ... 240
    r1 = _mm512_shuffle_i32x4(t1, t9, 0x88); //   1  17  33  49  66  81  97 113 ... 241
    r2 = _mm512_shuffle_i32x4(t2, ta, 0x88); //   2  18  34  50  67  82  98 114 ... 242
    r3 = _mm512_shuffle_i32x4(t3, tb, 0x88); //   3  19  35  51  68  83  99 115 ... 243
    r4 = _mm512_shuffle_i32x4(t4, tc, 0x88); //   4 ...
    r5 = _mm512_shuffle_i32x4(t5, td, 0x88); //   5 ...
    r6 = _mm512_shuffle_i32x4(t6, te, 0x88); //   6 ...
    r7 = _mm512_shuffle_i32x4(t7, tf, 0x88); //   7 ...
    r8 = _mm512_shuffle_i32x4(t0, t8, 0xdd); //   8 ...
    r9 = _mm512_shuffle_i32x4(t1, t9, 0xdd); //   9 ...
    ra = _mm512_shuffle_i32x4(t2, ta, 0xdd); //  10 ...
    rb = _mm512_shuffle_i32x4(t3, tb, 0xdd); //  11 ...
    rc = _mm512_shuffle_i32x4(t4, tc, 0xdd); //  12 ...
    rd = _mm512_shuffle_i32x4(t5, td, 0xdd); //  13 ...
    re = _mm512_shuffle_i32x4(t6, te, 0xdd); //  14 ...
    rf = _mm512_shuffle_i32x4(t7, tf, 0xdd); //  15  31  47  63  79  96 111 127 ... 255

У меня возникла идея использовать _mm512_shufflei32x4, посмотрев на транспонирование матрицы 4x4 с помощью _mm_shuffle_ps (это то, что MSVC использует в _MM_TRANSPOSE4_PS, но не GCC и ICC).

__m128 tmp0 ,tmp1, tmp2, tmp3;
tmp0 = _mm_shuffle_ps(row0, row1, 0x88); // 0 2 4 6
tmp1 = _mm_shuffle_ps(row0, row1, 0xdd); // 1 3 5 7
tmp2 = _mm_shuffle_ps(row2, row3, 0x88); // 8 a c e
tmp3 = _mm_shuffle_ps(row2, row3, 0xdd); // 9 b d f

row0 = _mm_shuffle_ps(tmp0, tmp2, 0x88); // 0 4 8 c 
row1 = _mm_shuffle_ps(tmp1, tmp3, 0x88); // 1 5 9 d
row2 = _mm_shuffle_ps(tmp0, tmp2, 0xdd); // 2 6 a e 
row3 = _mm_shuffle_ps(tmp1, tmp3, 0xdd); // 3 7 b f

та же идея применима к _mm512_shuffle_i32x4, но теперь полосы 128-битные вместо 32-битных и 16 строк вместо 4.

Наконец, для сравнения со скалярными операциями я модифицировал пример 9.5a из руководства по оптимизации C ++ Агнера Фога.

#define SIZE 16
void transpose(int a[SIZE][SIZE]) { // function to transpose matrix
    // define a macro to swap two array elements:
    #define swapd(x,y) {temp=x; x=y; y=temp;}
    int r, c; int temp;
    for (r = 1; r < SIZE; r++) {
        for (c = 0; c < r; c++) {
            swapd(a[r][c], a[c][r]);
        }
    }
}

это делает n*(n-1)/2 свопы (потому что диагональ менять местами не нужно). Свапы из сборки на 16х16 выглядят так

mov     r8d, DWORD PTR [rax+68]
mov     r9d, DWORD PTR [rdx+68]
mov     DWORD PTR [rax+68], r9d
mov     DWORD PTR [rdx+68], r8d

поэтому количество операций чтения / записи с использованием скалярных регистров составляет 2*n*(n-1).

person Z boson    schedule 12.04.2015
comment
1, как бы уродливо это ни выглядело, это, вероятно, все же будет быстрее, чем использование 16 сборочных нагрузок. - person Mysticial; 12.04.2015
comment
@Mysticial, правда ли, что только процессоры xeon и Skylake для рабочих станций будут иметь AVX512? Если это так, то что за # @ $! есть смысл Skylake ??? Это очень неутешительные новости, если это правда. Что делает Skylake крутым без AVX512? - person Z boson; 01.06.2015
comment
Да, я не осознавал, что все было так плохо, до недавней утечки информации о Перли. Похоже, что это будет Knights Landing в первом квартале 2016 года и Skylake Xeon с AVX512 в (конце?) 2017 года. Процессоры Intel обычно делятся на ноутбуки / недорогие настольные компьютеры (сокет 115x) и серверы / высокопроизводительные настольные компьютеры. (сокет 2011-х) линии. Похоже, что AVX512 для Skylake будет только в линейке серверов / высокопроизводительных настольных компьютеров для Skylake. Это потенциально позже, чем Cannonlake для ноутбуков / недорогих настольных компьютеров. - person Mysticial; 01.06.2015
comment
Конечно, я делаю эти предположения, основываясь на недавних утечках, а также на моих (ограниченных) знаниях о линейках продуктов Intel. Так что я определенно мог ошибаться. Вероятно, будет Xeon Skylake для сокета 1151, который выйдет в третьем квартале 2015 года. Но, вероятно, это просто прославленный процессор для настольных ПК, поэтому я не уверен, что у него будет AVX512. - person Mysticial; 01.06.2015
comment
@ Таинственный, ладно, думаю, я понимаю. Ток в Skylake должен быть для IGP. Но у IGP, похоже, нет цикличного цикла. Он значительно улучшается с каждым поколением. Мне действительно не нравится эта тенденция, когда все больше и больше потребительских процессоров Intel становятся IGP (кремний моего процессора Haswell более чем на 50% состоит из графического процессора). Графические процессоры имеют закрытый исходный код и обычно не программируются напрямую, а через драйверы устройств. Не только это, но двойная точность - отстой. Вот что происходит, когда AMD не предлагает конкуренции. - person Z boson; 01.06.2015
comment
С точки зрения видимых для пользователя функций, Skylake не так уж и много. С архитектурной точки зрения, это определенно огромный редизайн конвейера, чтобы освободить место для AVX512, которого они пока не поддерживают. AVX512 не так прост, как расширение исполнительных блоков. Такие вещи, как регистры маски и мопы зависимостей с 3 входами, должны быть встроены в конвейер. - person Mysticial; 01.06.2015
comment
Что касается того, почему они еще не поддерживают AVX512, у меня есть несколько предположений: 1) На 14-нм рабочем столе не хватает кремниевой области. 2) Занимает такую ​​большую площадь, что влияет на урожайность. 3) Intel может намеренно затягивать, чтобы дать AMD время наверстать упущенное или быть разорванным антимонопольным комитетом. - person Mysticial; 01.06.2015
comment
Клиент @Zboson Skylake - это новый uarch, что и означает tock (intel.com/content/www/us/en/silicon-innovations/). Произошел значительный рост производительности графического процессора (cdn.arstechnica.net /wp-content/uploads/2015/08/02.png), что гораздо важнее в клиентском пространстве, чем AVX-512. AVX-512 в первую очередь представляет интерес для специалистов по HPC, которые собираются покупать процессоры Xeon или Xeon Phi. - person Jeff Hammond; 23.12.2015
comment
@Jeff, AVX512 меня интересует, потому что я написал трассировщик лучей в реальном времени (а не только трассировщик лучей) в OpenCL, и он работает как минимум в пять раз быстрее на моем шестилетнем GTX 580, чем на всех процессорах Intel, которые я использовал до сих пор. - person Z boson; 23.12.2015
comment
Ради интереса, кто-то должен сделать транспонирование размером 64x64 байта. Это точно выполнимо с AVX512-VBMI. Может быть, и с AVX512-BW. Не то чтобы я пробовал, да и пользы от этого нет. - person Mysticial; 17.08.2016
comment
@Mysticial, у кого есть оборудование AVX512? Еще есть в наличии? - person Z boson; 17.08.2016
comment
@Mysticial, вы можете найти это интересно. - person Z boson; 17.08.2016
comment
Я уже видел ту статью. Ага, интересно. KNL доступен, но AFAICT, только через платформу разработчика Ninja для потребителей. Я не вижу этого нигде, кроме университетов и тому подобного. - person Mysticial; 17.08.2016
comment
@Mysticial Я заметил, что вы обновили тег AVX512. Возможно, вы захотите упомянуть, что будут некоторые расширения инструкций по глубокому обучению для AVX512, а также AVX512_4VNNIW и AVX512_4FMAPS lemire.me/blog/2016/10/14/. Я только что устроился на работу в сфере машинного обучения / искусственного интеллекта и больших данных. - person Z boson; 18.10.2016
comment
@Zboson Woah wtf. Для меня это в новинку. И они безумны? Их рот теперь на 3 поколения опережает их действия. - person Mysticial; 18.10.2016
comment
Похоже, они еще не опубликовали никаких подробностей реальных инструкций. - person Mysticial; 18.10.2016
comment
Кстати, перестановка / перемешивание Knights Landing, которые извлекают из двух векторов вместо одного, имеют половину пропускной способности. У меня нет оборудования для его тестирования, но я думаю, что может быть быстрее использовать альтернативу, например: _mm512_unpacklo_epi64(a, b) -> _mm512_mask_permutex_epi64(a, 0xaa, b, 177) или _mm512_shuffle_i64x2(a, b, 68) -> _mm512_inserti64x4(a, _mm512_castsi512_si256(b), 1) - person Mysticial; 20.01.2017
comment
Этот первый пример, в частности, имеет недостаток, заключающийся в том, что он требует маски и является деструктивным, поэтому вам понадобится ход reg-reg для каждой пары. Но он распространяется на множество различных комбинаций перемешивания. Можно построить транспонирование 8x8 или 16x16 только с пермутацией / перемешиванием одной ширины, которые имеют пропускную способность 1 / цикл на KNL. Потенциальное улучшение составляет около 25%, поскольку соотношение работы падает с 4 до 3, если у вас есть все предварительно загруженные векторы маски и перестановки. - person Mysticial; 20.01.2017
comment
@Mysticial, я сделал нечто подобное в моем последнем обновлении здесь в конце, но я использовал permutexvar_ вместо permutex_. Я попробую ваши предложения, если у меня будет время на этой неделе. - person Z boson; 23.01.2017
comment
Как ваши формулы масштабируются, когда размер строки больше ширины вектора? Например, транспонирование 8x8, но только с использованием AVX2, а не AVX-512. Или большие транспозиции, где даже любой будущий набор инструкций не будет достаточно широким, чтобы удерживать строку (я не задерживаю дыхание для AVX-65536)? - person BeeOnRope; 18.10.2019

Недавно я получил доступ к оборудованию Xeon Phi Knights Landing с AVX512. В частности, я использую процессор Intel (R) Xeon Phi (TM) 7250 @ 1,40 ГГц (http://ark.intel.com/products/94035/Intel-Xeon-Phi-Processor-7250-16GB-1_40-GHz-68-core). Это не вспомогательная карта. Xeon Phi - это главный компьютер.

Я протестировал инструкции по сбору AVX512 по сравнению с моим методом здесь https://stackoverflow.com/a/29587984/2542702 и это похоже, что сборка все еще медленнее. Мой код в этом ответе сработал с первой попытки без ошибок.

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

Я тестировал только ICC 17.0.0, потому что в настоящее время установленная ОС - это только CentOS 7.2 с ядром Linux 3.10 и GCC 4.8.5, а GCC 4.8 не поддерживает AVX512. Я могу убедить группу HPC на моей работе обновиться.

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

//icc -O3 -xCOMMON-AVX512 tran.c -fopenmp
#include <stdio.h>
#include <x86intrin.h>
#include <omp.h>    

void tran(int* mat, int* matT) {
    int i,j;

    __m512i t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, ta, tb, tc, td, te, tf;
    __m512i r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, ra, rb, rc, rd, re, rf;

    r0 = _mm512_load_epi32(&mat[ 0*16]);
    r1 = _mm512_load_epi32(&mat[ 1*16]);
    r2 = _mm512_load_epi32(&mat[ 2*16]);
    r3 = _mm512_load_epi32(&mat[ 3*16]);
    r4 = _mm512_load_epi32(&mat[ 4*16]);
    r5 = _mm512_load_epi32(&mat[ 5*16]);
    r6 = _mm512_load_epi32(&mat[ 6*16]);
    r7 = _mm512_load_epi32(&mat[ 7*16]);
    r8 = _mm512_load_epi32(&mat[ 8*16]);
    r9 = _mm512_load_epi32(&mat[ 9*16]);
    ra = _mm512_load_epi32(&mat[10*16]);
    rb = _mm512_load_epi32(&mat[11*16]);
    rc = _mm512_load_epi32(&mat[12*16]);
    rd = _mm512_load_epi32(&mat[13*16]);
    re = _mm512_load_epi32(&mat[14*16]);
    rf = _mm512_load_epi32(&mat[15*16]);

    t0 = _mm512_unpacklo_epi32(r0,r1); //   0  16   1  17   4  20   5  21   8  24   9  25  12  28  13  29 
    t1 = _mm512_unpackhi_epi32(r0,r1); //   2  18   3  19   6  22   7  23  10  26  11  27  14  30  15  31
    t2 = _mm512_unpacklo_epi32(r2,r3); //  32  48  33  49 ...
    t3 = _mm512_unpackhi_epi32(r2,r3); //  34  50  35  51 ...
    t4 = _mm512_unpacklo_epi32(r4,r5); //  64  80  65  81 ...  
    t5 = _mm512_unpackhi_epi32(r4,r5); //  66  82  67  83 ...
    t6 = _mm512_unpacklo_epi32(r6,r7); //  96 112  97 113 ...
    t7 = _mm512_unpackhi_epi32(r6,r7); //  98 114  99 115 ...
    t8 = _mm512_unpacklo_epi32(r8,r9); // 128 ...
    t9 = _mm512_unpackhi_epi32(r8,r9); // 130 ...
    ta = _mm512_unpacklo_epi32(ra,rb); // 160 ...
    tb = _mm512_unpackhi_epi32(ra,rb); // 162 ...
    tc = _mm512_unpacklo_epi32(rc,rd); // 196 ...
    td = _mm512_unpackhi_epi32(rc,rd); // 198 ...
    te = _mm512_unpacklo_epi32(re,rf); // 228 ...
    tf = _mm512_unpackhi_epi32(re,rf); // 230 ...

    r0 = _mm512_unpacklo_epi64(t0,t2); //   0  16  32  48 ...
    r1 = _mm512_unpackhi_epi64(t0,t2); //   1  17  33  49 ...
    r2 = _mm512_unpacklo_epi64(t1,t3); //   2  18  34  49 ...
    r3 = _mm512_unpackhi_epi64(t1,t3); //   3  19  35  51 ...
    r4 = _mm512_unpacklo_epi64(t4,t6); //  64  80  96 112 ...  
    r5 = _mm512_unpackhi_epi64(t4,t6); //  65  81  97 114 ...
    r6 = _mm512_unpacklo_epi64(t5,t7); //  66  82  98 113 ...
    r7 = _mm512_unpackhi_epi64(t5,t7); //  67  83  99 115 ...
    r8 = _mm512_unpacklo_epi64(t8,ta); // 128 144 160 176 ...  
    r9 = _mm512_unpackhi_epi64(t8,ta); // 129 145 161 178 ...
    ra = _mm512_unpacklo_epi64(t9,tb); // 130 146 162 177 ... 
    rb = _mm512_unpackhi_epi64(t9,tb); // 131 147 163 179 ...
    rc = _mm512_unpacklo_epi64(tc,te); // 192 208 228 240 ... 
    rd = _mm512_unpackhi_epi64(tc,te); // 193 209 229 241 ...
    re = _mm512_unpacklo_epi64(td,tf); // 194 210 230 242 ...
    rf = _mm512_unpackhi_epi64(td,tf); // 195 211 231 243 ...

    t0 = _mm512_shuffle_i32x4(r0, r4, 0x88); //   0  16  32  48   8  24  40  56  64  80  96  112 ...
    t1 = _mm512_shuffle_i32x4(r1, r5, 0x88); //   1  17  33  49 ...
    t2 = _mm512_shuffle_i32x4(r2, r6, 0x88); //   2  18  34  50 ...
    t3 = _mm512_shuffle_i32x4(r3, r7, 0x88); //   3  19  35  51 ...
    t4 = _mm512_shuffle_i32x4(r0, r4, 0xdd); //   4  20  36  52 ...
    t5 = _mm512_shuffle_i32x4(r1, r5, 0xdd); //   5  21  37  53 ...
    t6 = _mm512_shuffle_i32x4(r2, r6, 0xdd); //   6  22  38  54 ...
    t7 = _mm512_shuffle_i32x4(r3, r7, 0xdd); //   7  23  39  55 ...
    t8 = _mm512_shuffle_i32x4(r8, rc, 0x88); // 128 144 160 176 ...
    t9 = _mm512_shuffle_i32x4(r9, rd, 0x88); // 129 145 161 177 ...
    ta = _mm512_shuffle_i32x4(ra, re, 0x88); // 130 146 162 178 ...
    tb = _mm512_shuffle_i32x4(rb, rf, 0x88); // 131 147 163 179 ...
    tc = _mm512_shuffle_i32x4(r8, rc, 0xdd); // 132 148 164 180 ...
    td = _mm512_shuffle_i32x4(r9, rd, 0xdd); // 133 149 165 181 ...
    te = _mm512_shuffle_i32x4(ra, re, 0xdd); // 134 150 166 182 ...
    tf = _mm512_shuffle_i32x4(rb, rf, 0xdd); // 135 151 167 183 ...

    r0 = _mm512_shuffle_i32x4(t0, t8, 0x88); //   0  16  32  48  64  80  96 112 ... 240
    r1 = _mm512_shuffle_i32x4(t1, t9, 0x88); //   1  17  33  49  66  81  97 113 ... 241
    r2 = _mm512_shuffle_i32x4(t2, ta, 0x88); //   2  18  34  50  67  82  98 114 ... 242
    r3 = _mm512_shuffle_i32x4(t3, tb, 0x88); //   3  19  35  51  68  83  99 115 ... 243
    r4 = _mm512_shuffle_i32x4(t4, tc, 0x88); //   4 ...
    r5 = _mm512_shuffle_i32x4(t5, td, 0x88); //   5 ...
    r6 = _mm512_shuffle_i32x4(t6, te, 0x88); //   6 ...
    r7 = _mm512_shuffle_i32x4(t7, tf, 0x88); //   7 ...
    r8 = _mm512_shuffle_i32x4(t0, t8, 0xdd); //   8 ...
    r9 = _mm512_shuffle_i32x4(t1, t9, 0xdd); //   9 ...
    ra = _mm512_shuffle_i32x4(t2, ta, 0xdd); //  10 ...
    rb = _mm512_shuffle_i32x4(t3, tb, 0xdd); //  11 ...
    rc = _mm512_shuffle_i32x4(t4, tc, 0xdd); //  12 ...
    rd = _mm512_shuffle_i32x4(t5, td, 0xdd); //  13 ...
    re = _mm512_shuffle_i32x4(t6, te, 0xdd); //  14 ...
    rf = _mm512_shuffle_i32x4(t7, tf, 0xdd); //  15  31  47  63  79  96 111 127 ... 255

    _mm512_store_epi32(&matT[ 0*16], r0);
    _mm512_store_epi32(&matT[ 1*16], r1);
    _mm512_store_epi32(&matT[ 2*16], r2);
    _mm512_store_epi32(&matT[ 3*16], r3);
    _mm512_store_epi32(&matT[ 4*16], r4);
    _mm512_store_epi32(&matT[ 5*16], r5);
    _mm512_store_epi32(&matT[ 6*16], r6);
    _mm512_store_epi32(&matT[ 7*16], r7);
    _mm512_store_epi32(&matT[ 8*16], r8);
    _mm512_store_epi32(&matT[ 9*16], r9);
    _mm512_store_epi32(&matT[10*16], ra);
    _mm512_store_epi32(&matT[11*16], rb);
    _mm512_store_epi32(&matT[12*16], rc);
    _mm512_store_epi32(&matT[13*16], rd);
    _mm512_store_epi32(&matT[14*16], re);
    _mm512_store_epi32(&matT[15*16], rf);
}

void gather(int *mat, int *matT) {
    int i,j;
    int index[16] __attribute__((aligned(64)));

    __m512i vindex;

    for(i=0; i<16; i++) index[i] = 16*i;
    for(i=0; i<256; i++) mat[i] = i;
    vindex = _mm512_load_epi32(index);

    for(i=0; i<16; i++) 
    _mm512_store_epi32(&matT[16*i], _mm512_i32gather_epi32(vindex, &mat[i], 4));
}

int verify(int *mat) {
    int i,j;
    int error = 0;
    for(i=0; i<16; i++) {
      for(j=0; j<16; j++) {
        if(mat[j*16+i] != i*16+j) error++;
      }
    }
    return error;
}

void print_mat(int *mat) {
    int i,j;
    for(i=0; i<16; i++) {
      for(j=0; j<16; j++) printf("%2X ", mat[i*16+j]);
      puts("");
    }
    puts("");
}

int main(void) {
    int i,j, rep;
    int mat[256] __attribute__((aligned(64)));
    int matT[256] __attribute__((aligned(64)));
    double dtime;

    rep = 10000000;
    for(i=0; i<256; i++) mat[i] = i;
    print_mat(mat);

    gather(mat, matT);
    for(i=0; i<256; i++) mat[i] = i;
    dtime = -omp_get_wtime();
    for(i=0; i<rep; i++) gather(mat, matT);
    dtime += omp_get_wtime();
    printf("errors %d\n", verify(matT));
    printf("dtime %f\n", dtime);
    print_mat(matT);

    tran(mat,matT);
    dtime = -omp_get_wtime();
    for(i=0; i<rep; i++) tran(mat, matT);
    dtime += omp_get_wtime();
    printf("errors %d\n", verify(matT));
    printf("dtime %f\n", dtime);
    print_mat(matT);
}

Функция gather в этом случае занимает 1,5 с, а функция tran 1,15 с. Если кто-нибудь увидит ошибку или у вас есть предложения по моему тесту, дайте мне знать. Я только начинаю знакомиться с AVX512 и Knights Landing.


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

#include <stdio.h>
#include <x86intrin.h>
#include <omp.h>   

void tran(int* mat, int* matT, int rep) {
    int i;

    __m512i t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, ta, tb, tc, td, te, tf;
    __m512i r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, ra, rb, rc, rd, re, rf;

    for(i=0; i<rep; i++) {

    r0 = _mm512_load_epi32(&mat[ 0*16]);
    r1 = _mm512_load_epi32(&mat[ 1*16]);
    r2 = _mm512_load_epi32(&mat[ 2*16]);
    r3 = _mm512_load_epi32(&mat[ 3*16]);
    r4 = _mm512_load_epi32(&mat[ 4*16]);
    r5 = _mm512_load_epi32(&mat[ 5*16]);
    r6 = _mm512_load_epi32(&mat[ 6*16]);
    r7 = _mm512_load_epi32(&mat[ 7*16]);
    r8 = _mm512_load_epi32(&mat[ 8*16]);
    r9 = _mm512_load_epi32(&mat[ 9*16]);
    ra = _mm512_load_epi32(&mat[10*16]);
    rb = _mm512_load_epi32(&mat[11*16]);
    rc = _mm512_load_epi32(&mat[12*16]);
    rd = _mm512_load_epi32(&mat[13*16]);
    re = _mm512_load_epi32(&mat[14*16]);
    rf = _mm512_load_epi32(&mat[15*16]);

    t0 = _mm512_unpacklo_epi32(r0,r1); //   0  16   1  17   4  20   5  21   8  24   9  25  12  28  13  29 
    t1 = _mm512_unpackhi_epi32(r0,r1); //   2  18   3  19   6  22   7  23  10  26  11  27  14  30  15  31
    t2 = _mm512_unpacklo_epi32(r2,r3); //  32  48  33  49 ...
    t3 = _mm512_unpackhi_epi32(r2,r3); //  34  50  35  51 ...
    t4 = _mm512_unpacklo_epi32(r4,r5); //  64  80  65  81 ...  
    t5 = _mm512_unpackhi_epi32(r4,r5); //  66  82  67  83 ...
    t6 = _mm512_unpacklo_epi32(r6,r7); //  96 112  97 113 ...
    t7 = _mm512_unpackhi_epi32(r6,r7); //  98 114  99 115 ...
    t8 = _mm512_unpacklo_epi32(r8,r9); // 128 ...
    t9 = _mm512_unpackhi_epi32(r8,r9); // 130 ...
    ta = _mm512_unpacklo_epi32(ra,rb); // 160 ...
    tb = _mm512_unpackhi_epi32(ra,rb); // 162 ...
    tc = _mm512_unpacklo_epi32(rc,rd); // 196 ...
    td = _mm512_unpackhi_epi32(rc,rd); // 198 ...
    te = _mm512_unpacklo_epi32(re,rf); // 228 ...
    tf = _mm512_unpackhi_epi32(re,rf); // 230 ...

    r0 = _mm512_unpacklo_epi64(t0,t2); //   0  16  32  48 ...
    r1 = _mm512_unpackhi_epi64(t0,t2); //   1  17  33  49 ...
    r2 = _mm512_unpacklo_epi64(t1,t3); //   2  18  34  49 ...
    r3 = _mm512_unpackhi_epi64(t1,t3); //   3  19  35  51 ...
    r4 = _mm512_unpacklo_epi64(t4,t6); //  64  80  96 112 ...  
    r5 = _mm512_unpackhi_epi64(t4,t6); //  65  81  97 114 ...
    r6 = _mm512_unpacklo_epi64(t5,t7); //  66  82  98 113 ...
    r7 = _mm512_unpackhi_epi64(t5,t7); //  67  83  99 115 ...
    r8 = _mm512_unpacklo_epi64(t8,ta); // 128 144 160 176 ...  
    r9 = _mm512_unpackhi_epi64(t8,ta); // 129 145 161 178 ...
    ra = _mm512_unpacklo_epi64(t9,tb); // 130 146 162 177 ... 
    rb = _mm512_unpackhi_epi64(t9,tb); // 131 147 163 179 ...
    rc = _mm512_unpacklo_epi64(tc,te); // 192 208 228 240 ... 
    rd = _mm512_unpackhi_epi64(tc,te); // 193 209 229 241 ...
    re = _mm512_unpacklo_epi64(td,tf); // 194 210 230 242 ...
    rf = _mm512_unpackhi_epi64(td,tf); // 195 211 231 243 ...

    t0 = _mm512_shuffle_i32x4(r0, r4, 0x88); //   0  16  32  48   8  24  40  56  64  80  96  112 ...
    t1 = _mm512_shuffle_i32x4(r1, r5, 0x88); //   1  17  33  49 ...
    t2 = _mm512_shuffle_i32x4(r2, r6, 0x88); //   2  18  34  50 ...
    t3 = _mm512_shuffle_i32x4(r3, r7, 0x88); //   3  19  35  51 ...
    t4 = _mm512_shuffle_i32x4(r0, r4, 0xdd); //   4  20  36  52 ...
    t5 = _mm512_shuffle_i32x4(r1, r5, 0xdd); //   5  21  37  53 ...
    t6 = _mm512_shuffle_i32x4(r2, r6, 0xdd); //   6  22  38  54 ...
    t7 = _mm512_shuffle_i32x4(r3, r7, 0xdd); //   7  23  39  55 ...
    t8 = _mm512_shuffle_i32x4(r8, rc, 0x88); // 128 144 160 176 ...
    t9 = _mm512_shuffle_i32x4(r9, rd, 0x88); // 129 145 161 177 ...
    ta = _mm512_shuffle_i32x4(ra, re, 0x88); // 130 146 162 178 ...
    tb = _mm512_shuffle_i32x4(rb, rf, 0x88); // 131 147 163 179 ...
    tc = _mm512_shuffle_i32x4(r8, rc, 0xdd); // 132 148 164 180 ...
    td = _mm512_shuffle_i32x4(r9, rd, 0xdd); // 133 149 165 181 ...
    te = _mm512_shuffle_i32x4(ra, re, 0xdd); // 134 150 166 182 ...
    tf = _mm512_shuffle_i32x4(rb, rf, 0xdd); // 135 151 167 183 ...

    r0 = _mm512_shuffle_i32x4(t0, t8, 0x88); //   0  16  32  48  64  80  96 112 ... 240
    r1 = _mm512_shuffle_i32x4(t1, t9, 0x88); //   1  17  33  49  66  81  97 113 ... 241
    r2 = _mm512_shuffle_i32x4(t2, ta, 0x88); //   2  18  34  50  67  82  98 114 ... 242
    r3 = _mm512_shuffle_i32x4(t3, tb, 0x88); //   3  19  35  51  68  83  99 115 ... 243
    r4 = _mm512_shuffle_i32x4(t4, tc, 0x88); //   4 ...
    r5 = _mm512_shuffle_i32x4(t5, td, 0x88); //   5 ...
    r6 = _mm512_shuffle_i32x4(t6, te, 0x88); //   6 ...
    r7 = _mm512_shuffle_i32x4(t7, tf, 0x88); //   7 ...
    r8 = _mm512_shuffle_i32x4(t0, t8, 0xdd); //   8 ...
    r9 = _mm512_shuffle_i32x4(t1, t9, 0xdd); //   9 ...
    ra = _mm512_shuffle_i32x4(t2, ta, 0xdd); //  10 ...
    rb = _mm512_shuffle_i32x4(t3, tb, 0xdd); //  11 ...
    rc = _mm512_shuffle_i32x4(t4, tc, 0xdd); //  12 ...
    rd = _mm512_shuffle_i32x4(t5, td, 0xdd); //  13 ...
    re = _mm512_shuffle_i32x4(t6, te, 0xdd); //  14 ...
    rf = _mm512_shuffle_i32x4(t7, tf, 0xdd); //  15  31  47  63  79  96 111 127 ... 255

    _mm512_store_epi32(&matT[ 0*16], r0);
    _mm512_store_epi32(&matT[ 1*16], r1);
    _mm512_store_epi32(&matT[ 2*16], r2);
    _mm512_store_epi32(&matT[ 3*16], r3);
    _mm512_store_epi32(&matT[ 4*16], r4);
    _mm512_store_epi32(&matT[ 5*16], r5);
    _mm512_store_epi32(&matT[ 6*16], r6);
    _mm512_store_epi32(&matT[ 7*16], r7);
    _mm512_store_epi32(&matT[ 8*16], r8);
    _mm512_store_epi32(&matT[ 9*16], r9);
    _mm512_store_epi32(&matT[10*16], ra);
    _mm512_store_epi32(&matT[11*16], rb);
    _mm512_store_epi32(&matT[12*16], rc);
    _mm512_store_epi32(&matT[13*16], rd);
    _mm512_store_epi32(&matT[14*16], re);
    _mm512_store_epi32(&matT[15*16], rf);   
    }
}

void gather(int *mat, int *matT, int rep) {
    int i,j;
    int index[16] __attribute__((aligned(64)));

    __m512i vindex;

    for(i=0; i<16; i++) index[i] = 16*i;
    for(i=0; i<256; i++) mat[i] = i;
    vindex = _mm512_load_epi32(index);

    for(i=0; i<rep; i++) {
        _mm512_store_epi32(&matT[ 0*16], _mm512_i32gather_epi32(vindex, &mat[ 0], 4));
        _mm512_store_epi32(&matT[ 1*16], _mm512_i32gather_epi32(vindex, &mat[ 1], 4));
        _mm512_store_epi32(&matT[ 2*16], _mm512_i32gather_epi32(vindex, &mat[ 2], 4));
        _mm512_store_epi32(&matT[ 3*16], _mm512_i32gather_epi32(vindex, &mat[ 3], 4));
        _mm512_store_epi32(&matT[ 4*16], _mm512_i32gather_epi32(vindex, &mat[ 4], 4));
        _mm512_store_epi32(&matT[ 5*16], _mm512_i32gather_epi32(vindex, &mat[ 5], 4));
        _mm512_store_epi32(&matT[ 6*16], _mm512_i32gather_epi32(vindex, &mat[ 6], 4));
        _mm512_store_epi32(&matT[ 7*16], _mm512_i32gather_epi32(vindex, &mat[ 7], 4));
        _mm512_store_epi32(&matT[ 8*16], _mm512_i32gather_epi32(vindex, &mat[ 8], 4));
        _mm512_store_epi32(&matT[ 9*16], _mm512_i32gather_epi32(vindex, &mat[ 9], 4));
        _mm512_store_epi32(&matT[10*16], _mm512_i32gather_epi32(vindex, &mat[10], 4));
        _mm512_store_epi32(&matT[11*16], _mm512_i32gather_epi32(vindex, &mat[11], 4));
        _mm512_store_epi32(&matT[12*16], _mm512_i32gather_epi32(vindex, &mat[12], 4));
        _mm512_store_epi32(&matT[13*16], _mm512_i32gather_epi32(vindex, &mat[13], 4));
        _mm512_store_epi32(&matT[14*16], _mm512_i32gather_epi32(vindex, &mat[14], 4));
        _mm512_store_epi32(&matT[15*16], _mm512_i32gather_epi32(vindex, &mat[15], 4));
    }
}

int verify(int *mat) {
    int i,j;
    int error = 0;
    for(i=0; i<16; i++) {
      for(j=0; j<16; j++) {
        if(mat[j*16+i] != i*16+j) error++;
      }
    }
    return error;
}

void print_mat(int *mat) {
    int i,j;
    for(i=0; i<16; i++) {
      for(j=0; j<16; j++) printf("%2X ", mat[i*16+j]);
      puts("");
    }
    puts("");
}

int main(void) {
    int i,j, rep;
    int mat[256] __attribute__((aligned(64)));
    int matT[256] __attribute__((aligned(64)));
    double dtime;

    rep = 10000000;
    for(i=0; i<256; i++) mat[i] = i;
    print_mat(mat);

    gather(mat, matT,1);
    for(i=0; i<256; i++) mat[i] = i;
    dtime = -omp_get_wtime();
    gather(mat, matT, rep);
    dtime += omp_get_wtime();
    printf("errors %d\n", verify(matT));
    printf("dtime %f\n", dtime);
    print_mat(matT);

    tran(mat,matT,1);
    dtime = -omp_get_wtime();
    tran(mat, matT, rep);
    dtime += omp_get_wtime();
    printf("errors %d\n", verify(matT));
    printf("dtime %f\n", dtime);
    print_mat(matT);
}

Функция gather заняла 1,13 с, а функция tran 0,8 с.


Согласно руководству Agner Fog по микроархитектуре, инструкции перемешивания и перестановки имеют низкую производительность с KNL. Инструкции по перемешиванию и распаковке, использованные в моем исходном ответе https://stackoverflow.com/a/29587984/2542702, имеют обратный пропускная способность 2. Мне удалось значительно улучшить производительность, используя вместо этого vpermq, который имеет обратную пропускную способность 1. Кроме того, я улучшил первую 1/4 транспонирования, используя vinserti64x4 (см. tran_new2 ниже). Вот таблица времен. Функция tran занимает 0,8 секунды, а функция tran_new2 - 0,46 с.

void tran_new2(int* mat, int* matT, int rep) {
  __m512i t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, ta, tb, tc, td, te, tf;
  __m512i r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, ra, rb, rc, rd, re, rf;

  int mask;
  int64_t idx1[8] __attribute__((aligned(64))) = {2, 3, 0, 1, 6, 7, 4, 5}; 
  int64_t idx2[8] __attribute__((aligned(64))) = {1, 0, 3, 2, 5, 4, 7, 6}; 
  int32_t idx3[16] __attribute__((aligned(64))) = {1, 0, 3, 2, 5 ,4 ,7 ,6 ,9 ,8 , 11, 10, 13, 12 ,15, 14};
  __m512i vidx1 = _mm512_load_epi64(idx1);
  __m512i vidx2 = _mm512_load_epi64(idx2);
  __m512i vidx3 = _mm512_load_epi32(idx3);

  int i;

  for(i=0; i<rep; i++) {

  t0 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 0*16+0])), _mm256_load_si256((__m256i*)&mat[ 8*16+0]), 1);
  t1 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 1*16+0])), _mm256_load_si256((__m256i*)&mat[ 9*16+0]), 1);
  t2 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 2*16+0])), _mm256_load_si256((__m256i*)&mat[10*16+0]), 1);
  t3 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 3*16+0])), _mm256_load_si256((__m256i*)&mat[11*16+0]), 1);
  t4 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 4*16+0])), _mm256_load_si256((__m256i*)&mat[12*16+0]), 1);
  t5 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 5*16+0])), _mm256_load_si256((__m256i*)&mat[13*16+0]), 1);
  t6 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 6*16+0])), _mm256_load_si256((__m256i*)&mat[14*16+0]), 1);
  t7 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 7*16+0])), _mm256_load_si256((__m256i*)&mat[15*16+0]), 1);

  t8 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 0*16+8])), _mm256_load_si256((__m256i*)&mat[ 8*16+8]), 1);
  t9 = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 1*16+8])), _mm256_load_si256((__m256i*)&mat[ 9*16+8]), 1);
  ta = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 2*16+8])), _mm256_load_si256((__m256i*)&mat[10*16+8]), 1);
  tb = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 3*16+8])), _mm256_load_si256((__m256i*)&mat[11*16+8]), 1);
  tc = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 4*16+8])), _mm256_load_si256((__m256i*)&mat[12*16+8]), 1);
  td = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 5*16+8])), _mm256_load_si256((__m256i*)&mat[13*16+8]), 1);
  te = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 6*16+8])), _mm256_load_si256((__m256i*)&mat[14*16+8]), 1);
  tf = _mm512_inserti64x4(_mm512_castsi256_si512(_mm256_load_si256((__m256i*)&mat[ 7*16+8])), _mm256_load_si256((__m256i*)&mat[15*16+8]), 1);

  mask= 0xcc;
  r0 = _mm512_mask_permutexvar_epi64(t0, (__mmask8)mask, vidx1, t4);
  r1 = _mm512_mask_permutexvar_epi64(t1, (__mmask8)mask, vidx1, t5);
  r2 = _mm512_mask_permutexvar_epi64(t2, (__mmask8)mask, vidx1, t6);
  r3 = _mm512_mask_permutexvar_epi64(t3, (__mmask8)mask, vidx1, t7);
  r8 = _mm512_mask_permutexvar_epi64(t8, (__mmask8)mask, vidx1, tc);
  r9 = _mm512_mask_permutexvar_epi64(t9, (__mmask8)mask, vidx1, td);
  ra = _mm512_mask_permutexvar_epi64(ta, (__mmask8)mask, vidx1, te);
  rb = _mm512_mask_permutexvar_epi64(tb, (__mmask8)mask, vidx1, tf);

  mask= 0x33;
  r4 = _mm512_mask_permutexvar_epi64(t4, (__mmask8)mask, vidx1, t0);
  r5 = _mm512_mask_permutexvar_epi64(t5, (__mmask8)mask, vidx1, t1);
  r6 = _mm512_mask_permutexvar_epi64(t6, (__mmask8)mask, vidx1, t2);
  r7 = _mm512_mask_permutexvar_epi64(t7, (__mmask8)mask, vidx1, t3);
  rc = _mm512_mask_permutexvar_epi64(tc, (__mmask8)mask, vidx1, t8);
  rd = _mm512_mask_permutexvar_epi64(td, (__mmask8)mask, vidx1, t9);
  re = _mm512_mask_permutexvar_epi64(te, (__mmask8)mask, vidx1, ta);
  rf = _mm512_mask_permutexvar_epi64(tf, (__mmask8)mask, vidx1, tb);

  mask = 0xaa;
  t0 = _mm512_mask_permutexvar_epi64(r0, (__mmask8)mask, vidx2, r2);
  t1 = _mm512_mask_permutexvar_epi64(r1, (__mmask8)mask, vidx2, r3);
  t4 = _mm512_mask_permutexvar_epi64(r4, (__mmask8)mask, vidx2, r6);
  t5 = _mm512_mask_permutexvar_epi64(r5, (__mmask8)mask, vidx2, r7);
  t8 = _mm512_mask_permutexvar_epi64(r8, (__mmask8)mask, vidx2, ra);
  t9 = _mm512_mask_permutexvar_epi64(r9, (__mmask8)mask, vidx2, rb);
  tc = _mm512_mask_permutexvar_epi64(rc, (__mmask8)mask, vidx2, re);
  td = _mm512_mask_permutexvar_epi64(rd, (__mmask8)mask, vidx2, rf);

  mask = 0x55;
  t2 = _mm512_mask_permutexvar_epi64(r2, (__mmask8)mask, vidx2, r0);
  t3 = _mm512_mask_permutexvar_epi64(r3, (__mmask8)mask, vidx2, r1);
  t6 = _mm512_mask_permutexvar_epi64(r6, (__mmask8)mask, vidx2, r4);
  t7 = _mm512_mask_permutexvar_epi64(r7, (__mmask8)mask, vidx2, r5);
  ta = _mm512_mask_permutexvar_epi64(ra, (__mmask8)mask, vidx2, r8);
  tb = _mm512_mask_permutexvar_epi64(rb, (__mmask8)mask, vidx2, r9);
  te = _mm512_mask_permutexvar_epi64(re, (__mmask8)mask, vidx2, rc);
  tf = _mm512_mask_permutexvar_epi64(rf, (__mmask8)mask, vidx2, rd);

  mask = 0xaaaa;
  r0 = _mm512_mask_permutexvar_epi32(t0, (__mmask16)mask, vidx3, t1);
  r2 = _mm512_mask_permutexvar_epi32(t2, (__mmask16)mask, vidx3, t3);
  r4 = _mm512_mask_permutexvar_epi32(t4, (__mmask16)mask, vidx3, t5);
  r6 = _mm512_mask_permutexvar_epi32(t6, (__mmask16)mask, vidx3, t7);
  r8 = _mm512_mask_permutexvar_epi32(t8, (__mmask16)mask, vidx3, t9);
  ra = _mm512_mask_permutexvar_epi32(ta, (__mmask16)mask, vidx3, tb);
  rc = _mm512_mask_permutexvar_epi32(tc, (__mmask16)mask, vidx3, td);
  re = _mm512_mask_permutexvar_epi32(te, (__mmask16)mask, vidx3, tf);    

  mask = 0x5555;
  r1 = _mm512_mask_permutexvar_epi32(t1, (__mmask16)mask, vidx3, t0);
  r3 = _mm512_mask_permutexvar_epi32(t3, (__mmask16)mask, vidx3, t2);
  r5 = _mm512_mask_permutexvar_epi32(t5, (__mmask16)mask, vidx3, t4);
  r7 = _mm512_mask_permutexvar_epi32(t7, (__mmask16)mask, vidx3, t6);
  r9 = _mm512_mask_permutexvar_epi32(t9, (__mmask16)mask, vidx3, t8);  
  rb = _mm512_mask_permutexvar_epi32(tb, (__mmask16)mask, vidx3, ta);  
  rd = _mm512_mask_permutexvar_epi32(td, (__mmask16)mask, vidx3, tc);
  rf = _mm512_mask_permutexvar_epi32(tf, (__mmask16)mask, vidx3, te);

  _mm512_store_epi32(&matT[ 0*16], r0);
  _mm512_store_epi32(&matT[ 1*16], r1);
  _mm512_store_epi32(&matT[ 2*16], r2);
  _mm512_store_epi32(&matT[ 3*16], r3);
  _mm512_store_epi32(&matT[ 4*16], r4);
  _mm512_store_epi32(&matT[ 5*16], r5);
  _mm512_store_epi32(&matT[ 6*16], r6);
  _mm512_store_epi32(&matT[ 7*16], r7);
  _mm512_store_epi32(&matT[ 8*16], r8);
  _mm512_store_epi32(&matT[ 9*16], r9);
  _mm512_store_epi32(&matT[10*16], ra);
  _mm512_store_epi32(&matT[11*16], rb);
  _mm512_store_epi32(&matT[12*16], rc);
  _mm512_store_epi32(&matT[13*16], rd);
  _mm512_store_epi32(&matT[14*16], re);
  _mm512_store_epi32(&matT[15*16], rf);
  int* tmp = mat;
  mat = matT;
  matT = tmp;
  }
}
person Z boson    schedule 21.12.2016
comment
Отлично! В предыдущем ответе вы написали, что транспонирование 8x8 + r / w использует 40 инструкций. То есть: 8 загрузок, 24 перемешивания в порту выполнения 5 и 8 хранилищах. В параграфе 11.11.2 документа Intel 64-ia-32-architecture-optimisation-manual они заменяют 8 из этих перетасовок на 8 vinsertf128 инструкций с оперантом памяти. Это приводит к меньшему давлению на порт 5: 16 инструкций на порт 5. Фактически огромная полоса пропускания L1 используется для уменьшения узкого места на порту 5. Результат - более быстрый алгоритм. Как вы думаете, можно ли использовать здесь аналогичную идею для ускорения транспонирования 16x16? - person wim; 21.12.2016
comment
Недостаточно места в моем предыдущем комментарии: ссылка на документ: здесь - person wim; 21.12.2016
comment
@wim огромное спасибо за ссылку! Я быстро посмотрел на него. Когда я создавал ответ 8x8, я думал не о давлении порта, а только о количестве инструкций. Мне придется разобраться в этом и вернуться к вам. - person Z boson; 21.12.2016
comment
@wim: Хорошая идея. Но, судя по таблицам Агнера Фога, я думаю, что vinsert KNL с источником памяти все еще нуждается в блоке тасования. Он основан на Silvermont, сильно отличается от Haswell. В таблицах Agner Fog не указан порт для vinsertf128 или его вариантов AVX512, но, как и в Haswell, похоже, что есть только один блок перемешивания. Это на FP0. vinsertf32x4 z,z,m128/m256 имеют пропускную способность один на такт, а не один на 0,5 с, как нагрузки, поэтому они могут по-прежнему использовать устройство перемешивания. Широковещательные сообщения полностью обрабатываются портом загрузки, поэтому vbroadcastf64x4 z,m256 имеет один порт на 0,5 с. - person Peter Cordes; 21.12.2016
comment
Для сравнения, SKL vinsertf128 y,y,m128,i составляет 2 мопа: нагрузка и p015. Он рассчитан на пропускную способность 0,5 с. Версия с регистровым источником имеет размер 1 мупа, но работает только на порте 5. Возможно, версия для памяти использует моп широковещательной загрузки, поэтому мап ALU может быть смешанным, а не перемешиваемым. Смеси с большой степенью детализации могут работать на каждом порте ALU. В большинстве случаев это делает его лучше, чем однократное перемешивание с микроплавкой. - person Peter Cordes; 21.12.2016
comment
@PeterCordes Действительно, на KNL нет порта 5. Перемешивание переходит к блоку FP0. Из руководства Агнера неясно, какие ресурсы vinsertf64x4 использует. Но, по крайней мере, мы можем имитировать vinsertf64x4 KNL путем vbroadcastf6x4 загрузки из памяти плюс vblendmpd, которые, по словам Агнера Фога, имеют пропускную способность один на 0,5 с. vblendmpd работает на FP0 или FP1. Итак, насколько я могу судить (я совершенно не знаком с KNL, я только начал читать таблицы инструкций Агнера на KNL), за два цикла мы можем сделать 2 перемешивания на FP0 и эмулированный vinsertf64x4 на порту памяти и на FP1. - person wim; 22.12.2016
comment
@wim: да, хорошо, что мы можем транслировать и смешивать отдельно. По-видимому, они не реализовали vinsert таким образом внутри, потому что даже двухэлементные инструкции микрокодированы на KNL и имеют ужасную пропускную способность внешнего интерфейса. (один на 7 часов). Все, что я знаю о микрооптимизации KNL, исходит из руководства + таблиц Agner Fog по микроархитектуре, так что мы находимся в одной лодке с этим. - person Peter Cordes; 22.12.2016
comment
@PeterCordes Верно! Кстати, я только что заметил, что пропускная способность vshufi32x4 составляет всего одну инструкцию за два цикла вместо моего предыдущего предположения об одной инструкции за цикл. - person wim; 22.12.2016
comment
@wim: да, похоже, что некоторые перемешивания не полностью конвейеризированы. Согласно документам Агнера, интерфейс, вероятно, по-прежнему будет основным узким местом во многих случаях (до двух однократных инструкций за такт из 16 байтов кода. 7c для декодирования многопозиционных инструкций и некоторые другие возможные задержки ). Это делает разумным не сходить с ума с бюджетом транзисторов при построении сверхбыстрых устройств перетасовки в KNL. Декодирование вряд ли сможет обеспечить питание ООО ядра, если в него примешаны какие-либо не-векторные инструкции ALU (чистая загрузка, чистая память или целочисленный материал). - person Peter Cordes; 22.12.2016
comment
Так что IDK, если имитировать 1-uop vinsert с последовательностью из 2 инструкций, того стоит. Если этот код является узким местом во внешнем интерфейсе, то, вероятно, нет. Но если это не так, то он может выиграть. (И это может не быть узким местом во внешнем интерфейсе, если это в основном перетасовка на 1 мупа, а некоторые из них не полностью конвейеризованы). Окно не по порядку намного меньше, чем в SKL, поэтому, надеюсь, компилятор хорошо справляется с планированием этих загрузок, например перемешивание в некотором перемешивании перед отправкой всех 16 загрузок. - person Peter Cordes; 22.12.2016
comment
@PeterCordes, я думаю, было бы правильнее сказать, что KNL основан на Airmont. Я не уверен на 100%, но думаю, что Airmont - это больше, чем просто усадка Silvermont (с 22 до 14 нм). Агнер говорит, что Knights Landing имеет лучший конвейер вне очереди, чем ... Silvermont, который не может выполнять команды с плавающей запятой и векторные инструкции не по порядку. Может быть, у Airmont тоже есть такая функция? - person Z boson; 22.12.2016
comment
@PeterCordes, согласно этой ссылке, это всего лишь уменьшение размера кристалла и обновление графического процессора. . - person Z boson; 22.12.2016
comment
@PeterCordes, хорошо, чем больше я читал о KNL, я понял, что это собственная арка. Например. у него нет буфера цикла, в отличие от Silvermont. Таким образом, он больше всего похож на Silvermont, но, тем не менее, имеет важные отличия. - person Z boson; 22.12.2016
comment
Интересно, что они сбросили буфер петли, я не заметил, что у Сильвермонта он есть. Интересно, было ли это компромиссом с отказом от экономии энергии ради экономии площади кристалла для большего количества ядер на чип? Может быть, они считают, что большинство векторных петель будут слишком большими? Или, может быть, формат uop пришлось изменить для поддержки AVX512 и сделать поддержку loop buf неудобной? А может быть, с 4-х сторонним SMT было бы отстойно? Похоже, что было бы действительно хорошо иметь буфер цикла для случая без гиперпоточности, поскольку узкие места внешнего интерфейса, вероятно, обычны с накладными расходами цикла + загрузками и сохранениями. - person Peter Cordes; 22.12.2016
comment
@PeterCordes Если мы пренебрегаем загрузкой и сохранением, то есть 64 перемешивания, которые занимают как минимум: (64 перемешивания * 2 цикла / перемешивание на FP0 * 10e6 повторов) / 1,6e9 Гц = 0,8 секунды (частота процессора Turbo Boost = 1,6 ГГц. ). Случайно (?) tran также занимает 0,8 секунды, что указывает на то, что на FP0 может быть узкое место в случайном порядке. Все остальные инструкции, вероятно, выполняются на других портах выполнения. - person wim; 22.12.2016
comment
Насколько я могу судить, все приведенные здесь инструкции представляют собой одинарные инструкции. Интерфейс может выполнять не более двух инструкций по 1 шагу за такт. В скомпилированном коде (функция tran в компиляторе Godbold, icc 17 -O3) я насчитал 97 инструкций с шагом 1 на цикл (некоторые адреса памяти читаются дважды с помощью vpunpck, а не с помощью каких-либо vmovups). Общее количество циклов на tran составляет 0,8 сек. * 1,6 ГГц = 1,28e9 циклов. За эти циклы выполняется около 97 * 1e7 инструкций. 97e7 / 1,28e9 = 0,76 инструкций за цикл, что намного меньше, чем ограничение внешнего интерфейса, равное 2. - person wim; 22.12.2016
comment
Таким образом, веб-интерфейс вряд ли будет здесь узким местом. Перетасовка на KNL относительно дорога. Я все еще думаю, что можно было бы немного ускорить tran, заменив (например) 16 перетасовок (из 64) на 16 vinsertf64x4 или, если это не сработает, на 16 vbroadcastf64x4 + 16 vblendmpd. - person wim; 22.12.2016
comment
@wim, завтра попробую что-нибудь реализовать с vinsertf64x4 и vbroadcastf64x4 + vblendmpd. - person Z boson; 22.12.2016
comment
Это, конечно, было бы интересно, если бы у вас было время на такое тестирование. По крайней мере, приведенное выше обсуждение заставляет нас копаться в документах KNL, которые кое-что узнают о KNL. :-) - person wim; 22.12.2016
comment
Я нашел значительно более быструю версию, использующую _mm512_mask_permutexvar_epi64 (and epi32) с пропускной способностью 1 (вместо 2 с перетасовкой). Я также использую _mm512_mask_broadcast_i64x4 в первых восьми строках с пропускной способностью 0,5. Я отправлю код позже, после еще нескольких тестов. Я не был уверен, что это будет быстрее, потому что мне нужен векторный индекс (на самом деле три), и я уже использую 32 регистра AVX512, что означает некоторую утечку, но, тем не менее, это намного быстрее. - person Z boson; 26.12.2016
comment
Что ж, у нас в лаборатории есть Xeon Phi 7210. Я только что быстро протестировал ваш первый фрагмент кода, используя те же флаги компилятора, которые вы указали. Результат - 1,78 с против 1,37 с. - person lei_z; 27.12.2016
comment
@ lei.april, это звучит правильно, я думаю. Частота AVX у 7210 составляет 1,1 ГГц, а частота AVX у 7250 - 1,5 ГГц, поэтому 7210 примерно на 36% медленнее, чем 7250 для кода AVX intel.com/content/www/us/en/processors/xeon/ - person Z boson; 27.12.2016
comment
@ lei.april У меня есть более оптимизированная версия транспонирования 16x16 для KNL, которая занимает около 0,5 с по сравнению с 0,8 с с исходным кодом. - person Z boson; 27.12.2016
comment
@wim, я забыл отметить ваше имя в своем комментарии выше, чтобы вы его пропустили. Скоро выложу обновление кода. Смена сюжета. Я рассмотрел пример 11-19 в руководстве Intel, транспонируя матрицу 8x8. Они используют больше инструкций, чем я. Они используют 12 перемешиваний, которые они преобразуют в 4 перемешивания и восемь сочетаний. Использую только 8 перетасовок. Я очень сомневаюсь, что восемь смесей и четыре перетасовки лучше восьми перетасовок. Intel явно не понимала, что это можно сделать за 24 операции. - person Z boson; 27.12.2016
comment
@wim, я обновил свой ответ транспонирования 8x8 на основе примеров в руководстве Intel stackoverflow.com/a/25627536/2542702 - person Z boson; 27.12.2016
comment
@wim, хорошо, я обновил этот ответ, используя vinsertf64x4. - person Z boson; 27.12.2016
comment
@Zboson Это отличное обновление! Меня также удивляет, что _mm512_mask_permutexvar_epi64 намного быстрее, чем перемешивание. От 0,8 до 0,46 с - довольно много. - person wim; 28.12.2016
comment
О, я только что понял, что эта ветка комментариев охватывает почти то же самое, что я только что опубликовал под первым ответом. ржу не могу - person Mysticial; 20.01.2017
comment
Чтобы расширить мои комментарии под другим ответом, вам не нужно переставлять векторы для всего, что остается в пределах 256-битной половины. (спасибо _mm512_mask_shuffle_epi32() и _mm512_mask_permutex_epi64) И поскольку вы выполняете часть шириной 256 с использованием нагрузок, вам вообще не нужны никакие перестановочные векторы. Также можно объединить некоторые константы маски. Например, 0x3333 можно использовать для покрытия как 128-битной, так и 64-битной степени детализации, если вы выполняете 128-битную гранулярность с использованием _mm512_mask_permutex_epi64() и 64-битную с использованием _mm512_mask_shuffle_epi32. - person Mysticial; 20.01.2017