Optimizing Matrix Transpose
Matrix Transpose 還算一個(gè)常見(jiàn)的問(wèn)題. cache lab handout 也寫(xiě)出了最常見(jiàn)和 easy 的解法，但這并不是最memory friendly的。
先來(lái)讀懂文檔：
Matrices are stored in memory in a row major order.
If the entire matrix can’t fit in the cache, then after the cache is full with all the elements it can load. The next elements will evict the existing elements of the cache.
Example:- 4x4 Matrix of integers and a 32 byte cache.
The third row will evict the first row!
前兩個(gè)row 4 * 2 * 4 = 32 byte， int 4 byte，所以滿(mǎn)了自然就會(huì)擠出去。
The first row of Matrix A evicts the first row of Matrix B
Caches are memory aligned.
Matrix A and B are stored in memory at addresses such that both the first elements align to the same place in cache!
Diagonal elements evict each other.
這看起來(lái)很難的樣子，同時(shí)也解答了為什么在給的例子中，我們這樣來(lái)做計(jì)算：
void trans(int M, int N, int A[N][M], int B[M][N])
{
int i, j, tmp;
for (i = 0; i < N; i++) {
for (j = 0; j < M; j++) {
tmp = A[i][j];
B[j][i] = tmp;
}
}
}
就是引入一個(gè)temp。 A 是B 么？對(duì)于 M = N 的（此處存疑）因?yàn)檎f(shuō)他們?cè)趦?nèi)存中相同的位置？
再看 cache 的規(guī)格：
E = 1, S = 32, B = 32, total = 1 kilobyte
然后針對(duì)于 M = 32, N = 32, 一開(kāi)始 load 的狀況是：
A[0][0] ... A[0][7]
A[0][8] ... A[0][15]
A[0][16]... A[0][23]
A[0][24]... A[0][31]
.
.
A[7][24] ... A[7][31]
然后比如我們加載進(jìn) B 的哪一行之后，它會(huì) evict 其中一行。
然后再看提示：
- 最多12個(gè)int variable
- array A 我們雖然不能改變，但是 array B 我們可以隨意
然我們從最簡(jiǎn)單的實(shí)驗(yàn)開(kāi)始吧，先就按照寫(xiě)的例子來(lái)填 trans_submit，結(jié)果：
Function 0 (2 total)
Step 1: Validating and generating memory traces
Step 2: Evaluating performance (s=5, E=1, b=5)
func 0 (Transpose submission): hits:870, misses:1183, evictions:1151
Function 1 (2 total)
Step 1: Validating and generating memory traces
Step 2: Evaluating performance (s=5, E=1, b=5)
func 1 (Simple row-wise scan transpose): hits:870, misses:1183, evictions:1151
Summary for official submission (func 0): correctness=1 misses=1183
TEST_TRANS_RESULTS=1:1183
先來(lái)測(cè)試次簡(jiǎn)單的：既然我們?cè)试S12個(gè)local variable，每個(gè)block可以存8個(gè)int，那么我們每次把8個(gè)記住，然后來(lái)填B， ??，然而這樣寫(xiě)出來(lái)的代碼為什么misses比之前更多了呢？？？
思路二：
t0 = A[0][0];
t1 = A[0][1];
t2 = A[0][2];
t3 = A[0][3];
t4 = A[1][1];
t5 = A[1][2];
t6 = A[1][3];
t7 = A[2][2];
t8 = A[2][3];
t9 = A[3][3];
B[0][0] = t0;
for(j = 1; j < 8; j++)
B[0][j] = A[j][0];
B[1][0] = t1;
B[1][1] = t4;
for(j = 2; j < 8; j++)
B[1][j] = A[j][1];
B[2][0] = t2;
B[2][1] = t5;
B[2][2] = t7;
for(j = 3; j < 8; j++)
B[2][j] = A[j][2];
B[3][0] = t3;
B[3][1] = t6;
B[3][2] = t8;
B[3][3] = t9;
for(j = 4; j < 8; j++)
B[3][j] = A[j][3];
這樣操作下來(lái)之后明顯會(huì)減少一些miss。 我們 A[0~3] 被替換為 B[0~3], 剩余 A[4~7]，注意A[7]就一部分在cache里面,這里還有疑問(wèn)，就是 local variable 也是在cache里面么？知道它一定在stack上，但它也是占據(jù)了cache的空間么？（存疑
這個(gè)時(shí)候如果我們?cè)賮?lái)加載 A[8]， 是否我們會(huì)將B evict，然后我們可以用類(lèi)似的方法來(lái)處理B[0][8]~B[0][15], B[1][8]~B[1][15], B[2][8]~B[2][15]...
不過(guò)即使這樣做，A[4] ~ A[7]的填也麻煩，按照這樣的思路拼湊出來(lái)的代碼可以做到 misses 700 some what， 但是這個(gè)完全達(dá)不到要求，而且很僵化，我只考慮了 32 * 32
所以來(lái)讀一讀關(guān)于matrix transpose的 Q & A.
首先，發(fā)現(xiàn)有這個(gè)：
A Cache Efficient Matrix Transpose Program?
在 wikipedia 的這篇 Cache-oblivious algorithm 也提到了 matrix transpose的例子，甚至還畫(huà)了示意圖。
先看問(wèn)題中的代碼：
for( int i = 0; i < n; i++ )
for( int j = 0; j < n; j++ )
destination[j+i*n] = source[i+j*n];
這里假設(shè) n = 4, 那么我們依次：
- i = 0
- dst[0] = src[0]
- dst[1] = src[4]
- dst[2] = src[8]
- dst[3] = src[12]
- i = 1
- dst[4] = src[1]
- dst[5] = src[5]
- dst[6] = src[9]
- dst[7] = src[13]
無(wú)需繼續(xù)驗(yàn)證，如果 dst 和 src 都被裝入了 cache 中，那么這將是非常cache 友好的。
如果不是方形矩陣，維度是 M * N ,我們也可以寫(xiě)出類(lèi)似代碼：
for (int i = 0 ; i < N; i++)
for (int j = 0; j < M; j++)
dst[j + i * N] = src[i + j * M]
然后來(lái)嘗試明白這個(gè)top voted answer給的算法。所以其實(shí)還是要先回到課上給的matrix 相乘的例子。
for (int i = 0; i < n; i += blocksize) {
for (int j = 0; j < n; j += blocksize) {
// transpose the block beginning at [i,j]
for (int k = i; k < i + blocksize; ++k) {
for (int l = j; l < j + blocksize; ++l) {
dst[k + l*n] = src[l + k*n];
}
}
}
}
根據(jù)**改寫(xiě)的代碼：
void transpose_submit(int M, int N, int A[N][M], int B[M][N]){
int blocksize = 8;
int k, l;
int i, j;
for (i = 0; i < M; i += blocksize)
{
for (j = 0; j < N; j += blocksize)
{
for (k = i; k < i + blocksize && k < M; ++k)
{
for (l = j; l < j + blocksize && l < N; ++l)
{
B[k][l] = A[l][k];
}
}
}
}
}
結(jié)果：
Function 0 (2 total)
Step 1: Validating and generating memory traces
Step 2: Evaluating performance (s=5, E=1, b=5)
func 0 (Transpose submission): hits:1710, misses:343, evictions:311
Function 1 (2 total)
Step 1: Validating and generating memory traces
Step 2: Evaluating performance (s=5, E=1, b=5)
func 1 (Simple row-wise scan transpose): hits:870, misses:1183, evictions:1151
Summary for official submission (func 0): correctness=1 misses=343
TEST_TRANS_RESULTS=1:343
比 1183 好很多了。但是沒(méi)有達(dá)到 doc 的要求， 同時(shí)，發(fā)現(xiàn) blocksize = 8 對(duì)于 M = 64, N = 64 的時(shí)候完全沒(méi)有改善
- 32×32: 8 points if m<300,0 points if m>600
- 64×64:8 points if m<1,300,0 points if m>2,000
- 61×67: 10 points if m<2,000,0 points if m>3,000
blocksize = 4 32* 32 ：misses:487
blocksize = 4 64 * 64 ：misses:1891
blocksize = 4 61 * 67 ： misses:2425
當(dāng)然這里還有一些trick，如果是硬朝著減少miss方向努力的話(huà)，針對(duì)不同的 M 和 N 調(diào)整blocksize，然后如果是對(duì)角線(xiàn)我們就不用去改變它。
可以參考cachelab
實(shí)際上看這段代碼，情況會(huì)更清楚，突然，一頁(yè)slides make sense了：
假設(shè) block size 是 8 byte
Matrix A
[1] [2] 3 4
5 6 7 8
9 10 11 12
13 14 15 16
此時(shí) matrix B :
[1]
[2]
接下來(lái)處理 5 和 6 更好，因?yàn)?1/2 旁邊已經(jīng)加載進(jìn)去， 如果處理 3/4的話(huà)我們都要重新加載，更多 miss， 這也是我們?yōu)槭裁从胋lock來(lái)處理的原因吧。
實(shí)際上我們看別人的代碼，也更清晰，它的兩層 loop 針對(duì) A 是先走第一列的 block， 然后走第二列的block：
1 5 . .
2 6 . .
3 7 . .
4 8 . .
對(duì)于block 1 的內(nèi)部，就是先橫著走啦，因?yàn)檫@個(gè)block已經(jīng)被加載入了。。。
按照github 上的代碼，32 x 32 和 61 x 67 能達(dá)到要求，即使已經(jīng)對(duì) 64 x 64 有了這么多改善，也還是 misses:1635。
這可真不容易啊
