程序性能优化探讨（6）——矩阵乘法优化之分块矩阵
来源：互联网 发布：mac放电影有杂音 编辑：程序博客网 时间：2024/05/16 17:11
有一种性格叫做偏执，有一种矩阵优化运算叫做分块。实话说，也许我这辈子也用不上这种随牛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（我只分析了前一个），结果如下：
接下来是统计图：
从结果看，b_rck，b_rkc的性能，确实随着数组元素个数的增长而区域性能的稳定，but，这里居然发现，rkc的版本也基本围绕b_rck，b_rkc来走，而且周期数还有超越这两个分块版本的可能！这个和教科书上的完全不一样啊！呵呵，不过，我也不知道这是为什么o(╯□╰)o
- 程序性能优化探讨（6）——矩阵乘法优化之分块矩阵
- 程序性能优化探讨（5）——高速缓存、存储器山与矩阵乘法优化
- 矩阵乘法——CUDA 优化记录
- 程序碎片- 矩阵乘法优化(dp,递归)
- 程序碎片- 矩阵乘法优化(dp,循环)
- 矩阵乘法的优化
- CUDA 矩阵乘法优化
- 矩阵乘法的优化
- CUDA: 矩阵乘法优化
- 矩阵乘法的优化
- 矩阵乘法优化
- CUDA: 矩阵乘法优化
- 矩阵乘法cache优化
- 矩阵乘法的优化
- 矩阵乘法并行优化
- 矩阵乘法优化算法
- 矩阵乘法优化DP
- 矩阵乘法优化DP
- 【Foundation】17-字符串NSString和NSMutableString
- Arthur and Questions - CodeForces 518 E
- 刚无意间看到一个网站，页面挺简洁，就是网站里面没啥内容
- JAVASE--集合、I/O流
- java学习之File类--2015-2-26
- 程序性能优化探讨（6）——矩阵乘法优化之分块矩阵
- 谈谈微信红包海量运营--发10亿个红包难在哪里？
- CreateProcess的前两个参数究竟怎么用
- 线性回归概率解释(Linear Regression)
- SDUT 3145 Integer division 1
- C++学习起步
- pat1021 Deepest Root
- 【HDU 1512】Monkey King
- C语言操作Excel表格
