# 2017q1 Homework5(matrix) contributed by < `0xff07` > ``` Architecture: x86_64 CPU 作業模式： 32-bit, 64-bit Byte Order: Little Endian CPU(s): 4 On-line CPU(s) list: 0-3 每核心執行緒數：2 每通訊端核心數：2 Socket(s): 1 NUMA 節點： 1 供應商識別號： GenuineIntel CPU 家族： 6 型號： 69 Model name: Intel(R) Core(TM) i5-4258U CPU @ 2.40GHz 製程： 1 CPU MHz： 2899.951 CPU max MHz: 2900.0000 CPU min MHz: 800.0000 BogoMIPS: 4799.89 虛擬： VT-x L1d 快取： 32K L1i 快取： 32K L2 快取： 256K L3 快取： 3072K NUMA node0 CPU(s): 0-3 Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc aperfmperf eagerfpu pni pclmulqdq dtes64 monitor ds_cpl vmx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm epb tpr_shadow vnmi flexpriority ept vpid fsgsbase tsc_adjust bmi1 avx2 smep bmi2 erms invpcid xsaveopt dtherm ida arat pln pts ``` ## Code 這裡 assign 預設是用 4 X 4 陣列，然後把它包成 `struct` ，這樣就可以用 designatied initializer 很漂亮的做一些事，比說在測試的程式碼裡面： ```C algo->assign(&n, (Mat4x4) { .values = { { 1, 2, 3, 4, }, { 5, 6, 7, 8, }, { 1, 2, 3, 4, }, { 5, 6, 7, 8, }, }, }); ``` 不過如果接受的是大小不會事先知道的的矩陣，不知道可不可以用類似的方法。 還有包裝各種方法的手法： ```C MatrixAlgo NaiveMatrixProvider = { .assign = assign, .equal = equal, .mul = mul, }; ``` 然後再測試程式當中就可以這樣使用： ```C MatrixAlgo *algo = matrix_providers[0]; Matrix dst, m, n, fixed; algo->assign(&m, (Mat4x4) { .values = { { 1, 2, 3, 4, }, { 5, 6, 7, 8, }, { 1, 2, 3, 4, }, { 5, 6, 7, 8, }, }, }); ``` 所以只要把不同實作塞進 `matrix_providers` 陣列裡面，測試程式就可以一個迴圈跑完所有的版本。真是太酷了。 缺點大概就是所有東西都是寫死的。 ## 解除大小限制： 這裡希望可以自由的宣告不同大小的矩陣。因為希望這個資料結構可以使用如 `a[i][j]` 這種方法存取，但是如果只有 `malloc` 的話，要存取 $a_{ij}$ 就要用如 `a[col * i + j]` 這樣的方法。這樣看起來不是很方便，所以用下面這樣的資料結構去紀錄二維陣列： ```C float** malltrix (int row, int col) { int idx_size = row * sizeof(float*); int data_size = row * col * sizeof(float*); void* ret = malloc(idx_size + data_size); /* 前 idx_size 個 byte 用來存指向每一列的指摽 */ float** idx = (float**)ret; /* 後面才是放資料 */ float* data = (float*)((char*)ret + idx_size); for(int i = 0; i < row; i++) idx[i] = data + col * i; return ret; } ``` 這樣把回傳的東西轉型成 `float**` 就可以繼續使用 `a[i][j]` 這種方式存取。 ## 計畫：矩陣運算轉成向量運算 想到一種對 Cache 友善的計算方法。假設現在執行矩陣乘法： $$A_{m\times n}B_{n \times l} = C_{ml}$$ 這裡把它寫成這樣： $$ \begin{bmatrix} a_{11} & a_{12} & ... &a_{1n} \\ . & . & . & . \\ . & . & . & . \\ a_{m1} & a_{m2} & . & a_{mn} \end{bmatrix} \begin{bmatrix} \vec{b}_{1} \\ \vec{b}_{2}\\ . \\ \vec{b}_{n} \end{bmatrix}= \begin{bmatrix} \vec{c}_{1} \\ \vec{c}_{2}\\ . \\ \vec{c}_{m} \end{bmatrix} $$ 其中： $$ \vec{c}_{1} = a_{11}\cdot\vec{b}_{1} + a_{12}\cdot\vec{b}_{2}+ a_{13}\cdot\vec{b}_{3} + ... a_{1n}\cdot\vec{b}_{n} $$ <br> 這樣<s>參訪</s> 存取記憶體的方向是： :::warning access 在台灣慣用的翻譯是「存取」，在對岸則是「訪問」 --jserv ::: <br> | $$A$$ | $$B$$ | |--|--| | ![](https://i.imgur.com/MxRwowv.png)|![](https://i.imgur.com/XNdztQk.png)| <br> 可以發現這樣對 cache 應該會很有幫助。而且不管對 prefetch 或是平行化都很方便。 不過在這之前要先寫出評估效能的程式。 ## 中場笑話：~~奇妙~~的 DEBUG 過程？ 本來以為測試沒有 bug, 後來發現測試寫錯，不小心讓矩陣全部都是 0 。在改寫了一下之後，大概花了4個小時在找錯。 找 bug 的過程中， 一直看到 gdb 跑出如下面的內容： ``` 136 for (int k = upper; k < r_col; k++) { (gdb) 137 PRIV(dst)[i][k] += coeff * PRIV(r)[j][k]; (gdb) 136 for (int k = upper; k < r_col; k++) { (gdb) p dump(r) value has been optimized out ``` 所以我以為我誤會了編譯器的什麼地方。而且 avx 的部份一直沒執行到。所以我就做了： 1. 檢查編譯選項 2. 思考我是不是漏了什麼編譯條件 2. 關掉 -O2 3. 再度思考是不是少了什麼編譯條件 4. 插幾個我覺得像 memory barrier 的東西再看看 3. 換個作業系統 (macOS) 跑跑看 4. 檢查編譯器版本 5. 思考是不是我沒有用對編譯器 X3 2. 裝 rr 看看 3. 參考 matrix-multiplication 的實作，看漏了什麼 4. 改用 `_mm256_fmadd_ps()` 指令 5. 用 objdump 檢查東西是不是真的有編進去。很不幸的發現裡面有 AVX 的指令。 後來發現是回圈的邊界少一個 `=` (然後就噴掉一個晚上了 囧) >> 其實會一直懷疑編譯器是因為之前在看 lock free 的影片時也看到了這個影片： {%youtube nXaxk27zwlk %} >> 15:20 的時候有他要評估的程式(push_back()) >> 18:00 的時候有人會說記憶體配置比較花時間，所以他把建立 vector 跟 push_back 一起 benchmark >>然後就發現快了 10 倍 >> 大約 39:20 的時候他會看一下編出來的東西長怎樣，然後發現要 benchmark 的東西通通被編譯器優化砍掉了(囧)所以很名正言順的跑起來快了十倍(這一段超好笑)。 >> 不過裡面看到一些之前沒看過的 perf 的功能還真是精美~ ## stopwatch 其實用起來滿直覺的，只要看下面 API 給了什麼，然後再大概看看實作長怎樣，就大致上可以掌握該怎麼用。 整體概念其實跟之前的用法很類似，按下去之後按停就會算時間。 ## Row Major 實作 採用上面 row major 的實作如下： ```C static bool row_major_mul(Matrix *dst, const Matrix *l, const Matrix *r) { /* FIXME: error hanlding */ int l_row = l->row, l_col = l->col; int r_row = r->row, r_col = r->col; assert(l_col == r->row); dst->priv = matrix_alloc(l_row, r_col); for (int i = 0; i < l_row; i++) { for (int j = 0; j < l_col; j++) { float coeff = PRIV(l)[i][j]; for (int k = 0; k < r_col; k++) PRIV(dst)[i][k] += coeff * PRIV(r)[j][k]; } } return true; } ``` 其實只是換個順序而已。雖然看起來幾乎一樣，但是後面看執行時間時會發現快很多。 ## Row Major + AVX 一樣是用上面的概念，只是加了 AVX： ```C static bool avx_mul(Matrix *dst, const Matrix *l, const Matrix *r) { /* FIXME: error hanlding */ int l_row = l->row, l_col = l->col; int r_row = r->row, r_col = r->col; assert(l_col == r->row); free(dst->priv); dst->priv = matrix_alloc(l_row, r_col); for (int i = 0; i < l_row; i++) { for (int j = 0; j < l_col; j++) { float coeff = PRIV(l)[i][j]; int lower = r_col & 7; int upper = r_col - lower; for (int k = 0; k + 8 <= upper; k+=8) { __m256 ymm0; __m256 ymm1; __m256 ymm2; __m256 ymm3; ymm0 = _mm256_loadu_ps(&(PRIV(r)[j][k])); ymm1 = _mm256_set1_ps(coeff); ymm2 = _mm256_loadu_ps(&(PRIV(dst)[i][k])); ymm3 = _mm256_fmadd_ps(ymm0, ymm1, ymm2); _mm256_storeu_ps((float *)(&(PRIV(dst)[i][k])), ymm3); ![](https://i.imgur.com/nJLMQId.png) } for (int k = upper; k < r_col; k++) { PRIV(dst)[i][k] += coeff * PRIV(r)[j][k]; } } } return true; } ``` 主要是拆成 8 的倍數跟餘數, 然後就可以知道要做幾次向量運算。 這裡 `fmadd` 可以一口氣做完一次乘法跟一次加法。 剩下的部份就乖乖用本來的作法。 >> 也許 mask 的方式載入，應該可以讓程式更好看 & 更快？ ## 比較 對於 800 X 800 的方陣互乘，重複 100 次並比較每次花費的時間： ![](https://i.imgur.com/LootA1E.png) 可以發現幾件有趣的事： 1. naive 波動很大(不知道為什麼？) 2. row major 的運算時間大約是 naive 的 1/3(感覺就像什麼都沒做就變快了) 3. 再加上 avx 就又更快了，大約只有本來的1/6。 ## Perf 的分析 如果把 3 個東西放在同一個程式，然後用 `perf record` 來觀察程式熱點，大致會得到下面的結果： ``` Samples: 15K of event 'cycles:ppp', Event count (approx.): 11360111538 Overhead Command Shared Object Symbol 70.62% test-matrix test-matrix [.] naive_mul 17.91% test-matrix test-matrix [.] row_major_mul 9.56% test-matrix test-matrix [.] avx_mul 0.27% test-matrix [unknown] [k] 0x000000000040119b 0.24% test-matrix [unknown] [k] 0x00000000004011a0 ``` 可以發現 naive 的 overhead 超大，然後 row_major 雖然跟 naive 的程式幾乎一樣，但是 overhead 大幅變少。 然後我想試試看上面那個研討會影片裡的 perf 功能。 ### Call Relation 首先是 perf 也可以紀錄 call graph 資料。不過要先再編譯時加入`-fno-omit-frame-pointer` 這個 flag。然後使用 perf record : ```BAS $ perf record -g ./test-matrix $ perf report -g 'graph,0,caller' ``` >> 中間那個數字是決定要不要過濾某些不重要的函數，愈小的話過濾的愈少。不過我還沒仔細研究。 然後就會發現 perf 裏面多出一些東西： ```BAS Children Self Command Shared Object Symbol ◆ + 99.97% 0.00% test-matrix test-matrix [.] main ▒ + 99.97% 0.00% test-matrix libc-2.23.so [.] __libc_start_main ▒ + 99.97% 0.00% test-matrix [unknown] [k] 0x087e258d4c544155 ▒ + 69.82% 69.64% test-matrix test-matrix [.] naive_mul ▒ + 18.94% 17.86% test-matrix test-matrix [.] row_major_mul ▒ + 10.31% 10.18% test-matrix test-matrix [.] avx_mul ``` 就是左邊的`+` 。如果把他往下拉（按 enter），比如說我拉開 main 的可以看到： ```bas - 99.97% 0.00% test-matrix test-matrix [.] main ▒ - main ▒ + 69.82% naive_mul ▒ + 18.94% row_major_mul ▒ + 10.31% avx_mul ▒ + 0.24% __random_r ▒ + 0.21% fill ▒ + 0.16% __random ▒ + 0.16% create ▒ 0.08% equal ▒ 0.03% rand ▒ 0.01% native_irq_return_iret ▒ + 0.01% __GI___munmap ▒ + 99.97% 0.00% test-matrix libc-2.23.so [.] __libc_start_main ▒ ``` 然後就發現 **他呼叫了什麼函數都跟你說了阿阿阿阿阿阿阿阿！** 這真是超級潮的。不過這還不是更神的，更厲害的在下面。 ### 追蹤 assembly 如果把這東西拉開，找到某個函數, ```bas Children Self Command Shared Object Symbol ▒ - 99.97% 0.00% test-matrix test-matrix [.] main ◆ + main ▒ + 99.97% 0.00% test-matrix libc-2.23.so [.] __libc_start_main ▒ + 99.97% 0.00% test-matrix [unknown] [k] 0x087e258d4c544155 ▒ + 69.82% 69.64% test-matrix test-matrix [.] naive_mul ▒ + 18.94% 17.86% test-matrix test-matrix [.] row_major_mul ▒ - 10.31% 10.18% test-matrix test-matrix [.] avx_mul ▒ + 10.18% 0x87e258d4c544155 ▒ - 0.13% avx_mul ▒ 0.12% avx_mul ▒ 0.01% __memset_avx2 ▒ + 0.34% 0.34% test-matrix [unknown] [k] 0x00000000004011c0 ``` 比如說我移到 `avx_mul` 上面，然後按 `a`，這時候神奇的事就發生了： ![](https://i.imgur.com/rSapBz4.png) 他把 assembly 跟對應的程式碼標出來，而且還耗費資源多的地方上好色了。這真是太太太太厲害了！ (因為用文字很難表現他好用的地方，所以我就貼截圖) 更有趣的是如果移到 branch 類的指令上： ![](https://i.imgur.com/K86hLuq.png) 他就會提示這個 branch 會跳去哪裡。像這裡本來是一個 `assert` ，他就會提示 `assert` 會跳到哪裡去。不過這個 jump 有點長，如果要一直往下拉有點麻煩。但是如果按鍵盤上的右鍵： ![](https://i.imgur.com/FF2T9IM.png) **他就會幫你跳過去了**。這真是太潮了，一個 perf 有這麼多神奇的功能。 --- 另外有幾件事可以繼續討論： - [ ] 如果加上 prefetch 會有什麼效應？這時候 prefetch distance 要怎麼決定 ？ - [ ] 上面的計算方式也可以平行化。如果不用向量運算，改用平行化，哪個效能會比較好？ - [ ] 如果同時做平行，又用向量運算，這樣 YMM 暫存器算不算 shared resoure ? 如果是的話，要怎麼樣做才可以避免 race condition ? - [ ] loop unrolling