c++ - 用于平铺矩阵乘法的 AVX 内在函数

Question

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

矩阵结构附在以下代码部分中。该函数需要两个矩阵指针 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

Answer 1

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

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

存在一些性能错误，

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

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

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

更高级的性能问题，

• 不重新包装瓷砖。这意味着图 block 将跨越必要的更多页面，这会损害 TLB 的有效性。您很容易遇到这样的情况:您的切片适合缓存，但分布在太多页面上，无法适合 TLB。 TLB 颠簸不好。

使用不对称的切片大小，您可以安排 m1 切片或 m2 切片以支持 TLB，但不能同时安排两者。