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矩阵，矩阵乘法

回答

相关问题