matlab 多重循环在最外层加断点_循环优化之循环分块（loop tiling）
引言
編譯器里的循環優化有兩個重要的目標，一是提高局部性，二是提高并行性，loop tiling是提高數據局部性最重要的優化之一，是傳統編譯器和深度編譯器考慮的重中之重，我們今天來看看如何做loop tiling（循環分塊）（aka cache blocking）。
loop tiling的basic idea很簡單，就是一個cache line現在被用過以后，后面什么時候還會被用，但是按循環默認的執行方式，可能到下次再被用到的時候已經被evict了。于是我們就把循環怎么重排一下，使得一個cache line在被evict之前就被再次使用。
在該文中我們討論一個簡單的tiling方法，以及通過計算cache miss來估算tiling的作用。
例程一
我們從一個簡單的例子看起，為了更好的可讀性，我們使用Fortran語法（column major）：
DO I = 1, NDO J = 1, MA(I) = A(I) + B(J) END DO END DO首先我們分析一下該循環的重用模式（reuse pattern），可以看到A(I)在每一輪J循環中被重用、B(J)在每一輪I循環中會被重用。當前由于I是外層循環，所以當B(J)在J層循環中的跨度太大時（無法fit in the cache），則在被下一次I重用時數據已被清出緩存。譬如我們假設M和N是很大的數（大數組），所以對于每一輪I循環，等到B(M)被訪問時，B(1)、B(2)等已經被清出緩存了。
分析完了重用模式，我們現在計算cache miss。假設每個cache line能容納b個數組元素，associativity是fully associative，則可計算出A的cache miss次數是N/b，B的cache miss次數是N*M/b。
現在我們來看看什么是tiling，以及如何通過做tiling減少cache miss。
我們先考慮tile J層循環，將B的訪問模式變成如下，令tile size為T，T一般是b的倍數，也足夠小，可以認為N,M >> b,T，T的取值應該使得當B(T)被訪問時B(1)也依然還在緩存中：
I=1: B(1), B(2), B(3) ... B(T) I=2: B(1), B(2), B(3) ... B(T) I=3: B(1), B(2), B(3) ... B(T) ... I=N: B(1), B(2), B(3) ... B(T)I=1: B(T+1), B(T+2), B(T+3) ... B(2T) I=2: B(T+1), B(T+2), B(T+3) ... B(2T) I=3: B(T+1), B(T+2), B(T+3) ... B(2T) ... I=N: B(T+1), B(T+2), B(T+3) ... B(2T)小燈泡：tile的本質當總的數據量無法fit in cache的時候，把數據分成一個一個tile，讓每一個tile的數據fit in the cache。所以tiling一般從最內層循環開始（tile外層循環的話一個tile就包括整個內層循環，這樣的tile太大）。此時循環變成：
DO J = 1, M, TDO I = 1, NDO jj = J, min(J+T-1, M)A(I) = A(I) + B(jj)END DOEND DO END DO注意當我們tile J層循環的時候，tile以后J層循環就跑到了外層（loop interchange），這樣又會影響A的緩存命中率。該變換又稱為strip-mine-and-interchange。我們可以計算此時A的cache miss次數是(M/T) * (N/b)，B的cache miss次數是T/b * M/T = M/b（每一輪J層循環miss一次），所以總共的cache miss是MN/(bT) + M/b，tile之前是N/b+N*M/b。所以在M和N的大小也相當的情況下做tiling大約能將cache miss縮小T倍。
思考題1：為什么分塊以后還要交換循環I和J？循環交換是否總是可行的？那要是我們選擇tile I層循環呢？我們看看會發生什么。要tile I層循環首先將I、J循環調換：
DO J = 1, MDO I = 1, NA(I) = A(I) + B(J) END DO END DO依然假設tile size=T，訪問模式將變為：
J=1: A(1) B(1), A(2) B(1), A(3) B(1) ... A(T) B(1) J=2: A(1) B(2), A(2) B(2), A(3) B(2) ... A(T) B(2) ... J=M: A(1) B(M), A(2) B(M), A(3) B(M) ... A(T) B(M)J=1: A(T+1) B(1), A(T+2) B(1), A(T+3) B(1) ... A(2T) B(1) J=2: A(T+1) B(2), A(T+2) B(2), A(T+3) B(2) ... A(2T) B(2) ... J=M:A(T+1) B(M), A(T+2) B(M), A(T+3) B(M) ... A(2T) B(M) ...此時循環變成：
DO I = 1, N, TDO J = 1, MDO ii = I, min(I+T-1, N)A(ii) = A(ii) + B(J)END DOEND DO END DO此時A的cache miss次數是T/b * N/T，B的cache miss次數是M/b*(N/T)，加起來就是(MN)/(bT)+N/b，而tile J層循環得到的cache miss是(MN)* (bT)+M/b，所以對于該循環到底應該tile I層循環還是J層循環更好其實取決于M和N的相當大小。
思考題2：是否可以既tile I層循環也tile J層循環？如果可以，這樣做是否具有更好的局部性？例程二
DO J = 1, MDO I = 1, ND(I) = D(I) + B(I, J)END DO END DO注：fortran是column major layout，B(I, J)是訪問二維數組，相當于C語言中的B[I][J]。不管是column major還是row major不影響tiling的原理。注：對于多維數組我們在計算cache miss時假設數據的存儲是連續的。
該循環的cache miss是M*N/b + M*N/b（D和B的miss都是M*N/b）。
還是用之前的strip-mine-and-interchange方法來tile內層循環I：
DO I = 1, N, TDO J = 1, MDO ii = I, min(I+T-1, N)D(ii) = D(ii) + B(ii, J)END DOEND DO END DOtiling loop I 之后的cache miss是N/T * T/b + T/b * M * N/T（J層循環每一輪一開始都會miss B）= N/b + NM/b。
再考慮tile J層循環（需要首先調換I、J層循環）：
DO J = 1, M, TDO I = 1, NDO jj = J, min(J+T-1, M)D(I) = D(I) + B(I, jj)END DOEND DO END DO此時的cache miss是N/b * M/T + T*(N/b)*M/T = NM/(bT) + MN/b。跟tiling I循環相比，由于M一般遠大于T，所以會有更高的cache miss。造成這個區別的本質原因在于D(I)在J層的重用沒有被exploited。
例程三
我們再來看一個更有趣的例子：矩陣轉置。
DO J = 1, NDO I = 1, NA(I, J) = B(J, I)ENDDO ENDDO這個例子之所以有趣是因為，無論是否做交換I、J循環（交換是合法的），cache miss都是N*N/b + N*N。下面是交換I、J循環后的情況：
DO I = 1, NDO J = 1, NA(I, J) = B(J, I)ENDDO ENDDO如果I是內層循環，則A的miss是N*N/b，而B每次都是miss，所以一共miss N*N次。如果J是內層循環，則是A每次都是miss，而B miss N*N/b次。不管哪種情況，都有一個沒有achieve locality。
那有什么辦法讓A和B都能achieve locality嗎？答案是可以，請看如下的tiling：
DO I = 1, N, TDO J = 1, N, TDO ii = I, min(I+T-1, N)DO jj = J, min(J+T-1, N)A(ii, jj) = B(jj, ii)ENDDOENDDOENDDO ENDDO這里的精髓在于tiling不僅可以利用temporal locality，也能用上spatial locality。對于A(ii, jj)，雖然最內層循環是jj（沒有locality），但是當內層循環jj完了以后開始下一輪ii的時候，因為tiling使得A(ii, jj)、A(ii, jj+1)依然會在緩存，于是這時的A(ii+1, jj)、A(ii+1, jj+1)等就能命中緩存（spatial locality）。
下面是各層cache miss的breakdown（外層的cache miss包含了內層的）：
| jj | T | T/b |
| ii | T * T/b | T/b * T |
| J | T * T/b * N/T = NT/b | T/b * T * N/T = NT/b |
| I | NT/b * N/T = N*N/b | NT/b * N/T = N*N/b |
總結
循環分塊的本質是在外層循環挖掘內層的temporal和spatial reuse。基本的方法是strip-mine-and-interchange。以下是takeaway：
- 循環分塊是什么？如何做？
- 對于不同的循環層分塊有什么不同的結果？
- 如何計算cache miss？
想了解更復雜的循環如何做tiling（譬如經典的矩陣乘法）以及tiling的合法性請看本系列第二講：
https://zhuanlan.zhihu.com/p/301905385?zhuanlan.zhihu.com總結
以上是生活随笔為你收集整理的matlab 多重循环在最外层加断点_循环优化之循环分块（loop tiling）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： html细节
- 下一篇： Python小白的数学建模课-15.图论