SSE行列 - 行列乗算

Question

は、私がこれまでに得たものである：

#define N 1000 

void matmulSSE(int mat1[N][N], int mat2[N][N], int result[N][N]) { 
    int i, j, k; 
    __m128i vA, vB, vR; 

    for(i = 0; i < N; ++i) { 
    for(j = 0; j < N; ++j) { 
     vR = _mm_setzero_si128(); 
     for(k = 0; k < N; k += 4) { 
      //result[i][j] += mat1[i][k] * mat2[k][j]; 
      vA = _mm_loadu_si128((__m128i*)&mat1[i][k]); 
      vB = _mm_loadu_si128((__m128i*)&mat2[k][j]); //how well does the k += 4 work here? Should it be unrolled? 
      vR = _mm_add_epi32(vR, _mm_mul_epi32(vA, vB)); 
     } 
     vR = _mm_hadd_epi32(vR, vR); 
     vR = _mm_hadd_epi32(vR, vR); 
     result[i][j] += _mm_extract_epi32(vR, 0); 
    } 
    } 
} 

私はそれを作るように見えることはできません正しい結果を出してください。何か不足していますか？ そしてdosent検索は非常に役立つように見える - すべての結果は唯一の4x4行列をやって、マット-VECまたは非常に読みやすく、理解するのは難しいことではないthatsのいくつかの特別な魔法のどちらかです...

更新： Woho！私はついにそれを理解した。私のロジックのエラー（ヘルプPeter Cordesのおかげです）に加えて、_mm_mul_epi32（）の問題もありました。私は思ったとおりに動作しませんでした。代わりに_mm_mullo_epi32（）を使用していたはずです！

私はこれが最も効果的なコードではないことを知っていますが、最初に正しく動作するようになっています。今は最適化を進めることができます。

void matmulSSE(int mat1[N][N], int mat2[N][N], int result[N][N]) { 
    int i, j, k; 
    __m128i vA, vB, vR, vSum; 

    for(i = 0; i < N; ++i) { 
     for(j = 0; j < N; ++j) { 
      vR = _mm_setzero_si128(); 
      for(k = 0; k < N; k += 4) { 
       //result[i][j] += mat1[i][k] * mat2[k][j]; 
       vA = _mm_loadu_si128((__m128i*)&mat1[i][k]); 
       vB = _mm_insert_epi32(vB, mat2[k][j], 0); 
       vB = _mm_insert_epi32(vB, mat2[k + 1][j], 1); 
       vB = _mm_insert_epi32(vB, mat2[k + 2][j], 2); 
       vB = _mm_insert_epi32(vB, mat2[k + 3][j], 3); 
       vR = _mm_mullo_epi32(vA, vB); 
       vR = _mm_hadd_epi32(vR, vR); 
       vR = _mm_hadd_epi32(vR, vR); 
       result[i][j] += _mm_extract_epi32(vR, 0); 

       //DEBUG 
       //printf("vA: %d, %d, %d, %d\n", vA.m128i_i32[0], vA.m128i_i32[1], vA.m128i_i32[2], vA.m128i_i32[3]); 
       //printf("vB: %d, %d, %d, %d\n", vB.m128i_i32[0], vB.m128i_i32[1], vB.m128i_i32[2], vB.m128i_i32[3]); 
       //printf("vR: %d, %d, %d, %d\n", vR.m128i_i32[0], vR.m128i_i32[1], vR.m128i_i32[2], vR.m128i_i32[3]); 
       //printf("\n"); 
      } 
     } 
    } 
} 

アップデート2：がi-K-Jループ順バージョンにピーターズ例を変換します。 vRに余分な負荷が必要で、ストア内でインナーループに移動する必要がありますが、設定されたvAをループの上に移動できます。速くなった。

void matmulSSE_2(int mat1[N][N], int mat2[N][N], int result[N][N]) { 
    int i, j, k; 
    __m128i vA, vB, vR; 

    for(i = 0; i < N; ++i) { 
     for(k = 0; k < N; ++k) { 
      vA = _mm_set1_epi32(mat1[i][k]); 
      for(j = 0; j < N; j += 4) { 
       //result[i][j] += mat1[i][k] * mat2[k][j]; 
       vB = _mm_loadu_si128((__m128i*)&mat2[k][j]); 
       vR = _mm_loadu_si128((__m128i*)&result[i][j]); 
       vR = _mm_add_epi32(vR, _mm_mullo_epi32(vA, vB)); 
       _mm_storeu_si128((__m128i*)&result[i][j], vR); 

       //DEBUG 
       //printf("vA: %d, %d, %d, %d\n", vA.m128i_i32[0], vA.m128i_i32[1], vA.m128i_i32[2], vA.m128i_i32[3]); 
       //printf("vB: %d, %d, %d, %d\n", vB.m128i_i32[0], vB.m128i_i32[1], vB.m128i_i32[2], vB.m128i_i32[3]); 
       //printf("vR: %d, %d, %d, %d\n", vR.m128i_i32[0], vR.m128i_i32[1], vR.m128i_i32[2], vR.m128i_i32[3]); 

       //printf("\n"); 
      } 
     } 
    } 
} 

Answer 1

あなたが正しく、あなたのvBが問題です。あなたは4つの連続した整数をロードしているが、mat2[k+0..3][j]は連続していません。あなたは実際にmat2[k][j+0..3]を取得しています。

---

私はmatmulにとってうまくいくか忘れてしまいます。場合によっては、結果ごとに水平方向の合計を計算するのではなく、4つの結果を並列に生成することもできます。

入力行列の転置はうまく機能し、コストはO（N^2）です。これは、O（N^3）matmulがシーケンシャルアクセスを使用できることと、あなたの現在のループ構造がSIMDに優しいことになるためです。

小さなブロックを使用する直前に転置すると、L1キャッシュで再び熱くなります。ループ・タイリングとも呼ばれるキャッシュ・ブロッキングは、優れたマットル性能の鍵です。

行列乗算の最適化について、SIMDとキャッシュブロックを使用して詳しく説明しました。私はそれをグーグルに勧めます。ほとんどの場合、おそらくFPについて話しているのですが、すべてが整数にも適用されます。

（SSE/AVXのみFPのためにではなく、32ビット整数のFMAを有し、8および16ビット入力PMADD命令が水平を行うことを除いては、ペアの追加。）

---

実はあなたがここに並列にを4つの結果を生み出すことができると思います。

void matmulSSE(int mat1[N][N], int mat2[N][N], int result[N][N]) { 

    for(int i = 0; i < N; ++i) { 
    for(int j = 0; j < N; j+=4) { // vectorize over this loop 
     __m128i vR = _mm_setzero_si128(); 
     for(int k = 0; k < N; k++) { // not this loop 
      //result[i][j] += mat1[i][k] * mat2[k][j]; 
      __m128i vA = _mm_set1_epi32(mat1[i][k]); // load+broadcast is much cheaper than MOVD + 3 inserts (or especially 4x insert, which your new code is doing) 
      __m128i vB = _mm_loadu_si128((__m128i*)&mat2[k][j]); // mat2[k][j+0..3] 
      vR = _mm_add_epi32(vR, _mm_mullo_epi32(vA, vB)); 
     } 
     _mm_storeu_si128((__m128i*)&result[i][j], vR)); 
    } 
    } 
} 

（AVXなしまたは別の負荷+放送）の放送負荷まだ集めるよりもはるかに安いです。

現在のコードではなく、最初の要素のMOVDを使用して、前の反復の値に依存チェーンを壊すの4つのインサートとギャザーないので、それはさらに悪いです。しかし、4つの散在した要素の最高の集まりでさえ、負荷+ PSHUFDに比べてかなり悪いです。言ってみれば、あなたのコードにはSSE4.1が必要です。とにかく、PMULDQ (_mm_mul_epi32)の代わりに_mm_mullo_epi32のために。