c - 使用SSE vector 指令加速矩阵矩阵乘法

我在使用SSE矢量指令对某些C代码进行矢量化时遇到了一些困难。我必须战胜的准则是

#define N 1000
void matrix_mul(int mat1[N][N], int mat2[N][N], int result[N][N])
{
   int i, j, k;
   for (i = 0; i < N; ++i)
   {
      for (j = 0; j < N; ++j)
      {
         for (k = 0; k < N; ++k)
         {
              result[i][k] += mat1[i][j] * mat2[j][k];
         }
      }
   }
}

到目前为止我得到的是：

void  matrix_mul_sse(int mat1[N][N], int mat2[N][N], int result[N][N])
{
   int i, j, k; int* l;
   __m128i v1, v2, v3;
   v3 = _mm_setzero_si128();
   for (i = 0; i < N; ++i)
   {
       for (j = 0; j < N; j += 4)
       {

           for (k = 0; k < N; k += 4)
           {

               v1 = _mm_set1_epi32(mat1[i][j]);
               v2 = _mm_loadu_si128((__m128i*)&mat2[j][k]);
               v3 = _mm_add_epi32(v3, _mm_mul_epi32(v1, v2));
               _mm_storeu_si128((__m128i*)&result[i][k], v3);
               v3 = _mm_setzero_si128();
           }
       }
   }
}

处决后我得到了错误的结果。我知道原因是从内存加载到v2。我按行主顺序遍历mat1，因此需要加载mat2[0][0]、mat2[1][0]、mat2[2][0]、mat2[3][0]。。。。但实际加载的是mat2[0][0]、mat2[0][1]、mat2[0][2]、mat2[0][3]。。。因为mat2已经按行的主要顺序存储在内存中。我试图解决这个问题，但没有任何改进。
有人能帮帮我吗。

最佳答案

下面修复了您的实现：

void  matrix_mul_sse(int mat1[N][N], int mat2[N][N], int result[N][N])
{
   int i, j, k;
   __m128i v1, v2, v3, v4;
   for (i = 0; i < N; ++i)
   {
       for (j = 0; j < N; ++j) // 'j' must be incremented by 1
       {
           // read mat1 here because it does not use 'k' index
           v1 = _mm_set1_epi32(mat1[i][j]);
           for (k = 0; k < N; k += 4)
           {
               v2 = _mm_loadu_si128((const __m128i*)&mat2[j][k]);

               // read what's in the result array first as we will need to add it later to our calculations
               v3 = _mm_loadu_si128((const __m128i*)&result[i][k]);

               // use _mm_mullo_epi32 here instead _mm_mul_epi32 and add it to the previous result
               v4 = _mm_add_epi32(v3, _mm_mullo_epi32(v1, v2));

               // store the result
               _mm_storeu_si128((__m128i*)&result[i][k], v4);
           }
       }
   }
}

简而言之，_mm_mullo_epi32（需要SSE4.1）产生4个int32结果，而_mm_mul_epi32产生2个int64结果。如果您不能使用SSE4.1，请查看替代SSE2解决方案的答案here。
完整描述由Intel Intrinsic Guide：
_mm_mullo_epi32：将a和b中的压缩32位整数相乘，生成中间的64位整数，然后存储
dst中中间整数的低32位。
_mm_mul_epi32：将a和b中每个压缩64位元素的低32位整数相乘，并存储
有符号的64位结果是dst。

关于c - 使用SSE vector 指令加速矩阵矩阵乘法，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/53064185/