SSE矩陣，矩陣乘法

我無法做矩陣，矩陣乘法與上交所C.SSE矩陣，矩陣乘法

這是我走到這一步：

#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或一些特殊的魔力，那不是很有可讀性，很難理解......

更新： 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"); 
      } 
     } 
    } 
}

來源

2016-10-28 Erlisch

你有什麼麻煩？ –

@RushyPanchal我沒有得到正確的結果。對不起，我應該在我的帖子中指定... – Erlisch

對於你來說調用者是否爲''結果[]'？如果沒有，你應該先做！還要注意，在最內部循環內部進行水平求和是非常可怕的。如果你在同一個最內層循環中爲'result [i] [j]'做了所有的數學運算，只要做'result = hsum（vR）'，而不是'+ ='。 hsum是一個水平和函數，可以移植到非MSVC（如果重要的話），並且比編譯器爲您寫的內容可能產生的更少。見http://stackoverflow.com/questions/6996764/fastest-way-to-do-horizontal-float-vector-sum-on-x86，我的答案提到整數hsums。 –

你說得對，你的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只有FMA爲FP，而不是32位整數，並將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）仍然便宜得多一個聚集地。

您當前的代碼執行與收集4個插件，而不是通過使用第一元素的MOVD打破上一次迭代的價值依賴鏈，所以這是更糟。但是，與負載+ PSHUFD相比，即使是4個分散元素的最佳聚集也相當糟糕。更不用說那會讓你的代碼需要SSE4.1。儘管如此，對於_mm_mullo_epi32而不是擴大PMULDQ (_mm_mul_epi32)。

來源

2016-10-30 06:55:55

是的，這是一個很大的錯誤。除此之外，我還發現SSE中的乘法運算沒有像我想的那樣工作 - 請參閱我的編輯。感謝所有的幫助。 :) – Erlisch

@Erlisch：神聖的廢話，現在你在內部循環中有2個PHADDD指令，並且有一個標量'+ ='。你是否將其與標量代碼進行比較？它可能更慢。 –

@Erlisch：我剛剛注意到我們可以並行地生成'j + 0..3'的結果，這比通過集合對'k'進行向量化更有效。更新了我的答案。 –

SSE矩陣，矩陣乘法

回答

相關問題