我编写了以下 C 函数，用于使用平铺/分块和 AVX vector 将两个 NxN 矩阵相乘，以加快计算速度。现在，当我尝试将 AVX 内在函数与平铺结合起来时，我遇到了段错误。知道为什么会发生这种情况吗？
此外，矩阵 B 是否有更好的内存访问模式？也许先转置它，甚至改变 k 和 j 循环？因为现在，我正在按列遍历它，这在空间局部性和缓存行方面可能不是很有效。
1 void mmult(double A[SIZE_M][SIZE_N], double B[SIZE_N][SIZE_K], double C[SIZE_M][SIZE_K])
2 {
3 int i, j, k, i0, j0, k0;
4 // double sum;
5 __m256d sum;
6 for(i0 = 0; i0 < SIZE_M; i0 += BLOCKSIZE) {
7 for(k0 = 0; k0 < SIZE_N; k0 += BLOCKSIZE) {
8 for(j0 = 0; j0 < SIZE_K; j0 += BLOCKSIZE) {
9 for (i = i0; i < MIN(i0+BLOCKSIZE, SIZE_M); i++) {
10 for (j = j0; j < MIN(j0+BLOCKSIZE, SIZE_K); j++) {
11 // sum = C[i][j];
12 sum = _mm256_load_pd(&C[i][j]);
13 for (k = k0; k < MIN(k0+BLOCKSIZE, SIZE_N); k++) {
14 // sum += A[i][k] * B[k][j];
15 sum = _mm256_add_pd(sum, _mm256_mul_pd(_mm256_load_pd(&A[i][k]), _mm256_broadcast_sd(&B[k][j])));
16 }
17 // C[i][j] = sum;
18 _mm256_store_pd(&C[i][j], sum);
19 }
20 }
21 }
22 }
23 }
24 }
最佳答案
_mm256_load_pd 是需要对齐的加载，但您仅在最内层循环中按 k++ 步进，而不是 k+=4加载 4 个 double 的 32 字节 vector 。因此它会出现故障，因为每 4 个负载中有 3 个未对齐。
您不想进行重叠加载，您真正的错误是索引；如果您的输入指针是 32 字节对齐的，您应该能够继续使用 _mm256_load_pd 而不是 _mm256_loadu_pd。因此，使用 _mm256_load_pd 成功捕获了您的错误，而不是工作但给出了数字错误的结果。
向量化四个行*列点积的策略(以生成C[i][j+0..3] vector ) strong> 应该从 4 个不同的列加载 4 个连续的 double (B[k][j+0..3] 通过来自 B[k][j] 的 vector 加载)，并从 A[i][k] 广播 1 个 double 值。请记住，您需要并行 4 个点积。
另一种策略可能涉及到最后的水平总和到标量C[i][j] +=horizontal_add(__m256d)，但我认为这需要首先转置一个输入，以便两者对于一个点积，行 vector 和列 vector 位于连续内存中。但随后您需要在每个内循环末尾进行洗牌以获得水平总和。
您可能还希望使用至少 2 个 sum 变量，以便可以一次读取整个缓存行，并隐藏内部循环中的 FMA 延迟，并有望成为吞吐量瓶颈。 或者最好并行处理 4 或 8 个 vector 。因此，您可以将 C[i][j+0..15] 生成为 sum0, sum1、sum2、sum3。 (或者使用 __m256d 数组；编译器通常会完全展开 8 的循环并将数组优化到寄存器中。)
我认为你只需要 5 个嵌套循环来阻止行和列。虽然显然 6 个嵌套循环是一个有效的选项:请参阅 loop tiling/blocking for large dense matrix multiplication问题中有 5 个嵌套循环，但答案有 6 个嵌套循环。 (不过只是标量，而不是矢量化)。
除了这里的行*列点积策略之外，可能还有其他错误，我不确定。
如果您使用 AVX，您可能还想使用 FMA，除非您需要在 Sandbybridge/Ivybridge 和 AMD Bulldozer 上运行。 (打桩机及更高版本有 FMA3)。
其他 matmul 策略包括添加到内部循环内的目标中，以便您在内循环内加载 C 和 A，并从 B 加载举起。 (或者 B 和 A 交换了，我忘了。) What Every Programmer Should Know About Memory?附录中有一个针对 SSE2 __m128d vector 的矢量化缓存阻止示例，该示例以这种方式工作。 https://www.akkadia.org/drepper/cpumemory.pdf
关于c - 使用 AVX 的平铺矩阵乘法，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/59009628/