当行大小大于矢量宽度时，SIMD转置

Question

您可以找到很多good answers，用于转置随SIMD指令集的natural大小而下降的矩阵，尤其是其中一行的大小不超过矢量宽度的情况下。例如，在SSE中使用4x4 float转置，或者在AVX / AVX2中使用4x4 double或8x8 float转置（对于AVX-512，再次将其翻倍）。

但是，当矩阵大于该值时，有哪些选择呢？例如，使用AVX2的16x16 float矩阵？可以完全使用SIMD改组来加快处理速度，还是仅采用聚集+顺序写入？

Answer 1

如果所有矩阵尺寸都是数据包大小的倍数，则可以按块进行操作并根据需要交换块。使用SSE2的4x4双矩阵示例：

// transpose vectors i0 and i1 and store the result to addresses r0 and r1
void transpose2x2(double *r0, double* r1, __m128d i0, __m128d i1)
{
    __m128d t0 = _mm_unpacklo_pd(i0,i1);
    __m128d t1 = _mm_unpackhi_pd(i0,i1);
    _mm_storeu_pd(r0, t0);
    _mm_storeu_pd(r1, t1);
}


void transpose(double mat[4][4])
{
    // transpose [00]-block in-place
    transpose2x2(mat[0]+0, mat[1]+0,_mm_loadu_pd(mat[0]+0),_mm_loadu_pd(mat[1]+0));

    // load [20]-block
    __m128d t20 = _mm_loadu_pd(mat[2]+0), t30 = _mm_loadu_pd(mat[3]+0);
    // transpose [02]-block and store it to [20] position
    transpose2x2(mat[2]+0,mat[3]+0, _mm_loadu_pd(mat[0]+2),_mm_loadu_pd(mat[1]+2));
    // transpose temp-block and store it to [02] position
    transpose2x2(mat[0]+2,mat[1]+2, t20, t30);

    // transpose [22]-block in-place
    transpose2x2(mat[2]+2, mat[3]+2,_mm_loadu_pd(mat[2]+2),_mm_loadu_pd(mat[3]+2));
}

这应该相对容易地扩展到其他平方矩阵，其他标量类型和其他体系结构。不是数据包大小倍数的矩阵可能会更复杂（如果它们足够大，那么进行向量化的大部分工作，而只是手动完成最后几行/几列，可能是值得的）。

对于某些尺寸，例如3x4或3x8矩阵有特殊的算法[1]-如果您有1003x1003矩阵，则可以在最后一行/列中利用它（并且可能还有其他奇数大小的算法）。

[您也可以花一些力气为矩形矩阵编写代码（必须考虑如何避免一次只缓存一个以上的块，但是有可能）。

Godbolt演示：https://godbolt.org/z/tVk_Bc

[1] https://software.intel.com/en-us/articles/3d-vector-normalization-using-256-bit-intel-advanced-vector-extensions-intel-avx

Answer 2

也许可以将fortran TRANSPOSE固有地与ISO_C_BINDING一起使用，并将其与C链接为子例程或函数调用。