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

  • Z时代
  • 2024-01-10
  • 分类：问答

我有一个需要转置的矩阵（相对较大）。例如，假设我的矩阵是

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

最快的方法是什么？

回答：

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

首先，让我列出我用于移调的功能之一（ ）

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/zh-cn/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（

）。该代码适用于任何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);
                }
            }
        }
    }
}

以上是 用C ++转置矩阵最快的方法是什么？ 的全部内容， 来源链接： utcz.com/qa/417121.html

回到顶部