Ich habe eine Matrix (relativ groß), die ich transponieren muss. Nehmen wir zum Beispiel an, dass meine Matrix ist
a b c d e f
g h i j k l
m n o p q r
Ich möchte, dass das Ergebnis wie folgt lautet:
a g m
b h n
c I o
d j p
e k q
f l r
Was ist der schnellste Weg, dies zu tun?
Antworten:
Das ist eine gute Frage. Es gibt viele Gründe, warum Sie die Matrix tatsächlich im Speicher transponieren möchten, anstatt nur die Koordinaten zu tauschen, z. B. bei der Matrixmultiplikation und beim Gaußschen Verschmieren.
Lassen Sie mich zunächst eine der Funktionen auflisten, die ich für die Transponierung verwende ( BEARBEITEN: Bitte lesen Sie das Ende meiner Antwort, wo ich eine viel schnellere Lösung gefunden habe ).
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]; } }
Nun wollen wir sehen, warum die Transponierung nützlich ist. Betrachten Sie die Matrixmultiplikation C = A * B. Wir könnten es so machen.
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; } }
Auf diese Weise wird es jedoch viele Cache-Fehler geben. Eine viel schnellere Lösung besteht darin, zuerst die Transponierte von B zu nehmen
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);
Die Matrixmultiplikation ist O (n ^ 3) und die Transponierung ist O (n ^ 2), daher sollte die Transponierung einen vernachlässigbaren Einfluss auf die Rechenzeit haben (für große
n
). Bei der Matrixmultiplikationsschleife ist das Kacheln noch effektiver als das Transponieren, aber das ist viel komplizierter.Ich wünschte, ich wüsste einen schnelleren Weg, um die Transponierung durchzuführen ( Bearbeiten: Ich habe eine schnellere Lösung gefunden, siehe das Ende meiner Antwort ). Wenn Haswell / AVX2 in einigen Wochen herauskommt, hat es eine Sammelfunktion. Ich weiß nicht, ob das in diesem Fall hilfreich sein wird, aber ich könnte mir vorstellen, eine Spalte zu sammeln und eine Zeile zu schreiben. Vielleicht macht es die Transponierung unnötig.
Beim Gaußschen Verschmieren verschmieren Sie horizontal und dann vertikal. Vertikales Verschmieren hat jedoch das Cache-Problem. Sie tun dies also
Hier ist ein Artikel von Intel, der erklärt, dass http://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions
Schließlich mache ich bei der Matrixmultiplikation (und beim Gaußschen Verschmieren) nicht genau die Transponierung, sondern die Transponierung in Breiten einer bestimmten Vektorgröße (z. B. 4 oder 8 für SSE / AVX). Hier ist die Funktion, die ich benutze
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]; } }
BEARBEITEN:
Ich habe mehrere Funktionen ausprobiert, um die schnellste Transponierung für große Matrizen zu finden. Am Ende ist das schnellste Ergebnis die Verwendung der Schleifenblockierung mit
block_size=16
( Bearbeiten: Ich habe eine schnellere Lösung mit SSE und Schleifenblockierung gefunden - siehe unten ). Dieser Code funktioniert für jede NxM-Matrix (dh die Matrix muss nicht quadratisch sein).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); } } }
Die Werte
lda
undldb
sind die Breite der Matrix. Dies müssen Vielfache der Blockgröße sein. Um die Werte zu finden und den Speicher für zB eine 3000x1001-Matrix zuzuweisen, mache ich so etwas#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);
Für 3000x1001 gibt dies
ldb = 3008
und zurücklda = 1008
Bearbeiten:
Ich habe mit SSE Intrinsics eine noch schnellere Lösung gefunden:
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); } } } } }
quelle
Dies hängt von Ihrer Anwendung ab. Im Allgemeinen besteht der schnellste Weg, eine Matrix zu transponieren, darin, Ihre Koordinaten beim Nachschlagen zu invertieren. Dann müssen Sie keine Daten verschieben.
quelle
(i,j)
werden(j,i)
Einige Details zum Transponieren von 4x4 Square Float-Matrizen (ich werde später auf 32-Bit-Ganzzahlen eingehen) mit x86-Hardware. Es ist hilfreich, hier zu beginnen, um größere quadratische Matrizen wie 8x8 oder 16x16 zu transponieren.
_MM_TRANSPOSE4_PS(r0, r1, r2, r3)
wird von verschiedenen Compilern unterschiedlich implementiert. GCC und ICC (ich habe Clang nicht überprüft) verwenden,unpcklps, unpckhps, unpcklpd, unpckhpd
während MSVC nur verwendetshufps
. Wir können diese beiden Ansätze tatsächlich so miteinander kombinieren.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);
Eine interessante Beobachtung ist, dass zwei Shuffles wie folgt in ein Shuffle und zwei Blends (SSE4.1) umgewandelt werden können.
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);
Dies wandelte effektiv 4 Mischungen in 2 Mischungen und 4 Mischungen um. Dies verwendet zwei weitere Anweisungen als die Implementierung von GCC, ICC und MSVC. Der Vorteil besteht darin, dass der Anschlussdruck reduziert wird, was unter bestimmten Umständen von Vorteil sein kann. Derzeit können alle Shuffles und Unpacks nur an einen bestimmten Port gesendet werden, während die Mischungen an einen von zwei verschiedenen Ports gesendet werden können.
Ich habe versucht, 8 Shuffles wie MSVC zu verwenden und diese in 4 Shuffles + 8 Blends umzuwandeln, aber es hat nicht funktioniert. Ich musste noch 4 Auspackungen verwenden.
Ich habe dieselbe Technik für eine 8x8-Float-Transponierung verwendet (siehe gegen Ende dieser Antwort). https://stackoverflow.com/a/25627536/2542702 . In dieser Antwort musste ich noch 8 Auspackungen verwenden, aber ich habe die 8 Shuffles in 4 Shuffles und 8 Blends umgewandelt.
Für 32-Bit-Ganzzahlen gibt es nichts
shufps
Vergleichbares (außer für 128-Bit-Shuffles mit AVX512), sodass es nur mit Entpackungen implementiert werden kann, von denen ich glaube, dass sie nicht effizient in Blends konvertiert werden können. Mit AVX512 verhält es sichvshufi32x4
effektiv wie mitshufps
Ausnahme von 128-Bit-Lanes mit 4 Ganzzahlen anstelle von 32-Bit-Floats, sodass diese Technikvshufi32x4
in einigen Fällen möglicherweise angewendet werden kann . Mit Knights Landing sind Shuffles viermal langsamer (Durchsatz) als Mischungen.quelle
shufps
für ganzzahlige Daten verwenden. Wenn Sie viel mischen, lohnt es sich möglicherweise, alles in der FP-Domäne fürshufps
+ zu erledigenblendps
, insbesondere wenn Sie nicht über den ebenso effizienten AVX2 verfügenvpblendd
. Außerdem gibt es auf Hardware der Intel SnB-Familie keine zusätzliche Bypass-Verzögerung für die Verwendungshufps
zwischen ganzzahligen Anweisungen wiepaddd
. (Es gibt jedoch eine Bypass-Verzögerung für das Mischenblendps
mitpaddd
, laut Agner Fogs SnB-Test.)vinsertf64x4
in meiner 16x16-Transponierungsantwort immer wieder erwähnt habenvinserti64x4
. Wenn ich die Matrix lese und dann schreibe, spielt es sicherlich keine Rolle, ob ich die Gleitkommadomäne oder die Ganzzahldomäne verwende, da die Transponierung nur Daten verschiebt.Betrachten Sie jede Zeile als Spalte und jede Spalte als Zeile. Verwenden Sie j, i anstelle von 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; }
quelle
Transponieren ohne Overhead (Klasse nicht vollständig):
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); } }
kann wie folgt verwendet werden:
Matrix M(100,200); double x=M.get(17,45); M.transpose(); x=M.get(17,45); // = original M(45,17)
Natürlich habe ich mich hier nicht um die Speicherverwaltung gekümmert, was ein entscheidendes, aber anderes Thema ist.
quelle
Wenn die Größe der Arrays vorher bekannt ist, können wir die Union zu unserer Hilfe verwenden. So was-
#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; }
quelle
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]; } } }
quelle
Moderne lineare Algebra-Bibliotheken enthalten optimierte Versionen der gängigsten Operationen. Viele von ihnen enthalten einen dynamischen CPU-Versand, der die beste Implementierung für die Hardware zur Programmausführungszeit auswählt (ohne Kompromisse bei der Portabilität einzugehen).
Dies ist im Allgemeinen eine bessere Alternative zur manuellen Optimierung Ihrer Funktionen über intrinsische Funktionen für Vektorerweiterungen. Letzteres bindet Ihre Implementierung an einen bestimmten Hardwareanbieter und ein bestimmtes Hardwaremodell: Wenn Sie sich für einen Wechsel zu einem anderen Anbieter (z. B. Power, ARM) oder zu neueren Vektorerweiterungen (z. B. AVX512) entscheiden, müssen Sie diese erneut implementieren Holen Sie das Beste aus ihnen.
Die MKL-Transposition enthält beispielsweise die BLAS-Erweiterungsfunktion
imatcopy
. Sie finden es auch in anderen Implementierungen wie 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); }
Für ein C ++ - Projekt können Sie Armadillo C ++ verwenden:
#include <armadillo> void transpose( arma::mat &matrix ) { arma::inplace_trans(matrix); }
quelle
intel mkl schlägt In-Place- und Out-of-Place-Transpositions- / Kopiermatrizen vor. Hier ist der Link zur Dokumentation . Ich würde empfehlen, die Implementierung an Ort und Stelle zu versuchen, da die schnellere Installation vor Ort und in der Dokumentation der neuesten Version von mkl einige Fehler enthält.
quelle
Ich denke, dass der schnellste Weg nicht höher als O (n ^ 2) sein sollte, auch auf diese Weise können Sie nur O (1) -Raum verwenden:
Der Weg, dies zu tun, besteht darin, paarweise zu tauschen, denn wenn Sie eine Matrix transponieren, dann was Sie tun ist: M [i] [j] = M [j] [i], also speichere M [i] [j] in Temp, dann M [i] [j] = M [j] [i] und die letzter Schritt: M [j] [i] = temp. Dies könnte in einem Durchgang erfolgen, daher sollte O (n ^ 2) benötigt werden.
quelle
Meine Antwort ist aus 3x3 Matrix transponiert
#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; }
quelle