带你学习《深入理解计算机系统》程序性能优化探讨（5）——高速缓存、存储器山与矩阵乘法优化
? ? ? ? 這一節(jié)內(nèi)容將綜合（3）和（4），討論高速緩存相關(guān)的程序優(yōu)化。
?
一、牛B完了的存儲器山
? ? ? ? 一個程序從存儲系統(tǒng)中讀數(shù)據(jù)的速率被稱為讀吞吐量或讀帶寬。如果一個程序在s秒的時間段內(nèi)讀n個字節(jié)，那么讀吞吐量就是n/s，一般用MB/s作為單位。
? ? ? ? 第（4）節(jié)中討論的時間局部性和空間局部性。從緩存塊大小來權(quán)衡，時間局部性和空間局部性似乎剛好成反比，但當(dāng)我們?nèi)娴挠懻撜麄€存儲層次性能時會發(fā)現(xiàn)，其實兩種局部性可能各自有各自的規(guī)律。
? ? ? ? 如果我們把一次性讀取的數(shù)據(jù)量的大小稱為工作集，把上一次和下一次讀取數(shù)據(jù)的之間存儲距離稱為步長，如果我想獲得工作集和步長這兩種因素影響下程序的讀吞吐量，會是什么結(jié)果呢？教材中給出實現(xiàn)代碼，但我不打算在這里全部貼出來，因為詳細的代碼分析涉及到教材后面章節(jié)里的K次最優(yōu)測量方法←_←內(nèi)啥，我自己還沒看懂，準確的說是暫時沒耐著性子看下去，好像天書，等哪天有心情看懂了，再補上代碼實現(xiàn)也不遲，那我們就先把存儲器山的圖貼出來吧：
?
? ? ? ? 我們看到，存儲器山由三個坐標(biāo)確定點位，底面右邊的坐標(biāo)是工作集大小（wording set size），按字節(jié)計數(shù)；底面左邊的坐標(biāo)是步長（stride），按字計數(shù)。這里可以把字理解成4字節(jié)；垂直的坐標(biāo)就是讀吞吐量，根據(jù)吞吐量的大小判斷程序性能的優(yōu)劣。注意到這個存儲器山針對的CPU型號已經(jīng)寫在上面了，這個CPU里面分別有大小2*32kb、塊大小64字節(jié)的高速緩存L1，指令和數(shù)據(jù)緩存都是這么大。另外還有6M大小的高速緩存L2。
?
? ? ? ? 首先把步長作為常數(shù)，來看工作集32k大小時的情況。一次性處理32k字節(jié)的數(shù)據(jù)，由于L1有32k的指令緩存和32k數(shù)據(jù)緩存，因此，32k工作集的數(shù)據(jù)是完全可以存入L1緩存的，因此可以看到此時它的吞吐量處于最優(yōu)狀態(tài)。奇怪的是，當(dāng)工作集坐標(biāo)往右減小時，山峰沒有繼續(xù)升高，而是急速降低，你可以明顯的看到，從32k之后的下降陡坡，這是為什么？其實，這里需要了解測試程序的結(jié)構(gòu)：
#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));
}
? ? ? ? 上面這個是生成存儲器山的核心循環(huán)，第一層循環(huán)是遍歷工作集大小，從32MB開始兩倍變小，直到4kB為止；第二層循環(huán)是步進，從1到30字（4字節(jié)）。有了這兩個測試關(guān)鍵參數(shù)傳入存儲器山生成函數(shù)run里，就能將這座山遍歷出來。
? ? ? ? 好了，當(dāng)?shù)谝粚觭ize取值小于32k，而步長又大于20時，由于工作集size太小，因此每一個工作集處理時間會非常短，但正因為工作集太小，而步長又太大，數(shù)據(jù)集越不夠連續(xù)，命中的情況就越遭，因此懲罰也就越大，主要的時間就消耗在不命中處罰和循環(huán)本身的開銷上，因此山高峰背后那個明顯下降，其實并不能反映出L1高速緩存的真實性能。
? ? ? ? 我們再來看看這座山，很明顯，工作集越小或者步長越小，對應(yīng)的都是山峰越高，也就是吞吐量越好，從循環(huán)主題上看，無論是多大的步長，工作集越小，越可能在L1或者L2緩存中多停留，那么在第二層循環(huán)遍歷時，現(xiàn)成緩存數(shù)據(jù)被重復(fù)調(diào)用的可能性就越大， 所以時間局部性就越好；（比如工作集大小為1，第一次讀1，緩存123，第二次第三次就有現(xiàn)成的2、3用，但如果工作集大小為4，緩存裝不下4個數(shù)，就算裝下123下一次也用不上，因為下次讀的是5678，那么每次都沒現(xiàn)成的，都是冷不命中，并且還要往更慢的下層緩存要數(shù)據(jù)，當(dāng)然時間局部性就很差了）
? ? ? ? 另一方面，工作集不變時，你步長越小，在第二層遍歷內(nèi)部，上一次循環(huán)的數(shù)據(jù)被下一次循環(huán)共享的幾率就越大，而步長越大，共享機會當(dāng)然就越小，所以空間局部性就越差。（比如步長為2，那么第一次讀1、2，緩存了1234，第二次讀3、4時就有現(xiàn)成的；如果步長為4，那么每次幾乎都沒現(xiàn)成的了）?
? ? ? ? 正因為有上面的分析，我們看到，空間局部性的斜坡相對連續(xù)，而時間局部性的斜坡就涇渭分明，體現(xiàn)了L1、L2和內(nèi)存的讀吞吐量之間巨大的性能差別。
?
? ? ? ? 另外我們還注意到，L1數(shù)據(jù)緩存是32k，所以它在工作集超過32k時急速下降，符合預(yù)期。但L2大小是6M，那么我們預(yù)測的吞吐量也應(yīng)該在工作集為6M時下降才對，但從我圖上標(biāo)注的兩條線來看，下降處居然是工作集為4M的地方，why？按書上的解釋是，L2不像L1那樣把指令緩存和數(shù)據(jù)緩存獨立區(qū)分，L2的數(shù)據(jù)和指令統(tǒng)一緩存在一起，因此雖然大小是6M，但你數(shù)據(jù)可不能獨享這6M的空間，分2M給指令緩存也是可以理解的……
? ? ? ? 還有個有趣的現(xiàn)象，我們單獨來看L2區(qū)域的空間局部性。我們假定工作集大小固定在512kB，也就是L2區(qū)域的中間部位，我們看到，步長從0~16的部分，山體沿著步長的下坐標(biāo)的方向有明顯降趨勢，這也符合之前對空間局部性的預(yù)期——步長越長，空間局部性越差！步長越長，被L1命中的可能性越低。但是從步長16開始往后走，坡度消失，幾乎變成飛機場，這是為啥呢？原來，步長為16字，對應(yīng)的就是64字節(jié)，這剛好是L1的塊大小，也就是說，之前L1的一個塊里就可能緩存了幾個步進大小的數(shù)據(jù)，一次訪問可能導(dǎo)致后面的數(shù)據(jù)在L1緩存塊上命中，現(xiàn)在步進超過16，也就是超過64字節(jié)的步進，那么L1上的每一個塊數(shù)據(jù)都不可能再有命中的存在，都必須從L2得到服務(wù)，因此吞吐量完全取決于L2傳送到L1的性能，因此保持不變。
? ? ? ? 存儲器山是反映特定系統(tǒng)的時間和空間局部性的山，對于高級別的程序員，一旦有了這樣的信息，他就會盡可能的使得自己的程序中頻繁使用的字是從L1中存取，同時還要讓盡可能多的字從L1高速緩存中訪問到。這就分別利用時間局部性和空間局部性。
?
二、簡單應(yīng)用
? ? ? ? 作為關(guān)心性能的程序員，知道對存儲器層次結(jié)構(gòu)各個部分訪問時間的粗略估計值是很重要的。根據(jù)上面這匹山，估計下列這些位置讀出一個4字節(jié)所需的時間，以CPU周期為單位（T9300的頻率是2.5GHz）：
? ? ? ? 1、在芯片上的L1 d-cash；
? ? ? ? 2、在芯片上的L2 d-cash；
? ? ? ? 3、在主存上（工作集大小16M，步長=16），讀吞吐率為80MB/s
?
? ? ? ? 為什么講這看似簡單的題呢？也許從中能發(fā)現(xiàn)概念上的一個模糊點，至少我模糊得很，如果你也模糊，那恰好可以跟著理一理！
? ? ? ? 第一個問， L1上讀取4字節(jié)的時間，以CPU周期為單位。我們看L1的峰值吞吐率是10000MB/s，也就是10GB/s，而CPU的頻率是2.5GHz，因此4個字節(jié)的讀取時間就是（2.5/10）*4 = 1周期，也就是接近一周期。
? ? ? ? 首先明確下，關(guān)于G和M的參數(shù)描述中是（1GHZ=10^3MHZ=10^6KHZ= 10^9HZ）的關(guān)系，而不是1024。OK，關(guān)于上面這個式子，你完全理解透了沒？反正我剛開始是一頭霧水，這里詳細糾結(jié)下周期和頻率的概念。
? ? ? ? 下面是我面對糾結(jié)時的思考步驟：
? ? ? ? ? ? ? ? ①：CPU頻率是2.5GHz，說明一秒鐘運行2.5G次，那么一個周期就是1/2.5G秒，∵L1的峰值吞吐率是10GB每秒，現(xiàn)在CPU一個周期又是1/2.5G秒，那么一個周期的吞吐量就是兩者相乘：10/2.5=4B，也就是說，CPU一個周期內(nèi)的吞吐量是4字節(jié)——我去，這典型的是撞上答案！（當(dāng)然，如果算出來是8B，我會用8/4算出2周期答案）
? ? ? ? ? ? ? ? ②：上面是先算出單位周期的吞吐量，然后與4字節(jié)相除，得出周期數(shù)，現(xiàn)在從頻率出發(fā)來考慮：既然L1吞吐率是10GB每秒，那么每一個字節(jié)的耗時就是（1/10G）秒，4個字節(jié)耗時自然就是（4/10G）秒，既然算出了總耗時，接著就是把它轉(zhuǎn)化成CPU周期數(shù)就OK了。怎么轉(zhuǎn)換呢？考慮到CPU頻率是2.5GHz，說明一秒鐘運行2.5G次，現(xiàn)在我有（4/10G）秒，兩者相乘，就能算出總的運行次數(shù)：(2.5G/10G)*4 = 1次，這意味著什么？意味著（4/10G）秒的這個時間量，剛好是CPU運行一次所需的時間，也就是CPU的周期！因此得出答案1周期！（當(dāng)然，如果算出來是2次，那就說明答案是2周期！）
? ? ? ? ? ? ? ? ③：以上兩種方法都得出正確答案，但有存在兩個問題，其一，不可能每次算個時間周期數(shù)都要整這么麻煩；其二，感覺沒有徹徹底底的理解清楚什么吞吐率頻率和周期之間的關(guān)系，有云里霧里的錯覺。那么接下來我就好好的整理了一下概念。
? ? ? ? 周期T：每次循環(huán)的消耗多少時間、每次多少秒、Ts/次
? ? ? ? 頻率f： 每個單位時間多少次循環(huán)、每秒多少次 、f次/s
? ? ? ? 上面是我根據(jù)以往常識總結(jié)的不太嚴謹?shù)珘蛴玫母拍?#xff0c;根據(jù)這兩個概念可以得知，所謂的周期，它的單位是（秒/次），描述的是每次循環(huán)的耗時；所謂頻率，單位是（次/秒），描述每秒包含多少次循環(huán)，他們的單位和丈量標(biāo)準剛好是相反的倒數(shù)關(guān)系，因此有f=1/T，T=1/f的換算公式，具體怎么理解這兩個換算公式呢？我們先看第一個公式：
? ? ? ? f = 1/T，完整的寫就是f(次/s) = 1(s)/T(s/次)，我們發(fā)現(xiàn)，T表示的是每次多少秒，它雖然是個時間概念，但它被限制在僅僅一次循環(huán)范圍內(nèi)的耗時秒數(shù)，（s/次）可以叫做"單位次時間"。好了，有了這個"單位次時間"，現(xiàn)在我關(guān)心對于一個確定的周期T，在總共1s的時間內(nèi)會循環(huán)多少次呢？既然T衡量單位次（也就是一次）消耗的時間，現(xiàn)在用1s除以T，就相當(dāng)于這個1s總時間被“單位次時間”砍成一節(jié)一節(jié)的段（一個總時間被一個特殊的時間量除），每個段都代表1次循環(huán)的時間量，有多少段就代表有多少次循環(huán)，好，這1s被砍了多少刀，也就是1s內(nèi)有多少次循環(huán)，這恰恰就是每秒多少次——頻率的概念！
? ? ? ? 再來看1(s)/T(s/次)這個運算，“次”是分母T的單位中的分母，當(dāng)兩個s被約掉，“次”這個單位就會跑到分子上，最終答案就是(1/T)次，而我們的總時間是1s，因此也就是f(次/s)。注意，(1/T)計算結(jié)果本身的單位是“次”，只是因為分子剛好是1s，在定義頻率概念時，就有了f(次/s) =?1/T。
? ? ? ? 接下來是T = 1/f，完整的寫就是T(s/次) = 1次/f(次/s)，f表示的是一秒鐘有多少次，雖然它有"次"的概念，但它是被限制在僅僅1秒中內(nèi)的次數(shù)，（次/s）可以叫做“單位時間次數(shù)”，好了，有個了這個“單位時間次數(shù)”，對于一個確定的頻率f，在總共一次循環(huán)內(nèi)會消耗多少秒時間呢？既然f衡量單位時間（也就是一秒）循環(huán)的次數(shù)，現(xiàn)在用1次除以f，就相當(dāng)于這個1次的總循環(huán)次數(shù)被“單位時間次數(shù)”砍成一節(jié)一節(jié)的段（一個總次數(shù)被一個特殊的次數(shù)量除），每個段都代表1秒時間的次數(shù)，有多少段就代表有多少秒，好，這1次循環(huán)被砍了多少刀，也就是1次循環(huán)內(nèi)有多少秒，這恰恰就是每次多少秒——周期的概念！
? ? ? ? 同樣的，1(次)/f(次/s)這運算，結(jié)果是(1/f)s，只是由于分子剛好是1次，所以周期T的單位才是(s/次)。這里注意的是，如果不止一次，那單位就應(yīng)該是s了！
? ? ? ? 有個問題，1s對應(yīng)多少個周期呢？簡單，用總時間除以單位次時間，就能得出周期數(shù)，也就是1/T，哈哈，恰恰就是我們f!，也就是說，1s對應(yīng)的周期數(shù)是f，那2s對應(yīng)的周期數(shù)是2f，ns對應(yīng)的周期數(shù)就是nf！所以頻率還可以理解成單位秒周期數(shù)——方便算周期。同理，1次對應(yīng)多少時間呢？用總次數(shù)除以單位時間次，就能得出時間數(shù)，也就是1/f，恰恰就是周期T！也就是說，1次循環(huán)對應(yīng)時間是T，那2次對應(yīng)的時間就是2T，n次對應(yīng)的時間就是nT，所以周期還可以理解成單位周期秒數(shù)——方便算時間（秒數(shù)）……這不廢話么……
?
? ? ? ? 好了，有了上面冗長的推理，終于對所謂的單位秒，單位次，單位時間等概念徹底理清了，物理中的很多除法公式也是類似，下面來看題。很明顯，（2.5/10）*4 = 1周期這個式子是先計算出來一個字節(jié)的傳輸周期數(shù)，然后再計算4字節(jié)。怎么理解（2.5/10）計算的是傳輸1字節(jié)的周期數(shù)呢？有兩種思維模式。其一，既然吞吐率是10GB/s，說明傳輸每個字節(jié)的時間就是1/10G，根據(jù)上面的結(jié)論，有了時間n，周期數(shù)就是nf，于是就有(1/10G)*2.5G；其二，既然頻率是2.5GBHz，根據(jù)上面的結(jié)論，說明每秒鐘的周期數(shù)就是1*f也就是2.5G，既然1秒的周期數(shù)有2.5G，現(xiàn)在要處理1B，但每秒鐘的吞吐率是10GB，1B根本用不了1秒這么長的時間，因此用總周期數(shù)除以多余的10GB，2.5G/10G，得出的就是1B所需的周期數(shù)了……
? ? ? ? ?剩下斜字的部分屬于本人任性的部分，閑著沒事干可以細讀。
? ? ? ? 說到這，我再用自己的理解，來闡述下除法和分數(shù)的內(nèi)涵。除法的本質(zhì)到底是什么？比如8÷2 = 4，有兩種理解方式：
? ? ? ? 1、8本來是個整體，現(xiàn)在用長度為2的模具，從頭切到尾，發(fā)現(xiàn)剛好能切成4份長度為模具長度的塊，因此除法的結(jié)果是2——除法的結(jié)果，就是描述以模具長度標(biāo)準來切割被除數(shù)后，所切下的模具長度的塊的個數(shù)！那么3÷2呢？整體3被模具2切下一塊后，剩余的1就不夠切了，于是整數(shù)位還是模具長度的塊的個數(shù)1，而剩下的1仍然要用來描述模具長度的個數(shù)，不夠一塊，那就有了1/2的概念，于是1.5的答案仍然是描述切割后模具長度的塊的個數(shù)。
? ? ? ? 2、8本來是個整體，被均分成兩份，每份就是4這么多，因此除法的結(jié)果是4——除法的結(jié)果，就是描述以特定份數(shù)對目標(biāo)數(shù)進行均勻切割后，每份的數(shù)量。那么1÷2呢？均分無法用整數(shù)實現(xiàn)，于是1/2就是答案本身。
? ? ? ? 觀察兩種對除法的定義可以發(fā)現(xiàn)，定義1中的除法結(jié)果，是份，切割后剩余的份數(shù)；而定義2中的除法結(jié)果，是數(shù)量本身，和被除數(shù)一樣。
? ? ? ? 用定義1再來看1/2，我們發(fā)現(xiàn)，除法還可以理解成，將被除數(shù)這個整體，均分成除數(shù)份，比如1/2就是將1等分成兩份，有了這兩種理解方式？就可以解釋3/(1/2)了：以（1/2）為模具，切割3，我們發(fā)現(xiàn)這個模具比我們長度為1的模具還要小一半，因此切割出來的個數(shù)肯定更大，相當(dāng)于6個1/2的數(shù)量，因此答案就是6，這個例子相當(dāng)于再描述：先把模具切小，然后再用來切除數(shù)！
? ? ? ? 我們熟悉的路程S單位是m，時間T單位是s，速度v單位是m/s。這個速度單位明顯就是個除法，表示1m/1s。如果S是6m，v是2(m/s)，時間是多少呢？答案簡單是3s，這里面是有玄機的，為什么除法的結(jié)果是s這個單位呢？有兩種解釋，一種解釋就是純代數(shù)運算，6m/2m/s，可以換算成(6m/2m)*s，s作為分母的分母，自然可以轉(zhuǎn)換到分子上來，而剩下的m可以通過約分消掉，得出答案……關(guān)鍵是我們該怎么去理解這個運算的本質(zhì)呢？為啥一個m單位的量除以一個m/s單位的量，結(jié)果就是s？m/s到底是個啥東西？
? ? ? ? 可以觀察到，m/s是1m/1s這樣的一個除法式子，1m被1s這個模子切割，切得動么？當(dāng)然切不動，就像1/3也是切不動那樣，所以m/s就是最終結(jié)果。好關(guān)鍵在于這種復(fù)合單位如何做除法，它背后的意義是什么？要理解這個類似除法的式子m/s，我們不妨換一個角度：
? ? ? ? 一般的自然數(shù)1、2、3……如果代表數(shù)量，我們默認其單位為（1數(shù)量），比如7就代表7個（1數(shù)量）。現(xiàn)在來看這個復(fù)合除法：8/(6/3)，有兩種方式去理解：
? ? ? ? 1、把分母的6看成是單位為（1/3）的量，把分子8看成是單位為（1）的量，于是就有了8(1)/6(1/3)，由于除法結(jié)果的單位是在分子，因此有必要把分母的單位轉(zhuǎn)換成（1），那么依靠通分轉(zhuǎn)化：8*3(1)/6(1)——8(3)/6(1) = 4/3（3），我們得出4/3的單位為（3）！
? ? ? ? 2、先把分母約分成2，此時2的單位已經(jīng)是（1/3），我們再把分子單位進行轉(zhuǎn)換：8（1）——8/3（3），于是就有(8/3)（3）/2（1）= 4/3（3），結(jié)果相等。
? ? ? ? 為什么我一定要特別把答案構(gòu)造成（3）單位呢？這是為了體現(xiàn)（1/3）和（3）這兩個單位的在除法中的轉(zhuǎn)換關(guān)系。我們按照上面1的理解，分母6的單位是（1/3），經(jīng)過轉(zhuǎn)換，將分母單位轉(zhuǎn)換成（1），分子的單位變換成（3），出來的結(jié)果也是（3），這是不是很像路程和速度的計算公式(6m)/(2m)/s？原來分母的單位是（1/s），經(jīng)過計算，單位1/s轉(zhuǎn)換到分子變成單位s，最后計算出來的結(jié)果也是就是s了！
? ? ? ?
? ? ? ? 我們再回到上面的題目，周期的概念是s/次，頻率的概念是次/s，我們要計算耗時是多少周期，本質(zhì)就是用周期T這個特殊單位的衡量時間，就是要按照Ts/次，把實際時間分割成一段一段的Ts，多少份就是多少個周期，然而我們發(fā)現(xiàn)，Ts是每次循環(huán)的時間，如果事先已經(jīng)知道是多少次循環(huán)，那根本不需要關(guān)心實際時間。那到底是多少次呢？這道題的條件是頻率500M次/s，吞吐率1000MB/s，既然都是按“每秒”來衡量，每秒500M次，每秒吞吐量1000MB，那么要計算1B（字節(jié)）所執(zhí)行的次數(shù)，就自然有500M次/s/1000MB/s = (1/2)次/B，每個字節(jié)消耗1/2次，也就是1/2周期，那4個字節(jié)當(dāng)然就是2周期。事實上，如果按照另一個算法，每個字節(jié)所耗時間(1/1000)（s/B），每個周期（每次）耗時(1/500)（s/次），你把兩者一除會發(fā)現(xiàn)，答案的單位仍然是次/B。
?
三、矩陣乘法的優(yōu)化
? ? ? ? 這里我們討論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
……
?
? ? ? ? 一般來說，第一個下標(biāo)表示行，遞增方向是從上到下；第二個下標(biāo)表示列，遞增方向是從左到右。對應(yīng)C語言的二維數(shù)組，則訪問順序是先遍歷每行的各列，然后再跳到下一行。
? ? ? ? 總之，矩陣乘法的本質(zhì)就是，目標(biāo)矩陣的Exy的值，等于Ax(0~n)與B(0~n)y各自相乘后的加法結(jié)果。可以理解成A的x這行一排數(shù)，和B的y這一列數(shù)，拼成個十字架扔到Exy這個空空里去，而十字架的運算方式就是先乘再加，我們暫且稱它為“交叉乘加”。
? ? ? ? 英文里行是row，也就是常說的排，column是列，也就是常說的縱隊。因此要實現(xiàn)矩陣乘法，分別用r、c來表示行列，用k來表示遞增變量，出現(xiàn)六個版本，下面參看第一個版本：
?
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;
? ? }
}
? ? ? ? 這個版本應(yīng)該是最直觀明了的版本了。ABE都定義成數(shù)組，n是矩陣的位數(shù)。遇上多重循環(huán)，我本人習(xí)慣從最里面開始看起，sum += A[r][k]*B[k][c]，很明顯，是在計算上面的矩陣乘法，某個E的結(jié)果是多個乘數(shù)相加的結(jié)果，而sum += A[r][k]*B[k][c]就是在做其中一個乘法。k是遍歷n位，也就是說，把A確定的r行里的每一列數(shù)遍歷完，同理B[k][c]是在確定的列縱隊c中遍歷每一排。兩個十字交叉完成，于是對應(yīng)的E[r][c]位置也就有了它該有的值。
? ? ? ? 接下來是第二層循環(huán)，對于列c的遍歷，重復(fù)里層循環(huán)，那么A完全相同的一行數(shù)據(jù)與B不同的列的數(shù)據(jù)進行交叉乘加，得出不同的E[r][c]，好了，當(dāng)這層訓(xùn)循環(huán)結(jié)束時，B的所有列都和A的某一個行交叉乘加完成，此時會跳出到第一層循環(huán)，也就是遍歷r的步驟。當(dāng)A的每一行都重復(fù)之前同一行A數(shù)據(jù)交叉乘加B的每一列時，E[r][c]的每一個元素都被賦值，整個矩陣乘法結(jié)束。
?
? ? ? ? 接下來看第二個版本：
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;
? ? }
}
?
? ? ? ? 很明顯里層循環(huán)沒有動，外層循環(huán)顛倒了順序。具體分析方法和rck版本類似，先用A的每一行去交叉乘加B的某一個列，完成后調(diào)到第一層循環(huán)，變換B的列，然后重復(fù)里層循環(huán)，當(dāng)B的列也被遍歷完時，E[r][c]的每一個元素都被賦值，整個矩陣乘法結(jié)束。接著看第三個版本：
?
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時的，第二三層循環(huán)遍歷訪問到的ABE各元素區(qū)域：
?
| 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時，目標(biāo)矩陣計算的就是第一排數(shù)據(jù)E[0][c]，由于e00~e02三個元素都是由B的第一排a00~a02去交叉乘加B的三個列，因此a00~a02三個元素各自都會被調(diào)用三次，該算法就索性讓他們只調(diào)用一次。基本思想是讓e00~e02分成三次計算，每次只計算每個e0y式子里的其中一個乘法，由里層循環(huán)實現(xiàn)e00~e02的單次乘法，再由第二層循環(huán)變更a0y，重復(fù)里層運算，直到每個e0y交叉乘加的三個乘法加拼裝計算完成時，e00~e02就計算完成。根據(jù)上面的式子，相當(dāng)于在計算e同行元素結(jié)果時，把交叉乘加式子整體從左到右豎著遍歷。
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版本沒區(qū)別，認真看發(fā)現(xiàn)第二層和第一層的循環(huán)順序變了。那么這我們假設(shè)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時，第一層循環(huán)內(nèi)已經(jīng)把e00~e22所有元素過了一遍，由于受k=0的限制，在第二層循環(huán)內(nèi)A遍歷的是第一列，第三層循環(huán)內(nèi)B遍歷的是第一行。這里很明顯，k=0時，里層循環(huán)的每一步乘法運算，計算的都是exy是交叉乘加的第一個乘法，如上面的式子，相當(dāng)于把所有exy交叉乘加式子按E的行順序遍歷完，r每遞增一次就遍歷三個exy。當(dāng)?shù)谝粚友h(huán)k遞增時，A的列和B的行隨即遞增，計算上面式子的第二豎。
? ? ? ? 這個算法實質(zhì)上是把交叉乘加的這個十字架，徹底拆分，但基本思想和上面類似，都是要實現(xiàn)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調(diào)換位置，但實際上把cr調(diào)換的位置，假設(shè)3×3，k=0的情形：
?
?
| a00 | ? | ? | ? | b00 | b01 | b02 |
| a10 | ? | ? | ? | ? | ? | ? |
| a20 | ? | ? | ? | ? | ? | ? |
?
| e00 | e01 | e02 |
| e10 | e11 | e12 |
| e20 | e21 | e22 |
?
? ? ? ? 我們發(fā)現(xiàn)，kcr版本和krc版本在元素遍歷區(qū)域上幾乎一樣，只是里層循環(huán)遍歷的是a00~a20，第二層循環(huán)遍歷的是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再次交換順序，假設(shè)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
?
? ? ? ? 終于到了一輪外層循環(huán)就遍歷完所有A元素的算法了。我們看到，c每遞增一次，只能計算出E的一個列，而里層循環(huán)遍歷的是A的一個列，第二層循環(huán)遍歷的是b的一個列，因此從拆分交叉乘加的角度來看，里層循環(huán)拆分了十字架的橫杠，第二層循環(huán)拆分的是豎杠。在看上面的式子，按E的列順序，里層循環(huán)順序是豎著的，而第二層循環(huán)順序是橫著的。
?
? ? ? ??好了，六種算法全部展示，接著就該分析性能了。很明顯，最平凡調(diào)用的是里層循環(huán)，因此我們就以里層循環(huán)作為分析對象，統(tǒng)計不命中率：
?
| 版本 | 每次循環(huán)加載 | 每次循環(huán)的存儲 | 每次循環(huán)A不命中率 | 每次循環(huán)B不命中率 | 每次循環(huán)E不命中率 | 每次循環(huán)的總不命中率 |
| 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 |
?
?
????????注意到六個版本可以劃分成三個等價類，劃分原則是里層循環(huán)的內(nèi)容完全一致。我們這里有個前提：
????????1、數(shù)組元素是double類型，并且sizeof(double)=8，
????????2、只有一個高速緩存，塊大小為32字節(jié)（B=32）
????????3、n的實際值很大，使得矩陣的一行不能完全裝進L1高速緩存中“中的一塊或一行”
????????4、所有局部變量都被編譯器安排到寄存器中存儲，不存在為局部變量本身進行加載和存儲指令。
? ? “5”、
????????（引號部分是教材上漏寫的部分，很關(guān)鍵，很容易引起歧義）
? ? 好了，先來分析rck&crk,他們的里層循環(huán)部分都是
? ? for (k = 0; k < n; k++)
? ? ? ?sum += A[r][k]*B[k][c];
? ? 這里可以明顯看出，對A的訪問是線性的，利用了空間局部性，由于元素大小是8字節(jié)，L1緩存塊大小又是32字節(jié)，因此對A的訪問不命中率是0.25（具體原因在本章第(4)節(jié)有詳細講解）。而對于B的訪問，由于是跳著行訪問，是非線性的，而且n足夠大L1緩存不可能裝下B的整個一行，因此其不命中率就是1.0，所以rck&crk里層循環(huán)的總不命中率是1.25。
?
? ? 接下來分析ckr&kcr，他們的里層循環(huán)部分都是
? ? ? ? for (r = 0; r < n; r++)
? ? ? ? ? ? E[r][c] += A[r][k]*r;
? ? ? ? 顯然，這里每次對A的訪問仍然是跨行的。因此A的不命中率是1.0。同理，里層循環(huán)還用對E進行應(yīng)用，由于變化的是r，因此E也是跨行訪問的，命中率仍然是1.0，所以ckr&kcr的里層循環(huán)總不命中率是2.0
?
? ? ? ? 最后分析krc&rkc，他們的里層循環(huán)部分都是
? ? ? ? ? ? for (c = 0; c < n; c++)
? ? ? ? ? ? ? ? E[r][c] += m*B[k][c];
? ? ? ? 用同樣分分析方法得出，B和E的訪問不命中率都是0.25，因此里層循環(huán)總不命中率是0.5。
?
? ? ? ? 注意到“每次循環(huán)加載“這一項，其實就是類似A[r][k]或B[r][k]的讀取操作，當(dāng)然都是兩次。而”每次循環(huán)的存儲“則是對E[r][c]的寫入。由于rck&crk版本的里層循環(huán)都由臨時變量sum存儲加值，不會產(chǎn)生對E[r][c]的寫入操作（存儲），因此是0。
?
? ? ? ? 好了，光從不命中率上來看，勝負已分，性能關(guān)系是krc&rkc >?rck&crk >?ckr&kcr，實際情況是不是這樣呢？
? ? ? ? 下面我要進行實際測試數(shù)據(jù)，我的CPU是酷睿2 p8400，查詢參數(shù)：
? ? ? ? http://www.cpu-world.com/CPUs/Core_2/Intel-Core%202%20Duo%20Mobile%20P8400%20AV80577SH0513M.html
查詢到L1高速緩存的資料，如果只考慮單核，數(shù)據(jù)緩存總大小C=32KB，塊大小B=64字節(jié)，相聯(lián)度E=4，組數(shù)就是S = 32KB/(64B*4) = 2^15/2^8 = 2^7 =128，組數(shù)有128組，因此這款CPU的單核L1高速數(shù)據(jù)緩存參數(shù)（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 |
? ? ? ? 如果你不能根據(jù)這些資料熟練的推到并理解（S,E,B,m）= （128,4,64,32）這個結(jié)果，那還是回去看下本章的第（3）節(jié)內(nèi)容，那里有詳細分析。
? ? ? ? 總之我么得到B = 64，和我們之前假設(shè)的32有出入，好吧，那不命中率等于0.25的地方事實上就得換成0.125（理由懶得解釋了，有本章前面的知識鋪墊），那么三個版本的不命中率分別應(yīng)該是1.125、2和0.25。同時，為了增加n使得A、B的一行數(shù)據(jù)足夠大到不能完全裝進L1高速緩存"里的一塊或者說一行"好了，現(xiàn)在塊大小B = 64B，而數(shù)組元素sizeof(double) = 8B，那么n必須大于64B/8B = 8。
? ? 好吧，按照教材的實現(xiàn)方式，n從25到400，以25遞增，對6個版本進行性能測試，計算每次循環(huán)所消耗的周期數(shù)。運行結(jié)果如下：
?
[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?
?
?
? ? ?看著是不是有點暈？那我們就把這組數(shù)據(jù)轉(zhuǎn)換成性能統(tǒng)計圖來分析更加直觀：
?
? ? 這個圖直觀的反映了六個版本的執(zhí)行效率，縱坐標(biāo)表示每次里層循環(huán)所需的CPU周期數(shù)，越高說明耗的周期數(shù)越多，性能也就越差。當(dāng)n超過250時，我驚奇的發(fā)現(xiàn)，性能圖非常鮮明的把六個版本分成三個梯隊，krc&rkc性能最佳，rck&crk其次，ckr&kcr最糟糕，回過頭去看我們對6個版本的總不命中率分析，他們的結(jié)果竟然相同！
? ? 為啥我要說驚人和竟然呢？理論和實際結(jié)論相同很奇怪么？當(dāng)然奇怪了。因為教科書上得出的結(jié)論和我稍有不同，在教材里性能最優(yōu)的是rck&crk！這再次體現(xiàn)了盡信書不如無書。實際測試結(jié)果很有意思。當(dāng)教材作者得出rck&crk性能優(yōu)于krc&rkc的結(jié)論時給出了解釋，認為不命中率并不說明一切，krc&rkc雖然有最少的不命中率，卻有額外的存儲器訪問：krc&rkc的每次里層循環(huán)需要引用E[r][c]\A[r][k]\B[k][c]這三個存儲器，而rck&crk只有引用A[r][k]\B[k][c]兩個存儲器，因此存儲器引用影響了最終性能……而經(jīng)過我的實際測試結(jié)果，是不是可以得出結(jié)論：這個存儲器影響因素消失了？是不是由于CPU更高級，使得存儲器引用的效率得以大幅提升？
? ? 如果學(xué)過CPU內(nèi)部原理應(yīng)該有所了解，CPU在進行運算時，可能會經(jīng)歷取指、譯碼、執(zhí)行、訪存、寫回、更新PC等步驟，而其中的訪存階段可以將數(shù)據(jù)寫入存儲器，或者從存儲器讀取數(shù)據(jù)。而這里面的所謂存儲器，有可能是L1緩存，也有可能是L2緩存甚至有可能是L3、內(nèi)存、硬盤，因此訪存被認為是CPU執(zhí)行指令的過程中最可能耗費更多時間的步驟。因此L1~L3甚至內(nèi)存到硬盤的大小與讀取速度可能直接影響訪存的效率，我只能推斷，我運行測試時的電腦，從L1~L3到內(nèi)存到硬盤可能都全面由于作者的測試電腦，那我能不能嘗試下獲得整個我的電腦存儲器層次的性能對測試結(jié)果的影響呢？這至少能涉及到從L1到內(nèi)存的部分，好的，接下來再理一下我的虛擬機linux配置的配置：
L1：32kB
L2:3072kB
L3：p8400沒有L3
內(nèi)存：1024MB
? ? 也就是說，當(dāng)我對存儲器的引用增加到E\A\B三個時，n的大小將決定我的數(shù)據(jù)能緩存到那一級。比如，我要限制在L1中，那么n最大就不能超過對(32k/(8*3))開方那么多次，也就是n<=37；如果我要限制在L2中，那就是限制n不超過對（3072k/(8*3)）開方那么多次，也就是n<=362;如果同理，要限制在內(nèi)存中，n<=11585，內(nèi)啥，n大于2048時，我的小本本已經(jīng)快跑糊了，要不是果斷shutdown估計可憐p8400就要報銷了，因此我只能分別測試L1和L2，然后呢？n<=37時結(jié)論已經(jīng)有了，我們發(fā)現(xiàn)krc&rkc的性能仍然完爆其他對手，而在75~275之間，rck&crk昂首挺胸；也就是說，在數(shù)據(jù)借助到L2緩存時，L2的本身的性能劣勢使得訪存耗時增加。而當(dāng)大于275時，krc&rkc幾乎已不可戰(zhàn)勝，即便當(dāng)n>362，緩存鐵定到了內(nèi)存級別時，也不能改變這個趨勢。
? ? 好了，為什么分界點不在362而是在275呢？我猜呢，注意哈，是我猜，沒有充足依據(jù)，應(yīng)該是數(shù)據(jù)緩存還要處理其他數(shù)據(jù)，不只是三個矩陣獨享……
?
?
?
?
?
?
?
?
?
總結(jié)
以上是生活随笔為你收集整理的带你学习《深入理解计算机系统》程序性能优化探讨（5）——高速缓存、存储器山与矩阵乘法优化的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： java gettime_Java Ut
- 下一篇： 【系统分析师之路】 第八章 复盘软件测试