程序性能优化探讨（6）——矩阵乘法优化之分块矩阵

        有一种性格叫做偏执，有一种矩阵优化运算叫做分块。实话说，也许我这辈子也用不上这种随牛B但很复杂的算法，有些版本的教材直接删除这个内容。但越是这样我越想不过，因此借写这篇博客，把分块矩阵乘法彻底分析清楚。

         把矩阵乘法进行分块优化，很奇妙的算法，假设我们要做11X11矩阵乘法，AXB = C原本是这样算的：

        比如我们要计算C1.2，就要把上图的两条射线，交叉乘加。把A1.0~A1.10和B0.2~B10.2遍历完，也就是说要调用121个元素，还要进行11次乘法和10次加法运算，才能得出答案！如果按照上图中分出四个粗线块来进行分批交叉乘加，就能提高效率。

        关键是，怎么理解这种算法的转变呢？如何证明两种算法是等价的呢？这里我肯定没办法做完整的论证。不过我可以简单做个定性的判断。既然矩阵乘法的本质就是交叉乘加，既然最外层都是加法，那么根据加法结合律，把射线分段进行乘加，如果射线拼起来刚好是原来整个射线，那么两种算法等价了。

                  （1）                    （2）                （3）

        

        像上图这样，把计算C1.2这个元素的大的交叉乘加，分成了三个段组。而什么时候计算哪组射线，这个顺序是没有影响的，只要小的交叉组配对正确（比如第一个横条对应第一个竖条），再把交叉乘加的结果都加到C1.2的临时值中，就能正确得出C1.2的最终结论。

        好了，下面把分块矩阵乘法的完整代码一段段贴出来：

vord brck(array A, array B, array C, int n, int bsize) 
{
    int r, c, k, kk, cc;
    double sum;
    int en = bsize * (n/bsize); /* Amount that frts evenly into blocks */
  
    for (r = 0; r < n; r++)
        for (c = 0; c < n; c++)
            C[r][c] = 0.0;




    for (kk = 0; kk < en; kk += bsize) { 
        for (cc = 0; cc < en; cc += bsize) {
            for (r = 0; r < n; r++) {
                for (c = cc; c < cc + bsize; c++) {
                    sum = C[r][c];
                    for (k = kk; k < kk + bsize; k++) {
                         sum += A[r][k]*B[k][c];
                    }
                    C[r][c] = sum;
                }
            }
        }

        好，第一段代码截止在这里。咋一看很复杂，我们来一条条分析。以我的习惯是从多层循环的最里层进行分析。

    sum = C[r][c];
    for (k = kk; k < kk + bsrze; k++) {
        sum += A[r][k]*B[k][c];
    }
    C[r][c] = sum;

        我们先分析kk=0的情况，而bsize是块的大小4，通过这个，把循环限制在某个分块里。代码中的C[r][c]存储的是Cr.c临时值，通过sum不断增加C[r][c]的累积加法结果。

        假设kk = 0，r = 1、c = 2。那么代码执行的就是下一段：

    for (k = 0; k < bsize; k++) {
        sum += A[1][k]*B[k][2];
    }

        这就相当于在做第一条红横线和第一条蓝竖线的交叉乘加！能理解到这里，那么我们就能理解更多的代码了。

        假设kk = 0，cc = 0，r = 0，那么以下代码：

for (c = 0; c <  bsize; c++) {
    sum = C[r][c];
    for (k = 0; k < bsize; k++) {
        sum += A[r][k]*B[k][c];
    }
    C[r][c] = sum;
}

        

        这相当于在做以下的交叉乘加遍历：

        

        

        如果把外层循环for (r = 0; r < n; r++) {}加上，那么他们的循环就是：

        

        我们看到r把所有行遍历完，但由于k被限制在bsize范围内，因此每行只处理四个列。

        好了，接下来剩下两个外层循环：

    for (kk = 0; kk < en; kk += bsize) { 
        for (cc = 0; cc < en; cc += bsize) {

        首先是cc，cc被限制在en里，而en的计算是限制在四个小分块中，好，那么我们先分析cc。我们看到，cc影响的是B[k][c]，因此能把蓝横竖线遍历的范围扩大：

        

        我们看到，由于循环参数的限制，for (cc = 0; cc < en; cc += bsize) 的循环只增加了四条小射线，真正扩大的因素在 for (kk = 0; kk < en; kk += bsize) 。关于这个kk的变化，会同时影响到A[r][k]和B[k][c]，而由于en本身是bsize的整数倍结果，那么kk的递增的扩展结果：

        

        好了，我们发现，第一段代码能遍历的范围似乎并不全，这时候我发现，在外层循环for (kk = 0; kk < en; kk += bsize)里，还有一段代码忘分析了：

/* $end bmm-rck */
/* Now finish off rest of c values  (not shown in textbook) */
    for (r = 0; r < n; r++) {
       for (c = en; c < n; c++) {
           sum = C[r][c];
           for (k = kk; k < kk + bsize; k++) {
                 sum += A[r][k]*B[k][c];
           }
           C[r][c] = sum;
       }
    }
/* $begin bmm-rck */
}

        我们到，这段代码的for (r = 0; r < n; r++) 和A[r][k]为遍历A的行提供了条件，而for (c = en; c < n; c++)和B[k][c]就是遍历A的列。

        这段代码相当于遍历矩阵B中分块en以外的那部分列，要是单独拿来分析，当kk = 0时的遍历范围就应该是：

        

        再加上外层循环for (kk = 0; kk < en; kk += bsize)，遍历范围应该是：

        

        经过对循环for (kk = 0; kk < en; kk += bsize)内的两段代码分析，我们找出了他们实现的交叉乘加遍历范围：

        

        剩余未纳入计算的区域，矩阵A和矩阵B中一目了然，接下来分析第三段代码：

    /* $end bmm-rck */
    /* Now finish remaining k values (not shown in textbook) */ 
    for (cc = 0; cc < en; cc += bsize) {
        for (r = 0; r < n; r++) {
            for (c = cc; c < cc + bsize; c++) {
                sum = C[r][c];
                for (k = en; k < n; k++) {
                    sum += A[r][k]*B[k][c];
                }
                C[r][c] = sum;
            }
        }
    }

        很明显，这段代码和第一段代码的唯一区别就是for (k = en; k < n; k++) ，因此交叉乘加的遍历范围很容易得出：

        

        好了，就只剩下B的最后9个元素没有被交叉乘加过，剩下的这点代码，和第三段代码唯一区别就是for (c = en; c < n; c++)：

    /* Now finish off rest of c values (not shown in textbook) */
    for (r = 0; r < n; r++) {
        for (c = en; c < n; c++) {
           sum = C[r][c];
           for (k = en; k < n; k++) {
               sum += A[r][k]*B[k][c];
           }
           C[r][c] = sum;
        }
   }
    /* $begin bmm-rck */
}
/* $end bmm-rck */

        那么这段代码实现的交叉乘加遍历范围是：

        

        至此，我们终于将所有交叉乘加遍历干净！

        

        费了那么多功夫，有个问题是，干嘛我们要用这种分块来实现算法优化呢？主要考虑在于，当数组大小增加时，时间局部性会明显降低，高速缓存中不命中数目增加。当使用分块技术是，时间局部性由分块大小来决定，而不是数组总大小来决定。另外，分块虽然能明显提高效率，但使得代码灰常难懂，因此一般用于编译器或者频繁执行的库例程……所以说嘛，反正我这辈子是用不上的……

         那么我们就具体来实践一下，把之前的六个版本，和分块的两个版本b_rck、b_rkc（我只分析了前一个），结果如下：

	rck	crk	ckr	kcr	krc	rkc	b_rck	b_rkc
25	0.12	0.02	0.01	0	0	0.31	12.38	14.48
50	0.09	0.01	0.01	0.1	0	0.02	0.01	0.07
75	0.02	0.01	0.02	0.05	0.01	0	0	0.15
100	0.07	0.07	3.71	1.78	4.48	0.11	1.37	1.6
125	0.67	1.88	6.13	4.23	0.67	0.66	0.66	2.25
150	1.06	1.44	9.8	8.16	1.77	1.79	1.78	2.47
175	1.66	1.66	13.05