当行大小大于向量宽度时，进行SIMD矩阵转置

Question

当行大小大于向量宽度时，进行SIMD矩阵转置

3

您可以在这里和这里找到很多关于转置矩阵的好答案，尤其是当一行的大小不超过向量宽度时，自然符合SIMD指令集的大小。例如，在SSE中进行4x4的float转置，在AVX / AVX2中进行4x4的double或8x8的float转置（为AVX-512再加倍）。

然而，如果矩阵大于这个大小呢？例如，使用AVX2的16x16的float矩阵？是否可以使用SIMD洗牌来加速，还是只能使用gather + sequential write的方式？

- BeeOnRope

2个回答

-1

也许可以使用Fortran的TRANSPOSE内置函数与ISO_C_BINDING一起使用，并将其链接为C子例程或函数调用。

在Fortran中，TRANSPOSE非常优化。

有时候混合语言技能是通用的。我甚至将F90与GO链接在一起。

- Holmz

2

好的，但是它最终运行的汇编代码是什么，使得它在任何给定的Fortran实现上都能快速运行？如果我们知道了这一点，我们也可以使用内置函数或手写汇编语言在C中实现它。 - Peter Cordes

2

如果矩阵是方阵且大小是四或八的倍数（在编译时已知），gfortran似乎会将“TRANSPOSE”优化为高效的SIMD代码。对于我查看的其他情况，我只发现它生成了一个微不足道的嵌套循环：https://godbolt.org/z/Uiv1iR -- 另外，我没有成功让它生成一个原地转置而不分配临时空间（但这可能是由于我的有限Fortran技能...） - chtz

在编程中使用易于使用的语言特性并不罕见。混合语言可能不是每个人都喜欢的，但它相当普遍。 - Holmz

对于一个小矩阵，比如4x4的双精度矩阵，调用另一个函数而不是内联几个SIMD指令的开销是相当大的。 - Peter Cordes

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

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

- chtz · Accepted Answer

如果您的矩阵维度都是数据包大小的倍数，那么您可以按块操作并根据需要交换块。以下是使用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