用c++最快的转置矩阵的方法是什么?

这个问题问得好。有很多原因你想要在内存中实际转置矩阵，而不仅仅是交换坐标，例如在矩阵乘法和高斯涂抹中。

首先让我列出我用于转置的函数之一(编辑:请参阅我的答案结尾，我发现了一个更快的解决方案)

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);
                }
            }
        }
    }
}

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

关于在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条指令。优点是它降低了端口压力，这在某些情况下可能有好处。目前所有的洗牌和解包只能去一个特定的端口，而混合可以去两个不同的端口之一。

我尝试使用8洗牌像MSVC和转换成4洗牌+ 8混合，但它没有工作。我仍然需要使用4个解包。

我对8x8浮点数的转置使用了相同的技术(参见答案的末尾)。https://stackoverflow.com/a/25627536/2542702。在这个答案中，我仍然需要使用8次解包，但我设法将8次洗牌转换为4次洗牌和8次混合。

对于32位整数，没有像shufps这样的东西(除了AVX512的128位shuffle)，所以它只能通过解包来实现，我认为解包不能(有效地)转换为混合。对于AVX512, vshufi32x4的行为有效地像shufps，除了128位的4个整数通道而不是32位的浮点数，所以在某些情况下，vshufi32x4也可能采用相同的技术。在骑士登陆中，洗牌比混合慢四倍。

如果事先知道数组的大小，则可以使用联合来帮助我们。像这样——

#include <bits/stdc++.h>
using namespace std;
union ua{
    int arr[2][3];
    int brr[3][2];
};
int main() {
    union ua uav;
    int karr[2][3] = {{1,2,3},{4,5,6}};
    memcpy(uav.arr,karr,sizeof(karr));
    for (int i=0;i<3;i++)
    {
        for (int j=0;j<2;j++)
            cout<<uav.brr[i][j]<<" ";
        cout<<'n';
    }
    return 0;
}

将每一行视为一列，每一列视为一行。使用j,i代替i,j

演示:http://ideone.com/lvsxKZ

#include <iostream> 
using namespace std;
int main ()
{
    char A [3][3] =
    {
        { 'a', 'b', 'c' },
        { 'd', 'e', 'f' },
        { 'g', 'h', 'i' }
    };
    cout << "A = " << endl << endl;
    // print matrix A
    for (int i=0; i<3; i++)
    {
        for (int j=0; j<3; j++) cout << A[i][j];
        cout << endl;
    }
    cout << endl << "A transpose = " << endl << endl;
    // print A transpose
    for (int i=0; i<3; i++)
    {
        for (int j=0; j<3; j++) cout << A[j][i];
        cout << endl;
    }
    return 0;
}

没有任何开销的转置(类未完成):

class Matrix{
   double *data; //suppose this will point to data
   double _get1(int i, int j){return data[i*M+j];} //used to access normally
   double _get2(int i, int j){return data[j*N+i];} //used when transposed
   public:
   int M, N; //dimensions
   double (*get_p)(int, int); //functor to access elements  
   Matrix(int _M,int _N):M(_M), N(_N){
     //allocate data
     get_p=&Matrix::_get1; // initialised with normal access 
     }
   double get(int i, int j){
     //there should be a way to directly use get_p to call. but i think even this
     //doesnt incur overhead because it is inline and the compiler should be intelligent
     //enough to remove the extra call
     return (this->*get_p)(i,j);
    }
   void transpose(){ //twice transpose gives the original
     if(get_p==&Matrix::get1) get_p=&Matrix::_get2;
     else get_p==&Matrix::_get1; 
     swap(M,N);
     }
}

可以这样使用:

Matrix M(100,200);
double x=M.get(17,45);
M.transpose();
x=M.get(17,45); // = original M(45,17)

当然，我没有在这里麻烦内存管理，这是至关重要的，但不同的主题。

现代线性代数库包括最常见操作的优化版本。其中许多包括动态CPU调度，它在程序执行时为硬件选择最佳实现(不影响可移植性)。

这通常是通过向量扩展内禀函数对函数进行手动优化的更好选择。后者会将你的实现绑定到特定的硬件供应商和模型:如果你决定切换到不同的供应商(例如Power, ARM)或更新的向量扩展(例如AVX512)，你将需要再次重新实现它以获得其中的大部分。

MKL转置包含BLAS扩展函数imatcopy。您也可以在其他实现中找到它，例如OpenBLAS:

#include <mkl.h>
void transpose( float* a, int n, int m ) {
    const char row_major = 'R';
    const char transpose = 'T';
    const float alpha = 1.0f;
    mkl_simatcopy (row_major, transpose, n, m, alpha, a, n, n);
}

对于c++项目，您可以使用Armadillo c++:

#include <armadillo>
void transpose( arma::mat &matrix ) {
    arma::inplace_trans(matrix);
}

intel mkl建议使用就地和非就地置换/复制矩阵。这里是文档的链接。我建议尝试异地实现，因为异地实现比异地实现更快，而且最新版本的mkl文档中包含一些错误。

template <class T>
void transpose( const std::vector< std::vector<T> > & a,
std::vector< std::vector<T> > & b,
int width, int height)
{
    for (int i = 0; i < width; i++)
    {
        for (int j = 0; j < height; j++)
        {
            b[j][i] = a[i][j];
        }
    }
} 

最快的转置是将留在缓存中等待下一个操作(将使用它)的转置。

例如，不要一次全部转置。只转置一个子矩阵。然后将它用于下一个算法中需要转置数据的部分。然后转置下一个子矩阵。然后计算。然后转置另一个子矩阵。重复，直到整个矩阵被转置。通过这种方式，数据在缓存中保持热。

如果您一次在具有2MB缓存的CPU上完全转置128MB矩阵，那么在操作结束时，只有矩阵的最新位在缓存中。然后，您最好从最新的位开始乘矩阵，以使用热的2MB数据。

但是当你把工作分成更小的部分时，比如用子矩阵做乘法，那么你可以简单地做一个像这样的惰性转置:

multiply:
  for all sub_matrices in mat1 row
  for all sub_matrices in mat2 column
    select sub_matrix1
    select sub_matrix2
    if sub_mat2 is not transposed
        transpose sub_mat2
    multiply sub_mat1 and sub_mat2 <---- data in cache!
    accumulate result

优势:

• L1/L2缓存带宽用于下一个操作
• 转置延迟隐藏在下一个操作
的计算之后。
• 工作在小缓存低至64kB，取决于块大小

我认为最快速的方法不应该大于O(n^2)也可以这样使用O(1)空间:
这样做的方法是成对交换，因为当你对一个矩阵进行转置时，你要做的是:M[i][j]=M[j][i]，所以把M[i][j] [j]存储在temp中，然后M[i][j]=M[j][i]，最后一步:M[j][i]=temp。这可以一次完成所以它需要O(n^2)

我的答案是3x3矩阵的转置

 #include<iostream.h>
#include<math.h>

main()
{
int a[3][3];
int b[3];
cout<<"You must give us an array 3x3 and then we will give you Transposed it "<<endl;
for(int i=0;i<3;i++)
{
    for(int j=0;j<3;j++)
{
cout<<"Enter a["<<i<<"]["<<j<<"]: ";
cin>>a[i][j];
}
}
cout<<"Matrix you entered is :"<<endl;
 for (int e = 0 ; e < 3 ; e++ )
{
    for ( int f = 0 ; f < 3 ; f++ )
        cout << a[e][f] << "t";

    cout << endl;
    }
 cout<<"nTransposed of matrix you entered is :"<<endl;
 for (int c = 0 ; c < 3 ; c++ )
{
    for ( int d = 0 ; d < 3 ; d++ )
        cout << a[d][c] << "t";
    cout << endl;
    }
return 0;
}