原文 https://csapp.cs.cmu.edu/public/waside/waside-blocking.pdf

介绍

有一个有趣的技术叫 Blocking(分块)，它可以改善内循环的时间局部性。blocking 的一般思想是将程序中的数据结构组织成大的片(chunks)，这种 chunks 被称作块(blocks)。在本文中，block 指的是应用程序级别的数据块，不是 cache block(缓存块)。程序的结构是这样的：它将一个 chunk 加载到 L1 缓存中，对该 chunk 执行所需的所有读写操作，然后丢弃该 chunk，加载下一个 chunk，依次类推。

与改善空间局部性的简单循环转换不同，blocking 使代码更难阅读和理解。尽管如此，研究和理解该技术还是很有趣的，因为它是一个通用概念，可以在某些系统上产生很大的性能提升。

矩阵乘法的块版本

矩阵乘法分块程序的原理是先将原始的矩阵分成多个子矩阵，然后利用数学原理，将子矩阵的运算变成标量运算。例如，我们要计算 C = AB，这里 A, B, C 是 8 x 8 矩阵。然后，我们可以将每个矩阵划分为 4 x 4 的子矩阵：

void bijk(array A, array B, array C, int n, int bsize)
{
    int i, j, k, kk, jj;
    double sum;
    /* Amount that fits evently into blocks */
    int en = bsize * (n/bsize);    
    
    for(i = 0; i < n; i++)
        for(j = 0; j < n; j++)
            C[i][j] = 0.0;
    
    for(kk = 0; kk < en; kk += bsize) {
        for(jj = 0; jj < en; jj += bsize) {
            for(i = 0; i < n; i++) {
                for(j = jj; j < jj + bsize; j++) {
                    sum = C[i][j];
                    for(k = kk; k < kk + bsize; k++) {
                        sum += A[i][k]*B[k][j];
                    }
                    C[i][j] = sum;
                }
            }
        }
    }
}

array size (n)，block size(bsize)。

当前这个版本的解法，我们称之为 b_ijk 版本。这段代码背后的基本思想是将 A 和 C 划分为 1 x bszie 的行片段(row slivers)，并将 B 划分为bsize x bszie 的块。最内层的 (j,k) 循环对将 A 的一个 1 x bszie 的行片段乘以 B 的一个 bsize x bszie 的块，并累积成 C 的一个 1 x bszie 的行片段。i 循环遍历 A 和 C 的 n 个行片段，使用 B 中相同块。

关键思想是，它将 B 中的一个块加载到缓存中，用完它，然后丢弃它。对 A 的引用具有良好的空间局部性，因为每个行片段迭代的步长为 1，同时也有很好的时间局部性，因为这个行片段被连续引用 bsize次。对 B 的引用具有良好的时间局部性，因为整个 bsize x bszie 块被连续访问 n 次。最后，对 C 的引用具有良好的空间局部性，因为行片段中的每个元素都是连续写入的。注意，对C 的引用没有很好的时间局部性，因为每个块只被访问一次。

Blocking 会使代码更难阅读，但它也能带来很大的性能红利。上图展示了两个版本的分块矩阵相乘在 Pentium III Xeon 系统（bsize=25) 上的性能。请注意，与最佳的非分块版本相比，分块将运行时间提高了两倍，从每次迭代大约 20 个周期降低到每次迭代大约 10 个周期。关于分块的另一个有趣现象是，每次迭代时间随着数组大小的增加几乎保持不变。对于较小的数组，分块版本中的额外开销导致其运行速度比未分块版本慢。在 n=100 处有一个交叉点，在这个交叉点之后，分块版本运行得更快。

注意，分块矩阵乘法并不能提供所有系统的性能。例如，在 Core i7 系统上，存在未分块版本的矩阵乘法，其性能与最佳分块版本相同的情况。