Jaki jest najszybszy sposób transpozycji macierzy w C ++?

Question 1

Mam matrycę (stosunkowo dużą), którą muszę przetransponować. Załóżmy na przykład, że moja macierz to

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

Chcę, aby wynik był następujący:

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

Jaki jest najszybszy sposób na zrobienie tego?

Question 2

To jest dobre pytanie. Istnieje wiele powodów, dla których chciałbyś faktycznie transponować macierz w pamięci, a nie tylko zamienić współrzędne, np. W mnożeniu macierzy i rozmazaniu Gaussa.

Najpierw pozwól mi wymienić jedną z funkcji, których używam do transpozycji ( EDYCJA: zobacz koniec mojej odpowiedzi, gdzie znalazłem znacznie szybsze rozwiązanie )

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

Zobaczmy teraz, dlaczego transpozycja jest przydatna. Rozważ mnożenie macierzy C = A * B. Moglibyśmy to zrobić w ten sposó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;
    }
}

W ten sposób będzie jednak dużo chybień w pamięci podręcznej. Znacznie szybszym rozwiązaniem jest wykonanie transpozycji 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);

Mnożenie macierzy to O (n ^ 3), a transpozycja to O (n ^ 2), więc wykonanie transpozycji powinno mieć znikomy wpływ na czas obliczeń (dla dużych n ). W przypadku pętli mnożenia macierzy jest jeszcze bardziej efektywne niż wykonanie transpozycji, ale jest to znacznie bardziej skomplikowane.

Żałuję, że nie znam szybszego sposobu wykonania transpozycji ( Edycja: znalazłem szybsze rozwiązanie, zobacz koniec mojej odpowiedzi ). Kiedy Haswell / AVX2 wyjdzie za kilka tygodni, będzie miał funkcję zbierającą. Nie wiem, czy to będzie pomocne w tym przypadku, ale mógłbym sobie wyobrazić zbieranie kolumny i pisanie wiersza. Może sprawi, że transpozycja stanie się niepotrzebna.

W przypadku rozmazywania Gaussa rozmazujesz w poziomie, a następnie w pionie. Ale smużenie w pionie ma problem z pamięcią podręczną, więc to, co robisz, jest

Smear image horizontally
transpose output 
Smear output horizontally
transpose output

Oto dokument firmy Intel wyjaśniający, że http://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions

Wreszcie, to, co faktycznie robię w mnożeniu macierzy (i rozmazaniu Gaussa) nie polega na dokładnej transpozycji, ale na transpozycji w szerokościach o określonym rozmiarze wektora (np. 4 lub 8 dla SSE / AVX). Oto funkcja, której używam

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

EDYTOWAĆ:

Wypróbowałem kilka funkcji, aby znaleźć najszybszą transpozycję dla dużych macierzy. Ostatecznie najszybszym rezultatem jest użycie blokowania pętli z block_size=16( Edycja: znalazłem szybsze rozwiązanie wykorzystujące SSE i blokowanie pętli - patrz poniżej ). Ten kod działa dla dowolnej macierzy NxM (tj. Macierz nie musi być kwadratowa).

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

Wartości ldai ldbsą szerokością macierzy. Muszą to być wielokrotności rozmiaru bloku. Aby znaleźć wartości i przydzielić pamięć np. Dla macierzy 3000x1001, robię coś takiego

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

Dla 3000x1001 zwraca ldb = 3008i lda = 1008

Edytować:

Znalazłem jeszcze szybsze rozwiązanie przy użyciu funkcji wewnętrznych 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);
                }
            }
        }
    }
}

Question 3

Będzie to zależeć od aplikacji, ale generalnie najszybszym sposobem transpozycji macierzy byłoby odwrócenie współrzędnych podczas patrzenia w górę, wtedy nie trzeba faktycznie przenosić żadnych danych.

Question 4

Kilka szczegółów na temat transpozycji macierzy typu float kwadratowych 4x4 (omówię później 32-bitowe liczby całkowite) na sprzęcie x86. Warto zacząć tutaj, aby transponować większe macierze kwadratowe, takie jak 8x8 lub 16x16.

_MM_TRANSPOSE4_PS(r0, r1, r2, r3)jest implementowany w różny sposób przez różne kompilatory. GCC i ICC (nie sprawdzałem Clang) używają, unpcklps, unpckhps, unpcklpd, unpckhpdpodczas gdy MSVC używa tylko shufps. Właściwie możemy połączyć te dwa podejścia razem w ten sposób.

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

Interesującą obserwacją jest to, że dwa tasowania można przekształcić w jedno tasowanie i dwa mieszanki (SSE4.1) w ten sposób.

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

To skutecznie przekształciło 4 tasowania w 2 tasowania i 4 mieszanki. Wykorzystuje to 2 instrukcje więcej niż implementacja GCC, ICC i MSVC. Zaletą jest to, że zmniejsza ciśnienie w porcie, co może być korzystne w niektórych okolicznościach. Obecnie wszystkie tasowania i rozpakowywania mogą trafiać tylko do jednego konkretnego portu, podczas gdy mieszanki mogą trafiać do jednego z dwóch różnych portów.

Próbowałem użyć 8 tasowań, takich jak MSVC i przekonwertować to na 4 tasowania + 8 mieszanek, ale to nie zadziałało. Nadal musiałem użyć 4 rozpakowań.

Użyłem tej samej techniki do transpozycji zmiennoprzecinkowej 8x8 (patrz pod koniec odpowiedzi). https://stackoverflow.com/a/25627536/2542702 . W tej odpowiedzi nadal musiałem użyć 8 unpaków, ale udało mi się zamienić 8 tasowań na 4 tasowania i 8 mieszanek.

Dla 32-bitowych liczb całkowitych nie ma nic podobnego shufps(z wyjątkiem 128-bitowego tasowania z AVX512), więc można go zaimplementować tylko z rozpakowaniami, których nie sądzę, aby można je było przekonwertować na blendy (wydajnie). Z AVX512 vshufi32x4działa skutecznie tak, jak shufpsz wyjątkiem 128-bitowych ścieżek 4 liczb całkowitych zamiast 32-bitowych liczb zmiennoprzecinkowych, więc ta sama technika może być vshufi32x4w niektórych przypadkach. W Knights Landing tasowanie jest cztery razy wolniejsze (przepustowość) niż mieszanki.

Question 5

Traktuj każdy wiersz jako kolumnę, a każdą kolumnę jako wiersz… użyj j, i zamiast i, j

demo: 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;
}

Question 6

transpozycja bez narzutów (klasa niekompletna):

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

można używać w ten sposób:

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

oczywiście nie zawracałem sobie głowy zarządzaniem pamięcią, co jest kluczowe, ale inny temat.

Question 7

Jeśli rozmiar tablic jest znany wcześniej, możemy użyć unii do naszej pomocy. Lubię to-

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

Question 8

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

Question 9

Nowoczesne biblioteki algebry liniowej zawierają zoptymalizowane wersje najpopularniejszych operacji. Wiele z nich obejmuje dynamiczną wysyłkę procesora, która wybiera najlepszą implementację dla sprzętu w czasie wykonywania programu (bez uszczerbku dla przenośności).

Jest to zwykle lepsza alternatywa dla wykonywania ręcznej optymalizacji twoich functinos poprzez wewnętrzne funkcje rozszerzeń wektorowych. Ta ostatnia wiąże twoją implementację z konkretnym dostawcą sprzętu i modelem: jeśli zdecydujesz się na zamianę na innego dostawcę (np. Power, ARM) lub na nowsze rozszerzenia wektorowe (np. AVX512), będziesz musiał ponownie zaimplementować go, aby uzyskać jak najwięcej z nich.

Na przykład transpozycja MKL obejmuje funkcję rozszerzeń BLAS imatcopy. Możesz go znaleźć również w innych implementacjach, takich jak 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);
}

W przypadku projektu C ++ możesz skorzystać z Armadillo C ++:

#include <armadillo>

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

Question 10

Intel mkl sugeruje macierze transpozycji / kopiowania w miejscu i poza miejscem. tutaj jest link do dokumentacji . Zalecałbym wypróbowanie implementacji nie na miejscu, ponieważ szybsza dziesiątka w miejscu i dokumentacja najnowszej wersji mkl zawiera błędy.

Question 11

Myślę, że najszybszy sposób nie powinien przyjmować wartości wyższej niż O (n ^ 2), również w ten sposób możesz użyć tylko O (1) przestrzeni:
sposobem na to jest zamiana parami, ponieważ kiedy transponujesz macierz, to co do to: M [i] [j] = M [j] [i], więc zapisz M [i] [j] w temp., a następnie M [i] [j] = M [j] [i], a ostatni krok: M [j] [i] = temp. można to zrobić jednym przebiegiem, więc powinno zająć O (n ^ 2)

Question 12

moja odpowiedź jest transponowana z macierzy 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;
}