用于平铺矩阵乘法的 AVX 内在函数答案

【问题标题】：AVX intrinsics for tiled matrix multiplication [closed]用于平铺矩阵乘法的 AVX 内在函数
【发布时间】：2020-01-29 08:03:03
【问题描述】：

我试图使用 AVX512 内在函数来矢量化我的矩阵乘法循环（平铺）。我使用 __mm256d 作为变量来存储中间结果并将它们存储在我的结果中。但是，不知何故，这会触发内存损坏。我不知道为什么会这样，因为非 AVX 版本可以正常工作。此外，另一个奇怪的事情是，瓷砖大小现在会以某种方式影响结果。

矩阵结构附在以下代码部分中。该函数采用两个矩阵指针 m1 和 m2 以及一个整数作为 tileSize。感谢@harold 的反馈，我现在已将矩阵 m1 的 _mm256_load_pd 替换为广播。但是，内存损坏问题仍然存在。我还在下面附上了内存损坏的输出


__m256d rResult rm1, rm2, rmult;


    for (int bi = 0; bi < result->row; bi += tileSize) {
         for (int bj = 0; bj < result->col; bj += tileSize) {
             for (int bk = 0; bk < m1->col; bk += tileSize) {
                 for (int i = 0; i < tileSize; i++ ) {
                     for (int j = 0; j < tileSize; j+=4) {
                         rResult = _mm256_setzero_pd();
                         for (int k = 0; k < tileSize; k++) {

                            //  result->val[bi+i][bj+j] += m1.val[bi+i][bk+k]*m2.val[bk+k][bj+j];


                             rm1 = _mm256_broadcast_pd((__m128d const *) &m1->val[bi+i][bk+k]);
                             rm2 = _mm256_load_pd(&m2->val[bk+k][bj+j]);
                             rmult = _mm256_mul_pd(rm1,rm2);
                             rResult = _mm256_add_pd(rResult,rmult);
                             _mm256_store_pd(&result->val[bi+i][bj+j],rResult);
                         }
                     }  
                 }
             }
         }
     }
return result;

*** Error in `./matrix': free(): invalid next size (fast): 0x0000000001880910 ***
======= Backtrace: =========
/lib64/libc.so.6(+0x81609)[0x2b04a26d0609]
./matrix[0x4016cc]
/lib64/libc.so.6(__libc_start_main+0xf5)[0x2b04a2671495]
./matrix[0x400e29]
======= Memory map: ========
00400000-0040c000 r-xp 00000000 00:2c 6981358608                         /home/matrix
0060b000-0060c000 r--p 0000b000 00:2c 6981358608                         /home/matrix
0060c000-0060d000 rw-p 0000c000 00:2c 6981358608                         /home/matrix
01880000-018a1000 rw-p 00000000 00:00 0                                  [heap]
2b04a1f13000-2b04a1f35000 r-xp 00000000 00:16 12900                      /usr/lib64/ld-2.17.so
2b04a1f35000-2b04a1f3a000 rw-p 00000000 00:00 0
2b04a1f4e000-2b04a1f52000 rw-p 00000000 00:00 0
2b04a2134000-2b04a2135000 r--p 00021000 00:16 12900                      /usr/lib64/ld-2.17.so
2b04a2135000-2b04a2136000 rw-p 00022000 00:16 12900                      /usr/lib64/ld-2.17.so
2b04a2136000-2b04a2137000 rw-p 00000000 00:00 0
2b04a2137000-2b04a2238000 r-xp 00000000 00:16 13188                      /usr/lib64/libm-2.17.so
2b04a2238000-2b04a2437000 ---p 00101000 00:16 13188                      /usr/lib64/libm-2.17.so
2b04a2437000-2b04a2438000 r--p 00100000 00:16 13188                      /usr/lib64/libm-2.17.so
2b04a2438000-2b04a2439000 rw-p 00101000 00:16 13188                      /usr/lib64/libm-2.17.so
2b04a2439000-2b04a244e000 r-xp 00000000 00:16 12867                      /usr/lib64/libgcc_s-4.8.5-20150702.so.1
2b04a244e000-2b04a264d000 ---p 00015000 00:16 12867                      /usr/lib64/libgcc_s-4.8.5-20150702.so.1
2b04a264d000-2b04a264e000 r--p 00014000 00:16 12867                      /usr/lib64/libgcc_s-4.8.5-20150702.so.1
2b04a264e000-2b04a264f000 rw-p 00015000 00:16 12867                      /usr/lib64/libgcc_s-4.8.5-20150702.so.1
2b04a264f000-2b04a2811000 r-xp 00000000 00:16 13172                      /usr/lib64/libc-2.17.so
2b04a2811000-2b04a2a11000 ---p 001c2000 00:16 13172                      /usr/lib64/libc-2.17.so
2b04a2a11000-2b04a2a15000 r--p 001c2000 00:16 13172                      /usr/lib64/libc-2.17.so
2b04a2a15000-2b04a2a17000 rw-p 001c6000 00:16 13172                      /usr/lib64/libc-2.17.so
2b04a2a17000-2b04a2a1c000 rw-p 00000000 00:00 0
2b04a2a1c000-2b04a2a1e000 r-xp 00000000 00:16 13184                      /usr/lib64/libdl-2.17.so
2b04a2a1e000-2b04a2c1e000 ---p 00002000 00:16 13184                      /usr/lib64/libdl-2.17.so
2b04a2c1e000-2b04a2c1f000 r--p 00002000 00:16 13184                      /usr/lib64/libdl-2.17.so
2b04a2c1f000-2b04a2c20000 rw-p 00003000 00:16 13184                      /usr/lib64/libdl-2.17.so
2b04a4000000-2b04a4021000 rw-p 00000000 00:00 0
2b04a4021000-2b04a8000000 ---p 00000000 00:00 0
7ffc8448e000-7ffc844b1000 rw-p 00000000 00:00 0                          [stack]
7ffc845ed000-7ffc845ef000 r-xp 00000000 00:00 0                          [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0                  [vsyscall]
Aborted

【问题讨论】：

除了不提供minimal reproducible example之外，您正在计算的内容只是与矩阵乘积模糊相关（例如，您从未将rResult重置为零，但您多次覆盖result->val[bi+i][bj+j]） -- 我建议首先实现一个标量变量进行比较。您的“内存损坏”的原因可能是您的内存未对齐。试试_mm256_loadu_pd 和_mm256_storeu_pd（或者更好，确保你的记忆是对齐的）。
为什么要广播doubles 的对，而不是_mm256_broadcast_sd 或_mm256_set1_pd()？没有多大意义。还有，你看我的回答了吗？我仍然认为它可能解释了内存损坏，除非您的矩阵大小是瓦片宽度的倍数。 Harold 的回答只是谈论数学正确性问题，与数组边界问题分开。
您发布的代码不包含对free() 的任何调用，也就是说，您的（当前）内存损坏的原因在其他地方。提供minimal reproducible example！阅读该链接以了解如何以及为什么！
@chtz 事实上，没有一个免费的功能。我猜内存损坏是指针在某处溢出。
您是否故意无视我的 MRE 请求？您是否知道 C++ 可能会调用 free()，例如通过delete 运算符或容器的析构函数？

标签： c++ matrix optimization intrinsics avx

【解决方案1】：

该代码从 m1 加载一个小行向量，从 m2 加载一个小行向量并将它们相乘，这不是矩阵乘法的工作方式，我假设它是相同标量循环的直接向量化。您可以使用来自 m1 的广播负载，这样与来自 m2 的行向量的乘积会产生一个方便的结果的行向量（相反，从 m2 广播，您会得到一个结果的列向量存储起来很棘手 - 当然，除非您使用以列为主的矩阵布局）。

从不重置rResult也是错误的，使用平铺时要格外小心，因为平铺意味着将单个结果放在一边，然后再重新拾取。实现C += A*B 很方便，因为这样您就不必区分第二次处理结果（从结果矩阵中加载rResult）和第一次处理结果（或者将累加器归零，或者如果你实现C += A*B，那么它也只是将它从结果中加载出来）。

存在一些性能错误，

使用一个蓄能器。这限制了内部循环在长期内每 4 个周期（Skylake）运行一次，因为通过加法（或 FMA）的循环携带依赖。每个周期应该有 2 个 FMA，但这样每 4 个周期就会有一个 FMA，速度为 1/8。
使用 2:1 的负载与 FMA 比率（假设 mul+add 已收缩），它需要为 1:1 或更好，以避免负载吞吐量成为瓶颈。 2:1 的比率限制为半速。

它们的解决方案是在内循环中将来自 m1 的小列向量与来自 m2 的小行向量相乘，求和成一个小的累加器矩阵，而不仅仅是其中一个。例如，如果您使用 3x16 区域（3x4 向量，向量长度为 4，向量对应于来自 m2 的加载，来自 m1 您将执行广播加载），则有 12 个累加器，因此有 12 个独立的依赖链：足以隐藏 FMA 的高延迟吞吐量产品（每个周期 2 个，但在 Skylake 上需要 4 个周期，因此您需要至少 8 个独立链，Haswell 上至少需要 10 个）。这也意味着内循环有7个负载和12个FMA，甚至比1：1更好，它甚至可以支持Turbo频率而无需超频缓存。

我还想指出，在每个维度上设置相同的图块大小不一定是最好的。也许是，但可能不是，维度的作用确实有点不同。

更高级的性能问题，

不重新包装瓷砖。这意味着图块将跨越不必要的页面，这会损害 TLB 的有效性。您很容易遇到这样的情况：您的切片适合缓存，但分布在太多页面上以适应 TLB。 TLB 抖动不好。

使用非对称图块大小，您可以安排 m1 分块或 m2 分块以对 TLB 友好，但不能同时安排两者。

【讨论】：

大声笑，我认为 OP 的代码看起来非常简单。我想知道他们的输入是否已经转置，但我只是发布了一个关于循环边界的快速答案，并没有认真研究它实际上在做什么。
@harold 感谢您的详细解释。我已经阅读了有关负载的更多信息，现在将 _mm256_load 替换为矩阵 m1。事实上，我还在考虑矩阵乘法中的标量加载和存储。我已经编辑了我的代码部分，请您再看一下并给我更多反馈吗？我遇到的一个问题是，随着瓦片大小的变化，我无法获得一致的结果，它会导致矩阵 free() 函数中的内存损坏，从而导致程序中止。

【解决方案2】：

如果您关心性能，通常您需要一块连续的内存，而不是一组指向行的指针。

无论如何，如果您的图块大小不是每个向量的 4 个双精度数的倍数，您可能正在读取行尾。或者，如果您的行或列不是 tile 大小的倍数，那么您需要在最后一个 full tile 之后停止，并为末尾编写清理代码。

例如bi < result->row - (tileSize-1) 用于外部循环

如果您的图块大小不是 4 的倍数，那么您还需要 i < tileSize-3。但希望您正在执行 2 次幂循环平铺/缓存阻塞。但是您需要一个size - 3 边界用于部分平铺中的矢量清理。然后可能对最后几个元素进行标量清理。（或者，如果您可以使用在行末尾结束的未对齐最终向量，则可以使用屏蔽加载/存储。但对于 matmul 比只进行一次传递的算法更棘手。）

【讨论】：

谢谢！尽管在我为外循环设置新边界后它仍然会导致内存损坏。我认为边界是正确的，因为它们适用于我的标量版本。
@HerrGünther：是的，向量化时很容易忘记修复循环边界。但要理解为什么，请记住向量化是一种展开：每次迭代无条件地执行 4 或 8 个或任何数组元素。您展开 4，然后使用一个 SIMD 向量进行这 4 次迭代。可以安全地开始展开迭代的最高 i 值是 i+3 是您应该触摸的最后一个元素。无论如何，如果这仍然不能修复错误，请仔细检查您的其他循环（例如列循环）以确保您没有错过任何一个。还有像 valgrind 这样的工具。