转置 8x8 64 位矩阵

Posted

技术标签:

【中文标题】转置 8x8 64 位矩阵【英文标题】:Transpose 8x8 64-bits matrix 【发布时间】:2021-03-23 15:47:09 【问题描述】:

针对 AVX2,转置包含 64 位整数(或双精度)的 8x8 矩阵的最快方法是什么?

我搜索了这个网站,发现了几种进行 8x8 转置的方法,但主要用于 32 位浮点数。所以我主要是问,因为我不确定使这些算法快速转换为 64 位的原理是否容易,其次,显然 AVX2 只有 16 个寄存器,因此仅加载所有值会占用所有寄存器。

一种方法是调用 2x2 _MM_TRANSPOSE4_PD,但我想知道这是否是最佳选择:

  #define _MM_TRANSPOSE4_PD(row0,row1,row2,row3)                \
                                                               \
            __m256d tmp3, tmp2, tmp1, tmp0;                     \
                                                                \
            tmp0 = _mm256_shuffle_pd((row0),(row1), 0x0);       \
            tmp2 = _mm256_shuffle_pd((row0),(row1), 0xF);       \
            tmp1 = _mm256_shuffle_pd((row2),(row3), 0x0);       \
            tmp3 = _mm256_shuffle_pd((row2),(row3), 0xF);       \
                                                                \
            (row0) = _mm256_permute2f128_pd(tmp0, tmp1, 0x20);  \
            (row1) = _mm256_permute2f128_pd(tmp2, tmp3, 0x20);  \
            (row2) = _mm256_permute2f128_pd(tmp0, tmp1, 0x31);  \
            (row3) = _mm256_permute2f128_pd(tmp2, tmp3, 0x31);  \
        

仍然假设 AVX2,转置 double[8][8]int64_t[8][8] 在原则上基本相同吗?

PS:只是好奇,拥有 AVX512 会大大改变事情,对吗?

【问题讨论】:

您的输入是来自内存还是寄存器?您需要内存或寄存器中的输出,还是只处理它们? (您可能根本不需要显式转置数据)。 根据您的需要,这可能是重复的:***.com/questions/58454741/… @chtz 他们来自记忆,应该去记忆。抱歉,忘记添加了。 【参考方案1】:

在 cmets 中经过一些思考和讨论后,我认为这是最有效的版本,至少当源和目标数据都在 RAM 中时。不需要AVX2,AVX1就够了。

主要思想是,与存储相比,现代 CPU 可以执行两倍的加载微操作,并且在许多 CPU 上,使用vinsertf128 将内容加载到较高一半的向量中的成本与常规的 16 字节加载相同。与您的宏相比,此版本不再需要这些相对昂贵的(大多数 CPU 上的 3 个延迟周期)vperm2f128 shuffle。

struct Matrix4x4

    __m256d r0, r1, r2, r3;
;

inline void loadTransposed( Matrix4x4& mat, const double* rsi, size_t stride = 8 )

    // Load top half of the matrix into low half of 4 registers
    __m256d t0 = _mm256_castpd128_pd256( _mm_loadu_pd( rsi ) );     // 00, 01
    __m256d t1 = _mm256_castpd128_pd256( _mm_loadu_pd( rsi + 2 ) ); // 02, 03
    rsi += stride;
    __m256d t2 = _mm256_castpd128_pd256( _mm_loadu_pd( rsi ) );     // 10, 11
    __m256d t3 = _mm256_castpd128_pd256( _mm_loadu_pd( rsi + 2 ) ); // 12, 13
    rsi += stride;
    // Load bottom half of the matrix into high half of these registers
    t0 = _mm256_insertf128_pd( t0, _mm_loadu_pd( rsi ), 1 );    // 00, 01, 20, 21
    t1 = _mm256_insertf128_pd( t1, _mm_loadu_pd( rsi + 2 ), 1 );// 02, 03, 22, 23
    rsi += stride;
    t2 = _mm256_insertf128_pd( t2, _mm_loadu_pd( rsi ), 1 );    // 10, 11, 30, 31
    t3 = _mm256_insertf128_pd( t3, _mm_loadu_pd( rsi + 2 ), 1 );// 12, 13, 32, 33

    // Transpose 2x2 blocks in registers.
    // Due to the tricky way we loaded stuff, that's enough to transpose the complete 4x4 matrix.
    mat.r0 = _mm256_unpacklo_pd( t0, t2 ); // 00, 10, 20, 30
    mat.r1 = _mm256_unpackhi_pd( t0, t2 ); // 01, 11, 21, 31
    mat.r2 = _mm256_unpacklo_pd( t1, t3 ); // 02, 12, 22, 32
    mat.r3 = _mm256_unpackhi_pd( t1, t3 ); // 03, 13, 23, 33


inline void store( const Matrix4x4& mat, double* rdi, size_t stride = 8 )

    _mm256_storeu_pd( rdi, mat.r0 );
    _mm256_storeu_pd( rdi + stride, mat.r1 );
    _mm256_storeu_pd( rdi + stride * 2, mat.r2 );
    _mm256_storeu_pd( rdi + stride * 3, mat.r3 );


// Transpose 8x8 matrix of double values
void transpose8x8( double* rdi, const double* rsi )

    Matrix4x4 block;
    // Top-left corner
    loadTransposed( block, rsi );
    store( block, rdi );

#if 1
    // Using another instance of the block to support in-place transpose
    Matrix4x4 block2;
    loadTransposed( block, rsi + 4 );       // top right block
    loadTransposed( block2, rsi + 8 * 4 ); // bottom left block

    store( block2, rdi + 4 );
    store( block, rdi + 8 * 4 );
#else
    // Flip the #if if you can guarantee ( rsi != rdi )
    // Performance is about the same, but this version uses 4 less vector registers,
    // slightly more efficient when some registers need to be backed up / restored.
    assert( rsi != rdi );
    loadTransposed( block, rsi + 4 );
    store( block, rdi + 8 * 4 );

    loadTransposed( block, rsi + 8 * 4 );
    store( block, rdi + 4 );
#endif
    // Bottom-right corner
    loadTransposed( block, rsi + 8 * 4 + 4 );
    store( block, rdi + 8 * 4 + 4 );


为了完整起见,这里有一个版本,它使用的代码与您的宏非常相似,加载次数减少了两倍,存储次数相同,并且洗牌次数更多。尚未进行基准测试,但我希望它会稍微慢一些。

struct Matrix4x4

    __m256d r0, r1, r2, r3;
;

inline void load( Matrix4x4& mat, const double* rsi, size_t stride = 8 )

    mat.r0 = _mm256_loadu_pd( rsi );
    mat.r1 = _mm256_loadu_pd( rsi + stride );
    mat.r2 = _mm256_loadu_pd( rsi + stride * 2 );
    mat.r3 = _mm256_loadu_pd( rsi + stride * 3 );


inline void store( const Matrix4x4& mat, double* rdi, size_t stride = 8 )

    _mm256_storeu_pd( rdi, mat.r0 );
    _mm256_storeu_pd( rdi + stride, mat.r1 );
    _mm256_storeu_pd( rdi + stride * 2, mat.r2 );
    _mm256_storeu_pd( rdi + stride * 3, mat.r3 );


inline void transpose( Matrix4x4& m4 )

    // These unpack instructions transpose lanes within 2x2 blocks of the matrix
    const __m256d t0 = _mm256_unpacklo_pd( m4.r0, m4.r1 );
    const __m256d t1 = _mm256_unpacklo_pd( m4.r2, m4.r3 );
    const __m256d t2 = _mm256_unpackhi_pd( m4.r0, m4.r1 );
    const __m256d t3 = _mm256_unpackhi_pd( m4.r2, m4.r3 );
    // Produce the transposed matrix by combining these blocks
    m4.r0 = _mm256_permute2f128_pd( t0, t1, 0x20 );
    m4.r1 = _mm256_permute2f128_pd( t2, t3, 0x20 );
    m4.r2 = _mm256_permute2f128_pd( t0, t1, 0x31 );
    m4.r3 = _mm256_permute2f128_pd( t2, t3, 0x31 );


// Transpose 8x8 matrix with double values
void transpose8x8( double* rdi, const double* rsi )

    Matrix4x4 block;
    // Top-left corner
    load( block, rsi );
    transpose( block );
    store( block, rdi );

    // Using another instance of the block to support in-place transpose, with very small overhead
    Matrix4x4 block2;
    load( block, rsi + 4 );     // top right block
    load( block2, rsi + 8 * 4 ); // bottom left block

    transpose( block2 );
    store( block2, rdi + 4 );
    transpose( block );
    store( block, rdi + 8 * 4 );

    // Bottom-right corner
    load( block, rsi + 8 * 4 + 4 );
    transpose( block );
    store( block, rdi + 8 * 4 + 4 );

【讨论】:

您能否使用 vmovdqu / vinsertf128 加载以设置 256 位存储,而不是 256 位加载和拆分存储?大多数现代 CPU 每个时钟可以执行 2 次加载,但每个时钟只能执行 1 个存储(Ice Lake 是每个 2 个)。在 Zen 2 上,vinsertf128 y,mem 是前端的单微指令。 (在 Intel 上它是 2,并且在后端需要一个 ALU uop 以及负载,但该 uop 可以在任何矢量 ALU 端口上运行,可能就像一个广播负载馈送混合)。 @PeterCordes 好主意,添加了另一个版本。 为什么不先显示可能更好的版本,然后将旧版本下移?您的回答似乎确实同意我的观点,即这是更好的版本,因此 IMO 它应该位于顶部,或者像其他 8x2 方式一样位于 4x4 方式之后,未来的读者将首先看到这两个选项。然后你可以展示更多商店的方式并讨论为什么它可能更糟。或者只是删除代码并提及将商店数量增加 2 倍会更糟,除非您的目的地无法与 32 对齐。 @PeterCordes 没有测试过而且很难做到(取决于 CPU 和周围的代码),但是是的,我倾向于同意。已更新。 @PeterCordes 我明白为什么随机写入比顺序写入慢,只是随机读取太糟糕了。除非缓存,否则很可能会发生管道停顿。当我处理各种 RAM 绑定代码时,我不记得有一次用随机读取换随机写入不是一个好主意。【参考方案2】:

对于单个 SIMD 向量中可以容纳超过 1 行的小型矩阵,AVX-512 具有非常好的 2 输入通道交叉混洗,具有 32 位或 64 位粒度,并带有向量控制。 (与 _mm512_unpacklo_pd 不同,它基本上是 4 个单独的 128 位随机播放。)

一个 4x4 double 矩阵“只有”128 个字节,两个 ZMM __m512d 向量,因此您只需要两个 vpermt2ps (_mm512_permutex2var_pd) 即可生成两个输出向量:一个每个输出向量的随机播放,加载和存储都是全宽的。不过,您确实需要控制向量常量。

使用 512 位向量指令有一些缺点(时钟速度和执行端口吞吐量),但如果您的程序可以在使用 512 位向量的代码中花费大量时间,那么通过使用更多的指令可能会显着提高吞吐量每条指令的数据,并具有更强大的随机播放。

对于 256 位向量,vpermt2pd ymm 可能对 4x4 没有用,因为对于每个 __m256d 输出行,您想要的 4 个元素中的每一个都来自不同的输入行。所以一个 2-input shuffle 不能产生你想要的输出。

我认为小于 128 位粒度的车道交叉洗牌没有用处,除非您的矩阵小到足以在一个 SIMD 向量中容纳多行。请参阅 How to transpose a 16x16 matrix using SIMD instructions? 了解一些算法复杂性关于 32 位元素的推理 - 使用 AVX1 的 32 位元素的 8x8 xpose 与使用 AVX-512 的 64 位元素的 8x8 大致相同,其中每个 SIMD 向量恰好包含一整行。

因此不需要向量常量,只需立即对 128 位块进行洗牌,以及unpacklo/hi


用 512 位向量(8 个双精度)转置一个 8x8 会遇到同样的问题:8 个双精度的每个输出行都需要 8 个输入向量中的每个向量的 1 个双精度。 所以最终我认为您需要与 Soonts 的 AVX 答案类似的策略,从 _mm512_insertf64x4(v, load, 1) 开始作为将 2 个输入行的前半部分放入一个向量的第一步。

(如果您关心 KNL / Xeon Phi,@ZBoson 在How to transpose a 16x16 matrix using SIMD instructions? 上的另一个答案显示了一些有趣的想法,即使用合并掩码与 vpermpdvpermq 等 1 输入随机播放,而不是像 2 输入随机播放vunpcklpdvpermt2pd)

使用更宽的向量意味着更少的加载和存储,甚至可能更少的总洗牌次数,因为每一个都结合了更多的数据。但是您还需要做更多的改组工作,将一行的所有 8 个元素放入一个向量中,而不是仅仅以一行一半大小的块加载和存储到不同的位置。不明显更好;如果我有时间实际编写代码,我会更新这个答案。

请注意,Ice Lake(第一个使用 AVX-512 的消费级 CPU)每个时钟可以执行 2 次加载和 2 次存储。对于 some shuffle,它比 Skylake-X 具有更好的 shuffle 吞吐量,但对于任何对此或 Soonts 的答案有用的都没有。 (对于 ymm 和 zmm 版本,所有vperm2f128vunpcklpdvpermt2pd 仅在端口 5 上运行。https://uops.info/。vinsertf64x4 zmm, mem, 1 是前端的 2 微指令,需要一个加载端口和一个uop 用于 p0/p5。(不是 p1,因为它是 512 位 uop,另请参阅 SIMD instructions lowering CPU frequency。)

【讨论】:

以上是关于转置 8x8 64 位矩阵的主要内容,如果未能解决你的问题,请参考以下文章

仅使用 avx 而不是 avx2 转置 64 位元素

nyoj 29-求转置矩阵问题 (行,列位置调换)

python 矩阵转置transpose

matlab转置矩阵命令

矩阵与转置

matlab转置矩阵?