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);
    }
  }
}

我似乎无法让它给出正确的结果。我错过了什么吗？ 搜索似乎有很大帮助 - 每个结果要么只做 4x4 矩阵、mat-vec，要么是一些不太可读且难以理解的特殊魔法......

你说得对，你的vB是问题所在。您正在加载 4 个连续的整数，但是mat2[k+0..3][j]不连续。你实际上得到了mat2[k][j+0..3].

我忘记什么对 matmul 有效。有时并行产生 4 个结果效果很好，而不是对每个结果进行水平求和。

转置输入矩阵之一是可行的，并且成本为 O(N^2)。这是值得的，因为这意味着 O(N^3) matmul 可以使用顺序访问，并且当前的循环结构变得 SIMD 友好。

还有更好的方法，例如在使用之前调换小块，这样当您再次读取它们时，它们在 L1 缓存中仍然是热的。或者循环遍历目标行并添加一个结果，而不是累积单个或一小组行*列点积的完整结果。缓存阻塞（又称为循环平铺）是良好的 matmul 性能的关键之一。也可以看看每个程序员都应该了解哪些关于内存的知识？附录中有一个缓存阻塞的 SIMD FP matmul 示例，没有转置。

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

（除了 SSE/AVX 仅具有用于 FP 的 FMA，而不具有用于 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 的单独加载+广播）仍然比聚集便宜得多。

您当前的代码通过 4 次插入进行收集，而不是通过对第一个元素使用 MOVD 来破坏上一次迭代值的依赖链，因此情况更糟。但与负载 + PUNPCKLDQ 相比，即使是 4 个分散元素的最佳聚集也相当糟糕。更不用说这使得您的代码需要 SSE4.1。

虽然无论如何它都需要SSE4.1_mm_mullo_epi32而不是扩大PMULDQ (_mm_mul_epi32).

请注意，整数乘法吞吐量通常比 FP 乘法差，尤其是在 Haswell 及更高版本上。 FP FMA 单元每个 32 位元素仅具有 24 位宽乘法器（对于 FP 尾数），因此将这些乘法器用于 32x32=>32 位整数需要拆分为两个微指令。