带你学习《深入理解计算机系统》程序性能优化探讨（5）——高速缓存、存储器山与矩阵乘法优化
? ? ? ? 這一節內容將綜合（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下一次也用不上，因為下次讀的是5678，那么每次都沒現成的，都是冷不命中，并且還要往更慢的下層緩存要數據，當然時間局部性就很差了）
? ? ? ? 另一方面，工作集不變時，你步長越小，在第二層遍歷內部，上一次循環的數據被下一次循環共享的幾率就越大，而步長越大，共享機會當然就越小，所以空間局部性就越差。（比如步長為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
? ? ? ? 上面是我根據以往常識總結的不太嚴謹但夠用的概念，根據這兩個概念可以得知，所謂的周期，它的單位是（秒/次），描述的是每次循環的耗時；所謂頻率，單位是（次/秒），描述每秒包含多少次循環，他們的單位和丈量標準剛好是相反的倒數關系，因此有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，于是就有(1/10G)*2.5G；其二，既然頻率是2.5GBHz，根據上面的結論，說明每秒鐘的周期數就是1*f也就是2.5G，既然1秒的周期數有2.5G，現在要處理1B，但每秒鐘的吞吐率是10GB，1B根本用不了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，我們要計算耗時是多少周期，本質就是用周期T這個特殊單位的衡量時間，就是要按照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 |
| 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個版本進行性能測試，計算每次循環所消耗的周期數。運行結果如下：
?
[root@localhost?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呢？我猜呢，注意哈，是我猜，沒有充足依據，應該是數據緩存還要處理其他數據，不只是三個矩陣獨享……
?
?
?
?
?
?
?
?
?
總結
以上是生活随笔為你收集整理的带你学习《深入理解计算机系统》程序性能优化探讨（5）——高速缓存、存储器山与矩阵乘法优化的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： java gettime_Java Ut
- 下一篇： 【系统分析师之路】 第八章 复盘软件测试