CS:APP3e, Bryant and O'Hallaron 可以参考这里
void bijk(array A, array B, array C, int n, int bsize) {int i, j, k, kk, jj;double sum;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;}}}}
}
先来谈一下参考资料内bijk函数中的blocking技术吧,和题目的两层嵌套循环不同, bijk函数是五层的嵌套 作为人类似乎很难去理解为啥我就处理个矩阵,要整它五层嵌套,完了它还对性能有好处.
你想象一下transpose函数中如果dim=9999999999999...时,这世界不会存在一个cache能存下这个数组,假设此时cache就只有bsize*bsize(bsize<dim)大小,写完一列bsize个dst后就开始写下一列dst这样只有第一列是不命中的,其他bsize-1列都是命中的.
for (k = kk; k < kk + bsize; k++)
bijk函数中 k<kk+bsize就是控制程序写完一个bsize 后就开始写下一列.
这就是blocking技术的核心了.
我们现在开始改transpose函数:
#include <stdio.h>void transpose(int *dst, int *src, int n, int bsize) { // n为数组大小(假设是方阵的边长), bsize为块大小, bsize宜接近高速缓存大小 int i, j, kk, jj; // 处理能够完整被块大小分割的部分 for (kk = 0; kk < n; kk += bsize) { // 注意这里应该使用n而不是en for (jj = 0; jj < n; jj += bsize) { // 同上 for (i = kk; i < kk + (kk + bsize < n ? bsize : n - kk); i++) { // 确保不越界 for (j = jj; j < jj + (jj + bsize < n ? bsize : n - jj); j++) { // 确保不越界 // 计算一维数组中的索引 int src_index = i * n + j; int dst_index = j * n + i; dst[dst_index] = src[src_index]; // 复制元素} } } }
}int main() {int dim=500;int src[dim][dim];int dst[dim][dim];int i, j;
//给数组赋值for (i = 0; i < dim; i++)for (j = 0; j < dim; j++)src[i][j] = i+j;
//转置transpose(dst,src,dim,500);
//检查转置后的结果for (i = 0; i < dim; i++)for (j = 0; j < dim; j++){if(src[i][j]!=dst[j][i])printf("转置出错\n");}return 0;
}