使用AVX/AVX2转置一个8x8的浮点数

Question

使用AVX/AVX2转置一个8x8的浮点数

simdavxavx2

22

通过制作四个4x4矩阵并对它们进行转置，可以实现8x8矩阵的转置。但这不是我要的。

在另一个问题中，一个答案提供了一个解决方案，只需要24条指令就能处理8x8矩阵。然而，这不适用于浮点数。

由于AVX2包含256位寄存器，每个寄存器可以容纳八个32位整数（浮点数）。但问题是：

如何使用AVX/AVX2转置一个8x8浮点矩阵，并且使用最少的指令？

- DavidS

1

嗨，David。我认为我已经给出了你问题的确切答案。你还想要什么？ - Z boson

@Zboson：这些天我非常忙，但是我已经将你的答案标记为被接受的。非常感谢！ - DavidS

5个回答

8

这里有一个AVX2解决方案，适用于8 x 8 32位整数。如果您想转置8 x 8浮点数，则可以将浮点向量转换为整数，然后再转回来。可能也可以只针对浮点数进行AVX版本（即不需要AVX2），但我尚未尝试过。

//
// tranpose_8_8_avx2.c
//

#include <stdio.h>

#include <immintrin.h>

#define V_ELEMS 8

static inline void _mm256_merge_epi32(const __m256i v0, const __m256i v1, __m256i *vl, __m256i *vh)
{
    __m256i va = _mm256_permute4x64_epi64(v0, _MM_SHUFFLE(3, 1, 2, 0));
    __m256i vb = _mm256_permute4x64_epi64(v1, _MM_SHUFFLE(3, 1, 2, 0));
    *vl = _mm256_unpacklo_epi32(va, vb);
    *vh = _mm256_unpackhi_epi32(va, vb);
}

static inline void _mm256_merge_epi64(const __m256i v0, const __m256i v1, __m256i *vl, __m256i *vh)
{
    __m256i va = _mm256_permute4x64_epi64(v0, _MM_SHUFFLE(3, 1, 2, 0));
    __m256i vb = _mm256_permute4x64_epi64(v1, _MM_SHUFFLE(3, 1, 2, 0));
    *vl = _mm256_unpacklo_epi64(va, vb);
    *vh = _mm256_unpackhi_epi64(va, vb);
}

static inline void _mm256_merge_si128(const __m256i v0, const __m256i v1, __m256i *vl, __m256i *vh)
{
    *vl = _mm256_permute2x128_si256(v0, v1, _MM_SHUFFLE(0, 2, 0, 0));
    *vh = _mm256_permute2x128_si256(v0, v1, _MM_SHUFFLE(0, 3, 0, 1));
}

//
// Transpose_8_8
//
// in place transpose of 8 x 8 int array
//

static void Transpose_8_8(
    __m256i *v0,
    __m256i *v1,
    __m256i *v2,
    __m256i *v3,
    __m256i *v4,
    __m256i *v5,
    __m256i *v6,
    __m256i *v7)
{
    __m256i w0, w1, w2, w3, w4, w5, w6, w7;
    __m256i x0, x1, x2, x3, x4, x5, x6, x7;

    _mm256_merge_epi32(*v0, *v1, &w0, &w1);
    _mm256_merge_epi32(*v2, *v3, &w2, &w3);
    _mm256_merge_epi32(*v4, *v5, &w4, &w5);
    _mm256_merge_epi32(*v6, *v7, &w6, &w7);

    _mm256_merge_epi64(w0, w2, &x0, &x1);
    _mm256_merge_epi64(w1, w3, &x2, &x3);
    _mm256_merge_epi64(w4, w6, &x4, &x5);
    _mm256_merge_epi64(w5, w7, &x6, &x7);

    _mm256_merge_si128(x0, x4, v0, v1);
    _mm256_merge_si128(x1, x5, v2, v3);
    _mm256_merge_si128(x2, x6, v4, v5);
    _mm256_merge_si128(x3, x7, v6, v7);
}

int main(void)
{
    int32_t buff[V_ELEMS][V_ELEMS] __attribute__ ((aligned(32)));
    int i, j;
    int k = 0;

    // init buff
    for (i = 0; i < V_ELEMS; ++i)
    {
        for (j = 0; j < V_ELEMS; ++j)
        {
            buff[i][j] = k++;
        }
    }

    // print buff
    printf("\nBEFORE:\n");
    for (i = 0; i < V_ELEMS; ++i)
    {
        for (j = 0; j < V_ELEMS; ++j)
        {
            printf("%4d", buff[i][j]);
        }
        printf("\n");
    }

    // transpose
    Transpose_8_8((__m256i *)buff[0], (__m256i *)buff[1], (__m256i *)buff[2], (__m256i *)buff[3], (__m256i *)buff[4], (__m256i *)buff[5], (__m256i *)buff[6], (__m256i *)buff[7]);

    // print buff
    printf("\nAFTER:\n");
    for (i = 0; i < V_ELEMS; ++i)
    {
        for (j = 0; j < V_ELEMS; ++j)
        {
            printf("%4d", buff[i][j]);
        }
        printf("\n");
    }

    // transpose
    Transpose_8_8((__m256i *)buff[0], (__m256i *)buff[1], (__m256i *)buff[2], (__m256i *)buff[3], (__m256i *)buff[4], (__m256i *)buff[5], (__m256i *)buff[6], (__m256i *)buff[7]);

    // print buff
    printf("\nAFTER x2:\n");
    for (i = 0; i < V_ELEMS; ++i)
    {
        for (j = 0; j < V_ELEMS; ++j)
        {
            printf("%4d", buff[i][j]);
        }
        printf("\n");
    }

    return 0;
}

< p >Transpose_8_8 在使用clang编译时，大约需要56条指令，包括加载和存储操作。我认为通过更多的努力，应该能够进一步优化。

编译并测试：

$ gcc -Wall -mavx2 -O3 transpose_8_8_avx2.c && ./a.out

BEFORE:
   0   1   2   3   4   5   6   7
   8   9  10  11  12  13  14  15
  16  17  18  19  20  21  22  23
  24  25  26  27  28  29  30  31
  32  33  34  35  36  37  38  39
  40  41  42  43  44  45  46  47
  48  49  50  51  52  53  54  55
  56  57  58  59  60  61  62  63

AFTER:
   0   8  16  24  32  40  48  56
   1   9  17  25  33  41  49  57
   2  10  18  26  34  42  50  58
   3  11  19  27  35  43  51  59
   4  12  20  28  36  44  52  60
   5  13  21  29  37  45  53  61
   6  14  22  30  38  46  54  62
   7  15  23  31  39  47  55  63

AFTER x2:
   0   1   2   3   4   5   6   7
   8   9  10  11  12  13  14  15
  16  17  18  19  20  21  22  23
  24  25  26  27  28  29  30  31
  32  33  34  35  36  37  38  39
  40  41  42  43  44  45  46  47
  48  49  50  51  52  53  54  55
  56  57  58  59  60  61  62  63

$

- Paul R

4

我决定对三种不同的例行程序进行全面测试，以苹果对苹果的比较方式进行。

  // transpose.cpp : 
  /*
  Transpose8x8Shuff 100,000,000 times
  (0.750000 seconds).
  Transpose8x8Permute 100,000,000 times
  (0.749000 seconds).
  Transpose8x8Insert 100,000,000 times
  (0.858000 seconds).
  */

  #include <stdio.h>
  #include <time.h>
  #include <thread>
  #include <chrono>
  #include <xmmintrin.h>
  #include <emmintrin.h>
  #include <tmmintrin.h>
  #include <immintrin.h>

  inline void Transpose8x8Shuff(unsigned long *in)
  {
     __m256 *inI = reinterpret_cast<__m256 *>(in);
     __m256 rI[8];
     rI[0] = _mm256_unpacklo_ps(inI[0], inI[1]); 
     rI[1] = _mm256_unpackhi_ps(inI[0], inI[1]); 
     rI[2] = _mm256_unpacklo_ps(inI[2], inI[3]); 
     rI[3] = _mm256_unpackhi_ps(inI[2], inI[3]); 
     rI[4] = _mm256_unpacklo_ps(inI[4], inI[5]); 
     rI[5] = _mm256_unpackhi_ps(inI[4], inI[5]); 
     rI[6] = _mm256_unpacklo_ps(inI[6], inI[7]); 
     rI[7] = _mm256_unpackhi_ps(inI[6], inI[7]); 

     __m256 rrF[8];
     __m256 *rF = reinterpret_cast<__m256 *>(rI);
     rrF[0] = _mm256_shuffle_ps(rF[0], rF[2], _MM_SHUFFLE(1,0,1,0));
     rrF[1] = _mm256_shuffle_ps(rF[0], rF[2], _MM_SHUFFLE(3,2,3,2));
     rrF[2] = _mm256_shuffle_ps(rF[1], rF[3], _MM_SHUFFLE(1,0,1,0)); 
     rrF[3] = _mm256_shuffle_ps(rF[1], rF[3], _MM_SHUFFLE(3,2,3,2));
     rrF[4] = _mm256_shuffle_ps(rF[4], rF[6], _MM_SHUFFLE(1,0,1,0));
     rrF[5] = _mm256_shuffle_ps(rF[4], rF[6], _MM_SHUFFLE(3,2,3,2));
     rrF[6] = _mm256_shuffle_ps(rF[5], rF[7], _MM_SHUFFLE(1,0,1,0));
     rrF[7] = _mm256_shuffle_ps(rF[5], rF[7], _MM_SHUFFLE(3,2,3,2));

     rF = reinterpret_cast<__m256 *>(in);
     rF[0] = _mm256_permute2f128_ps(rrF[0], rrF[4], 0x20);
     rF[1] = _mm256_permute2f128_ps(rrF[1], rrF[5], 0x20);
     rF[2] = _mm256_permute2f128_ps(rrF[2], rrF[6], 0x20);
     rF[3] = _mm256_permute2f128_ps(rrF[3], rrF[7], 0x20);
     rF[4] = _mm256_permute2f128_ps(rrF[0], rrF[4], 0x31);
     rF[5] = _mm256_permute2f128_ps(rrF[1], rrF[5], 0x31);
     rF[6] = _mm256_permute2f128_ps(rrF[2], rrF[6], 0x31);
     rF[7] = _mm256_permute2f128_ps(rrF[3], rrF[7], 0x31);
  }

  inline void Transpose8x8Permute(unsigned long *in)
  {
     __m256i *inI = reinterpret_cast<__m256i *>(in);
     __m256i rI[8];
     rI[0] = _mm256_permute2f128_si256(inI[0], inI[4], 0x20); 
     rI[1] = _mm256_permute2f128_si256(inI[0], inI[4], 0x31); 
     rI[2] = _mm256_permute2f128_si256(inI[1], inI[5], 0x20); 
     rI[3] = _mm256_permute2f128_si256(inI[1], inI[5], 0x31); 
     rI[4] = _mm256_permute2f128_si256(inI[2], inI[6], 0x20); 
     rI[5] = _mm256_permute2f128_si256(inI[2], inI[6], 0x31); 
     rI[6] = _mm256_permute2f128_si256(inI[3], inI[7], 0x20); 
     rI[7] = _mm256_permute2f128_si256(inI[3], inI[7], 0x31); 

     __m256 rrF[8];
     __m256 *rF = reinterpret_cast<__m256 *>(rI);
     rrF[0] = _mm256_unpacklo_ps(rF[0], rF[4]); 
     rrF[1] = _mm256_unpackhi_ps(rF[0], rF[4]); 
     rrF[2] = _mm256_unpacklo_ps(rF[1], rF[5]); 
     rrF[3] = _mm256_unpackhi_ps(rF[1], rF[5]); 
     rrF[4] = _mm256_unpacklo_ps(rF[2], rF[6]); 
     rrF[5] = _mm256_unpackhi_ps(rF[2], rF[6]); 
     rrF[6] = _mm256_unpacklo_ps(rF[3], rF[7]); 
     rrF[7] = _mm256_unpackhi_ps(rF[3], rF[7]); 

     rF = reinterpret_cast<__m256 *>(in);
     rF[0] = _mm256_unpacklo_ps(rrF[0], rrF[4]); 
     rF[1] = _mm256_unpackhi_ps(rrF[0], rrF[4]); 
     rF[2] = _mm256_unpacklo_ps(rrF[1], rrF[5]); 
     rF[3] = _mm256_unpackhi_ps(rrF[1], rrF[5]); 
     rF[4] = _mm256_unpacklo_ps(rrF[2], rrF[6]); 
     rF[5] = _mm256_unpackhi_ps(rrF[2], rrF[6]); 
     rF[6] = _mm256_unpacklo_ps(rrF[3], rrF[7]); 
     rF[7] = _mm256_unpackhi_ps(rrF[3], rrF[7]); 
  }

  inline void Transpose8x8Insert(unsigned long *in)
  {
     __m256i *inIa = reinterpret_cast<__m256i *>(in);
     __m256i *inIb = reinterpret_cast<__m256i *>(&reinterpret_cast<__m128i *>(in)[1]);
     __m128i *inI128 = reinterpret_cast<__m128i *>(in);
     __m256i rI[8];
     rI[0] = _mm256_insertf128_si256(inIa[0], inI128[8], 1); 
     rI[1] = _mm256_insertf128_si256(inIb[0], inI128[9], 1); 
     rI[2] = _mm256_insertf128_si256(inIa[1], inI128[10], 1); 
     rI[3] = _mm256_insertf128_si256(inIb[1], inI128[11], 1); 
     rI[4] = _mm256_insertf128_si256(inIa[2], inI128[12], 1); 
     rI[5] = _mm256_insertf128_si256(inIb[2], inI128[13], 1); 
     rI[6] = _mm256_insertf128_si256(inIa[3], inI128[14], 1); 
     rI[7] = _mm256_insertf128_si256(inIb[3], inI128[15], 1); 

     __m256 rrF[8];
     __m256 *rF = reinterpret_cast<__m256 *>(rI);
     rrF[0] = _mm256_unpacklo_ps(rF[0], rF[4]); 
     rrF[1] = _mm256_unpackhi_ps(rF[0], rF[4]); 
     rrF[2] = _mm256_unpacklo_ps(rF[1], rF[5]); 
     rrF[3] = _mm256_unpackhi_ps(rF[1], rF[5]); 
     rrF[4] = _mm256_unpacklo_ps(rF[2], rF[6]); 
     rrF[5] = _mm256_unpackhi_ps(rF[2], rF[6]); 
     rrF[6] = _mm256_unpacklo_ps(rF[3], rF[7]); 
     rrF[7] = _mm256_unpackhi_ps(rF[3], rF[7]); 

     rF = reinterpret_cast<__m256 *>(in);
     rF[0] = _mm256_unpacklo_ps(rrF[0], rrF[4]); 
     rF[1] = _mm256_unpackhi_ps(rrF[0], rrF[4]); 
     rF[2] = _mm256_unpacklo_ps(rrF[1], rrF[5]); 
     rF[3] = _mm256_unpackhi_ps(rrF[1], rrF[5]); 
     rF[4] = _mm256_unpacklo_ps(rrF[2], rrF[6]); 
     rF[5] = _mm256_unpackhi_ps(rrF[2], rrF[6]); 
     rF[6] = _mm256_unpacklo_ps(rrF[3], rrF[7]); 
     rF[7] = _mm256_unpackhi_ps(rrF[3], rrF[7]); 
  }

  int main()
  {
  #define dwordCount 64
     alignas(32) unsigned long mat[dwordCount];
     for (int i = 0; i < dwordCount; i++) {
        mat[i] = i;
     }
     clock_t t;
     printf ("Transpose8x8Shuff 100,000,000 times\n");
     Transpose8x8Shuff(mat);
     t = clock();
     int i = 100000000;
     do {
        Transpose8x8Shuff(mat);
     } while (--i >= 0);
     t = clock() - t;
     volatile int dummy = mat[2];
     printf ("(%f seconds).\n",((float)t)/CLOCKS_PER_SEC);
     printf ("Transpose8x8Permute 100,000,000 times\n");
     Transpose8x8Permute(mat);
     t = clock();
     i = 100000000;
     do {
        Transpose8x8Permute(mat);
     } while (--i >= 0);
     t = clock() - t;
     volatile int dummy = mat[2];
     printf ("(%f seconds).\n",((float)t)/CLOCKS_PER_SEC);
     printf ("Transpose8x8Insert 100,000,000 times\n");
     Transpose8x8Insert(mat);
     t = clock();
     i = 100000000;
     do {
        Transpose8x8Insert(mat);
     } while (--i >= 0);
     t = clock() - t;
     volatile int dummy = mat[2];
     printf ("(%f seconds).\n",((float)t)/CLOCKS_PER_SEC);
     char c = getchar(); 
     return 0;
  }

这是基准测试延迟，而不是吞吐量（因为一个转置的输出是下一个转置的输入），但很可能会在洗牌吞吐量上产生瓶颈。

Results on Skylake i7-6700k @ 3.9GHz for a modified version of the above code (see it on the Godbolt compiler explorer), fixing the following bugs:

在开始clock()之前，在定时区域之外使用printf
volatile dummy = in[2]在最后，以防止所有的转置被优化掉（否则gcc会这样做）。
C++11具有可移植性，不需要MSVC（使用alignas(32)代替__declspec，不要包含stdafx.h。）
删除了sleep，以使CPU在测试之间不会降到空闲速度。

我没有修复__m256i*/__m256*的不必要混合，也没有检查是否导致gcc或clang生成更差的代码。我也没有使用std::chrono高分辨率时钟，因为clock()对于这么多重复足够准确。

Arch Linux上的g++7.3 -O3 -march=native：Z Boson的版本是最快的

Transpose8x8Shuff 100,000,000 times
(0.637479 seconds).
Transpose8x8Permute 100,000,000 times
(0.792658 seconds).
Transpose8x8Insert 100,000,000 times
(0.846590 seconds).

clang++ 5.0.1 -O3 -march=native: 在优化级别为-O3和本地架构的情况下，8x8Permute被优化为比gcc更快的代码，但8x8Insert的性能却变得非常差。

Transpose8x8Shuff 100,000,000 times
(0.642185 seconds).
Transpose8x8Permute 100,000,000 times
(0.622157 seconds).
Transpose8x8Insert 100,000,000 times
(2.149958 seconds).

从源代码生成的汇编指令与内置函数不完全匹配：尤其是clang具有一个优化Shuffle的优化器，它确实像对整数的标量代码一样编译Shuffle。 Transpose8x8Insert 不应该比那慢太多，所以clang肯定选错了。

- Andrew Hallendorff

C内置函数和变量并不是指令或寄存器。加载和存储需要额外的指令。(AVX大多数情况下避免了额外的寄存器复制指令。)为什么要强制编译器发出存储指令来更新in[]，而不是使用局部变量呢？相同数量的存储指令，但在编译x86-64时，它们都可以适应寄存器，并且存储可以被优化掉。(此外，使用_mm256_permute2f128_ps而不是si256可以避免将指针转换为两种不同类型。) - Peter Cordes

没有使用__restrict，编译器不知道输入和输出是否重叠，因此它实际上会执行所有存储操作，包括对outI[]的存储。修复这个问题后，您的版本与Z Boson的版本（适应于读取和写入in/out数组而不是就地修改引用）在gcc -m32 -O3 -march=haswell下具有大致相同的指令计数：https://godbolt.org/g/bEwS4r。gcc通过32字节对齐堆栈，并在两个版本中将临时变量溢出到堆栈中。因此，您使用`in[]`来避免额外的溢出的想法行不通，只会导致汇编代码中执行更多的工作。 - Peter Cordes

gcc会将一些 _mm256_permute2f128_si256 优化为128位加载 + vinserti128，就像Z Boson在第二个版本中手动执行的那样。（是的，它使用了整数版本，因为您使用了 si256 而不是 ps）。 - Peter Cordes

我已经按照你的提示进行了一些编辑。不确定volatile dummy = in[2]应该放在哪里。清理一些注释并将其聚焦到有用的东西会很好。也许可以倾向于预处理数据，何时进行以及对某些编译器有什么影响。我还有更多问题/建议，但不想堵塞评论区。（顺便说一句：我有一台旧推土机，我可能会挖出来构建一个Linux盒子，看看它对这段代码的反应如何） - Andrew Hallendorff

随意根据您的需求修改我的代码。我更喜欢将我的变量声明为数组。谢谢。 - Andrew Hallendorff

显示剩余5条评论

2

针对之前的答案，这种情况下使用shuffleps就有些过度了，因为我们可以通过unpacklo/unpackhi来得到结果。Shuffle和unpack指令具有相同的延迟/吞吐量，但是Shuffle在机器代码操作中生成一个额外的字节（即Shuffle需要5个字节，而Unpack只需要4个字节）。

在某个时候，我们需要在通道之间进行8次置换。这是一种较慢的操作（延迟为3个周期），因此如果可能的话，我们希望尽早启动这些操作。假设transpose8f方法被内联（应该会被内联！），那么a->h参数所需的任何加载都应该融合到unpack指令中。

你可能会面临的唯一小问题是，因为你在这里使用了超过8个寄存器，所以你可能会溢出到YMM9及以上。这可能会导致VEX2操作生成为VEX3，每个操作会增加一个字节。

因此，稍微调整一下，你最终将得到以下结果：

typedef __m256 f256;
#define unpacklo8f _mm256_unpacklo_ps
#define unpackhi8f _mm256_unpackhi_ps

template<uint8_t X, uint8_t Y>
inline f256 permute128f(const f256 a, const f256 b)
{
  return _mm256_permute2f128_ps(a, b, X | (Y << 4)); 
}

inline void transpose8f(
  const f256 a, const f256 b, const f256 c, const f256 d, 
  const f256 e, const f256 f, const f256 g, const f256 h, 
  f256& s, f256& t, f256& u, f256& v,
  f256& x, f256& y, f256& z, f256& w)
{
  const f256 t00 = unpacklo8f(a, c);
  const f256 t02 = unpacklo8f(b, d);
  const f256 t20 = unpacklo8f(e, g);
  const f256 t22 = unpacklo8f(f, h);

  const f256 t10 = unpacklo8f(t00, t02);
  const f256 t30 = unpacklo8f(t20, t22);
  const f256 t11 = unpackhi8f(t00, t02);
  s = permute128f<0, 2>(t10, t30);
  const f256 t31 = unpackhi8f(t20, t22);
  x = permute128f<1, 3>(t10, t30);
  const f256 t01 = unpackhi8f(a, c);
  t = permute128f<0, 2>(t11, t31);
  const f256 t21 = unpackhi8f(e, g);
  y = permute128f<1, 3>(t11, t31);
  const f256 t03 = unpackhi8f(b, d);
  const f256 t23 = unpackhi8f(f, h);

  const f256 t12 = unpacklo8f(t01, t03);
  const f256 t13 = unpackhi8f(t01, t03);
  const f256 t32 = unpacklo8f(t21, t23);
  const f256 t33 = unpackhi8f(t21, t23);

  u = permute128f<0, 2>(t12, t32);
  z = permute128f<1, 3>(t12, t32);
  v = permute128f<0, 2>(t13, t33);
  w = permute128f<1, 3>(t13, t33);
}

您不会有所改进 （您可以先进行128位排列，然后进行解压缩，但它们最终将是相同的）。

- robthebloke

0

这是我的解决方案，指令更少，性能非常好，大约快了8倍。我已经在Fedora中使用ICC、GCC和Clang进行了测试。

#include <stdio.h>
#include <x86intrin.h>
#define MAX1    128
#define MAX2 MAX1

float __attribute__(( aligned(32))) a_tra[MAX2][MAX1], __attribute__(( aligned(32))) a[MAX1][MAX2] ;
int main()
{

    int i,j;//, ii=0,jj=0;
    // variables for vector section
    int vindexm [8]={0, MAX1, MAX1*2, MAX1*3, MAX1*4, MAX1*5, MAX1*6, MAX1*7 };
    __m256i vindex = _mm256_load_si256((__m256i *) &vindexm[0]);
    __m256 vec1, vec2, vec3, vec4, vec5, vec6, vec7, vec8;

        for(i=0; i<MAX1;  i+=8){            
            for(j=0; j<MAX2;  j+=8){
                //loading from columns
                vec1 = _mm256_i32gather_ps (&a[i][j+0],vindex,4);
                vec2 = _mm256_i32gather_ps (&a[i][j+1],vindex,4);
                vec3 = _mm256_i32gather_ps (&a[i][j+2],vindex,4);
                vec4 = _mm256_i32gather_ps (&a[i][j+3],vindex,4);
                vec5 = _mm256_i32gather_ps (&a[i][j+4],vindex,4);
                vec6 = _mm256_i32gather_ps (&a[i][j+5],vindex,4);
                vec7 = _mm256_i32gather_ps (&a[i][j+6],vindex,4);
                vec8 = _mm256_i32gather_ps (&a[i][j+7],vindex,4);

                //storing to the rows
                _mm256_store_ps(&a_tra[j+0][i], vec1);
                _mm256_store_ps(&a_tra[j+1][i], vec2);
                _mm256_store_ps(&a_tra[j+2][i], vec3);
                _mm256_store_ps(&a_tra[j+3][i], vec4);
                _mm256_store_ps(&a_tra[j+4][i], vec5);
                _mm256_store_ps(&a_tra[j+5][i], vec6);
                _mm256_store_ps(&a_tra[j+6][i], vec7);
                _mm256_store_ps(&a_tra[j+7][i], vec8);  
            }
        }
    return 0;
}

- Amiri

等等，你只是在将其与标量进行比较，而不是与这里的其他答案进行比较？ - Peter Cordes

@PeterCordes，我在一年前测试过这些答案，在某些情况下其他指令可以获得更好的性能，但例如9x vs 8x或8x vs 7.5x。 - Amiri

即使在Skylake中，如果您将gather指令与其他计算混合使用，性能也会下降。需要引用。如果您可以并行执行gather和其他有用的工作，则可以让更多的执行单元保持繁忙状态。与其他答案相比，“vpgatherdd”每个指令为4个uops，而其他答案中每个洗牌操作为1个uop。因此，uop缓存占用可能更大，即使L1i缓存占用可能更小。 - Peter Cordes

是的，我注意到如果我们以这种方式使用gather指令，我们将失去性能。它的加速（实际上是减速）为0.5倍。当我们在乘法之前对其进行转置时，我们分离了计算部分并重新组织了矩阵，使得内部循环中不需要gather或其他洗牌指令。 - Amiri

哦，你是在内部循环中转置了O(n^3)次，而不是保存/重复使用一个列向量的转置结果来处理所有行向量吗？是的，当然这样会更慢，但不是因为混合收集+计算，而是因为你进行了n倍的转置！如果您的大小大于8x8，则重复循环矩阵时的空间局部性更差。这就像缓存块一样，只不过您将转置结果保存到向量寄存器中，而不是保存到内存中。 - Peter Cordes

显示剩余15条评论

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Z boson · Accepted Answer

我已经回答了这个问题使用SSE、AVX和OpenMP进行快速内存转置。

让我重申一下使用AVX转置8x8浮点矩阵的解决方案。如果与使用4x4块和_MM_TRANSPOSE4_PS相比，这是否更快？我在一个较大的矩阵转置内核中使用它，该内核受限于内存，因此可能不是公平的测试。

inline void transpose8_ps(__m256 &row0, __m256 &row1, __m256 &row2, __m256 &row3, __m256 &row4, __m256 &row5, __m256 &row6, __m256 &row7) {
__m256 __t0, __t1, __t2, __t3, __t4, __t5, __t6, __t7;
__m256 __tt0, __tt1, __tt2, __tt3, __tt4, __tt5, __tt6, __tt7;
__t0 = _mm256_unpacklo_ps(row0, row1);
__t1 = _mm256_unpackhi_ps(row0, row1);
__t2 = _mm256_unpacklo_ps(row2, row3);
__t3 = _mm256_unpackhi_ps(row2, row3);
__t4 = _mm256_unpacklo_ps(row4, row5);
__t5 = _mm256_unpackhi_ps(row4, row5);
__t6 = _mm256_unpacklo_ps(row6, row7);
__t7 = _mm256_unpackhi_ps(row6, row7);
__tt0 = _mm256_shuffle_ps(__t0,__t2,_MM_SHUFFLE(1,0,1,0));
__tt1 = _mm256_shuffle_ps(__t0,__t2,_MM_SHUFFLE(3,2,3,2));
__tt2 = _mm256_shuffle_ps(__t1,__t3,_MM_SHUFFLE(1,0,1,0));
__tt3 = _mm256_shuffle_ps(__t1,__t3,_MM_SHUFFLE(3,2,3,2));
__tt4 = _mm256_shuffle_ps(__t4,__t6,_MM_SHUFFLE(1,0,1,0));
__tt5 = _mm256_shuffle_ps(__t4,__t6,_MM_SHUFFLE(3,2,3,2));
__tt6 = _mm256_shuffle_ps(__t5,__t7,_MM_SHUFFLE(1,0,1,0));
__tt7 = _mm256_shuffle_ps(__t5,__t7,_MM_SHUFFLE(3,2,3,2));
row0 = _mm256_permute2f128_ps(__tt0, __tt4, 0x20);
row1 = _mm256_permute2f128_ps(__tt1, __tt5, 0x20);
row2 = _mm256_permute2f128_ps(__tt2, __tt6, 0x20);
row3 = _mm256_permute2f128_ps(__tt3, __tt7, 0x20);
row4 = _mm256_permute2f128_ps(__tt0, __tt4, 0x31);
row5 = _mm256_permute2f128_ps(__tt1, __tt5, 0x31);
row6 = _mm256_permute2f128_ps(__tt2, __tt6, 0x31);
row7 = _mm256_permute2f128_ps(__tt3, __tt7, 0x31);
}

根据这个评论，我学到了更高效的8x8转置方法。请看Intel优化手册中11.11节下的范例11-19和11-20。范例11-19使用相同数量的指令，但通过使用同时送往端口0的混合器来减少对端口5的压力。我可能会在某个时候使用固有函数来实现它，但目前还没有这个需求。

我仔细研究了我之前提到的英特尔手册中的示例11-19和11-20。结果发现示例11-19使用了比必要多4个的Shuffle操作。它有8个Unpack、12个Shuffle和8个128位Permute。我的方法使用了4个更少的Shuffle。他们用Blend替换了其中的8个Shuffle。因此需要4个Shuffle和8个Blend。我怀疑这不如我的方法只使用8个Shuffle好。

然而，如果您需要从内存中加载矩阵，则示例11-20是一种改进方法。它使用8个Unpack、8个Insert、8个Shuffle、8个128位Load和8个Store。128位Load减轻了端口压力。我已经使用Intrinsics实现了这个方法。

//Example 11-20. 8x8 Matrix Transpose Using VINSERTF128 loads
void tran(float* mat, float* matT) {
  __m256  r0, r1, r2, r3, r4, r5, r6, r7;
  __m256  t0, t1, t2, t3, t4, t5, t6, t7;

  r0 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[0*8+0])), _mm_load_ps(&mat[4*8+0]), 1);
  r1 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[1*8+0])), _mm_load_ps(&mat[5*8+0]), 1);
  r2 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[2*8+0])), _mm_load_ps(&mat[6*8+0]), 1);
  r3 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[3*8+0])), _mm_load_ps(&mat[7*8+0]), 1);
  r4 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[0*8+4])), _mm_load_ps(&mat[4*8+4]), 1);
  r5 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[1*8+4])), _mm_load_ps(&mat[5*8+4]), 1);
  r6 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[2*8+4])), _mm_load_ps(&mat[6*8+4]), 1);
  r7 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[3*8+4])), _mm_load_ps(&mat[7*8+4]), 1);

  t0 = _mm256_unpacklo_ps(r0,r1);
  t1 = _mm256_unpackhi_ps(r0,r1);
  t2 = _mm256_unpacklo_ps(r2,r3);
  t3 = _mm256_unpackhi_ps(r2,r3);
  t4 = _mm256_unpacklo_ps(r4,r5);
  t5 = _mm256_unpackhi_ps(r4,r5);
  t6 = _mm256_unpacklo_ps(r6,r7);
  t7 = _mm256_unpackhi_ps(r6,r7);

  r0 = _mm256_shuffle_ps(t0,t2, 0x44);
  r1 = _mm256_shuffle_ps(t0,t2, 0xEE);
  r2 = _mm256_shuffle_ps(t1,t3, 0x44);
  r3 = _mm256_shuffle_ps(t1,t3, 0xEE);
  r4 = _mm256_shuffle_ps(t4,t6, 0x44);
  r5 = _mm256_shuffle_ps(t4,t6, 0xEE);
  r6 = _mm256_shuffle_ps(t5,t7, 0x44);
  r7 = _mm256_shuffle_ps(t5,t7, 0xEE);

  _mm256_store_ps(&matT[0*8], r0);
  _mm256_store_ps(&matT[1*8], r1);
  _mm256_store_ps(&matT[2*8], r2);
  _mm256_store_ps(&matT[3*8], r3);
  _mm256_store_ps(&matT[4*8], r4);
  _mm256_store_ps(&matT[5*8], r5);
  _mm256_store_ps(&matT[6*8], r6);
  _mm256_store_ps(&matT[7*8], r7);
}

我看了一下例子11-19。据我所知，基本思路是用两个shuffle指令（shufps）代替一个shuffle和两个blends。例如：

r0 = _mm256_shuffle_ps(t0,t2, 0x44);
r1 = _mm256_shuffle_ps(t0,t2, 0xEE);

可以替换为

v = _mm256_shuffle_ps(t0,t2, 0x4E);
r0 = _mm256_blend_ps(t0, v, 0xCC);
r1 = _mm256_blend_ps(t2, v, 0x33);

这就解释了为什么我的原始代码使用了8次洗牌，而Example 11-19只使用了4次洗牌和8次混合。

混合对吞吐量有好处，因为洗牌只能到达一个端口（在洗牌端口上创建瓶颈），但混合可以在多个端口上运行，因此不会竞争。但是哪种更好：8次洗牌还是4次洗牌和8次混合？

这必须进行测试，并且可能取决于周围的代码。如果您在循环中大部分需要端口5之外的其他uop时总体uop吞吐量受到瓶颈限制，您可能会选择纯洗牌版本。理想情况下，您应该在存储转置数据之前对其进行一些计算，因为它已经在寄存器中。请参见https://agner.org/optimize/以及the x86 tag wiki中的其他性能链接。

但是，我并没有看到用混合替换解包指令的方法。

这是一个完整的代码示例，结合了示例11-19将2个洗牌操作转换为1个洗牌操作和两个混合操作以及示例11-20使用vinsertf128加载的方法（在英特尔Haswell/Skylake CPU上是2个微操作：一个ALU用于任何端口，一个内存。不幸的是，它们不能微融合。vinsertf128带有所有寄存器操作数是1个洗牌端口的微操作，在英特尔上很好，因此编译器将加载折叠到vinsertf128的内存操作数中。）这样做的好处是只需要让源数据16字节对齐即可达到最大性能，避免任何缓存行分裂。

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

/*
void tran(float* mat, float* matT) {
  __m256  r0, r1, r2, r3, r4, r5, r6, r7;
  __m256  t0, t1, t2, t3, t4, t5, t6, t7;

  r0 = _mm256_load_ps(&mat[0*8]);
  r1 = _mm256_load_ps(&mat[1*8]);
  r2 = _mm256_load_ps(&mat[2*8]);
  r3 = _mm256_load_ps(&mat[3*8]);
  r4 = _mm256_load_ps(&mat[4*8]);
  r5 = _mm256_load_ps(&mat[5*8]);
  r6 = _mm256_load_ps(&mat[6*8]);
  r7 = _mm256_load_ps(&mat[7*8]);

  t0 = _mm256_unpacklo_ps(r0, r1);
  t1 = _mm256_unpackhi_ps(r0, r1);
  t2 = _mm256_unpacklo_ps(r2, r3);
  t3 = _mm256_unpackhi_ps(r2, r3);
  t4 = _mm256_unpacklo_ps(r4, r5);
  t5 = _mm256_unpackhi_ps(r4, r5);
  t6 = _mm256_unpacklo_ps(r6, r7);
  t7 = _mm256_unpackhi_ps(r6, r7);

  r0 = _mm256_shuffle_ps(t0,t2,_MM_SHUFFLE(1,0,1,0));  
  r1 = _mm256_shuffle_ps(t0,t2,_MM_SHUFFLE(3,2,3,2));
  r2 = _mm256_shuffle_ps(t1,t3,_MM_SHUFFLE(1,0,1,0));
  r3 = _mm256_shuffle_ps(t1,t3,_MM_SHUFFLE(3,2,3,2));
  r4 = _mm256_shuffle_ps(t4,t6,_MM_SHUFFLE(1,0,1,0));
  r5 = _mm256_shuffle_ps(t4,t6,_MM_SHUFFLE(3,2,3,2));
  r6 = _mm256_shuffle_ps(t5,t7,_MM_SHUFFLE(1,0,1,0));
  r7 = _mm256_shuffle_ps(t5,t7,_MM_SHUFFLE(3,2,3,2));

  t0 = _mm256_permute2f128_ps(r0, r4, 0x20);
  t1 = _mm256_permute2f128_ps(r1, r5, 0x20);
  t2 = _mm256_permute2f128_ps(r2, r6, 0x20);
  t3 = _mm256_permute2f128_ps(r3, r7, 0x20);
  t4 = _mm256_permute2f128_ps(r0, r4, 0x31);
  t5 = _mm256_permute2f128_ps(r1, r5, 0x31);
  t6 = _mm256_permute2f128_ps(r2, r6, 0x31);
  t7 = _mm256_permute2f128_ps(r3, r7, 0x31);

  _mm256_store_ps(&matT[0*8], t0);
  _mm256_store_ps(&matT[1*8], t1);
  _mm256_store_ps(&matT[2*8], t2);
  _mm256_store_ps(&matT[3*8], t3);
  _mm256_store_ps(&matT[4*8], t4);
  _mm256_store_ps(&matT[5*8], t5);
  _mm256_store_ps(&matT[6*8], t6);
  _mm256_store_ps(&matT[7*8], t7);
}
*/

void tran(float* mat, float* matT) {
  __m256  r0, r1, r2, r3, r4, r5, r6, r7;
  __m256  t0, t1, t2, t3, t4, t5, t6, t7;

  r0 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[0*8+0])), _mm_load_ps(&mat[4*8+0]), 1);
  r1 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[1*8+0])), _mm_load_ps(&mat[5*8+0]), 1);
  r2 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[2*8+0])), _mm_load_ps(&mat[6*8+0]), 1);
  r3 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[3*8+0])), _mm_load_ps(&mat[7*8+0]), 1);
  r4 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[0*8+4])), _mm_load_ps(&mat[4*8+4]), 1);
  r5 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[1*8+4])), _mm_load_ps(&mat[5*8+4]), 1);
  r6 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[2*8+4])), _mm_load_ps(&mat[6*8+4]), 1);
  r7 = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm_load_ps(&mat[3*8+4])), _mm_load_ps(&mat[7*8+4]), 1);

  t0 = _mm256_unpacklo_ps(r0,r1);
  t1 = _mm256_unpackhi_ps(r0,r1);
  t2 = _mm256_unpacklo_ps(r2,r3);
  t3 = _mm256_unpackhi_ps(r2,r3);
  t4 = _mm256_unpacklo_ps(r4,r5);
  t5 = _mm256_unpackhi_ps(r4,r5);
  t6 = _mm256_unpacklo_ps(r6,r7);
  t7 = _mm256_unpackhi_ps(r6,r7);

  __m256 v;

  //r0 = _mm256_shuffle_ps(t0,t2, 0x44);
  //r1 = _mm256_shuffle_ps(t0,t2, 0xEE);  
  v = _mm256_shuffle_ps(t0,t2, 0x4E);
  r0 = _mm256_blend_ps(t0, v, 0xCC);
  r1 = _mm256_blend_ps(t2, v, 0x33);

  //r2 = _mm256_shuffle_ps(t1,t3, 0x44);
  //r3 = _mm256_shuffle_ps(t1,t3, 0xEE);
  v = _mm256_shuffle_ps(t1,t3, 0x4E);
  r2 = _mm256_blend_ps(t1, v, 0xCC);
  r3 = _mm256_blend_ps(t3, v, 0x33);

  //r4 = _mm256_shuffle_ps(t4,t6, 0x44);
  //r5 = _mm256_shuffle_ps(t4,t6, 0xEE);
  v = _mm256_shuffle_ps(t4,t6, 0x4E);
  r4 = _mm256_blend_ps(t4, v, 0xCC);
  r5 = _mm256_blend_ps(t6, v, 0x33);

  //r6 = _mm256_shuffle_ps(t5,t7, 0x44);
  //r7 = _mm256_shuffle_ps(t5,t7, 0xEE);
  v = _mm256_shuffle_ps(t5,t7, 0x4E);
  r6 = _mm256_blend_ps(t5, v, 0xCC);
  r7 = _mm256_blend_ps(t7, v, 0x33);

  _mm256_store_ps(&matT[0*8], r0);
  _mm256_store_ps(&matT[1*8], r1);
  _mm256_store_ps(&matT[2*8], r2);
  _mm256_store_ps(&matT[3*8], r3);
  _mm256_store_ps(&matT[4*8], r4);
  _mm256_store_ps(&matT[5*8], r5);
  _mm256_store_ps(&matT[6*8], r6);
  _mm256_store_ps(&matT[7*8], r7);
}


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

void print_mat(float *mat) {
    int i,j;
    for(i=0; i<8; i++) {
      for(j=0; j<8; j++) printf("%2.0f ", mat[i*8+j]);
      puts("");
    }
    puts("");
}

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

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

    tran(mat,matT);
    //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);
}