C中實現矩陣乘法的一種高效的方法

如何計算矩陣乘法，這個大家都知道。通常情況下，我們都是用以下代碼實現的：

複製代碼 代碼如下:for(i=0;i<n;++i)
for(j=0;j<n;++j){
sum=0;
for(k=0;k<n;++k)
sum+=A[i][k]*B[k][j];
C[i][j]+=sum;
}

但是考慮了快取的問題後，其實有一種更好的實現方式：複製代碼 代碼如下:

for(i=0;i<n;++i)
for(k=0;k<n;++k){
r=A[i][k];
for(j=0;j<n;++j)
C[i][j]+=r*B[k][j];
}

細看一番就會發現這兩種實現語義是等價的，但是後者的實際運行效率卻比前者高。

那為什麼會如此呢？

那是因為CPU讀資料時，並不是直接存取記憶體，而是先查看緩衝中是否有資料，有的話直接從緩衝讀取。而從緩衝讀取資料比從記憶體讀資料快很多。

當資料不在緩衝中時，CPU會將包含資料在內的一個資料區塊讀到緩衝，如果程式具有良好空間局部性，那麼第一次cache miss後，之後的幾次資料訪問就可以直接在緩衝中完成。除了空間局部性（程式傾向於引用與當前資料鄰近的資料）之外，還有時間局部性（程式傾向於引用最近被引用過的資料）。

回到矩陣乘法。（我們只考慮內迴圈）

前者對矩陣A，有良好的空間局部性，假設一次能緩衝四個元素，則每次迭代對於A只有0.25次miss，但是對於B，則不然，因此B是按列訪問的，每次訪問都會miss，因此每次迭代總的miss數是1.25。

後者對於矩陣C和矩陣B都有良好的局部性，每次迭代都只有0.25詞miss，因此總的miss數是0.5。後者每次迭代多了一次儲存（對C[i][j]寫入），但是即便如此，後者的運行效率也比前者高。

總而言之，要想程式跑得快，就要在程式中多利用局部性，讓緩衝hold住你的資料，減少訪存次數。要知道CPU可以在3個刻度內訪問到L1 cache，10個刻度左右的時間訪問到L2 cache。訪問記憶體卻要上百個刻度，孰快孰慢，很清楚了吧？