SSE矩阵-矩阵乘法

Question

我在用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);
    }
  }
}

我似乎无法使它给出正确的结果。我想念什么吗？而且搜索搜索量似乎有很大帮助-每个结果要么只是做4x4矩阵，mat-vec要么是某些特殊的魔术，这些魔术不是很易读并且很难理解...

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缓存中仍然很热。或循环遍历目标行并添加一个结果，而不是针对单个或一小组行*列点积累加完整结果。高速缓存阻塞（又名循环切片）是获得良好的性能的关键之一。另请参见What Every Programmer Should Know About Memory?，该文件的附录中有一个缓存阻止的SIMD FP matmul示例，没有转置。

关于使用SIMD和缓存块优化矩阵乘法的文章很多。我建议你用谷歌搜索。如果可能是在谈论FP，则大多数情况下，但它也适用于整数。

---

实际上，我认为您可以在此处并行生成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来打破对先前迭代值的依赖链，因此更糟。但是与负载+ PUNPCKLDQ相比，即使是4个分散元素的最佳集合也非常糟糕。更不用说这使您的代码需要SSE4.1。
尽管对于_mm_mullo_epi32仍然需要SSE4.1，而不是加宽PMULDQ (_mm_mul_epi32)。
请注意，整数乘法吞吐量通常比FP乘法差，特别是在Haswell及更高版本上。 FP FMA单元的每个32位元素只有24位宽的乘法器（对于FP尾数），因此对于32x32 => 32位整数使用乘数需要分成两个微码。

Answer 2

这是由OP发布的，是对它不属于的问题的编辑。

---

更新：哇！我终于弄明白了。除了我的逻辑中的错误（感谢Peter Cordes的帮助）之外，还有_mm_mul_epi32不能按我认为的那样工作的问题-我应该一直使用_mm_mul_epi32()！

我知道这不是最有效的代码，但它是为了使它首先正常工作而设计的-现在，我可以继续对其进行优化。

_mm_mullo_epi32()

Update 2：

将Peters的示例转换为i-k-j循环顺序版本。 vR需要额外的负载，并且需要在存储中移动到内部循环，但是可以将设置vA向上移动一个循环。结果更快。

// editor's note: this is the most naive way to vectorize
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);     // false dependency on old vB
                vB = _mm_insert_epi32(vB, mat2[k + 1][j], 1);  // bad spatial locality
                vB = _mm_insert_epi32(vB, mat2[k + 2][j], 2);  // striding down a column
                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");
            }
        }
    }
}