如何进行正确的缓存阻塞矩阵转置？

Question

我正在尝试在 C 中执行缓存阻塞矩阵转置，但我在代码中遇到了一些问题。我的猜测是它与索引有关。你能告诉我哪里出错了吗？

我正在考虑我在网上找到的这两种算法：http://users.cecs.anu.edu.au/~Alistair.Rendell/papers/coa.pdf 和 http://iosrjen.org/Papers/vol3_issue11%20(part-4)/I031145055.pdf

但我还不知道如何正确编码。

for (i = 0; i < N; i += block) {
    for (j = 0; j < i; j += block ) {

        for(ii=i;ii<i+block;ii++){
            for(jj=j;jj<j+block;jj++){

                temp1[ii][jj] = A2[ii][jj];
                temp2[ii][jj] = A2[jj][ii];

                A2[ii][jj] = temp1[ii][jj];
                A2[ii][jj] = temp2[ii][jj];
            }     
        }                                   
    }                       
}   

temp1 和 temp2 是两个大小为块 x 块的矩阵，用零填充。 当我将值返回到A2（转置矩阵前后）时，我不确定是否需要再做一次for。

我也试过这个：

for (i = 0; i < N; i += block) {
    for (j = 0; j < N; j += block ) {
    ii = A2[i][j];
    jj = A2[j][i];

        A2[j][i] = ii;
    A2[i][j] = jj;
    }                       
}

我期待比“天真的”矩阵转置算法有更好的性能：

for (i = 1; i < N; i++) {
    for(j = 0; j < i; j++) {

        TEMP= A[i][j];
        A[i][j]=A[j][i];
        A[j][i]=TEMP;

    }
}

Answer 1

执行阻塞矩阵转置的正确方法不是程序中的内容。额外的 temp1 和 temp2 数组将无用地填充您的缓存。而且您的第二个版本不正确。更多你做了太多操作：元素被转置两次，对角线元素被“转置”。

但首先我们可以做一些简单的（和近似的）缓存行为分析。我假设您有一个双精度矩阵，并且缓存行是 64 字节（8 个双精度）。

如果缓存可以完全包含矩阵，则阻塞实现等效于幼稚实现。您只有强制缓存未命中来获取矩阵元素。处理 N×N 个元素的缓存未命中数将为 N×N/8，平均每个元素的未命中数为 1/8。

现在，对于简单的实现，看看你在缓存中处理了 1 行之后的情况。假设您的缓存足够大，您的缓存中将有：
* 完整的行 A[0][i]
* 矩阵 A[i][0..7] 每隔一行的前 8 个元素

这意味着，如果您的缓存足够大，您可以处理 7 个连续的行，而不会出现任何缓存未命中的情况，除了获取行的缓存。因此，如果您的矩阵是 N×N，如果缓存大小大于 ~2×N×8，则只有 8×N/8(lines)+N(cols)=2N 缓存未命中来处理 8×N 个元素，每个元素的平均未命中数为 1/4。从数值上看，如果 L1 缓存大小为 32k，则在 N

如果您有一个非常大的矩阵，则在第一行结束后，第二行的开头已从缓存中弹出。如果矩阵的一行完全填满缓存，就会发生这种情况。在这种情况下，缓存未命中的数量将更为重要。每行将有 N/8（获取行）+ N（获取列的第一个元素）缓存未命中，每个元素平均有 (9×N/8)/N& 大约 1 个未命中。

因此，您可以通过阻塞实现获得收益，但仅适用于大型矩阵。

这是矩阵转置的正确实现。它避免了元素 A[l][m] 的双重处理（当 i=l 和 j=m 或 i=m 和 j=l 时），不转置对角线元素并使用寄存器进行转置。

朴素版

for (i=0;i<N;i++)
  for (j=i+1;j<N;j++)
    {
       temp=A[i][j];
       A[i][j]=A[j][i];
       A[j][i]=temp;
    }

和阻塞版本（我们假设矩阵大小是块大小的倍数）

for (ii=0;ii<N;ii+=block)
  for (jj=0;jj<N;jj+=block)
    for (i=ii;i<ii+block;i++)
      for (j=jj+i+1;j<jj+block;j++)
        {
           temp=A[i][j];
           A[i][j]=A[j][i];
           A[j][i]=temp;
        }

Answer 2

我正在使用您的代码，但是当我将幼稚算法与阻塞算法进行比较时，我没有得到相同的答案。我把这个矩阵 A 和我得到的矩阵 At 如下：

一个

2 8 1 8
6 8 2 4
7 2 6 5
6 8 6 5

在

2 6 1 6
8 8 2 4
7 2 6 5
8 8 6 5

使用大小为 N=4 且块数为 2 的矩阵