Лучше матрица 8x8 байтов транспонировать с SSE?
Я нашел этот пост, который объясняет, как транспонировать матрицу 8x8 байтов с 24 операциями, и несколько прокручиваний позже есть код, который реализует транспонирование Тем не менее, этот метод не использует тот факт, что мы можем блокировать транспонирование 8x8 на четыре транспонирования 4x4, и каждое из них может быть выполнено только в одной команде тасования ( этот пост является справочным). Итак, я вышел с этим решением:
__m128i transpose4x4mask = _mm_set_epi8(15, 11, 7, 3, 14, 10, 6, 2, 13, 9, 5, 1, 12, 8, 4, 0);
__m128i shuffle8x8Mask = _mm_setr_epi8(0, 1, 2, 3, 8, 9, 10, 11, 4, 5, 6, 7, 12, 13, 14, 15);
void TransposeBlock8x8(uint8_t *src, uint8_t *dst, int srcStride, int dstStride) {
__m128i load0 = _mm_set_epi64x(*(uint64_t*)(src + 1 * srcStride), *(uint64_t*)(src + 0 * srcStride));
__m128i load1 = _mm_set_epi64x(*(uint64_t*)(src + 3 * srcStride), *(uint64_t*)(src + 2 * srcStride));
__m128i load2 = _mm_set_epi64x(*(uint64_t*)(src + 5 * srcStride), *(uint64_t*)(src + 4 * srcStride));
__m128i load3 = _mm_set_epi64x(*(uint64_t*)(src + 7 * srcStride), *(uint64_t*)(src + 6 * srcStride));
__m128i shuffle0 = _mm_shuffle_epi8(load0, shuffle8x8Mask);
__m128i shuffle1 = _mm_shuffle_epi8(load1, shuffle8x8Mask);
__m128i shuffle2 = _mm_shuffle_epi8(load2, shuffle8x8Mask);
__m128i shuffle3 = _mm_shuffle_epi8(load3, shuffle8x8Mask);
__m128i block0 = _mm_unpacklo_epi64(shuffle0, shuffle1);
__m128i block1 = _mm_unpackhi_epi64(shuffle0, shuffle1);
__m128i block2 = _mm_unpacklo_epi64(shuffle2, shuffle3);
__m128i block3 = _mm_unpackhi_epi64(shuffle2, shuffle3);
__m128i transposed0 = _mm_shuffle_epi8(block0, transpose4x4mask);
__m128i transposed1 = _mm_shuffle_epi8(block1, transpose4x4mask);
__m128i transposed2 = _mm_shuffle_epi8(block2, transpose4x4mask);
__m128i transposed3 = _mm_shuffle_epi8(block3, transpose4x4mask);
__m128i store0 = _mm_unpacklo_epi32(transposed0, transposed2);
__m128i store1 = _mm_unpackhi_epi32(transposed0, transposed2);
__m128i store2 = _mm_unpacklo_epi32(transposed1, transposed3);
__m128i store3 = _mm_unpackhi_epi32(transposed1, transposed3);
*((uint64_t*)(dst + 0 * dstStride)) = _mm_extract_epi64(store0, 0);
*((uint64_t*)(dst + 1 * dstStride)) = _mm_extract_epi64(store0, 1);
*((uint64_t*)(dst + 2 * dstStride)) = _mm_extract_epi64(store1, 0);
*((uint64_t*)(dst + 3 * dstStride)) = _mm_extract_epi64(store1, 1);
*((uint64_t*)(dst + 4 * dstStride)) = _mm_extract_epi64(store2, 0);
*((uint64_t*)(dst + 5 * dstStride)) = _mm_extract_epi64(store2, 1);
*((uint64_t*)(dst + 6 * dstStride)) = _mm_extract_epi64(store3, 0);
*((uint64_t*)(dst + 7 * dstStride)) = _mm_extract_epi64(store3, 1);
}
Исключая операции загрузки / сохранения, эта процедура состоит только из 16 инструкций вместо 24.
Что мне не хватает?
7 ответов
Помимо грузов, магазинов и pinsrq-s для чтения и записи в память, возможно, с шагом, не равным 8 байтам, вы можете выполнить транспонирование всего с 12 инструкциями (этот код можно легко использовать в сочетании с тестовым кодом Z-бозона):
void tran8x8b_SSE_v2(char *A, char *B) {
__m128i pshufbcnst = _mm_set_epi8(15,11,7,3, 14,10,6,2, 13,9,5,1, 12,8,4,0);
__m128i B0, B1, B2, B3, T0, T1, T2, T3;
B0 = _mm_loadu_si128((__m128i*)&A[ 0]);
B1 = _mm_loadu_si128((__m128i*)&A[16]);
B2 = _mm_loadu_si128((__m128i*)&A[32]);
B3 = _mm_loadu_si128((__m128i*)&A[48]);
T0 = _mm_castps_si128(_mm_shuffle_ps(_mm_castsi128_ps(B0),_mm_castsi128_ps(B1),0b10001000));
T1 = _mm_castps_si128(_mm_shuffle_ps(_mm_castsi128_ps(B2),_mm_castsi128_ps(B3),0b10001000));
T2 = _mm_castps_si128(_mm_shuffle_ps(_mm_castsi128_ps(B0),_mm_castsi128_ps(B1),0b11011101));
T3 = _mm_castps_si128(_mm_shuffle_ps(_mm_castsi128_ps(B2),_mm_castsi128_ps(B3),0b11011101));
B0 = _mm_shuffle_epi8(T0,pshufbcnst);
B1 = _mm_shuffle_epi8(T1,pshufbcnst);
B2 = _mm_shuffle_epi8(T2,pshufbcnst);
B3 = _mm_shuffle_epi8(T3,pshufbcnst);
T0 = _mm_unpacklo_epi32(B0,B1);
T1 = _mm_unpackhi_epi32(B0,B1);
T2 = _mm_unpacklo_epi32(B2,B3);
T3 = _mm_unpackhi_epi32(B2,B3);
_mm_storeu_si128((__m128i*)&B[ 0], T0);
_mm_storeu_si128((__m128i*)&B[16], T1);
_mm_storeu_si128((__m128i*)&B[32], T2);
_mm_storeu_si128((__m128i*)&B[48], T3);
}
Здесь мы используем 32-битное перемешивание с плавающей точкой, которое является более гибким, чем epi32 перетасовать. Приведения не генерируют дополнительных инструкций (код, созданный с помощью gcc 5.4):
tran8x8b_SSE_v2:
.LFB4885:
.cfi_startproc
vmovdqu 48(%rdi), %xmm5
vmovdqu 32(%rdi), %xmm2
vmovdqu 16(%rdi), %xmm0
vmovdqu (%rdi), %xmm1
vshufps $136, %xmm5, %xmm2, %xmm4
vshufps $221, %xmm5, %xmm2, %xmm2
vmovdqa .LC6(%rip), %xmm5
vshufps $136, %xmm0, %xmm1, %xmm3
vshufps $221, %xmm0, %xmm1, %xmm1
vpshufb %xmm5, %xmm3, %xmm3
vpshufb %xmm5, %xmm1, %xmm0
vpshufb %xmm5, %xmm4, %xmm4
vpshufb %xmm5, %xmm2, %xmm1
vpunpckldq %xmm4, %xmm3, %xmm5
vpunpckldq %xmm1, %xmm0, %xmm2
vpunpckhdq %xmm4, %xmm3, %xmm3
vpunpckhdq %xmm1, %xmm0, %xmm0
vmovups %xmm5, (%rsi)
vmovups %xmm3, 16(%rsi)
vmovups %xmm2, 32(%rsi)
vmovups %xmm0, 48(%rsi)
ret
.cfi_endproc
В некоторых, но не во всех старых процессорах может быть небольшая задержка обхода (между 0 и 2 циклами) для перемещения данных между целым числом и модулями с плавающей запятой. Это увеличивает задержку функции, но не обязательно влияет на пропускную способность кода.
Простой тест на задержку с перестановками 1e9:
for (int i=0;i<500000000;i++){
tran8x8b_SSE(A,C);
tran8x8b_SSE(C,A);
}
print8x8b(A);
Это займет около 5,5 секунд (19,7e9 циклов) с tran8x8b_SSE и 4,5 секунды (16,0e9 циклов) с tran8x8b_SSE_v2 (Intel Core i5-6500). Обратите внимание, что загрузка и хранилища не были устранены компилятором, хотя функции были встроены в цикл for.
Обновление: решение AVX2-128 / SSE 4.1 со смесями.
"Перемешивание" (распаковка, перемешивание) обрабатывается портом 5 с 1 инструкцией на цикл процессора на современном процессоре. Иногда окупается замена одного "шаффла" на две смеси. На Skylake 32-битные инструкции смешивания могут выполняться на портах 0, 1 или 5.
К несчастью, _mm_blend_epi32 только AVX2-128. Эффективной альтернативой SSE 4.1 является _mm_blend_ps в сочетании с несколькими кастами (которые обычно бесплатны). 12 "шаффлов" заменяются 8 шаффлами в сочетании с 8 блендами.
Простой тест на задержку теперь выполняется за 3,6 секунды (циклы 13e9 ЦП), что на 18 % быстрее, чем результаты с tran8x8b_SSE_v2,
Код:
/* AVX2-128 version, sse 4.1 version see ----------------> SSE 4.1 version of tran8x8b_AVX2_128() */
void tran8x8b_AVX2_128(char *A, char *B) { /* void tran8x8b_SSE4_1(char *A, char *B) { */
__m128i pshufbcnst_0 = _mm_set_epi8(15, 7,11, 3,
13, 5, 9, 1, 14, 6,10, 2, 12, 4, 8, 0); /* __m128i pshufbcnst_0 = _mm_set_epi8(15, 7,11, 3, 13, 5, 9, 1, 14, 6,10, 2, 12, 4, 8, 0); */
__m128i pshufbcnst_1 = _mm_set_epi8(13, 5, 9, 1,
15, 7,11, 3, 12, 4, 8, 0, 14, 6,10, 2); /* __m128i pshufbcnst_1 = _mm_set_epi8(13, 5, 9, 1, 15, 7,11, 3, 12, 4, 8, 0, 14, 6,10, 2); */
__m128i pshufbcnst_2 = _mm_set_epi8(11, 3,15, 7,
9, 1,13, 5, 10, 2,14, 6, 8, 0,12, 4); /* __m128i pshufbcnst_2 = _mm_set_epi8(11, 3,15, 7, 9, 1,13, 5, 10, 2,14, 6, 8, 0,12, 4); */
__m128i pshufbcnst_3 = _mm_set_epi8( 9, 1,13, 5,
11, 3,15, 7, 8, 0,12, 4, 10, 2,14, 6); /* __m128i pshufbcnst_3 = _mm_set_epi8( 9, 1,13, 5, 11, 3,15, 7, 8, 0,12, 4, 10, 2,14, 6); */
__m128i B0, B1, B2, B3, T0, T1, T2, T3; /* __m128 B0, B1, B2, B3, T0, T1, T2, T3; */
/* */
B0 = _mm_loadu_si128((__m128i*)&A[ 0]); /* B0 = _mm_loadu_ps((float*)&A[ 0]); */
B1 = _mm_loadu_si128((__m128i*)&A[16]); /* B1 = _mm_loadu_ps((float*)&A[16]); */
B2 = _mm_loadu_si128((__m128i*)&A[32]); /* B2 = _mm_loadu_ps((float*)&A[32]); */
B3 = _mm_loadu_si128((__m128i*)&A[48]); /* B3 = _mm_loadu_ps((float*)&A[48]); */
/* */
B1 = _mm_shuffle_epi32(B1,0b10110001); /* B1 = _mm_shuffle_ps(B1,B1,0b10110001); */
B3 = _mm_shuffle_epi32(B3,0b10110001); /* B3 = _mm_shuffle_ps(B3,B3,0b10110001); */
T0 = _mm_blend_epi32(B0,B1,0b1010); /* T0 = _mm_blend_ps(B0,B1,0b1010); */
T1 = _mm_blend_epi32(B2,B3,0b1010); /* T1 = _mm_blend_ps(B2,B3,0b1010); */
T2 = _mm_blend_epi32(B0,B1,0b0101); /* T2 = _mm_blend_ps(B0,B1,0b0101); */
T3 = _mm_blend_epi32(B2,B3,0b0101); /* T3 = _mm_blend_ps(B2,B3,0b0101); */
/* */
B0 = _mm_shuffle_epi8(T0,pshufbcnst_0); /* B0 = _mm_castsi128_ps(_mm_shuffle_epi8(_mm_castps_si128(T0),pshufbcnst_0)); */
B1 = _mm_shuffle_epi8(T1,pshufbcnst_1); /* B1 = _mm_castsi128_ps(_mm_shuffle_epi8(_mm_castps_si128(T1),pshufbcnst_1)); */
B2 = _mm_shuffle_epi8(T2,pshufbcnst_2); /* B2 = _mm_castsi128_ps(_mm_shuffle_epi8(_mm_castps_si128(T2),pshufbcnst_2)); */
B3 = _mm_shuffle_epi8(T3,pshufbcnst_3); /* B3 = _mm_castsi128_ps(_mm_shuffle_epi8(_mm_castps_si128(T3),pshufbcnst_3)); */
/* */
T0 = _mm_blend_epi32(B0,B1,0b1010); /* T0 = _mm_blend_ps(B0,B1,0b1010); */
T1 = _mm_blend_epi32(B0,B1,0b0101); /* T1 = _mm_blend_ps(B0,B1,0b0101); */
T2 = _mm_blend_epi32(B2,B3,0b1010); /* T2 = _mm_blend_ps(B2,B3,0b1010); */
T3 = _mm_blend_epi32(B2,B3,0b0101); /* T3 = _mm_blend_ps(B2,B3,0b0101); */
T1 = _mm_shuffle_epi32(T1,0b10110001); /* T1 = _mm_shuffle_ps(T1,T1,0b10110001); */
T3 = _mm_shuffle_epi32(T3,0b10110001); /* T3 = _mm_shuffle_ps(T3,T3,0b10110001); */
/* */
_mm_storeu_si128((__m128i*)&B[ 0], T0); /* _mm_storeu_ps((float*)&B[ 0], T0); */
_mm_storeu_si128((__m128i*)&B[16], T1); /* _mm_storeu_ps((float*)&B[16], T1); */
_mm_storeu_si128((__m128i*)&B[32], T2); /* _mm_storeu_ps((float*)&B[32], T2); */
_mm_storeu_si128((__m128i*)&B[48], T3); /* _mm_storeu_ps((float*)&B[48], T3); */
} /* } */
Выкладываю это как ответ. Я также собираюсь изменить название вопроса с "... с SSE" на "... с SIMD" из-за некоторых ответов и комментариев, полученных до сих пор.
Мне удалось транспонировать матрицу с AVX2 только в 8 инструкциях, 10 из которых включают загрузку / сохранение (исключая загрузки масок). РЕДАКТИРОВАТЬ: я нашел более короткую версию. Увидеть ниже. Это тот случай, когда все матрицы находятся в памяти, поэтому можно использовать прямую загрузку / сохранение.
Вот код C:
void tran8x8b_AVX2(char *src, char *dst) {
__m256i perm = _mm256_set_epi8(
0, 0, 0, 7,
0, 0, 0, 5,
0, 0, 0, 3,
0, 0, 0, 1,
0, 0, 0, 6,
0, 0, 0, 4,
0, 0, 0, 2,
0, 0, 0, 0
);
__m256i tm = _mm256_set_epi8(
15, 11, 7, 3,
14, 10, 6, 2,
13, 9, 5, 1,
12, 8, 4, 0,
15, 11, 7, 3,
14, 10, 6, 2,
13, 9, 5, 1,
12, 8, 4, 0
);
__m256i load0 = _mm256_loadu_si256((__m256i*)&src[ 0]);
__m256i load1 = _mm256_loadu_si256((__m256i*)&src[32]);
__m256i perm0 = _mm256_permutevar8x32_epi32(load0, perm);
__m256i perm1 = _mm256_permutevar8x32_epi32(load1, perm);
__m256i transpose0 = _mm256_shuffle_epi8(perm0, tm);
__m256i transpose1 = _mm256_shuffle_epi8(perm1, tm);
__m256i unpack0 = _mm256_unpacklo_epi32(transpose0, transpose1);
__m256i unpack1 = _mm256_unpackhi_epi32(transpose0, transpose1);
perm0 = _mm256_castps_si256(_mm256_permute2f128_ps(_mm256_castsi256_ps(unpack0), _mm256_castsi256_ps(unpack1), 32));
perm1 = _mm256_castps_si256(_mm256_permute2f128_ps(_mm256_castsi256_ps(unpack0), _mm256_castsi256_ps(unpack1), 49));
_mm256_storeu_si256((__m256i*)&dst[ 0], perm0);
_mm256_storeu_si256((__m256i*)&dst[32], perm1);
}
GCC был достаточно умен, чтобы выполнить перестановку во время загрузки AVX, сохранив две инструкции. Вот вывод компилятора:
tran8x8b_AVX2(char*, char*):
vmovdqa ymm1, YMMWORD PTR .LC0[rip]
vmovdqa ymm2, YMMWORD PTR .LC1[rip]
vpermd ymm0, ymm1, YMMWORD PTR [rdi]
vpermd ymm1, ymm1, YMMWORD PTR [rdi+32]
vpshufb ymm0, ymm0, ymm2
vpshufb ymm1, ymm1, ymm2
vpunpckldq ymm2, ymm0, ymm1
vpunpckhdq ymm0, ymm0, ymm1
vinsertf128 ymm1, ymm2, xmm0, 1
vperm2f128 ymm0, ymm2, ymm0, 49
vmovdqu YMMWORD PTR [rsi], ymm1
vmovdqu YMMWORD PTR [rsi+32], ymm0
vzeroupper
ret
Он испустил vzerupper инструкция с -O3, но снижение до -O1 удаляет это.
В случае моей первоначальной проблемы (большая матрица, и я увеличиваю ее до размера 8x8), обработка шагов уничтожает вывод довольно плохо:
void tran8x8b_AVX2(char *src, char *dst, int srcStride, int dstStride) {
__m256i load0 = _mm256_set_epi64x(*(uint64_t*)(src + 3 * srcStride), *(uint64_t*)(src + 2 * srcStride), *(uint64_t*)(src + 1 * srcStride), *(uint64_t*)(src + 0 * srcStride));
__m256i load1 = _mm256_set_epi64x(*(uint64_t*)(src + 7 * srcStride), *(uint64_t*)(src + 6 * srcStride), *(uint64_t*)(src + 5 * srcStride), *(uint64_t*)(src + 4 * srcStride));
// ... the same as before, however we can skip the final permutations because we need to handle the destination stride...
*((uint64_t*)(dst + 0 * dstStride)) = _mm256_extract_epi64(unpack0, 0);
*((uint64_t*)(dst + 1 * dstStride)) = _mm256_extract_epi64(unpack0, 1);
*((uint64_t*)(dst + 2 * dstStride)) = _mm256_extract_epi64(unpack1, 0);
*((uint64_t*)(dst + 3 * dstStride)) = _mm256_extract_epi64(unpack1, 1);
*((uint64_t*)(dst + 4 * dstStride)) = _mm256_extract_epi64(unpack0, 2);
*((uint64_t*)(dst + 5 * dstStride)) = _mm256_extract_epi64(unpack0, 3);
*((uint64_t*)(dst + 6 * dstStride)) = _mm256_extract_epi64(unpack1, 2);
*((uint64_t*)(dst + 7 * dstStride)) = _mm256_extract_epi64(unpack1, 3);
}
Вот вывод компилятора:
tran8x8b_AVX2(char*, char*, int, int):
movsx rdx, edx
vmovq xmm5, QWORD PTR [rdi]
lea r9, [rdi+rdx]
vmovdqa ymm3, YMMWORD PTR .LC0[rip]
movsx rcx, ecx
lea r11, [r9+rdx]
vpinsrq xmm0, xmm5, QWORD PTR [r9], 1
lea r10, [r11+rdx]
vmovq xmm4, QWORD PTR [r11]
vpinsrq xmm1, xmm4, QWORD PTR [r10], 1
lea r8, [r10+rdx]
lea rax, [r8+rdx]
vmovq xmm7, QWORD PTR [r8]
vmovq xmm6, QWORD PTR [rax+rdx]
vpinsrq xmm2, xmm7, QWORD PTR [rax], 1
vinserti128 ymm1, ymm0, xmm1, 0x1
vpinsrq xmm0, xmm6, QWORD PTR [rax+rdx*2], 1
lea rax, [rsi+rcx]
vpermd ymm1, ymm3, ymm1
vinserti128 ymm0, ymm2, xmm0, 0x1
vmovdqa ymm2, YMMWORD PTR .LC1[rip]
vpshufb ymm1, ymm1, ymm2
vpermd ymm0, ymm3, ymm0
vpshufb ymm0, ymm0, ymm2
vpunpckldq ymm2, ymm1, ymm0
vpunpckhdq ymm0, ymm1, ymm0
vmovdqa xmm1, xmm2
vmovq QWORD PTR [rsi], xmm1
vpextrq QWORD PTR [rax], xmm1, 1
vmovdqa xmm1, xmm0
add rax, rcx
vextracti128 xmm0, ymm0, 0x1
vmovq QWORD PTR [rax], xmm1
add rax, rcx
vpextrq QWORD PTR [rax], xmm1, 1
add rax, rcx
vextracti128 xmm1, ymm2, 0x1
vmovq QWORD PTR [rax], xmm1
add rax, rcx
vpextrq QWORD PTR [rax], xmm1, 1
vmovq QWORD PTR [rax+rcx], xmm0
vpextrq QWORD PTR [rax+rcx*2], xmm0, 1
vzeroupper
ret
Тем не менее, это не имеет большого значения, если сравнивать с выводом моего исходного кода.
РЕДАКТИРОВАТЬ: я нашел более короткую версию. Всего 4 инструкции, 8 считая оба груза / хранилища. Это возможно, потому что я читаю матрицу по-другому, скрывая некоторые "перемешивания" в инструкции "собирать" во время загрузки. Также обратите внимание, что окончательная перестановка необходима для выполнения сохранения, потому что AVX2 не имеет инструкции "разброса". Наличие инструкции разброса сводит все к двум инструкциям. Кроме того, обратите внимание, что я могу без проблем справиться с src, изменив содержимое vindex вектор.
К сожалению, этот AVX_v2 кажется медленнее, чем предыдущий. Вот код:
void tran8x8b_AVX2_v2(char *src1, char *dst1) {
__m256i tm = _mm256_set_epi8(
15, 11, 7, 3,
14, 10, 6, 2,
13, 9, 5, 1,
12, 8, 4, 0,
15, 11, 7, 3,
14, 10, 6, 2,
13, 9, 5, 1,
12, 8, 4, 0
);
__m256i vindex = _mm256_setr_epi32(0, 8, 16, 24, 32, 40, 48, 56);
__m256i perm = _mm256_setr_epi32(0, 4, 1, 5, 2, 6, 3, 7);
__m256i load0 = _mm256_i32gather_epi32((int*)src1, vindex, 1);
__m256i load1 = _mm256_i32gather_epi32((int*)(src1 + 4), vindex, 1);
__m256i transpose0 = _mm256_shuffle_epi8(load0, tm);
__m256i transpose1 = _mm256_shuffle_epi8(load1, tm);
__m256i final0 = _mm256_permutevar8x32_epi32(transpose0, perm);
__m256i final1 = _mm256_permutevar8x32_epi32(transpose1, perm);
_mm256_storeu_si256((__m256i*)&dst1[ 0], final0);
_mm256_storeu_si256((__m256i*)&dst1[32], final1);
}
И вот вывод компилятора:
tran8x8b_AVX2_v2(char*, char*):
vpcmpeqd ymm3, ymm3, ymm3
vmovdqa ymm2, YMMWORD PTR .LC0[rip]
vmovdqa ymm4, ymm3
vpgatherdd ymm0, DWORD PTR [rdi+4+ymm2*8], ymm3
vpgatherdd ymm1, DWORD PTR [rdi+ymm2*8], ymm4
vmovdqa ymm2, YMMWORD PTR .LC1[rip]
vpshufb ymm1, ymm1, ymm2
vpshufb ymm0, ymm0, ymm2
vmovdqa ymm2, YMMWORD PTR .LC2[rip]
vpermd ymm1, ymm2, ymm1
vpermd ymm0, ymm2, ymm0
vmovdqu YMMWORD PTR [rsi], ymm1
vmovdqu YMMWORD PTR [rsi+32], ymm0
vzeroupper
ret
Упрощенный
void tp128_8x8(char *A, char *B) {
__m128i sv = _mm_set_epi8(15, 7, 14, 6, 13, 5, 12, 4, 11, 3, 10, 2, 9, 1, 8, 0);
__m128i iv[4], ov[4];
ov[0] = _mm_shuffle_epi8(_mm_loadu_si128((__m128i*)A), sv);
ov[1] = _mm_shuffle_epi8(_mm_loadu_si128((__m128i*)(A+16)), sv);
ov[2] = _mm_shuffle_epi8(_mm_loadu_si128((__m128i*)(A+32)), sv);
ov[3] = _mm_shuffle_epi8(_mm_loadu_si128((__m128i*)(A+48)), sv);
iv[0] = _mm_unpacklo_epi16(ov[0], ov[1]);
iv[1] = _mm_unpackhi_epi16(ov[0], ov[1]);
iv[2] = _mm_unpacklo_epi16(ov[2], ov[3]);
iv[3] = _mm_unpackhi_epi16(ov[2], ov[3]);
_mm_storeu_si128((__m128i*)B, _mm_unpacklo_epi32(iv[0], iv[2]));
_mm_storeu_si128((__m128i*)(B+16), _mm_unpackhi_epi32(iv[0], iv[2]));
_mm_storeu_si128((__m128i*)(B+32), _mm_unpacklo_epi32(iv[1], iv[3]));
_mm_storeu_si128((__m128i*)(B+48), _mm_unpackhi_epi32(iv[1], iv[3]));
}
Benchmark:i5-5300U 2.3GHz (cycles per byte)
tran8x8b : 2.140
tran8x8b_SSE : 1.602
tran8x8b_SSE_v2 : 1.551
tp128_8x8 : 1.535
tran8x8b_AVX2 : 1.563
tran8x8b_AVX2_v2 : 1.731
Мне это было действительно интересно, и я хотел сделать именно это, но по разным причинам мне пришлось делать это в Go вместо C, и у меня не было векторных встроенных функций, поэтому я подумал: "ну, Я просто напишу что-нибудь и посмотрю, как будет ".
Мои сообщения о времени на CPU ~ 3,6 ГГц составляют около 28 нс на 64-байтовый блок, транспонированный для наивной реализации, и около 19 нс для каждого блока, выполненного с использованием битовых сдвигов. Я использовал perf для подтверждения чисел, что мне показалось маловероятным, и они, похоже, сходятся. Реализация причудливого сдвига битов составляет чуть более 250 инструкций и дает около 3,6 инструкций за цикл, так что получается около 69-70 циклов на операцию.
Это Go, но, честно говоря, реализовать его несложно; он просто обрабатывает входной массив из 64 байтов как 8 uint64_t.
Вы можете получить еще одну наносекунду или около того, объявив некоторые из этих вещей как новые переменные, которые будут намекать распределителю регистров.
import (
"unsafe"
)
const (
hi16 = uint64(0xFFFF0000FFFF0000)
lo16 = uint64(0x0000FFFF0000FFFF)
hi8 = uint64(0xFF00FF00FF00FF00)
lo8 = uint64(0x00FF00FF00FF00FF)
)
// Okay, this might take some explaining. We are working on a logical
// 8x8 matrix of bytes, which we get as a 64-byte array. We want to transpose
// it (row/column).
//
// start:
// [[00 08 16 24 32 40 48 56]
// [01 09 17 25 33 41 49 57]
// [02 10 18 26 34 42 50 58]
// [03 11 19 27 35 43 51 59]
// [04 12 20 28 36 44 52 60]
// [05 13 21 29 37 45 53 61]
// [06 14 22 30 38 46 54 62]
// [07 15 23 31 39 47 55 63]]
//
// First, let's make sure everything under 32 is in the top four rows,
// and everything over 32 is in the bottom four rows. We do this by
// swapping pairs of 32-bit words.
// swap32:
// [[00 08 16 24 04 12 20 28]
// [01 09 17 25 05 13 21 29]
// [02 10 18 26 06 14 22 30]
// [03 11 19 27 07 15 23 31]
// [32 40 48 56 36 44 52 60]
// [33 41 49 57 37 45 53 61]
// [34 42 50 58 38 46 54 62]
// [35 43 51 59 39 47 55 63]]
//
// Next, let's make sure everything over 16 or 48 is in the bottom two
// rows of the two four-row sections, and everything under 16 or 48 is
// in the top two rows of the section. We do this by swapping masked
// pairs in much the same way:
// swap16:
// [[00 08 02 10 04 12 06 14]
// [01 09 03 11 05 13 07 15]
// [16 24 18 26 20 28 22 30]
// [17 25 19 27 21 29 23 31]
// [32 40 34 42 36 44 38 46]
// [33 41 35 43 37 45 39 47]
// [48 56 50 58 52 60 54 62]
// [49 57 51 59 53 61 55 63]]
//
// Now, we will do the same thing to each pair -- but because of
// clever choices in the specific arrange ment leading up to this, that's
// just one more byte swap, where each 2x2 block has its upper right
// and lower left corners swapped, and that turns out to be an easy
// shift and mask.
func UnswizzleLazy(m *[64]uint8) {
// m32 treats the 8x8 array as a 2x8 array, because
// it turns out we only need to swap a handful of the
// bits...
m32 := (*[16]uint32)(unsafe.Pointer(&m[0]))
m32[1], m32[8] = m32[8], m32[1]
m32[3], m32[10] = m32[10], m32[3]
m32[5], m32[12] = m32[12], m32[5]
m32[7], m32[14] = m32[14], m32[7]
m64 := (*[8]uint64)(unsafe.Pointer(&m[0]))
// we're now at the state described above as "swap32"
tmp0, tmp1, tmp2, tmp3 :=
(m64[0]&lo16)|(m64[2]&lo16)<<16,
(m64[1]&lo16)|(m64[3]&lo16)<<16,
(m64[0]&hi16)>>16|(m64[2]&hi16),
(m64[1]&hi16)>>16|(m64[3]&hi16)
tmp4, tmp5, tmp6, tmp7 :=
(m64[4]&lo16)|(m64[6]&lo16)<<16,
(m64[5]&lo16)|(m64[7]&lo16)<<16,
(m64[4]&hi16)>>16|(m64[6]&hi16),
(m64[5]&hi16)>>16|(m64[7]&hi16)
// now we're at "swap16".
lo8 := lo8
hi8 := hi8
m64[0], m64[1] = (tmp0&lo8)|(tmp1&lo8)<<8, (tmp0&hi8)>>8|tmp1&hi8
m64[2], m64[3] = (tmp2&lo8)|(tmp3&lo8)<<8, (tmp2&hi8)>>8|tmp3&hi8
m64[4], m64[5] = (tmp4&lo8)|(tmp5&lo8)<<8, (tmp4&hi8)>>8|tmp5&hi8
m64[6], m64[7] = (tmp6&lo8)|(tmp7&lo8)<<8, (tmp6&hi8)>>8|tmp7&hi8
}
Что это делает, я надеюсь, достаточно очевидно: перемешайте полуслова так, чтобы первые четыре слова имели все значения, которые им принадлежат, а последние четыре - все значения, которые им принадлежат. Затем проделайте то же самое с каждым набором из четырех слов, чтобы в итоге получить то, что входит в два верхних слова в два верхних и т. Д.
Я не собирался комментировать, пока не понял, что, если числа циклов / байтов выше, это на самом деле превосходит решение тасования / распаковки.
(Обратите внимание, что это транспонирование на месте, но на промежуточных этапах легко использовать temps, а последнее сохранить где-то еще. На самом деле это, вероятно, быстрее.)
ОБНОВЛЕНИЕ: я изначально немного неправильно описал свой алгоритм, затем я понял, что действительно могу делать то, что описал. Это работает около 65,7 цикла на 64 бита.
РЕДАКТИРОВАТЬ № 2: Пробовал одну из вышеуказанных версий AVX на этой машине. На моем оборудовании (Xeon E3-1505M, номинально 3 ГГц) я получаю чуть более 10 циклов на 64-байтовый блок, то есть около 6 байтов на цикл. Мне это кажется более разумным, чем 1,5 цикла на байт.
РЕДАКТИРОВАТЬ № 3: Я пошел немного дальше, примерно 45 циклов на 64 бита, просто записав первую часть в виде сдвигов и масок на uint64, вместо того, чтобы пытаться быть "умным" и просто переместить 32 бита, о которых я заботился.
Обычно, когда команды загрузки и сохранения не учитываются, это происходит потому, что код работает с матрицей в регистре, например, выполняет несколько операций в дополнение к транспонированию в цикле. Нагрузки и накопители в этом случае не учитываются, поскольку они не являются частью основного цикла.
Но в вашем коде загрузки и хранилища (или, скорее, наборы и выдержки) выполняют часть транспонирования.
GCC реализует _mm_set_epi64x для SSE4.1 в вашем коде с _mm_insert_epi64 а также _mm_loadl_epi64, Инструкция вставки выполняет часть транспонирования, т.е. транспонирование начинается с load0,1,2,3 не в shuffle0,1,2,3, И тогда ваш финал store0,1,2,3 значения также не содержат транспонирования. Вы должны использовать восемь _mm_extract_epi64 инструкции для завершения транспонирования в памяти. Так что на самом деле не имеет смысла не считать набор и извлекать внутренности.
В любом случае оказывается, что вы можете выполнить транспонирование из регистра только с 16 инструкциями, используя только SSSE3, например:
//__m128i B0, __m128i B1, __m128i B2, __m128i B3
__m128i mask = _mm_setr_epi8(0x0,0x04,0x01,0x05, 0x02,0x06,0x03,0x07, 0x08,0x0c,0x09,0x0d, 0x0a,0x0e,0x0b,0x0f);
__m128i T0, T1, T2, T3;
T0 = _mm_unpacklo_epi8(B0,B1);
T1 = _mm_unpackhi_epi8(B0,B1);
T2 = _mm_unpacklo_epi8(B2,B3);
T3 = _mm_unpackhi_epi8(B2,B3);
B0 = _mm_unpacklo_epi16(T0,T2);
B1 = _mm_unpackhi_epi16(T0,T2);
B2 = _mm_unpacklo_epi16(T1,T3);
B3 = _mm_unpackhi_epi16(T1,T3);
T0 = _mm_unpacklo_epi32(B0,B2);
T1 = _mm_unpackhi_epi32(B0,B2);
T2 = _mm_unpacklo_epi32(B1,B3);
T3 = _mm_unpackhi_epi32(B1,B3);
B0 = _mm_shuffle_epi8(T0,mask);
B1 = _mm_shuffle_epi8(T1,mask);
B2 = _mm_shuffle_epi8(T2,mask);
B3 = _mm_shuffle_epi8(T3,mask);
Я не уверен, имеет ли смысл исключать загрузки и хранить их здесь, потому что я не уверен, насколько удобно работать с матрицей 8x8 байтов в четырех 128-битных регистрах.
Вот тестирование кода это:
#include <stdio.h>
#include <x86intrin.h>
void print8x8b(char *A) {
for(int i=0; i<8; i++) {
for(int j=0; j<8; j++) {
printf("%2d ", A[i*8+j]);
} puts("");
} puts("");
}
void tran8x8b(char *A, char *B) {
for(int i=0; i<8; i++) {
for(int j=0; j<8; j++) {
B[j*8+i] = A[i*8+j];
}
}
}
void tran8x8b_SSE(char *A, char *B) {
__m128i mask = _mm_setr_epi8(0x0,0x04,0x01,0x05, 0x02,0x06,0x03,0x07, 0x08,0x0c,0x09,0x0d, 0x0a,0x0e,0x0b,0x0f);
__m128i B0, B1, B2, B3, T0, T1, T2, T3;
B0 = _mm_loadu_si128((__m128i*)&A[ 0]);
B1 = _mm_loadu_si128((__m128i*)&A[16]);
B2 = _mm_loadu_si128((__m128i*)&A[32]);
B3 = _mm_loadu_si128((__m128i*)&A[48]);
T0 = _mm_unpacklo_epi8(B0,B1);
T1 = _mm_unpackhi_epi8(B0,B1);
T2 = _mm_unpacklo_epi8(B2,B3);
T3 = _mm_unpackhi_epi8(B2,B3);
B0 = _mm_unpacklo_epi16(T0,T2);
B1 = _mm_unpackhi_epi16(T0,T2);
B2 = _mm_unpacklo_epi16(T1,T3);
B3 = _mm_unpackhi_epi16(T1,T3);
T0 = _mm_unpacklo_epi32(B0,B2);
T1 = _mm_unpackhi_epi32(B0,B2);
T2 = _mm_unpacklo_epi32(B1,B3);
T3 = _mm_unpackhi_epi32(B1,B3);
B0 = _mm_shuffle_epi8(T0,mask);
B1 = _mm_shuffle_epi8(T1,mask);
B2 = _mm_shuffle_epi8(T2,mask);
B3 = _mm_shuffle_epi8(T3,mask);
_mm_storeu_si128((__m128i*)&B[ 0], B0);
_mm_storeu_si128((__m128i*)&B[16], B1);
_mm_storeu_si128((__m128i*)&B[32], B2);
_mm_storeu_si128((__m128i*)&B[48], B3);
}
int main(void) {
char A[64], B[64], C[64];
for(int i=0; i<64; i++) A[i] = i;
print8x8b(A);
tran8x8b(A,B);
print8x8b(B);
tran8x8b_SSE(A,C);
print8x8b(C);
}
В большинстве ответов здесь используется комбинация перетасовок и перестановок разного размера с использованием
_mm_shuffle_epi8, который доступен только в SSSE3 и выше.
Чистая реализация SSE2 с ядром из 12* инструкций может быть сформирована путем чередования первых 32 элементов с последними 32 элементами три раза подряд:
void systolic_kernel(__m128i a[4]) {
__m128i a0 = _mm_unpacklo_epi8(a[0], a[2]);
__m128i a1 = _mm_unpackhi_epi8(a[0], a[2]);
__m128i a2 = _mm_unpacklo_epi8(a[1], a[3]);
__m128i a3 = _mm_unpackhi_epi8(a[1], a[3]);
a[0] = a0;
a[1] = a1;
a[2] = a2;
a[3] = a3;
}
void transpose(__m128i a[4]) {
systolic_kernel(a);
systolic_kernel(a);
systolic_kernel(a);
}
*без кодирования VEX (для инструкций с тремя операндами) будет 6 потенциально нулевых затрат
movdqaдобавлены инструкции.
Та же стратегия может быть легче применена к транспонированию 4x4, 16x16 и т. д., поскольку вычисление индексов, подлежащих перестановке, и размеров блоков не учитывается в уравнении.
AVX512VBMI (Каскадное озеро / Ледяное озеро)
AVX512VBMI представляет vpermb, тасование с пересечением полосы 64 байтами с байтовой детализацией._mm512_permutexvar_epi8( __m512i idx, __m512i a);
Существующие процессоры, которые его поддерживают, запускают его как единый uop с пропускной способностью 1/ такт. (https://www.uops.info/html-tp/CNL/VPERMB_ZMM_ZMM_M512-Measurements.html)
Это упрощает задачу, делая возможным с помощью 1 инструкции (по крайней мере, для случая stride=8, когда весь блок 8x8 является смежным). В противном случае вам следует посмотреть наvpermt2bперемешать байты из 2 источников. Но это 3 мопа на Кэннон-Лейк.
// TODO: strided loads / stores somehow for stride != 8
// AVX512VBMI
void TransposeBlock8x8_contiguous(uint8_t *src, uint8_t *dst)
{
const __m512i trans8x8shuf = _mm512_set_epi8(
63, 63-8*1, 63-8*2, 63-8*3, 63-8*4, ...
...
57, 49, 41, 33, 25, 17, 9, 1,
56, 48, 40, 32, 24, 16, 8, 0
);
__m512i vsrc = _mm512_loadu_si512(src);
__m512i shuffled = _mm512_permutexvar_epi8(trans8x8shuf, vsrc);
_mm512_storeu_si512(dst, shuffled);
}
По-видимому _mm512_setr_epi8 не существует для gcc/clang (только версии 256 и 128), поэтому вам нужно определить константу в порядке от последнего к первому, в противоположность порядку инициализатора массива C.
vpermbдаже работает с данными в качестве операнда источника памяти, поэтому он может загружать + перемешивать в одной инструкции. Но согласно https://uops.info/, на CannonLake он не срабатывает: в отличие от vpermd zmm, zmm, [r14] который декодируется в 1 uop с объединенным доменом (обратите внимание на "retire_slots: 1.0")
vpermd zmm, zmm, [r14] декодирует до 2 отдельных мопов для внешнего / объединенного домена: "retire_slots: 2.0"). Это результат экспериментального тестирования с использованием счетчиков производительности на реальном процессоре CannonLake. На uops.info пока нет доступных Cascade Lake или Ice Lake, так что, возможно, там будет еще эффективнее.
Таблицы uops.info бесполезно подсчитывают общее количество uops неиспользуемых доменов, поэтому вам нужно щелкнуть инструкцию, чтобы увидеть, срабатывают ли она микропредохранителями или нет.
Исходные данные или данные с постепенным переходом, не смежные, данные
Я предполагаю, что вы захотите выполнить qword (8-байтовую) загрузку в регистры XMM и перемешать пары входных данных или объединить их с movhps или pinsrq. Возможно, стоит использовать загрузку qword-gather с индексированными индексами, но часто это того не стоит.
Я не уверен, стоит ли комбинировать регистры YMM, не говоря уже о ZMM, или лучше всего получить только регистры XMM, чтобы мы могли эффективно разбрасывать qword обратно в память вручную с помощью vmovq а также vmovhps(что не требует перетасовки, просто магазин на процессорах Intel). Если dst является непрерывным, объединение несмежных src с последовательным разделением имеет гораздо больший смысл.
AVX512VBMI vpermt2b ymm похоже, это было бы полезно для перемешивания + слияния, как при пересечении полосы движения punpcklbw, выбирая любые 32 байта из конкатенации двух других 32-байтовых регистров YMM. (Или 64 из 2x 64-байтовых регистров для версии ZMM). Но, к сожалению, на CannonLake он стоит 3 мопса, напримерvpermt2w на Skylake-X и Cannon Lake.
Если мы можем позже побеспокоиться о байтах, vpermt2dэффективен на процессорах, которые его поддерживают (одиночный uop)! Skylake-X и новее.
Ice Lake имеет пропускную способность за 2 цикла для vpermt2b(instlat), возможно, потому что у него есть дополнительный модуль перемешивания, который может запускать некоторые (но не все) команды перемешивания. Обратите внимание, например, чтоvpshufb xmm а также ymm составляет 0,5 с, но vpshufb zmmпропускная способность 1с. Ноvpermb всегда пропускная способность всего 1с.
Интересно, можем ли мы воспользоваться маскированием слияния? Может бытьvpmovzxbqобнулить входные байты до qwords. (одна 8-байтовая строка -> 64-байтовый регистр ZMM). Тогда, может быть, двойной сдвиг влево с маскированием слияния в другой регистр? Нет, это не помогает, полезные данные находятся в одних и тех же элементах двойного слова для обоих входов, если только вы сначала не сделаете что-то с одним регистром, что противоречит цели.
Перекрывающиеся хранилища с байтовой маской (vmovdqu8 [rdi + 0..7]{k1}, zmm0..7) из vpmovzxbqрезультаты загрузки также возможны, но, вероятно, неэффективны. Все, кроме одного, в лучшем случае будут смещены. Однако аппаратный буфер хранилища и / или кеш-память может эффективно фиксировать 8-кратное замаскированное хранилище.
Может быть интересна гибридная стратегия с перемещением по регистрам и некоторым замаскированным хранилищам, чтобы сбалансировать перемешивание / смешивание и работу с хранилищами для смежных dst. Особенно, если все хранилища могут быть выровнены, что потребует перемещения данных в каждом векторе, чтобы они находились в нужном месте.
Ice Lake имеет 2 магазина исполнения. (IDK, если фиксация кэша L1d может поспевать за этим, или если слияние в буфере хранилища обычно помогает, или если это просто помогает с всплесками работы.)