這一節內容將綜合（3）和（4），討論快取記憶體相關的程式優化。

一、牛B完了的儲存器山

        一個程式從儲存系統中讀資料的速率被稱為讀吞吐量或讀頻寬。如果一個程式在s秒的時間段內讀n個位元組，那麼讀吞吐量就是n/s，一般用MB/s作為單位。

        第（4）節中討論的時間區域性性和空間區域性性。從快取塊大小來權衡，時間區域性性和空間區域性性似乎剛好成反比，但當我們全面的討論整個儲存層次效能時會發現，其他兩種區域性性可能各自有各自的規律。

        如果我們把一次性讀取的資料量的大小稱為工作集，把上一次和下一次讀取資料的之間儲存距離稱為步長，如果我想獲得工作集和步長這兩種因素影響下程式的讀吞吐量，會是什麼結果呢？教材中給出實現程式碼，但我不打算在這裡全部貼出來，因為詳細的程式碼分析涉及到教材後面章節裡的K次最優測量方法←_←內啥，我自己還沒看懂，準確的說是暫時沒耐著性子看下去，好像天書，等哪天有心情看懂了，再補上程式碼實現也不遲，那我們就先把儲存器山的圖貼出來吧：

        我們看到，儲存器山由三個座標確定點位，底面右邊的座標是工作集大小（wording set size），按位元組計數；底面左邊的座標是步長（stride），按字計數。這裡可以把字理解成4位元組；垂直的座標就是讀吞吐量，根據吞吐量的大小判斷程式效能的優劣。注意到這個儲存器山針對的CPU型號已經寫在上面了，這個CPU裡面分別有大小2*32kb、塊大小64位元組的快取記憶體L1，指令和資料快取都是這麼大。另外還有6M大小的快取記憶體L2。

        首先把步長作為常數，來看工作集32k大小時的情況。一次性處理32k位元組的資料，由於L1有32k的指令快取和32k資料快取，因此，32k工作集的資料是完全可以存入L1快取的，因此可以看到此時它的吞吐量處於最優狀態。奇怪的是，當工作集座標往右減小時，山峰沒有繼續升高，而是急速降低，你可以明顯的看到，從32k之後的下降陡坡，這是為什麼？其實，這裡需要了解測試程式的結構：

#define MINBYTES (1 << 12)  /* 最小工作集 ranges from 4 KB */

#define MAXBYTES (1 << 25)  /*最大工作集 ... up to 32 MB */
#define MAXSTRIDE 30           /*最大步長*/

    for (size = MAXBYTES; size >= MINBYTES; size >>= 1) {
for (stride = 1; stride <= MAXSTRIDE; stride++) {
   printf("%.1f\t", run(size, stride, Mhz));
}

        上面這個是生成儲存器山的核心迴圈，第一層迴圈是遍歷工作集大小，從32MB開始兩倍變小，直到4kB為止；第二層迴圈是步進，從1到30字（4位元組）。有了這兩個測試關鍵引數傳入儲存器山生成函式run裡，就能將這座山遍歷出來。

        好了，當第一層size取值小於32k，而步長又大於20時，由於工作集size太小，因此每一個工作集處理時間會非常短，但正因為工作集太小，而步長又太大，資料集越不夠連續，命中的情況就越遭，因此懲罰也就越大，主要的時間就消耗在不命中處罰和迴圈本身的開銷上，因此山高峰背後那個明顯下降，其實並不能反映出L1快取記憶體的真實效能。

        我們再來看看這座山，很明顯，工作集越小或者步長越小，對應的都是山峰越高，也就是吞吐量越好，從迴圈主題上看，無論是多大的步長，工作集越小，越可能在L1或者L2快取中多停留，那麼在第二層迴圈遍歷時，現成快取資料被重複呼叫的可能性就越大， 所以時間區域性性就越好；（比如工作集大小為1，第一次讀1，快取123，第二次第三次就有現成的2、3用，但如果工作集大小為4，快取裝不下4個數，就算裝下123下一次也用不上，那麼每次都沒現成的，並且還要往更慢的下層快取要資料，當然時間區域性性就很差了）

        另一方面，工作集不變時，你步長越小，在第二層遍歷內部，上一次迴圈的資料被下一次迴圈共享的機率就越大，而步長越大，共享機會當然就越小，所以空間區域性性就越好。（比如步長為2，那麼第一次讀1、2，快取了1234，第二次讀3、4時就有現成的；如果步長為4，那麼每次幾乎都沒現成的了） 

        正因為有上面的分析，我們看到，空間區域性性的斜坡相對連續，而時間區域性性的斜坡就涇渭分明，體現了L1、L2和記憶體的讀吞吐量之間巨大的效能差別。

        另外我們還注意到，L1資料快取是32k，所以它在工作集超過32k時急速下降，符合預期。但L2大小是6M，那麼我們預測的吞吐量也應該在工作集為6M時下降才對，但從我圖上標註的兩條線來看，下降處居然是工作集為4M的地方，why？按書上的解釋是，L2不像L1那樣把指令快取和資料快取獨立區分，L2的資料和指令統一快取在一起，因此雖然大小是6M，但你資料可不能獨享這6M的空間，分2M給指令快取也是可以理解的……

        還有個有趣的現象，我們單獨來看L2區域的空間區域性性。我們假定工作集大小固定在512kB，也就是L2區域的中間部位，我們看到，步長從0~16的部分，山體沿著步長座標的方向有明顯的下降趨勢，這也符合之前對空間區域性性的預期——步長越長，空間區域性性越差！步長越長，被L1命中的可能性越低。但是從步長16開始往後走，坡度消失，幾乎變成飛機場，這是為啥呢？原來，步長為16字，對應的就是64位元組，這剛好是L1的塊大小，也就是說，之前L1的一個塊裡就可能快取了幾個步進大小的資料，一次訪問可能導致後面的資料在L1快取塊上命中，現在步進達到16，也就是64位元組的步進，那麼L1的每一個塊都不可能再有命中的存在，都必須從L2得到服務，因此吞吐量完全取決於L2傳送到L1的效能，因此保持不變。

        儲存器山是反映特定系統的時間和空間區域性性的山，對於高級別的程式設計師，一旦有了這樣的資訊，他就會盡可能的使得自己的程式中頻繁使用的字是從L1中存取，同時還要讓儘可能多的字從L1快取記憶體中訪問到。這就分別利用時間區域性性和空間區域性性。

二、簡單應用

        作為關心效能的程式設計師，知道對儲存器層次結構各個部分訪問時間的粗略估計值是很重要的。根據上面這匹山，估計下列這些位置讀出一個4位元組所需的時間，以CPU週期為單位（T9300的頻率是2.5GHz）：

        1、在晶片上的L1 d-cash；

        2、在晶片上的L2 d-cash；

        3、在主存上（工作集大小16M，步長=16），讀吞吐率為80MB/s

        為什麼要講這個看似簡單的題呢？也許從中能發現概念上的一個模糊點，至少我模糊得很，如果你也模糊，那恰好可以跟著理一理！

        第一個問， L1上讀取4位元組的時間，以CPU週期為單位。我們看L1的峰值吞吐率是10000MB/s，也就是10GB/s，而CPU的頻率是2.5GHz，因此4個位元組的讀取時間就是（2.5/10）*4 = 1週期，也就是接近一週期。

        首先明確下，關於G和M的引數描述中是（1GHZ=10^3MHZ=10^6KHZ= 10^9HZ）的關係，而不是1024。OK，關於上面這個式子，你完全理解透了沒？反正我剛開始是一頭霧水，這裡詳細糾結下週期和頻率的概念。

        下面是我面對糾結時的思考步驟：

                ①：CPU頻率是2.5GHz，說明一秒鐘執行2.5G次，那麼一個週期就是1/2.5G秒，∵L1的峰值吞吐率是10GB每秒，現在CPU一個週期又是1/2.5G秒，那麼一個週期的吞吐量就是兩者相乘：10/2.5=4B，也就是說，CPU一個週期內的吞吐量是4位元組——我去，這典型的是撞上答案！（當然，如果算出來是8B，我會用8/4算出2週期答案）

                ②：上面是先算出單位週期的吞吐量，然後與4位元組相除，得出週期數，現在從頻率出發來考慮：既然L1吞吐率是10GB每秒，那麼每一個位元組的耗時就是（1/10G）秒，4個位元組耗時自然就是（4/10G）秒，既然算出了總耗時，接著就是把它轉化成CPU週期數就OK了。怎麼轉換呢？考慮到CPU頻率是2.5GHz，說明一秒鐘執行2.5G次，現在我有（4/10G）秒，兩者相乘，就能算出總的執行次數：(2.5G/10G)*4 = 1次，這意味著什麼？意味著（4/10G）秒的這個時間量，剛好是CPU執行一次所需的時間，也就是CPU的週期！因此得出答案1週期！（當然，如果算出來是2次，那就說明答案是2週期！）

                ③：以上兩種方法都得出正確答案，但有存在兩個問題，其一，不可能每次算個時間週期數都要整這麼麻煩；其二，感覺沒有徹徹底底的理解清楚什麼吞吐率頻率和週期之間的關係，有云裡霧裡的錯覺。那麼接下來我就好好的整理了一下概念。

        週期T：每次迴圈的消耗多少時間、每次多少秒、Ts/次

        頻率f： 每個單位時間多少次迴圈、每秒多少次 、f次/s

        上面是我根據以往常識總結的不太嚴謹但夠用的概念，根據這兩個概念可以得知，所謂的週期，它的單位是（秒/次），描述的是每次迴圈的耗時；所謂頻率，單位是（次/s），描述每秒包含多少次迴圈，他們的單位和丈量標準剛好是相反的倒數關係，因此有f=1/T，T=1/f的換算公式，具體怎麼理解這兩個換算公式呢？我們先看第一個公式：

        f = 1/T，完整的寫就是f(次/s) = 1(s)/T(s/次)，我們發現，T表示的是每次多少秒，它雖然是個時間概念，但它被限制在僅僅一次迴圈範圍內的耗時秒數，（s/次）可以叫做"單位次時間"。好了，有了這個"單位次時間"，現在我關心對於一個確定的週期T，在總共1s的時間內會迴圈多少次呢？既然T衡量單位次（也就是一次）消耗的時間，現在用1s除以T，就相當於這個1s總時間被“單位次時間”砍成一節一節的段（一個總時間被一個特殊的時間量除），每個段都代表1次迴圈的時間量，有多少段就代表有多少次迴圈，好，這1s被砍了多少刀，也就是1s內有多少次迴圈，這恰恰就是每秒多少次——頻率的概念！

        再來看1(s)/T(s/次)這個運算，“次”是分母T的單位中的分母，當兩個s被約掉，“次”這個單位就會跑到分子上，最終答案就是(1/T)次，而我們的總時間是1s，因此也就是f(次/s)。注意，(1/T)計算結果本身的單位是“次”，只是因為分子剛好是1s，在定義頻率概念時，就有了f(次/s) = 1/T。

        接下來是T = 1/f，完整的寫就是T(s/次) = 1次/f(次/s)，f表示的是一秒鐘有多少次，雖然它有"次"的概念，但它是被限制在僅僅1秒中內的次數，（次/s）可以叫做“單位時間次數”，好了，有個了這個“單位時間次數”，對於一個確定的頻率f，在總共一次迴圈內會消耗多少秒時間呢？既然f衡量單位時間（也就是一秒）迴圈的次數，現在用1次除以f，就相當於這個1次的總迴圈次數被“單位時間次數”砍成一節一節的段（一個總次數被一個特殊的次數量除），每個段都代表1秒時間的次數，有多少段就代表有多少秒，好，這1次迴圈被砍了多少刀，也就是1次迴圈內有多少秒，這恰恰就是每次多少秒——週期的概念！

        同樣的，1(次)/f(次/s)這運算，結果是(1/f)s，只是由於分子剛好是1次，所以週期T的單位才是(s/次)。這裡注意的是，如果不止一次，那單位就應該是s了！

        有個問題，1s對應多少個週期呢？簡單，用總時間除以單位次時間，就能得出週期數，也就是1/T，哈哈，恰恰就是我們f!，也就是說，1s對應的週期數是f，那2s對應的週期數是2f，ns對應的週期數就是nf！所以頻率還可以理解成單位秒週期數。同理，1次對應多少時間呢？用總次數除以單位時間次，就能得出時間數，也就是1/f，恰恰就是週期T！也就是說，1次迴圈對應時間是T，那2次對應的時間就是2T，n對應的時間就是nT，所以週期還可以理解成單位週期時間數……這不廢話麼……

        好了，有了上面冗長的推理，終於對所謂的單位秒，單位次，單位時間等概念徹底理清了，物理中的很多除法公式也是類似，下面來看題。很明顯，（2.5/10）*4 = 1週期這個式子是先計算出來一個位元組的傳輸週期數，然後再計算4位元組。怎麼理解（2.5/10）計算的是傳輸1位元組的週期數呢？有兩種思維模式。其一，既然吞吐率是10GB/s，說明傳輸每個位元組的時間就是1/10G，根據上面的結論，有了時間n，週期數就是nf，於是就有2.5G*(1/10G)；其二，既然頻率是2.5GBHz，根據上面的結論，說明每秒鐘的週期數就是1*f也就是2.5G，既然1秒的週期數有2.5G，現在要處理1B，但每秒鐘的吞吐率是10GB，根本用不了1秒這麼長的時間，因此用總週期數除以多餘的10GB，2.5G/10G，得出的就是1B所需的週期數了……

         剩下斜字的部分屬於本人任性的部分，閒著沒事幹可以細讀。

        說到這，我再用自己的理解，來闡述下除法和分數的內涵。除法的本質到底是什麼？比如8÷2 = 4，有兩種理解方式：

        1、8本來是個整體，現在用長度為2的模具，從頭切到尾，發現剛好能切成4份長度為模具長度的塊，因此除法的結果是2——除法的結果，就是描述以模具長度標準來切割被除數後，所切下的模具長度的塊的個數！那麼3÷2呢？整體3被模具2切下一塊後，剩餘的1就不夠切了，於是整數位還是模具長度的塊的個數1，而剩下的1仍然要用來描述模具長度的個數，不夠一塊，那就有了1/2的概念，於是1.5的答案仍然是描述切割後模具長度的塊的個數。

        2、8本來是個整體，被均分成兩份，每份就是4這麼多，因此除法的結果是4——除法的結果，就是描述以特定份數對目標數進行均勻切割後，每份的數量。那麼1÷2呢？均分無法用整數實現，於是1/2就是答案本身。

        觀察兩種對除法的定義可以發現，定義1中的除法結果，是份，切割後剩餘的份數；而定義2中的除法結果，是數量本身，和被除數一樣。

        用定義1再來看1/2，我們發現，除法還可以理解成，將被除數這個整體，均分成除數份，比如1/2就是將1等分成兩份，有了這兩種理解方式？就可以解釋3/(1/2)了：以（1/2）為模具，切割3，我們發現這個模具比我們長度為1的模具還要小一半，因此切割出來的個數肯定更大，相當於6個1/2的數量，因此答案就是6，這個例子相當於再描述：先把模具切小，然後再用來切除數！

        我們熟悉的路程S單位是m，時間T單位是s，速度v單位是m/s。這個速度單位明顯就是個除法，表示1m/1s。如果S是6m，v是2(m/s)，時間是多少呢？答案簡單是3s，這裡面是有玄機的，為什麼除法的結果是s這個單位呢？有兩種解釋，一種解釋就是純代數運算，6m/2m/s，可以換算成(6m/2m)*s，s作為分母的分母，自然可以轉換到分子上來，而剩下的m可以通過約分消掉，得出答案……關鍵是我們該怎麼去理解這個運算的本質呢？為啥一個m單位的量除以一個m/s單位的量，結果就是s？m/s到底是個啥東西？

        可以觀察到，m/s是1m/1s這樣的一個除法式子，1m被1s這個模子切割，切得動麼？當然切不動，就像1/3也是切不動那樣，所以m/s就是最終結果。好關鍵在於這種複合單位如何做除法，它背後的意義是什麼？要理解這個類似除法的式子m/s，我們不妨換一個角度：

        一般的自然數1、2、3……如果代表數量，我們預設其單位為（1數量），比如7就代表7個（1數量）。現在來看這個複合除法：8/(6/3)，有兩種方式去理解：

        1、把分母的6看成是單位為（1/3）的量，把分子8看成是單位為（1）的量，於是就有了8(1)/6(1/3)，由於除法結果的單位是在分子，因此有必要把分母的單位轉換成（1），那麼依靠通分轉化：8*3(1)/6(1)——8(3)/6(1) = 4/3（3），我們得出4/3的單位為（3）！

        2、先把分母約分成2，此時2的單位已經是（1/3），我們再把分子單位進行轉換：8（1）——8/3（3），於是就有(8/3)（3）/2（1）= 4/3（3），結果相等。

        為什麼我一定要特別把答案構造成（3）單位呢？這是為了體現（1/3）和（3）這兩個單位的在除法中的轉換關係。我們按照上面1的理解，分母6的單位是（1/3），經過轉換，將分母單位轉換成（1），分子的單位變換成（3），出來的結果也是（3），這是不是很像路程和速度的計算公式(6m)/(2m)/s？原來分母的單位是（1/s），經過計算，單位1/s轉換到分子變成單位s，最後計算出來的結果也是就是s了！

        我們再回到上面的題目，週期的概念是s/次，頻率的概念是次/s，我們要計算耗時是多少週期，就是要按照Ts/次，把實際時間分割分割成一段一段的Ts，多少份就是多少個週期，然而我們發現，Ts是每次迴圈的時間，如果事先已經知道是多少次迴圈，那根本不需要關心實際時間。那到底是多少次呢？這道題的條件是頻率500M次/s，吞吐率1000MB/s，既然都是按“每秒”來衡量，每秒500M次，每秒吞吐量1000MB，那麼要計算1B（位元組）所執行的次數，就自然有500M次/s/1000MB/s = (1/2)次/B，每個位元組消耗1/2次，也就是1/2週期，那4個位元組當然就是2週期。事實上，如果按照另一個演算法，每個位元組所耗時間(1/1000)（s/B），每個週期（每次）耗時(1/500)（s/次），你把兩者一除會發現，答案的單位仍然是次/B。

三、矩陣乘法的優化

        這裡我們討論n×n矩陣的乘法：E=AB，比如n=3時，那麼：

e00	e01	e02
e10	e11	e12
e20	e21	e22

=  

a00	a01	a02
a10	a11	a12
a20	a21	a22

×

b00	b01	b02
b10	b11	b12
b20	b21	b22

e00 = a00×b00 + a01×b10 + a02×b20

e11 = a00×b01 + a01×b11 + a02×b21

……

        一般來說，第一個下標表示行，遞增方向是從上到下；第二個下標表示列，遞增方向是從左到右。對應C語言的二維陣列，則訪問順序是先遍歷每行的各列，然後再跳到下一行。

        總之，矩陣乘法的本質就是，目標矩陣的Exy的值，等於Ax(0~n)與B(0~n)y各自相乘後的加法結果。可以理解成A的x這行一排數，和B的y這一列數，拼成個十字架扔到Exy這個空空裡去，而十字架的運算方式就是先乘再加，我們暫且稱它為“交叉乘加”。

        英文裡行是row，也就是常說的排，column是列，也就是常說的縱隊。因此要實現矩陣乘法，分別用r、c來表示行列，用k來表示遞增變數，出現六個版本，下面參看第一個版本：

typedef double array[MAXN][MAXN]

void rck(array A, array B, array E, int n) 
{
    int r, c, k;
    double sum;

for (r = 0; r < n; r++) 
    for (c = 0; c < n; c++) {
sum = 0.0;
for (k = 0; k < n; k++)
   sum += A[r][k]*B[k][c];
E[r][c] += sum;
    }

}

        這個版本應該是最直觀明瞭的版本了。ABE都定義成陣列，n是矩陣的位數。遇上多重迴圈，我本人習慣從最裡面開始看起，sum += A[r][k]*B[k][c]，很明顯，是在計算上面的矩陣乘法，某個E的結果是多個乘數相加的結果，而sum += A[r][k]*B[k][c]就是在做其中一個乘法。k是遍歷n位，也就是說，把A確定的r行裡的每一列數遍歷完，同理B[k][c]是在確定的列縱隊c中遍歷每一排。兩個十字交叉完成，於是對應的E[r][c]位置也就有了它該有的值。

        接下來是第二層迴圈，對於列c的遍歷，重複裡層迴圈，那麼A完全相同的一行資料與B不同的列的資料進行交叉乘加，得出不同的E[r][c]，好了，當這層訓迴圈結束時，B的所有列都和A的某一個行交叉乘加完成，此時會跳出到第一層迴圈，也就是遍歷r的步驟。當A的每一行都重複之前同一行A資料交叉乘加B的每一列時，E[r][c]的每一個元素都被賦值，整個矩陣乘法結束。

        接下來看第二個版本：

vord crk(array A, array B, array E, int n) 
{
    int r, c, k;
    double sum;

for (c = 0; c < n; c++) 
    for (r = 0; r < n; r++) {
sum = 0.0;
for (k = 0; k < n; k++)
   sum += A[r][k]*B[k][c];
E[r][c] += sum;
    }
}

        很明顯裡層迴圈沒有動，外層迴圈顛倒了順序。具體分析方法和rck版本類似，先用A的每一行去交叉乘加B的某一個列，完成後調到第一層迴圈，變換B的列，然後重複裡層迴圈，當B的列也被遍歷完時，E[r][c]的每一個元素都被賦值，整個矩陣乘法結束。接著看第三個版本：

vord rkc(array A, array B, array E, int n) 
{
    int r, c, k;
    double m;
    
    for (r = 0; r < n; r++)
        for (k = 0; k < n; k++) {
            m = A[r][k];
            for (c = 0; c < n; c++)
                E[r][c] += m*B[k][c];
        }
}

        這個版本的演算法有了明顯的不同，我們先分析3×3矩陣，下圖描繪的是r=0時的，第二三層迴圈遍歷訪問到的ABE各元素區域：

a00	a01	a02	b00	b01	b02
b10	b11	b12
b20	b21	b22

e00

e01

e02

        e00 = a00×b00 + a01×b10 + a02×b20

        e01 = a00×b01 + a01×b11 + a02×b21

        e02 = a00×b02 + a01×b12 + a02×b22

        e10 = a10×b00 + a11×b10 + a12×b20

        r=0時，目標矩陣計算的就是第一排資料E[0][c]，由於e00~e02三個元素都是由B的第一排a00~a02去交叉乘加B的三個列，因此a00~a02三個元素各自都會被呼叫三次，該演算法就索性讓他們只調用一次。基本思想是讓e00~e02分成三次計算，每次只計算每個e0y式子裡的其中一個乘法，由裡層迴圈實現e00~e02的單次乘法，再由第二層迴圈變更a0y，重複裡層運算，直到每個e0y交叉乘加的三個乘法加拼裝計算完成時，e00~e02就計算完成。根據上面的式子，相當於在計算e同行元素結果時，把交叉乘加式子整體從左到右豎著遍歷。

e的其餘行以此類推，無非就是遍歷a1y和a2y。接著看krc版本：

void krc(array A, array B, array E, int n)
{
    int r, c, k;
    double m;
    for (k = 0; k < n; k++)
       for (r = 0; r < n; r++) {
    m = A[r][k];
    for (c = 0; c < n; c++)
        E[r][c] += m*B[k][c];
       }
}

        咋一看和rkc版本沒區別，認真看發現第二層和第一層的迴圈順序變了。那麼這我們假設3×3，k=0的情形：

a00	b00	b01	b02
a10
a20

e00	e01	e02
e10	e11	e12
e20	e21	e22

        e00 = a00×b00 + a01×b10 + a02×b20

        e01 = a00×b01 + a01×b11 + a02×b21

        e02 = a00×b02 + a01×b12 + a02×b22

        e10 = a10×b00 + a11×b10 + a12×b20

        e11 = a10×b01 + a11×b11 + a12×b21

        e12 = a10×b02 + a11×b12 + a12×b22

        ……

        e22 = a20×b02 + a21×b12 + a22×b22

        k=0時，第一層迴圈內已經把e00~e22所有元素過了一遍，由於受k=0的限制，在第二層迴圈內A遍歷的是第一列，第三層迴圈內B遍歷的是第一行。這裡很明顯，k=0時，裡層迴圈的每一步乘法運算，計算的都是exy是交叉乘加的第一個乘法，如上面的式子，相當於把所有exy交叉乘加式子按E的行順序遍歷完，r每遞增一次就遍歷三個exy。當第一層迴圈k遞增時，A的列和B的行隨即遞增，計算上面式子的第二豎。

        這個演算法實質上是把交叉乘加的這個十字架，徹底拆分，但基本思想和上面類似，都是要實現axy的多次利用，只是遍歷的方向不同罷了。



void kcr(array A, array B, array E, int n)
{
    int r, c, k;
    double m;

    for (k = 0; k < n; k++)
        for (c = 0; c < n; c++) {
    m = B[k][c];
    for (r = 0; r < n; r++)
        E[r][c] += A[r][k]*m;
    }
}

        咋一看就是把AB調換位置，但實際上把cr調換的位置，假設3×3，k=0的情形：

a00	b00	b01	b02
a10
a20

e00	e01	e02
e10	e11	e12
e20	e21	e22

        我們發現，kcr版本和krc版本在元素遍歷區域上幾乎一樣，只是裡層迴圈遍歷的是a00~a20，第二層迴圈遍歷的是b00~b02，光從對交叉乘加的拆分思想上來看，幾乎和kcr完全一樣。只是對E而言你，遍歷順序變成了e00 、e10 、e20 、e01、e11……也就是按E的列順序遍歷。

void ckr(array A, array B, array E, int n)
{
    int r, c, k;
    double m;

    for (c = 0; c < n; c++)
        for (k = 0; k < n; k++) {
   m = B[k][c];
    for (r = 0; r < n; r++)
        E[r][c] += A[r][k]*m;
        }
}

        ck再次交換順序，假設3×3，c=0的情形：

a00	a01	a02	b00
a10	a11	a12	b10
a20	a21	a22	b20

e00

e10

e20

        e00 = a00×b00 + a01×b10 + a02×b20

        e10 = a10×b00 + a11×b10 + a12×b20

        e20 = a20×b00 + a21×b10 + a22×b20

        e01 = a00×b01 + a01×b11 + a02×b21

        e11 = a10×b01 + a11×b11 + a12×b21

        e21 = a20×b01 + a21×b11 + a22×b21

        終於到了一輪外層迴圈就遍歷完所有A元素的演算法了。我們看到，c每遞增一次，只能計算出E的一個列，而裡層迴圈遍歷的是A的一個列，第二層迴圈遍歷的是b的一個列，因此從拆分交叉乘加的角度來看，裡層迴圈拆分了十字架的橫槓，第二層迴圈拆分的是豎槓。在看上面的式子，按E的列順序，裡層迴圈順序是豎著的，而第二層迴圈順序是橫著的。

        好了，六種演算法全部展示，接著就該分析效能了。很明顯，最平凡呼叫的是裡層迴圈，因此我們就以裡層迴圈作為分析物件，統計不命中率：

版本	每次迴圈載入	每次迴圈的儲存	每次迴圈A不命中率	每次迴圈B不命中率	每次迴圈E不命中率	每次迴圈的總不命中率
rck&crk	2	0	0.25	1.00	0.00	1.25
ckr&kcr	2	1	1.00	0.00	1.00	2.00
krc&rkc	2	1	0.00	0.25	0.25	0.50

        注意到六個版本可以劃分成三個等價類，劃分原則是裡層迴圈的內容完全一致。我們這裡有個前提：

        1、陣列元素是double型別，並且sizeof(double)=8，

        2、只有一個快取記憶體，塊大小為32位元組（B=32）

        3、n的實際值很大，使得矩陣的一行不能完全裝進L1快取記憶體中“中的一塊或一行”

        4、所有區域性變數都被編譯器安排到暫存器中儲存，不存在為區域性變數本身進行載入和儲存指令。

    “5”、

        （引號部分是教材上漏寫的部分，很關鍵，很容易引起歧義）

    好了，先來分析rck&crk,他們的裡層迴圈部分都是

    for (k = 0; k < n; k++)
       sum += A[r][k]*B[k][c];

    這裡可以明顯看出，對A的訪問是線性的，利用了空間區域性性，由於元素大小是8位元組，L1快取塊大小又是32位元組，因此對A的訪問不命中率是0.25（具體原因在本章第(4)節有詳細講解）。而對於B的訪問，由於是跳著行訪問，是非線性的，而且n足夠大L1快取不可能裝下B的整個一行，因此其不命中率就是1.0，所以rck&crk裡層迴圈的總不命中率是1.25。

    接下來分析ckr&kcr，他們的裡層迴圈部分都是

        for (r = 0; r < n; r++)
            E[r][c] += A[r][k]*r;

        顯然，這裡每次對A的訪問仍然是跨行的。因此A的不命中率是1.0。同理，裡層迴圈還用對E進行應用，由於變化的是r，因此E也是跨行訪問的，命中率仍然是1.0，所以ckr&kcr的裡層迴圈總不命中率是2.0

        最後分析krc&rkc，他們的裡層迴圈部分都是

            for (c = 0; c < n; c++)
                E[r][c] += m*B[k][c];

        用同樣分分析方法得出，B和E的訪問不命中率都是0.25，因此裡層迴圈總不命中率是0.5。

        注意到“每次迴圈載入“這一項，其實就是類似A[r][k]或B[r][k]的讀取操作，當然都是兩次。而”每次迴圈的儲存“則是對E[r][c]的寫入。由於rck&crk版本的裡層迴圈都由臨時變數sum儲存加值，不會產生對E[r][c]的寫入操作（儲存），因此是0。

        好了，光從不命中率上來看，勝負已分，效能關係是krc&rkc > rck&crk > ckr&kcr，實際情況是不是這樣呢？

        下面我要進行實際測試資料，我的CPU是酷睿2 p8400，查詢引數：

        http://www.cpu-world.com/CPUs/Core_2/Intel-Core%202%20Duo%20Mobile%20P8400%20AV80577SH0513M.html

查詢到L1快取記憶體的資料，如果只考慮單核，資料快取總大小C=32KB，塊大小B=64位元組，相聯度E=4，組數就是S = 32KB/(64B*4) = 2^15/2^8 = 2^7 =128，組數有128組，因此這款CPU的單核L1高速資料快取引數（S,E,B,m）= （128,4,64,32）。

Level 1 cache size  ?

2 x 32 KB 8-way set associative instruction caches 2 x 32 KB 8-way set associative write-back data caches

Cache:	L1 data	L1 instruction	L2
Size:	2 x 32 KB	2 x 32 KB	3 MB
Associativity:	8-way set associative	8-way set associative	12-way set associative
Line size:	64 bytes	64 bytes	64 bytes
Comments:	Direct-mapped	Direct-mapped	Non-inclusive Direct-mapped Shared between all cores

        如果你不能根據這些資料熟練的推到並理解（S,E,B,m）= （128,4,64,32）這個結果，那還是回去看下本章的第（3）節內容，那裡有詳細分析。

        總之我麼得到B = 64，和我們之前假設的32有出入，好吧，那不命中率等於0.25的地方事實上就得換成0.125（理由懶得解釋了，有本章前面的知識鋪墊），那麼三個版本的不命中率分別應該是1.125、2和0.25。同時，為了增加n使得A、B的一行資料足夠大到不能完全裝進L1快取記憶體"裡的一塊或者說一行"好了，現在塊大小B = 64B，而陣列元素sizeof(double) = 8B，那麼n必須大於64B/8B = 8。

    好吧，按照教材的實現方式，n從25到400，以25遞增，對6個版本進行效能測試，計算每次迴圈所消耗的週期數。執行結果如下：

[[email protected] matmult]# ./mm

matmult cycles/loop iteration

  n   rck   crk   ckr   kcr   krc   rkc

 25  0.13  0.46  0.23  0.02  0.01  0.00 

 50  0.08  0.06  0.05  0.01  0.04  0.01 

 75  0.04  0.01  0.03  0.02  0.04  0.03 

100  2.22  0.07  1.44  2.89  0.09  4.46 

125  0.72  2.30  3.04  5.77  0.71  1.91 

150  1.88  1.34  4.57  6.54  1.78  1.97 

175  1.56  1.66  9.69  8.43  1.95  2.36 

200  2.81  1.86 13.45 11.29  2.50  2.53 

225  2.70  2.70 18.04 14.86  2.70  2.90 

250  2.31  3.17 19.28 16.67  2.84  2.26 

275  3.73  3.55 19.64 17.28  3.15  2.91 

300  4.76  5.26 20.42 19.08  2.65  2.84 

325  5.00  5.25 20.40 19.26  3.16  3.16 

350  5.65  5.44 20.48 19.44  3.32  2.79 

375  7.28  7.47 20.61 17.47  3.26  3.16 

400  7.97  8.06 20.92 17.76  3.03  2.74 

     看著是不是有點暈？那我們就把這組資料轉換成效能統計圖來分析更加直觀：

    這個圖直觀的反映了六個版本的執行效率，縱座標表示每次裡層迴圈所需的CPU週期數，越高說明耗的週期數越多，效能也就越差。當n超過250時，我驚奇的發現，效能圖非常鮮明的把六個版本分成三個梯隊，krc&rkc效能最佳，rck&crk其次，ckr&kcr最糟糕，回過頭去看我們對6個版本的總不命中率分析，他們的結果竟然相同！

    為啥我要說驚人和竟然呢？理論和實際結論相同很奇怪麼？當然奇怪了。因為教科書上得出的結論和我稍有不同，在教材裡效能最優的是rck&crk！這再次體現了盡信書不如無書。實際測試結果很有意思。當教材作者得出rck&crk效能優於krc&rkc的結論時給出瞭解釋，認為不命中率並不說明一切，krc&rkc雖然有最少的不命中率，卻有額外的儲存器訪問：krc&rkc的每次裡層迴圈需要引用E[r][c]\A[r][k]\B[k][c]這三個儲存器，而rck&crk只有引用A[r][k]\B[k][c]兩個儲存器，因此儲存器引用影響了最終效能……而經過我的實際測試結果，是不是可以得出結論：這個儲存器影響因素消失了？是不是由於CPU更高階，使得儲存器引用的效率得以大幅提升？

    如果學過CPU內部原理應該有所瞭解，CPU在進行運算時，可能會經歷取指、譯碼、執行、訪存、寫回、更新PC等步驟，而其中的訪存階段可以將資料寫入儲存器，或者從儲存器讀取資料。而這裡面的所謂儲存器，有可能是L1快取，也有可能是L2快取甚至有可能是L3、記憶體、硬碟，因此訪存被認為是CPU執行指令的過程中最可能耗費更多時間的步驟。因此L1~L3甚至記憶體到硬碟的大小與讀取速度可能直接影響訪存的效率，我只能推斷，我執行測試時的電腦，從L1~L3到記憶體到硬碟可能都全面由於作者的測試電腦，那我能不能嘗試下獲得整個我的電腦儲存器層次的效能對測試結果的影響呢？這至少能涉及到從L1到記憶體的部分，好的，接下來再理一下我的虛擬機器linux配置的配置：

L1：32kB

L2:3072kB

L3：p8400沒有L3

記憶體：1024MB

    也就是說，當我對儲存器的引用增加到E\A\B三個時，n的大小將決定我的資料能快取到那一級。比如，我要限制在L1中，那麼n最大就不能超過對(32k/(8*3))開方那麼多次，也就是n<=37；如果我要限制在L2中，那就是限制n不超過對（3072k/(8*3)）開方那麼多次，也就是n<=362;如果同理，要限制在記憶體中，n<=11585，內啥，n大於2048時，我的小本本已經快跑糊了，要不是果斷shutdown估計可憐p8400就要報銷了，因此我只能分別測試L1和L2，然後呢？n<=37時結論已經有了，我們發現krc&rkc的效能仍然完爆其他對手，而在75~275之間，rck&crk昂首挺胸；也就是說，在資料藉助到L2快取時，L2的本身的效能劣勢使得訪存耗時增加。而當大於275時，krc&rkc幾乎已不可戰勝，即便當n>362，快取鐵定到了記憶體級別時，也不能改變這個趨勢。

    好了，為什麼分界點不在362而是在275呢？我猜呢，注意哈，是我猜，沒有充足依據，應該是資料快取還要處理其他資料，不只是三個矩陣獨享……