Wie kann ich Transpositionsmatrizen effizienter multiplizieren?

Question

Ich prügele mich immer mit dem Kopf über den Kopf. Ich habe einen SSE-basierten Algorithmus zum Multiplizieren der Matrix A mit der Matrix B. Ich muss auch die Operationen implementieren, für die A, B oder beide transponiert sind. Ich habe eine naive Implementierung davon gemacht, den 4x4-Matrix-Code, der unten dargestellt ist (was ziemlich normale SSE-Operationen sind, glaube ich), aber die Operation A*B^T dauert ungefähr doppelt so lang wie A*B. Die ATLAS-Implementierung gibt ähnliche Werte für A*B zurück, und fast identische Ergebnisse für die Multiplikation mit einer Transponierung, was für mich einen effizienten Weg zu dieser Vorgehensweise nahelegt.

MM-Multiplikation:

m1 = (mat1.m_>>2)<<2; 
n2 = (mat2.n_>>2)<<2; 
n = (mat1.n_>>2)<<2; 

for (k=0; k<n; k+=4) { 
    for (i=0; i<m1; i+=4) { 
    // fetch: get 4x4 matrix from mat1 
    // row-major storage, so get 4 rows 
    Float* a0 = mat1.el_[i]+k; 
    Float* a1 = mat1.el_[i+1]+k; 
    Float* a2 = mat1.el_[i+2]+k; 
    Float* a3 = mat1.el_[i+3]+k; 

    for (j=0; j<n2; j+=4) { 
     // fetch: get 4x4 matrix from mat2 
     // row-major storage, so get 4 rows 
     Float* b0 = mat2.el_[k]+j; 
     Float* b1 = mat2.el_[k+1]+j; 
     Float* b2 = mat2.el_[k+2]+j; 
     Float* b3 = mat2.el_[k+3]+j; 

     __m128 b0r = _mm_loadu_ps(b0); 
     __m128 b1r = _mm_loadu_ps(b1); 
     __m128 b2r = _mm_loadu_ps(b2); 
     __m128 b3r = _mm_loadu_ps(b3); 

     { // first row of result += first row of mat1 * 4x4 of mat2 
     __m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a0+0), b0r), _mm_mul_ps(_mm_load_ps1(a0+1), b1r)); 
     __m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a0+2), b2r), _mm_mul_ps(_mm_load_ps1(a0+3), b3r)); 
     Float* c0 = this->el_[i]+j; 
     _mm_storeu_ps(c0, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c0))); 
     } 

     { // second row of result += second row of mat1 * 4x4 of mat2 
     __m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a1+0), b0r), _mm_mul_ps(_mm_load_ps1(a1+1), b1r)); 
     __m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a1+2), b2r), _mm_mul_ps(_mm_load_ps1(a1+3), b3r)); 
     Float* c1 = this->el_[i+1]+j; 
     _mm_storeu_ps(c1, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c1))); 
     } 

     { // third row of result += third row of mat1 * 4x4 of mat2 
     __m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a2+0), b0r), _mm_mul_ps(_mm_load_ps1(a2+1), b1r)); 
     __m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a2+2), b2r), _mm_mul_ps(_mm_load_ps1(a2+3), b3r)); 
     Float* c2 = this->el_[i+2]+j; 
     _mm_storeu_ps(c2, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c2))); 
     } 

     { // fourth row of result += fourth row of mat1 * 4x4 of mat2 
     __m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a3+0), b0r), _mm_mul_ps(_mm_load_ps1(a3+1), b1r)); 
     __m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a3+2), b2r), _mm_mul_ps(_mm_load_ps1(a3+3), b3r)); 
     Float* c3 = this->el_[i+3]+j; 
     _mm_storeu_ps(c3, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c3))); 
     } 
    } 
// Code omitted to handle remaining rows and columns 
} 

für das MT Multiplikation (Matrix durch transponierten Matrix multipliziert), anstelle I B0R gespeichert mit den folgenden Befehlen B3R und verändert die Schleifenvariablen entsprechend:

__m128 b0r = _mm_set_ps(b3[0], b2[0], b1[0], b0[0]); 
__m128 b1r = _mm_set_ps(b3[1], b2[1], b1[1], b0[1]); 
__m128 b2r = _mm_set_ps(b3[2], b2[2], b1[2], b0[2]); 
__m128 b3r = _mm_set_ps(b3[3], b2[3], b1[3], b0[3]); 

Ich vermute, dass die Verlangsamung teilweise auf den Unterschied zwischen dem Ziehen einer Reihe zu einer Zeit und dem Speichern von 4 Werten jedes Mal zurückzuführen ist, um die Spalte zu erhalten, aber ich fühle mich wie der andere Weg, um in Reihen von B zu ziehen und dann mit dem co multiplizieren lumn of As, verschiebt die Kosten einfach auf 4 Spalten mit Ergebnissen.

Ich habe auch versucht, B-Zeilen als Zeilen ziehen und dann _MM_TRANSPOSE4_PS(b0r, b1r, b2r, b3r); verwenden, um die Transposition (ich dachte, es könnte einige zusätzliche Optimierungen in diesem Makro sein), aber es gibt keine echte Verbesserung.

Auf der Oberfläche, ich fühle mich wie dies schneller sein sollte ... die Punktprodukte eine Reihe durch eine Reihe wäre beteiligt, die von Natur aus effizienten scheint, aber versucht, die Punktprodukte gerade nach oben gerade ergibt, die zu tun Machen Sie dasselbe, um die Ergebnisse zu speichern.

Was fehlt mir hier?

Hinzugefügt: Nur um zu verdeutlichen, versuche ich nicht die Matrizen zu transponieren. Ich würde es vorziehen, entlang ihnen zu iterieren. Das Problem, so gut ich es beurteilen kann, ist, dass der Befehl _mm_set_ps viel langsamer ist als _mm_load_ps.

Ich habe auch eine Variation, wo ich die 4 Zeilen der Matrix A gespeichert und dann ersetzt die 4 Curly-klammert Segmente, die 1 Last, 4 Multiplikationen und 2 addiert mit 4 Multiplikationsbefehle und 3 hadds, aber mit wenig Erfolg . Der Zeitpunkt bleibt gleich (und ja, ich versuchte es mit einer Debug-Anweisung, um sicherzustellen, dass der Code in meinem Test-Kompilierung Said Debug-Anweisung vor dem Profilieren entfernt wurde geändert hatte, natürlich.):

{ // first row of result += first row of mat1 * 4x4 of mat2 
     __m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a0r, b0r), _mm_mul_ps(a0r, b1r)); 
     __m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a0r, b2r), _mm_mul_ps(a0r, b3r)); 
     Float* c0 = this->el_[i]+j; 
     _mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0))); 
    } 

    { // second row of result += second row of mat1 * 4x4 of mat2 
     __m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a1r, b0r), _mm_mul_ps(a1r, b1r)); 
     __m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a1r, b2r), _mm_mul_ps(a1r, b3r)); 
     Float* c0 = this->el_[i+1]+j; 
     _mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0))); 
    } 

    { // third row of result += third row of mat1 * 4x4 of mat2 
     __m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a2r, b0r), _mm_mul_ps(a2r, b1r)); 
     __m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a2r, b2r), _mm_mul_ps(a2r, b3r)); 
     Float* c0 = this->el_[i+2]+j; 
     _mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0))); 
    } 

    { // fourth row of result += fourth row of mat1 * 4x4 of mat2 
     __m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a3r, b0r), _mm_mul_ps(a3r, b1r)); 
     __m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a3r, b2r), _mm_mul_ps(a3r, b3r)); 
     Float* c0 = this->el_[i+3]+j; 
     _mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0))); 
    } 

Update: Richtig, und das Laden der Zeilen von a0r zu a3r in die geschweiften Klammern in einem Versuch, Register Thrashing zu vermeiden fehlgeschlagen.

Answer 1

Ich denke, das ist ein paar Fälle, wo horizontale hinzufügen ist nützlich. Sie wollen C = A B^T, aber B wird nicht als Transponierung gespeichert. Das ist das Problem. Es ist Speicher wie ein AoS anstelle eines SoA. In diesem Fall ist die Transponierung von B und das vertikale Addieren langsamer als die Verwendung horizontaler Adds, denke ich. Dies gilt zumindest für die Matrix Vektor Efficient 4x4 matrix vector multiplication with SSE: horizontal add and dot product - what's the point?. In dem Code unter der Funktion m4x4 ist nicht SSE 4x4 Matrix-Produkt, m4x4_vec verwendet SSE, m4x4T nicht = C A B^T ohne SSE, und tut m4x4T_vec C A B^T usisg SSE =. Der letzte ist der, den du willst, denke ich.

Hinweis: Für größere Matrizen würde ich diese Methode nicht verwenden. In diesem Fall ist es schneller, die Transponierung zuerst zu nehmen und vertikale Add zu verwenden (mit SSE/AVX macht man etwas komplizierteres, man transponiert Streifen mit der SSE/AVX-Breite). Das liegt daran, dass die Transponierte als O (n^2) und das Matrixprodukt als O (n^3) gilt. Für große Matrizen ist die Transponierung also unbedeutend. Für 4x4 ist die Transponierung jedoch signifikant, so dass die horizontale Addition gewinnt.

Edit: Ich missverstanden, was Sie wollten. Sie wollen C = (A B)^T. Das sollte nur so schnell wie (A B) und der Code ist fast das gleiche Sie im Grunde tauschen nur die Rollen von A und B.
Wir können die Mathematik schreiben wie folgt:

C = A*B in Einstein notation is C_i,j = A_i,k * B_k,j. 
Since (A*B)^T = B^T*A^T we can write 
C = (A*B)^T in Einstein notation is C_i,j = B^T_i,k * A^T_k,j = A_j,k * B_k,i 

Wenn Sie vergleichen Die zwei die einzige Sache, die sich ändert, ist, wir tauschen die Rollen von j und ich aus. Ich habe am Ende dieser Antwort einen Code dafür eingegeben.

#include "stdio.h" 
#include <nmmintrin.h>  

void m4x4(const float *A, const float *B, float *C) { 
    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      float sum = 0.0f; 
      for(int k=0; k<4; k++) { 
       sum += A[i*4+k]*B[k*4+j]; 
      } 
      C[i*4 + j] = sum; 
     } 
    } 
} 

void m4x4T(const float *A, const float *B, float *C) { 
    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      float sum = 0.0f; 
      for(int k=0; k<4; k++) { 
       sum += A[i*4+k]*B[j*4+k]; 
      } 
      C[i*4 + j] = sum; 
     } 
    } 
} 

void m4x4_vec(const float *A, const float *B, float *C) { 
    __m128 Brow[4], Mrow[4]; 
    for(int i=0; i<4; i++) { 
     Brow[i] = _mm_load_ps(&B[4*i]); 
    } 

    for(int i=0; i<4; i++) { 
     Mrow[i] = _mm_set1_ps(0.0f); 
     for(int j=0; j<4; j++) { 
      __m128 a = _mm_set1_ps(A[4*i +j]); 
      Mrow[i] = _mm_add_ps(Mrow[i], _mm_mul_ps(a, Brow[j])); 
     } 
    } 
    for(int i=0; i<4; i++) { 
     _mm_store_ps(&C[4*i], Mrow[i]); 
    } 
} 

void m4x4T_vec(const float *A, const float *B, float *C) { 
    __m128 Arow[4], Brow[4], Mrow[4]; 
    for(int i=0; i<4; i++) { 
     Arow[i] = _mm_load_ps(&A[4*i]); 
     Brow[i] = _mm_load_ps(&B[4*i]); 
    } 

    for(int i=0; i<4; i++) { 
     __m128 prod[4]; 
     for(int j=0; j<4; j++) { 
      prod[j] = _mm_mul_ps(Arow[i], Brow[j]); 
     } 
     Mrow[i] = _mm_hadd_ps(_mm_hadd_ps(prod[0], prod[1]), _mm_hadd_ps(prod[2], prod[3]));  
    } 
    for(int i=0; i<4; i++) { 
     _mm_store_ps(&C[4*i], Mrow[i]); 
    } 

} 

float compare_4x4(const float* A, const float*B) { 
    float diff = 0.0f; 
    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      diff += A[i*4 +j] - B[i*4+j]; 
      printf("A %f, B %f\n", A[i*4 +j], B[i*4 +j]); 
     } 
    } 
    return diff;  
} 

int main() { 
    float *A = (float*)_mm_malloc(sizeof(float)*16,16); 
    float *B = (float*)_mm_malloc(sizeof(float)*16,16); 
    float *C1 = (float*)_mm_malloc(sizeof(float)*16,16); 
    float *C2 = (float*)_mm_malloc(sizeof(float)*16,16); 

    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      A[i*4 +j] = i*4+j; 
      B[i*4 +j] = i*4+j; 
      C1[i*4 +j] = 0.0f; 
      C2[i*4 +j] = 0.0f; 
     } 
    } 
    m4x4T(A, B, C1); 
    m4x4T_vec(A, B, C2); 
    printf("compare %f\n", compare_4x4(C1,C2)); 

} 

Edit:

Hier sind die skalare und SSE-Funktion, C = (A B)^T tun. Sie sollten genauso schnell wie ihre A B-Versionen sein.

void m4x4TT(const float *A, const float *B, float *C) { 
    for(int i=0; i<4; i++) { 
     for(int j=0; j<4; j++) { 
      float sum = 0.0f; 
      for(int k=0; k<4; k++) { 
       sum += A[j*4+k]*B[k*4+i]; 
      } 
      C[i*4 + j] = sum; 
     } 
    } 
} 

void m4x4TT_vec(const float *A, const float *B, float *C) { 
    __m128 Arow[4], Crow[4]; 
    for(int i=0; i<4; i++) { 
     Arow[i] = _mm_load_ps(&A[4*i]); 
    } 

    for(int i=0; i<4; i++) { 
     Crow[i] = _mm_set1_ps(0.0f); 
     for(int j=0; j<4; j++) { 
      __m128 a = _mm_set1_ps(B[4*i +j]); 
      Crow[i] = _mm_add_ps(Crow[i], _mm_mul_ps(a, Arow[j])); 
     } 
    } 

    for(int i=0; i<4; i++) { 
     _mm_store_ps(&C[4*i], Crow[i]); 
    } 
} 

Answer 2

Einige Empfehlungen, die helfen können:

• nicht unaligned Speicher Verwenden Sie (die _mm_loadu * sind langsam).
• Sie greifen nicht sequenziell auf den Speicher zu, wodurch der Cache beendet wird. Versuchen Sie, die Matrix vor dem eigentlichen Zugriff auf diesen Speicher zu transponieren, damit die CPU den Cache so oft wie möglich abrufen und verwenden kann. Auf diese Weise wird das folgende __m128 b0r = _mm_set_ps(b3[0], b2[0], b1[0], b0[0]); // and b1r, etc.. nicht benötigt.Die Idee wäre, die ganzen 4 Komponenten nacheinander zu greifen. Wenn Sie Ihren Speicher reorganisieren müssen, bevor der SSE-Code aufgerufen wird, tun Sie dies.
• Sie laden in der inneren Schleife dies: _mm_load_ps1(a0+0)(das gleiche für a1, a2 und a3) aber ist für alle Iterationen in der inneren Schleife konstant. Sie können diese Werte nach draußen laden und einige Zyklen speichern. Behalten Sie im Auge, was Sie aus der vorherigen Iteration wiederverwenden können.
• Profil. Verwenden Sie Intel VTune oder etwas Ähnliches, es wird Ihnen sagen, wo der Flaschenhals ist.