C ++에서 행렬을 전치하는 가장 빠른 방법은 무엇입니까?

C ++에서 행렬을 전치하는 가장 빠른 방법은 무엇입니까?

---

전치해야 할 행렬 (상대적으로 큼)이 있습니다. 예를 들어 내 행렬이

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)이므로 전치를 수행하면 계산 시간에 무시할 수있는 영향을 미칩니다 (large 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를 설명하는 Intel의 문서입니다 .

마지막으로 매트릭스 곱셈 (및 가우스 스미어 링)에서 실제로 수행하는 작업은 정확히 전치하는 것이 아니라 특정 벡터 크기 (예 : 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 개의 명령어를 더 사용합니다. 장점은 일부 상황에서 이점을 가질 수있는 포트 압력을 감소 시킨다는 것입니다. 현재 모든 셔플 및 압축 풀기는 하나의 특정 포트로만 이동할 수 있지만 블렌드는 두 개의 다른 포트 중 하나로 이동할 수 있습니다.

MSVC와 같은 8 개의 셔플을 사용하여 4 개의 셔플 + 8 개의 블렌드로 변환하려고 시도했지만 작동하지 않았습니다. 나는 여전히 4 번의 포장을 풀어야했다.

나는이 같은 기술을 8x8 float transpose에 사용했습니다 (해답의 끝 부분 참조). https://stackoverflow.com/a/25627536/2542702 . 그 대답에서 나는 여전히 8 개의 압축 풀기를 사용해야했지만 8 개의 셔플을 4 개의 셔플과 8 개의 블렌드로 변환하도록 관리했습니다.

32 비트 정수의 경우 shufps(AVX512를 사용한 128 비트 셔플을 제외하고) 같은 것이 없으므로 블렌드로 변환 할 수 없다고 생각하는 언팩으로 만 구현할 수 있습니다 (효율적으로). AVX512 vshufi32x4는 shufps32 비트 부동 소수점 대신 4 개의 정수로 구성된 128 비트 레인을 제외하고 효과적으로 작동 하므로이 동일한 기술이 경우에 따라 가능할 수 있습니다 vshufi32x4. Knights Landing의 경우 셔플은 혼합보다 4 배 더 느립니다 (처리량).

---

template <class T>
void transpose( 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];
        }
    }
} 

---

각 행을 열로, 각 열을 행으로 간주 .. i, j 대신 j, i 사용

데모 : 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)

물론 중요하지만 다른 주제 인 메모리 관리에 대해서는 신경 쓰지 않았습니다.

---

배열의 크기를 미리 알고있는 경우 통합을 사용하여 도움을받을 수 있습니다. 이렇게-

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

---

최신 선형 대수 라이브러리에는 가장 일반적인 연산의 최적화 된 버전이 포함되어 있습니다. 이들 중 다수에는 동적 CPU 디스패치가 포함되어 프로그램 실행시 하드웨어에 대한 최상의 구현을 선택합니다 (이동성 저하없이).

This is commonly a better alternative to performing manual optimization of your functinos via vector extensions intrinsic functions. The latter will tie your implementation to a particular hardware vendor and model: if you decide to swap to a different vendor (e.g. Power, ARM) or to a newer vector extensions (e.g. AVX512), you will need to re-implement it again to get the most of them.

MKL transposition, for example, includes the BLAS extensions function imatcopy. You can find it in other implementations such as OpenBLAS as well:

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

For a C++ project, you can make use of the Armadillo C++:

#include <armadillo>

void transpose( arma::mat &matrix ) {
    arma::inplace_trans(matrix);
}

---

I think that most fast way should not taking higher than O(n^2) also in this way you can use just O(1) space :
the way to do that is to swap in pairs because when you transpose a matrix then what you do is: M[i][j]=M[j][i] , so store M[i][j] in temp, then M[i][j]=M[j][i],and the last step : M[j][i]=temp. this could be done by one pass so it should take O(n^2)

---

my answer is transposed of 3x3 matrix

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

참고URL : https://stackoverflow.com/questions/16737298/what-is-the-fastest-way-to-transpose-a-matrix-in-c