如何进行适当的缓存阻塞矩阵转换？

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块的两个矩阵。当我将值返回到for（转置矩阵之前和之后）时，我不确定是否必须再做另一个A2。

我也试过这个：

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;
    }                       
}

我期待比“天真”Matrix Transposition算法更好的性能：

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 / 8以处理N×N个元素，每个元素的平均未命中数为1/8。

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

这意味着，如果缓存足够大，则可以处理7个连续的行，而不会有任何缓存未命中，而不是获取行。因此，如果你的矩阵是N×N，如果缓存大小大于~2×N×8，你将只有8×N / 8（行）+ N（cols）= 2N缓存未命中来处理8×N个元素，每个元素的平均未命中数是1/4。在数值上，如果L1高速缓存大小为32k，则如果N <2k则会发生这种情况。如果L2缓存为256k，如果N <16k，数据将保留在缓存L2中。由于现代处理器中的预取非常有效，我认为L1中的数据与L2中的数据之间的区别不会真正可见。

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

因此，您可以获得阻止的实现，但仅适用于大型矩阵。

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

天真的版本

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

我正在使用你的代码，但是当我将naive与阻塞算法进行比较时，我得不到同样的答案。我把这个矩阵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的矩阵