如何使用 SIMD 指令转置 16x16 矩阵？答案

【问题标题】：How to transpose a 16x16 matrix using SIMD instructions?如何使用 SIMD 指令转置 16x16 矩阵？
【发布时间】：2015-06-13 16:20:40
【问题描述】：

我目前正在编写一些针对英特尔即将推出的支持 512 位操作的 AVX-512 SIMD 指令的代码。

现在假设有一个由 16 个 SIMD 寄存器表示的矩阵，每个寄存器包含 16 个 32 位整数（对应于一行），我如何用纯 SIMD 指令转置矩阵？

已经有解决方案分别用 SSE 和 AVX2 转置 4x4 或 8x8 矩阵。但我不知道如何使用 AVX-512 将其扩展到 16x16。

有什么想法吗？

【问题讨论】：

通常最快的做某事的方法是什么都不做——基本上，给每个矩阵一个“转置”标志，然后反转那个标志。当然，这意味着您需要检查“转置”标志并在任何其他可能处理转置矩阵的代码中交换列索引和行索引。例如。如果您有添加 2 个矩阵的函数，您最终可能会遇到 3 种情况（既没有转置，也没有转置，都转置），其中加法的结果始终是未转置的矩阵。
出于好奇，您能解释一下您为什么对 16x16 转置感兴趣吗？例如。这是用于更大转置的内核吗？读/写重要吗？还是生成的数据？
@Zboson 这是我们尝试使用 AVX512 优化的加密算法的一部分。事实上，我们可以在从内存加载时使用gather指令来转置矩阵。但是当没有收集/分散指令时，我们设法用 SSE/AVX2 做到了这一点，所以我只是好奇我们如何用 AVX512 做同样的事情，即寄存器内转置。
@Zboson KNL 提供了一些粗略的延迟/吞吐量数据。正如预期的那样，聚集/分散仍然很慢。 2 个元素/周期加载，1 个/周期存储。所以 8 个周期/float-gather 和 16 个周期/float-scatter。 IOW，收集/分散指令仍在为每个元素分解为单独的微指令并进入其适当的端口。它比前几代产品更有效率，前几代产品有大量其他额外的微指令。
@Mysticial 工作中的 HPC 小组给了我一个使用 AVX512 的 Knights Landing 卡帐户。我尝试了我的代码，它第一次尝试就成功了。很高兴知道。我还没有做任何性能测试。我在大约 30 分钟前获得了帐户。

标签： assembly matrix intel simd avx512

【解决方案1】：

对于使用 SIMD 的两条操作数指令，您可以证明转置 nxn 矩阵所需的操作数为 n*log_2(n)，而使用标量操作则为 O(n^2)。事实上，稍后我将展示使用标量寄存器的读写操作数为2*n*(n-1)。下表显示了与标量操作相比，使用 SSE、AVX、AVX512 和 AVX1024 转置 4x4、8x8、16x16 和 32x32 矩阵的操作数

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 次。

I already showed that 8x8 can be done with 24 operations with AVX。所以问题是如何在 64 次操作中使用 AVX512 为 16x16 执行此操作？一般算法是：

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

我通过查看使用_mm_shuffle_ps 转置4x4 矩阵（这是MSVC 在_MM_TRANSPOSE4_PS 而不是GCC 和ICC 中使用的），得到了使用_mm512_shufflei32x4 的想法。

__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 的 optimizing C++ manual 的示例 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)。

【讨论】：

＋1，虽然丑，但可能还是比使用 16 个收集负载要快。
@Mysticial, is it true that only xeon and workstation Skylake processors will have AVX512?如果是这种情况，那么 #@$!是Skylake的重点？？？如果这是真的，这是非常令人失望的消息。是什么让 Skylake 成为没有 AVX512 的“tock”？
是的，直到最近有关 Purley 的泄密事件，我才意识到情况如此糟糕。似乎它将是 2016 年第一季度的 Knights Landing 和 2017 年（后期？）带有 AVX512 的 Skylake Xeon。英特尔处理器通常分为笔记本/低端台式机（插槽 115x）和服务器/高端台式机（套接字 2011-x）线。 Skylake 的 AVX512 似乎只会出现在 Skylake 的服务器/高端桌面线上。这可能比笔记本/低端台式机的 Cannonlake 晚。
当然，我是根据最近的泄漏以及我对英特尔产品线的（有限）了解做出这些猜测的。所以我肯定是错的。可能会在 2015 年第三季度推出用于 1151 插槽的“Xeon Skylake”。但它可能只是一款经过美化的桌面处理器，所以我不确定它是否会有 AVX512。
顺便说一句，从两个向量而不是一个向量中提取的 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)

【解决方案2】：

我最近获得了具有 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-core)。这不是辅助卡。 Xeon Phi 是主计算机。

与https://stackoverflow.com/a/29587984/2542702 此处的方法相比，我测试了 AVX512 收集指令，但收集似乎仍然较慢。我在那个答案中的代码第一次尝试没有错误。

我已经有大约 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 的 Micro-architecture 手动 shuffle 和 permute 指令，KNL 的性能很差。在我的原始答案https://stackoverflow.com/a/29587984/2542702 中使用的 shuffle 和 unpack 指令的倒数吞吐量为 2。我设法通过使用vpermq 来大大提高性能，而不是倒数吞吐量为 1。此外，我使用 vinserti64x4 改进了转置的前 1/4（参见下面的 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;
  }
}

【讨论】：

不错！在您之前的回答中，您写道 8x8 转置 +r/w 使用 40 条指令。即：8 次加载，执行端口 5 和 8 次存储上的 24 次随机播放。在英特尔的文档 64-ia-32-architectures-optimization-manual 的第 11.11.2 段中，他们用内存操作符的 8 个 vinsertf128 指令替换了其中的 8 个随机播放。这导致端口 5 的压力较小：端口 5 上有 16 条指令。事实上，大量的 L1 带宽用于减少端口 5 上的瓶颈。结果是一个更快的算法。你认为你可以在这里使用类似的想法来加速 16x16 转置吗？
@wim 非常感谢您提供的链接！我连忙看了看。当我创建 8x8 答案时，我并没有考虑端口压力只是指令数量。我将不得不对此进行调查并回复您。
@wim：好主意。但是根据 Agner Fog 的表格，我认为 KNL 的 vinsert 带有内存源仍然需要随机播放单元。它基于 Silvermont，与 Haswell 非常不同。 Agner Fog 的表格没有列出 vinsertf128 或其 AVX512 变体的端口，但与 Haswell 一样，似乎只有一个 shuffle 单元。它在 FP0 上。 vinsertf32x4 z,z,m128/m256 都是每时钟一个吞吐量，而不是每 0.5c 负载一个，因此它们可能仍在使用随机播放单元。广播完全由加载端口处理，因此vbroadcastf64x4 z,m256 每 0.5c 吞吐量有一个。
@PeterCordes 实际上，KNL 上没有端口 5。随机播放到 FP0 单元。从 Agner 的手册中不清楚 vinsertf64x4 使用哪些资源。但至少我们可以通过从内存加载vbroadcastf6x4 加上vblendmpd 来模拟KNL 的vinsertf64x4，根据Agner Fog 的说法，这两者的吞吐量都是每0.5c 一个。 vblendmpd 在 FP0 或 FP1 上运行。所以，据我所知（我对 KNL 很不熟悉，我刚开始阅读 Agner 在 KNL 上的指令表），在两个周期内，我们可以在 FP0 上进行 2 次随机播放，并在内存端口和 FP1 上模拟vinsertf64x4 .
因此，前端不太可能成为这里的瓶颈。洗牌在 KNL 上相对昂贵。我仍然认为可以通过用 16 vinsertf64x4 替换（例如）16 个随机播放（共 64 个）来稍微加速 tran，或者，如果这不起作用，16 vbroadcastf64x4+ 16 @987654350 @.