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

Question

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

我正在考虑我在网上找到的这两种算法：http://users.cecs.anu.edu.au/~Alistair.Rendell/papers/coa.pdf and 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;
    }                       
}

我期待比 "naive" 矩阵转置算法更好的性能：

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 数组将无用地填充您的缓存。你的第二个版本不正确。更多操作太多：元素转置两次，对角线元素为 "tranposed".

但首先我们可以做一些简单的（和近似的）缓存行为分析。我假设您有一个双精度矩阵并且缓存行是 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(lines)+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[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 的矩阵

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

How to do a proper Cache Blocked Matrix Transposition?

c

caching

transpose

block

matrix