在C++中转置矩阵的最快方法是什么？

Question

我有一个矩阵(相对较大),我需要转置.例如假设我的矩阵是

a b c d e f
g h i j k l
m n o p q r 

我希望结果如下:

a g m
b h n
c I o
d j p
e k q
f l r

最快的方法是什么？

Answer 1

这是一个很好的问题.有很多原因你想要在内存中实际转置矩阵而不是交换坐标,例如在矩阵乘法和高斯拖尾中.

首先让我列出我用于转置的一个功能(编辑:请参阅我的答案的结尾,我找到了一个更快的解决方案)

void transpose(float *src, float *dst, const int N, const int M) {
    #pragma omp parallel for
    for(int n = 0; n<N*M; n++) {
        int i = n/N;
        int j = n%N;
        dst[n] = src[M*j + i];
    }
}

现在让我们看看为什么转置是有用的.考虑矩阵乘法C = A*B. 我们可以这样做.

for(int i=0; i<N; i++) {
    for(int j=0; j<K; j++) {
        float tmp = 0;
        for(int l=0; l<M; l++) {
            tmp += A[M*i+l]*B[K*l+j];
        }
        C[K*i + j] = tmp;
    }
}

但是,这样会有很多缓存未命中.更快的解决方案是首先采用B的转置

transpose(B);
for(int i=0; i<N; i++) {
    for(int j=0; j<K; j++) {
        float tmp = 0;
        for(int l=0; l<M; l++) {
            tmp += A[M*i+l]*B[K*j+l];
        }
        C[K*i + j] = tmp;
    }
}
transpose(B);

矩阵乘法是O(n ^ 3)并且转置是O(n ^ 2),因此对转置应该对计算时间(对于大n)具有可忽略的影响.在矩阵乘法循环中,平铺甚至比采用转置更有效,但这要复杂得多.

我希望我知道更快的方式进行转置(编辑:我发现了一个更快的解决方案,请参阅我的答案的结尾).当Haswell/AVX2在几周内发布时,它将具有聚集功能.我不知道在这种情况下这是否有用,但我可以用图像收集一列并写出一行.也许它会使转置变得不必要.

对于高斯涂抹,你要做的是水平涂抹,然后垂直涂抹.但垂直涂抹有缓存问题,所以你要做的就是

Smear image horizontally
transpose output 
Smear output horizontally
transpose output

以下是英特尔的一篇论文,解释说 http://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions

最后,我在矩阵乘法(以及高斯拖尾)中实际做的并不是完全采用转置,而是以某个矢量大小的宽度(例如,SSE/AVX为4或8)进行转置.这是我使用的功能

void reorder_matrix(const float* A, float* B, const int N, const int M, const int vec_size) {
    #pragma omp parallel for
    for(int n=0; n<M*N; n++) {
        int k = vec_size*(n/N/vec_size);
        int i = (n/vec_size)%N;
        int j = n%vec_size;
        B[n] = A[M*i + k + j];
    }
}

编辑:

我尝试了几个函数来找到大矩阵的最快转置.最后,最快的结果是使用循环阻塞block_size=16(编辑:我发现使用SSE和循环阻塞的更快的解决方案 - 见下文).此代码适用于任何NxM矩阵(即矩阵不必是正方形).

inline void transpose_scalar_block(float *A, float *B, const int lda, const int ldb, const int block_size) {
    #pragma omp parallel for
    for(int i=0; i<block_size; i++) {
        for(int j=0; j<block_size; j++) {
            B[j*ldb + i] = A[i*lda +j];
        }
    }
}

inline void transpose_block(float *A, float *B, const int n, const int m, const int lda, const int ldb, const int block_size) {
    #pragma omp parallel for
    for(int i=0; i<n; i+=block_size) {
        for(int j=0; j<m; j+=block_size) {
            transpose_scalar_block(&A[i*lda +j], &B[j*ldb + i], lda, ldb, block_size);
        }
    }
}

值lda和ldb矩阵的宽度.这些需要是块大小的倍数.为了找到值并为例如3000x1001矩阵分配内存,我做了类似的事情

#define ROUND_UP(x, s) (((x)+((s)-1)) & -(s))
const int n = 3000;
const int m = 1001;
int lda = ROUND_UP(m, 16);
int ldb = ROUND_UP(n, 16);

float *A = (float*)_mm_malloc(sizeof(float)*lda*ldb, 64);
float *B = (float*)_mm_malloc(sizeof(float)*lda*ldb, 64);

对于3000x1001,这将返回 ldb = 3008和 lda = 1008

编辑:

我找到了使用SSE内在函数的更快的解决方案:

inline void transpose4x4_SSE(float *A, float *B, const int lda, const int ldb) {
    __m128 row1 = _mm_load_ps(&A[0*lda]);
    __m128 row2 = _mm_load_ps(&A[1*lda]);
    __m128 row3 = _mm_load_ps(&A[2*lda]);
    __m128 row4 = _mm_load_ps(&A[3*lda]);
     _MM_TRANSPOSE4_PS(row1, row2, row3, row4);
     _mm_store_ps(&B[0*ldb], row1);
     _mm_store_ps(&B[1*ldb], row2);
     _mm_store_ps(&B[2*ldb], row3);
     _mm_store_ps(&B[3*ldb], row4);
}

inline void transpose_block_SSE4x4(float *A, float *B, const int n, const int m, const int lda, const int ldb ,const int block_size) {
    #pragma omp parallel for
    for(int i=0; i<n; i+=block_size) {
        for(int j=0; j<m; j+=block_size) {
            int max_i2 = i+block_size < n ? i + block_size : n;
            int max_j2 = j+block_size < m ? j + block_size : m;
            for(int i2=i; i2<max_i2; i2+=4) {
                for(int j2=j; j2<max_j2; j2+=4) {
                    transpose4x4_SSE(&A[i2*lda +j2], &B[j2*ldb + i2], lda, ldb);
                }
            }
        }
    }
}

Answer 2

这将取决于您的应用程序,但一般来说,转置矩阵的最快方法是在查找时反转坐标,然后您不必实际移动任何数据.

Answer 3

有关使用x86硬件转置4x4正方形浮点数（稍后我将讨论32位整数）矩阵的一些详细信息。从此处开始以便转置较大的正方形矩阵（例如8x8或16x16）会很有帮助。

_MM_TRANSPOSE4_PS(r0, r1, r2, r3)由不同的编译器以不同的方式实现。GCC和ICC（我尚未检查Clang）的使用，unpcklps, unpckhps, unpcklpd, unpckhpd而MSVC仅使用shufps。我们实际上可以像这样将这两种方法结合在一起。

t0 = _mm_unpacklo_ps(r0, r1);
t1 = _mm_unpackhi_ps(r0, r1);
t2 = _mm_unpacklo_ps(r2, r3);
t3 = _mm_unpackhi_ps(r2, r3);

r0 = _mm_shuffle_ps(t0,t2, 0x44);
r1 = _mm_shuffle_ps(t0,t2, 0xEE);
r2 = _mm_shuffle_ps(t1,t3, 0x44);
r3 = _mm_shuffle_ps(t1,t3, 0xEE);

一个有趣的发现是，可以将两种混洗转换为一种混洗和两种混合（SSE4.1）。

t0 = _mm_unpacklo_ps(r0, r1);
t1 = _mm_unpackhi_ps(r0, r1);
t2 = _mm_unpacklo_ps(r2, r3);
t3 = _mm_unpackhi_ps(r2, r3);

v  = _mm_shuffle_ps(t0,t2, 0x4E);
r0 = _mm_blend_ps(t0,v, 0xC);
r1 = _mm_blend_ps(t2,v, 0x3);
v  = _mm_shuffle_ps(t1,t3, 0x4E);
r2 = _mm_blend_ps(t1,v, 0xC);
r3 = _mm_blend_ps(t3,v, 0x3);

这有效地将4个混洗转换为2个混洗和4个混合。这比GCC，ICC和MSVC的实现多使用2条指令。优点是它可以降低端口压力，这在某些情况下可能会有所帮助。当前，所有混洗和拆包只能进入一个特定的港口，而混合酒可以进入两个不同的港口之一。

我尝试使用MSVC等8种混洗并将其转换为4种混洗+ 8种混合，但没有用。我仍然必须使用4个解压缩包。

我对8x8浮点转置使用了相同的技术（请参见该答案的结尾）。 /sf/answers/1793927551/。在该答案中，我仍然必须使用8个解包程序，但我设法将8个混洗转换为4个混洗和8种混合处理。

对于32位整数，没有什么比shufps（AVX512的128位改组除外），因此只能使用解压缩来实现，我认为这些解压缩不能（有效地）转换为混合。使用AVX512时，其vshufi32x4行为类似于shufps4个整数的128位通道而不是32位浮点数，因此vshufi32x4在某些情况下可能也可以使用相同的技术。使用Knights Landing，混洗的速度（吞吐量）比混合慢4倍。