我已经编写了以下C函数，用于使用平铺/分块和AVX向量将两个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倍的32字节向量。因此出现故障，因为每4个负载中有3个未对齐。
您不想做重叠的加载，您真正的错误是索引；如果输入指针是32字节对齐的，则应该可以继续使用_mm256_load_pd而不是_mm256_loadu_pd。因此，使用_mm256_load_pd成功捕获了您的错误，而不是正常工作，但给出了错误的数字结果。
您的矢量化四个row*column点乘积（以生成C[i][j+0..3]矢量的策略）应该从4个不同的列中加载4个连续的double（通过B[k][j+0..3]的矢量加载来加载B[k][j]），并进行广播A[i][k]的1倍。请记住，您需要并行4个点积。
另一种策略可能涉及到最后将水平和降低到标量C[i][j] += horizontal_add(__m256d)，但是我认为这需要先移置一个输入，以便行和列向量都在一个点积的连续存储器中。但是然后，您需要在每个内部循环的末尾对水平和进行混洗。
您可能还希望至少使用2个sum变量，以便您可以一次读取整个高速缓存行，并在内部循环中隐藏FMA延迟，并希望对吞吐量造成瓶颈。 或者最好并行处理4或8个向量。因此，将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上运行。 （Piledriver及更高版本具有FMA3）。
[其他matmul策略包括在内部循环中添加到目标，以便在内部循环中加载C和A，并从B上加载一个负载。 （或者我将B和A互换了。）What Every Programmer Should Know About Memory?有一个矢量化的缓存块示例，该示例在附录中以这种方式工作，适用于SSE2 __m128d矢量。 https://www.akkadia.org/drepper/cpumemory.pdf