一種高效的矩陣乘法實現

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

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。訪問內存卻要上百個時鍾周期，孰快孰慢，很清楚了吧？

 

程序摘自CSAPP 2nd edition