logo标签列表

SSE矩阵乘法

cssematrix-multiplication
3

我在使用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矩阵，矩阵-向量操作，或者一些特殊的魔法操作，这些都不容易读懂和理解...
1
你遇到了什么问题？ - Rushy Panchal
@RushyPanchal 我没有得到正确的结果。抱歉，我应该在我的帖子中指明这一点... - Erlisch
1
调用者是否为您清零了result[]？如果没有，您应该先这样做！还要注意，在最内层循环中进行水平求和是可怕的。如果您在同一最内层循环中为result[i][j]执行所有数学运算，则只需执行result = hsum(vR)而不是+=。其中hsum是一个水平求和函数，可移植到非MSVC（如果有必要）且比您编写的编译器产生的更好。请参见https://dev59.com/g2w05IYBdhLWcg3w11Qk，我的答案提到了整数hsums。 - Peter Cordes
是的，它被清零了。+= 主要用于基本测试。使用 hadd 进行更新 - 这是你所说的吗？ - Erlisch
糟糕，我正在进行一些测试并将其注释掉了——我想我在一个不适当的时间复制了它。在进一步优化之前，我希望先让它正常工作起来。我本来希望一次计算出4个结果，但现在我意识到这并不像我希望的那样有效。然而，它仍然没有正常工作……我非常怀疑问题出在我如何使用vB上。 - Erlisch
显示剩余2条评论
2个回答
3
你说得对，你的存在问题。 你正在加载4个连续整数，但是不是连续的。 你实际上得到了。
我忘记了哪种方法适用于矩阵乘法（matmul）。 有时候，每次获得4个结果比为每个结果进行水平求和更好。
转置其中一个输入矩阵的成本是O(N^2)。这是值得的，因为它意味着O(N^3)的matmul可以使用顺序访问，并且您当前的循环结构变得友好于SIMD。
甚至有更好的方法，例如在使用之前立即转置小块，以便在再次读取它们时仍然在L1缓存中处于热状态。或者遍历目标行并添加一个结果，而不是累加单个或小组行*列点积的完整结果。 缓存阻塞，也称为循环平铺，是获得良好matmul性能的关键。还请参见每个程序员都应该知道的有关内存的知识？ ，其中附录中有一个缓存阻塞的SIMD FP matmul示例（无需转置）。
已经编写了许多有关优化矩阵乘法的文章，包括使用SIMD和缓存阻塞。我建议你搜索一下。大多数情况下，它可能在讨论FP，但它同样适用于整数。
（除了SSE / AVX仅针对FP具有FMA，而不是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个插入执行gather，而不是通过使用MOVD来打破先前迭代值上的依赖链，因此情况更糟。但是，即使是最好的4个分散元素的收集，与load + PUNPCKLDQ相比也非常糟糕。更不用说这会使您的代码需要SSE4.1。

尽管无论如何都需要SSE4.1来使用_mm_mullo_epi32而不是扩展PMULDQ (_mm_mul_epi32)。

请注意，整数乘法吞吐量通常比FP乘法差，特别是在Haswell及以后的处理器上。 FP FMA单元每个32位元素（用于FP尾数）只有24位宽的乘法器，因此将其用于32x32 => 32位整数需要拆分为两个微操作。

是的，这是一个大错误。另外我也发现SSE中的乘法并不像我想象的那样工作 - 请看我的修改。感谢所有的帮助。 :) - Erlisch
@Erlisch：天啊，现在你在内部循环中有2个PHADDD指令和一个标量+=。 你有没有将其与标量代码进行基准测试？它可能会更慢。 - Peter Cordes
@Erlisch：我刚刚注意到我们可以并行地为j+0..3生成结果，这比使用gather在k上进行向量化要高效得多。我更新了我的答案。 - Peter Cordes
是的，它慢了很多，但正如我所说，我首先制作了一个可工作的示例。你的版本将时间缩短了一半以上。它看起来也很容易改为i-k-j循环顺序，以使向量化循环成为最内层循环。 - Erlisch
@Erlisch：是的，但是你需要在内部循环中加载/添加/存储到&result[i][j]，而不是在寄存器中累积结果。最好从内部循环中并行生成4个结果。您可以通过并行使用4个j值向量来提高缓存利用率，以获得更多的广播加载重用，因此内部循环使用mat2中每个缓存行的所有64B（假设mat2是64B对齐的）。当前版本仅使用每个缓存行的16B，并具有跨距访问模式。 - Peter Cordes
我对同时处理4个向量以获得更多广播负载重用的做法感到似曾相识；我很确定我之前为了解决同样的问题想过这个想法，或者至少是非常相似的。 - Peter Cordes
0

这个第一版是由OP作为问题的编辑发布的，但它并不属于那里。出于纪念的目的，将其移动到社区维基答案中。

那个第一版在性能方面非常糟糕，是向量化的最差方式，内部循环中进行标量hsum，并使用insert_epi32进行手动收集，甚至没有4x4转置。

更新： 哇喔！我终于搞清楚了。除了我的逻辑错误（感谢Peter Cordes的帮助），还有一个问题是_mm_mul_epi32()并不像我想象的那样工作 - 我应该使用_mm_mullo_epi32()代替！

我知道这不是最有效的代码，但它是为了让它正常工作而制作的 - 现在我可以开始优化它了。

（注意，不要使用这个，它非常非常慢）

// editor's note: this is the most naive and horrible way to vectorize
void matmulSSE_inefficient(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]);
                          // less braindead would be to start vB with movd, avoiding a false dep and one shuffle uop.
                          // vB = _mm_cvtsi32_si128(mat2[k][j]);   // but this manual gather is still very bad
                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);                // very slow inside the inner loop
                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可以向上移动一个循环。结果更快。

// this is significantly better but doesn't do any cache-blocking
void matmulSSE_v2(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");
            }
        }
    }
}

这些假设N是向量宽度的倍数

如果不是这种情况，通常更容易仍然将数组存储填充到向量宽度的倍数，这样每行末尾都有填充，你可以使用简单的 j < N; j += 4 循环条件。你将想要跟踪实际的 N 大小，它与行跨度一起是 4 或 8 的倍数而分开存放。

否则，你需要一个类似于 j < N-3; j += 4` 的循环条件，并进行标量清除以结束一行。

或者进行掩蔽操作或保持最后一个完整向量在寄存器中，这样就可以使用可能重叠行末的最终向量来执行 _mm_alignr_epi8，并可能进行矢量存储。这在 AVX 或特别是 AVX512 掩蔽下更容易实现。

你的代码是正确的，但只能用于边长是4的倍数的图像。能否修改代码使其适用于任意尺寸的图像？ - rojan sudev
@rojansudev: 是的，如果它不是4的倍数，则可以对行中最后的N＆3个元素进行标量清理。或者更好的方法是填充您的行。已更新答案并提供了一些提示。 - Peter Cordes
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接