Z b*_*son 20
对于使用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).
可以在n*log_2(n)操作中完成SIMD转置的原因是该算法是:
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可以通过AVX进行24次操作.那么问题是如何16x16在64次操作中使用AVX512呢？一般算法是:
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行.
最后,为了与标量操作进行比较,我从Agner Fog的优化C++手册中修改了例9.5a
#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交换(因为对角线不需要交换).16x16的装配交换看起来像
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).
我最近获得了配备 AVX512 的 Xeon Phi Knights Landing 硬件。具体来说，我使用的硬件是 Intel(R) Xeon Phi(TM) CPU 7250 @ 1.40GHz ( http://ark.intel.com/products/94035/Intel-Xeon-Phi-Processor-7250-16GB-1_40 -GHz-68 核）。这不是一张辅助卡。Xeon Phi 是主计算机。
与我的方法相比，我测试了 AVX512 收集指令/sf/answers/2071158911/，看来收集仍然较慢。我在该答案中的代码第一次尝试就成功了，没有错误。
我大约 3 个月没有编写内在函数，也没有在这段时间考虑太多优化问题，所以也许我的测试不够稳健。当然存在一些开销，但尽管如此，我相信结果清楚地表明在这种情况下聚集速度较慢。
我只使用 ICC 17.0.0 进行测试，因为当前安装的操作系统仅为 CentOS 7.2，Linux Kernel 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 中的性能较差。我原来的答案中使用的 shuffle 和 unpack 指令/sf/answers/2071158911/的倒数吞吐量为 2。我设法通过使用vpermq吞吐量倒数为 1 的替代方案来极大地提高性能。另外我使用 using 改进了转置的前 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_cas
-
好的！在您之前的回答中，您写道 8x8 转置 +r/w 使用 40 条指令。即：执行端口 5 上进行 8 次加载、24 次洗牌和 8 次存储。在英特尔文档 64-ia-32-architectures-optimization-manual 的第 11.11.2 段中，他们用带有内存操作的 8 个“vinsertf128”指令替换了其中 8 个混洗。这导致端口 5 的压力较小：端口 5 上有 16 条指令。实际上，大量的 L1 带宽用于减少端口 5 上的瓶颈。结果是更快的算法。你认为你可以在这里使用类似的想法来加速 16x16 转置吗？ (2认同)
-
因此，前端不太可能成为这里的瓶颈。KNL 上的洗牌费用相对较高。我仍然认为，通过用 16 个“vinsertf64x4”替换（例如）16 个随机播放（共 64 个），或者如果这不起作用，用 16 个“vbroadcastf64x4”+ 16 个“vblendmpd”，可以稍微加快“tran”的速度`。 (2认同)
| 归档时间： |
|
| 查看次数： |
4722 次 |
| 最近记录： |