転置行列をより効率的に乗算するにはどうすればよいですか？

Question

私はこれで頭の中で自分自身を倒し続けます。行列Aに行列Bを乗ずるためのSSEベースのアルゴリズムがあります。私は、A、B、またはその両方が転置される場所の操作も実装する必要があります。私はそれを素朴な実装、4x4マトリックスコード（これはかなり標準的なSSE操作ですが、私は思う）を示しましたが、A*B^T操作は約の2倍の時間がかかります。 ATLASの実装では、A*Bの類似の値が返されます。これは、効率的な方法があることを示唆しています。

MM-乗算：MT乗算（転置行列を乗じた行列）について

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 
} 

は、私は、次のコマンドを使用してB3Rするb0rに格納され、適宜ループ変数を変更代え：

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

私は、減速は部分的に一度に行を引っ張ることと、列を取得するたびに4つの値を格納しなければならないという違いが原因だと思っていますが、私はBの行を引っ張って、その後、 Asの列は、結果の4列の格納にコストをシフトするだけです。

また、Bの行を行として引っ張ってから、_MM_TRANSPOSE4_PS(b0r, b1r, b2r, b3r);を使用して転置しましたが（そのマクロにはいくつかの追加の最適化があると思いましたが）、実際の改善はありません。

表面上、私はこれがより速くなければならないと感じています...関与するドットプロダクトは、本質的に効率的なように見える行ですが、ドットプロダクトをまっすぐにしようとすると同じ結果を保存する。

私はここで何が欠けていますか？

投稿日：私は行列を転置しないようにしようとしています。私は彼らに沿って繰り返すことを好むでしょう。問題は、私が知る限り、_mm_set_psコマンドは_mm_load_psよりもはるかに遅いということです。

また、私はA行列の4行を格納し、1ロード、4乗、および4乗算命令と3を含む2つの中括弧で囲まれたセグメントを置き換えました。 。タイミングは同じまま（はい、私はコードは私のテストコンパイルに変更したことを確認するために、デバッグ文でそれを試してみました。前記デバッグ文はもちろん、プロファイリングする前に削除されました。）：

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

更新日： 右に移動して、a0rの行のロードを中括弧に入れてa3rにすると、レジスタのスラッシングを回避することができません。

Answer 1

これは、水平方向の追加が便利な場合が1つあると思います。あなたはC = A B^Tを望んでいますが、Bは転置としてメモリに格納されていません。それが問題です。これは、SoAではなくAoSのような店舗です。この場合、Bの転置を取って垂直の加算を行うのは、水平の加算を使うよりも遅いと思います。マトリックスベクトルEfficient 4x4 matrix vector multiplication with SSE: horizontal add and dot product - what's the point?の場合、これは少なくとも真です。機能以下のコードm4x4非SSEの4×4マトリクス生成物である、m4x4_vecはSSEを使用し、 m4x4TでCは、SSEなしB^Tを=、及びm4x4T_vec Cは B^T usisg SSEを=いません。最後の1つは、私が思っているものです。

注：大規模な行列の場合は、このメソッドを使用しません。その場合は、最初にトランスポーズを取って垂直方向の追加（SSE/AVXを使用すると複雑な作業を行う場合、ストリップをSSE/AVX幅に置き換える方が速い）です。これは、転置がO（n^2）となり、行列積がO（n^3）となるためです。大きな行列の場合、転置は重要ではありません。しかし、4x4の場合、転置は重要なので水平方向の加算が勝つ。

編集： 私はあなたが望むものを誤解しました。あなたはC =（A B）^ Tが欲しいです。それは（ B）と同じくらい高速である必要がありますし、コードを使用すると、基本的には次のように私たちは数学を書くことができますAとB
の役割を入れ替えるほぼ同じである：

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 

あなたが比較した場合2つの唯一の変化は、jとiの役割を交換することです。私はこの答えの最後にこれを行うためのいくつかのコードを書いています。

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

} 

編集：ここ

はC =（B）^ T行うスカラー及びSSEの関数です。それらは、A Bバージョンと同じくらい速くなければなりません。

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

助けるかもしれないいくつかの提言：

• （* _mm_loaduそれらが遅い）非整列のメモリを使用しないでください。
• メモリを順番にアクセスしていないため、キャッシュが強制終了されます。そのメモリへの実際のアクセスの前に行列を転置してみてください。これにより、CPUがキャッシュをできるだけフェッチして使用できるようになります。このようにして、次の__m128 b0r = _mm_set_ps(b3[0], b2[0], b1[0], b0[0]); // and b1r, etc..は必要ありません。考え方は、4つのコンポーネント全体を順番に取得することです。 SSEコードが呼び出される前にメモリを再編成する必要がある場合は、そうしてください。
• 内部ループ内にロードしています。_mm_load_ps1(a0+0)（a1、a2、a3と同じ）ですが、内側ループのすべての反復で定数です。これらの値を外部にロードして、いくつかのサイクルを保存することができます。以前の繰り返しから再利用できるものを見守ってください。
• プロファイル。 Intel VTuneなどを使用すると、ボトルネックがどこにあるかがわかります。