# 2017q1 Homework5 (matrix) contributed by < `tina0405` > >>中英文字間記得以空白隔開 >>[name=課程助教][color=red] >>不好意思，已更正。 >>[name=tina0405] ## 開發環境 ~~~ tina@tina-X550VB:~$ lscpu 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 型號： 58 Model name: Intel(R) Core(TM) i5-3230M CPU @ 2.60GHz 製程： 9 CPU MHz： 1270.750 CPU max MHz: 3200.0000 CPU min MHz: 1200.0000 BogoMIPS: 5188.47 虛擬： VT-x L1d 快取： 32K L1i 快取： 32K L2 快取： 256K L3 快取： 3072K NUMA node0 CPU(s): 0-3 ~~~ ## 作業要求1 * ### 詳細閱讀 [Matrix Multiplication using SIMD](https://hackmd.io/s/Hk-llHEyx)： * 搭配[影片](https://www.youtube.com/watch?v=3rU6BX7w8Tk&index=1&list=PLKT8ER2pEV3umVSMwd06LY_eSIX-DnU6A)講解的整理： * 一般 C 語言是 row-major 所以再讀取 row 的第一個元素時，會主動去抓第二個元素放進 cache ，就不會有 cache misses 的問題。 >> 此時我們要解決的就是 column 方面的問題，因為它並不會主動抓取下一個元素進去 cache ， 解決 cache misses。 * 把大矩陣拆成 4x4 的小矩陣運算，有加快做用，因為它讓column 方面的 cache 變得有作用一點。 * 將 4x4 的小矩陣做轉置再相乘相加，就有向量 dot 的效果。 * 在 AVX 版本上的實作把原本操做的 128 改為 256 bit，因此改為 8x8 的矩陣加快速度。 * 參考[Strassen algorithm](https://en.wikipedia.org/wiki/Strassen_algorithm) 一般二維矩陣運算需要 8 次乘，4 次加： ![](https://i.imgur.com/AIIfeC9.png) 而經過 strassen algorithm 轉換，做了7次的乘法，18 次加法: ![](https://i.imgur.com/M3W6m0J.png) * 在 SSE 下 malloc 和 calloc 會在內建型態下對齊好。 * 加減搭配 SIMD 利用並行和對齊等特性，證明效能增加。 * 亂數夠亂，構均勻嗎？利用亂數獨立度驗證 (RUNS-UP-AND-DOWN TEST) ## 作業要求2 * ### 嘗試重現實驗的 [github](https://github.com/tina0405/Matrix-multiplication)： 一開始先把原版 naive 複製下來，參考[內容所附github](https://github.com/0140454/matrix-multiplication/commit/4115be2de09dc95f823b41a90bd789c8cd9dbeed) ~~~clike= void naive_multiply(int *src1, int *src2, int *dst, int src1_w, int src1_h, int src2_w, int src2_h) { for (int i = 0; i < src1_h; ++i) { for (int j = 0; j < src2_w; ++j) { dst[i * src2_w + j] = 0; for (int k = 0; k < src2_h; ++k) { dst[i * src2_w + j] += src1[i * src1_w + k] * src2[k * src2_w + j]; } } } } ~~~ * naive 版本： * 是先設定出大小為 16 的矩陣，然後我們將他想成 4X4 做矩陣的乘法，以 i=0， j=0的情況下， k 從 1 到 4 的過程： >>就是把 scr1[0]*scr2[0] 存回 dst[0] 再加 scr1[1]*scr2[4] 再加 scr1[2]*scr2[8] 再加 scr1[3]*scr2[12] 再加 scr1[4]*scr2[16] 如此一來就得到dst[0]的值，其他以公式類推： ![](https://i.imgur.com/xNtueiN.png) * 先將矩陣擴充成 1024 就是 32x32，但這邊的 rand() 我想考慮一下之後替換，因為這篇實驗最後有提到亂數夠不夠亂的問題： * #### [程式碼](https://github.com/tina0405/Matrix-multiplication/commit/50bb1566ee58f784862e7)： ~~~ clike= #define TEST_W 1024 #define TEST_H 1024 ... int test_src1[TEST_W]; int test_src2[TEST_H]; for(int w=0;w<TEST_W;w++) { test_src1[w]=rand(); test_src2[w]=rand(); } int testout[TEST_H]; naive_multiply(test_src1, test_src2, testout, 32, 32, 32, 32); for (int y = 0; y <32 ; y++) { for (int x = 0; x < 32; x++) printf(" %2d", testout[y * 32 + x]); printf("\n"); } ... ~~~ * 感覺我這想不太好用，應該要可以讓使用者選擇 nxn 之後寫 makefile 再加好了，如此一來用來驗證矩陣乘法的 4x4 矩陣也可以一直存在。 * 還有我發現一件事，我根本沒計算時間，參考之前的 [prefetcher](https://github.com/tina0405/prefetcher) 計算時間的方式，還有這篇[計算程式碼時間](https://starforcefield.wordpress.com/2012/08/06/c%E8%AA%9E%E8%A8%80%EF%BC%9A%E4%BD%BF%E7%94%A8clock_gettime%E8%A8%88%E7%AE%97%E7%A8%8B%E5%BC%8F%E7%A2%BC%E7%9A%84%E6%99%82%E9%96%93%E9%9C%80%E6%B1%82/)。 clock_gettime函式原型為： ~~~ int clock_gettime(clockid_t clk_id, struct timespec *tp); ~~~ struct timespec的宣告如下： ~~~ struct timespec { time_t tv_sec; /* seconds */ long tv_nsec; /* nanoseconds */ }; ~~~ 所以我們在計算時間時才可以使用類似這種： ~~~ t1.tv_sec t1.tv_nsec ~~~ * 跑完之後發現有一個問題，為什麼明明是 rand()，每次跑出來的執行檔矩陣乘法的答案都一樣呢？ 在這裡我去翻了[prefetcher](https://github.com/tina0405/prefetcher)裡面對亂數的處理，發現了我忘了使用 srand()，參考[亂數的使用](http://dhcp.tcgs.tc.edu.tw/c/p005.htm)。 * 原因： ~~~ 因為亂數都是由上一個決定下一個，而開始一樣過程也會一樣，因此我們就可以藉由 srand（）改變 seed 且由當下時間決定 time() ，如此一來就不會重複了喔！ ~~~ * 嘗試使用 SSE 去拆解 column（每4個一組） * 證明有效解決 row-major 問題： 我目前想設計的實驗是以之前，[轉置矩陣](https://www.randombit.net/bitbashing/2009/10/08/integer_matrix_transpose_in_sse2.html)為參考。 ![](https://i.imgur.com/tTKZ6OP.png) * 配合 [sse2函式] (http://blog.csdn.net/fengbingchun/article/details/18460199)下去做重組 * 在裡頭剛好有看到相乘的用法，或許我可以利用乘法 A[ ] X B[ ] 的一橫乘一直，把 B[ ] 矩陣做轉置，變成兩個矩陣橫向和橫向相乘。 * 後來在做的時候，發現根本沒必要做轉置，只要一個抓直的，一個抓橫的就好，目前以4X4做驗證，但是還沒找到錯誤，答案並不相符。 ~~~ clike= for (int x = 0; x < src2_h; x += 4) { for (int y = 0; y < src2_w; y += 4) { __m128i I0 = _mm_loadu_si128((__m128i *)(src2 + (y + 0) * src2_w + x)); __m128i I1 = _mm_loadu_si128((__m128i *)(src2 + (y + 1) * src2_w + x)); __m128i I2 = _mm_loadu_si128((__m128i *)(src2 + (y + 2) * src2_w + x)); __m128i I3 = _mm_loadu_si128((__m128i *)(src2 + (y + 3) * src2_w + x)); __m128i I4 = _mm_loadu_si128((__m128i *)(src1 + y * src1_h + x+0)); __m128i I5 = _mm_loadu_si128((__m128i *)(src1 + y * src1_h + x+1)); __m128i I6 = _mm_loadu_si128((__m128i *)(src1 + y * src1_h + x+2)); __m128i I7 = _mm_loadu_si128((__m128i *)(src1 + y * src1_h + x+3)); __m128i T4 = _mm_mullo_epi16(I4,I0); __m128i T5 = _mm_mullo_epi16(I5,I1); __m128i T6 = _mm_mullo_epi16(I6,I2); __m128i T7 = _mm_mullo_epi16(I7,I3); _mm_storeu_si128((__m128i *)(dst + ((x + 0) * src1_h) + y), T4); _mm_storeu_si128((__m128i *)(dst + ((x + 1) * src1_h) + y), T5); _mm_storeu_si128((__m128i *)(dst + ((x + 2) * src1_h) + y), T6); _mm_storeu_si128((__m128i *)(dst + ((x + 3) * src1_h) + y), T7); } } ~~~ * 想說不知道錯再哪就先全部印出來： ~~~ tina@tina-X550VB:~/Matrix-multiplication$ ./main //SCR1矩陣 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 //印出I4，I5，I6，I7 0 1 2 3 1 2 3 4 2 3 4 5 4 5 6 7 ~~~ * 根本完全不一樣，已修正 scr1 ~~~ tina@tina-X550VB:~/Matrix-multiplication$ ./main //SCR2 矩陣 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 //印出I0，I1，I2，I3 16 20 24 28 17 21 25 29 18 22 26 30 19 23 27 31 ~~~ * SCR2 轉置成功 * 修改過後的程式碼，但乘法目前還是失敗，感覺是我的 SSE 乘法函式用錯，等等繼續看： ~~~ clike= for (int x = 0; x < src2_h; x += 4) { for (int y = 0; y < src2_w; y += 4) { __m128i I0 = _mm_loadu_si128((__m128i *)(src2 + (y + 0) * src2_w + x)); __m128i I1 = _mm_loadu_si128((__m128i *)(src2 + (y + 1) * src2_w + x)); __m128i I2 = _mm_loadu_si128((__m128i *)(src2 + (y + 2) * src2_w + x)); __m128i I3 = _mm_loadu_si128((__m128i *)(src2 + (y + 3) * src2_w + x)); __m128i I4 = _mm_loadu_si128((__m128i *)(src1 +0*src1_w)); __m128i I5 = _mm_loadu_si128((__m128i *)(src1 +1*src1_w)); __m128i I6 = _mm_loadu_si128((__m128i *)(src1 +2*src1_w)); __m128i I7 = _mm_loadu_si128((__m128i *)(src1 +3*src1_w)); __m128i T0 = _mm_unpacklo_epi32(I0, I1); __m128i T1 = _mm_unpacklo_epi32(I2, I3); __m128i T2 = _mm_unpackhi_epi32(I0, I1); __m128i T3 = _mm_unpackhi_epi32(I2, I3); I0 = _mm_unpacklo_epi64(T0, T1); I1 = _mm_unpackhi_epi64(T0, T1); I2 = _mm_unpacklo_epi64(T2, T3); I3 = _mm_unpackhi_epi64(T2, T3); __m128i T4 = _mm_mul_epu32(I4,I0); __m128i T5 = _mm_mul_epu32(I5,I1); __m128i T6 = _mm_mul_epu32(I6,I2); __m128i T7 = _mm_mul_epu32(I7,I3); _mm_storeu_si128((__m128i *)(dst + ((x + 0) * 4) + y), T4); _mm_storeu_si128((__m128i *)(dst + ((x + 1) * 4) + y), T5); _mm_storeu_si128((__m128i *)(dst + ((x + 2) * 4) + y), T6); _mm_storeu_si128((__m128i *)(dst + ((x + 3) * 4) + y), T7); } ~~~ * 繼續印出來 T4，T5，T6，T7，就發現問題了！ ~~~ 0 0 48 0 68 0 150 0 144 0 260 0 228 0 378 0 ~~~ * 函式中：（所以只有做 _A0*_B0 和 _A2*_B2 的乘法） >>請使用繁體中文術語 >>[name=課程助教][color=red] ~~~ //回傳一個__m128i的暫存器，r0=_A0*_B0, r1=_A2*_B2 extern __m128i _mm_mul_epu32(__m128i _A, __m128i _B); ~~~ * 改使用 ~~~ //返回一个__m128i的寄存器，它含8个有符号或无符号16bit的整数，分别为_A和_B 对应位置的16bit //有符号或无符号整数相乘结果的低16bit数据，即ri=(_Ai*_Bi)[15:0] (r0=(_A0*_B0)[15:0], //r1=(_A1*_B1)[15:0] ... r7=(_A7*_B7)[15:0]) extern __m128i _mm_mullo_epi16(__m128i _A, __m128i _B); ~~~ * 結果，是我要的向量乘法但矩陣乘完還得全部加起來 ~~~ 0 20 48 84 68 105 150 203 144 198 260 330 228 299 378 465 ~~~ * 但 SSE 裡並沒有 A0+A1+A1+A4 該不會得一個個拿出來存入暫存吧！ * 後來想一想不如我就把乘好的矩陣再轉置就可以直的相加。 ~~~ 舉例： 4X4矩陣 0 20 48 84// 加起來放置 A[0] 0 21 50 87// 加起來放置 A[1] 0 22 52 90// 加起來放置 A[2] 0 23 54 93// 加起來放置 A[3] ... 直到 A[16] ~~~ * 結果就成功了，但在 4X4 矩陣中與 naive 的速度差不多。 * [SSE 程式碼](https://github.com/tina0405/Matrix-multiplication/commit/6279fe6774c8f1552da94787f513e49dd765b373)： ~~~ clike= for (int x = 0; x < src2_h; x += 4) { for (int y = 0; y < src2_w; y += 4) { __m128i I0 = _mm_loadu_si128((__m128i *)(src2 + (y + 0) * src2_w + x)); __m128i I1 = _mm_loadu_si128((__m128i *)(src2 + (y + 1) * src2_w + x)); __m128i I2 = _mm_loadu_si128((__m128i *)(src2 + (y + 2) * src2_w + x)); __m128i I3 = _mm_loadu_si128((__m128i *)(src2 + (y + 3) * src2_w + x)); __m128i I4 = _mm_loadu_si128((__m128i *)(src1 +0*src1_w)); __m128i I5 = _mm_loadu_si128((__m128i *)(src1 +1*src1_w)); __m128i I6 = _mm_loadu_si128((__m128i *)(src1 +2*src1_w)); __m128i I7 = _mm_loadu_si128((__m128i *)(src1 +3*src1_w)); __m128i T0 = _mm_unpacklo_epi32(I0, I1); __m128i T1 = _mm_unpacklo_epi32(I2, I3); __m128i T2 = _mm_unpackhi_epi32(I0, I1); __m128i T3 = _mm_unpackhi_epi32(I2, I3); I0 = _mm_unpacklo_epi64(T0, T1); I1 = _mm_unpackhi_epi64(T0, T1); I2 = _mm_unpacklo_epi64(T2, T3); I3 = _mm_unpackhi_epi64(T2, T3); /****** the first row ******/ __m128i T40 = _mm_mullo_epi16(I4,I0); __m128i T41 = _mm_mullo_epi16(I4,I1); __m128i T42 = _mm_mullo_epi16(I4,I2); __m128i T43 = _mm_mullo_epi16(I4,I3); __m128i T00 = _mm_unpacklo_epi32(T40, T41); __m128i T10 = _mm_unpacklo_epi32(T42, T43); __m128i T20 = _mm_unpackhi_epi32(T40, T41); __m128i T30 = _mm_unpackhi_epi32(T42, T43); T40 = _mm_unpacklo_epi64(T00, T10); T41 = _mm_unpackhi_epi64(T00, T10); T42 = _mm_unpacklo_epi64(T20, T30); T43 = _mm_unpackhi_epi64(T20, T30); __m128i A0 = _mm_add_epi64(T40, T41); __m128i A1 = _mm_add_epi64(T42, T43); __m128i AA = _mm_add_epi64(A0, A1); /****** the second row ******/ __m128i T50 = _mm_mullo_epi16(I5,I0); __m128i T51 = _mm_mullo_epi16(I5,I1); __m128i T52 = _mm_mullo_epi16(I5,I2); __m128i T53 = _mm_mullo_epi16(I5,I3); __m128i T01 = _mm_unpacklo_epi32(T50, T51); __m128i T11 = _mm_unpacklo_epi32(T52, T53); __m128i T21 = _mm_unpackhi_epi32(T50, T51); __m128i T31 = _mm_unpackhi_epi32(T52, T53); T50 = _mm_unpacklo_epi64(T01, T11); T51 = _mm_unpackhi_epi64(T01, T11); T52 = _mm_unpacklo_epi64(T21, T31); T53 = _mm_unpackhi_epi64(T21, T31); __m128i B0 = _mm_add_epi64(T50, T51); __m128i B1 = _mm_add_epi64(T52, T53); __m128i BB = _mm_add_epi64(B0, B1); /****** the third row *****/ __m128i T60 = _mm_mullo_epi16(I6,I0); __m128i T61 = _mm_mullo_epi16(I6,I1); __m128i T62 = _mm_mullo_epi16(I6,I2); __m128i T63 = _mm_mullo_epi16(I6,I3); __m128i T02 = _mm_unpacklo_epi32(T60, T61); __m128i T12 = _mm_unpacklo_epi32(T62, T63); __m128i T22 = _mm_unpackhi_epi32(T60, T61); __m128i T32 = _mm_unpackhi_epi32(T62, T63); T60 = _mm_unpacklo_epi64(T02, T12); T61 = _mm_unpackhi_epi64(T02, T12); T62 = _mm_unpacklo_epi64(T22, T32); T63 = _mm_unpackhi_epi64(T22, T32); __m128i C0 = _mm_add_epi64(T60, T61); __m128i C1 = _mm_add_epi64(T62, T63); __m128i CC = _mm_add_epi64(C0, C1); /****** the forth row ******/ __m128i T70 = _mm_mullo_epi16(I7,I0); __m128i T71 = _mm_mullo_epi16(I7,I1); __m128i T72 = _mm_mullo_epi16(I7,I2); __m128i T73 = _mm_mullo_epi16(I7,I3); __m128i T03 = _mm_unpacklo_epi32(T70, T71); __m128i T13 = _mm_unpacklo_epi32(T72, T73); __m128i T23 = _mm_unpackhi_epi32(T70, T71); __m128i T33 = _mm_unpackhi_epi32(T72, T73); T70 = _mm_unpacklo_epi64(T03, T13); T71 = _mm_unpackhi_epi64(T03, T13); T72 = _mm_unpacklo_epi64(T23, T33); T73 = _mm_unpackhi_epi64(T23, T33); __m128i D0 = _mm_add_epi64(T70, T71); __m128i D1 = _mm_add_epi64(T72, T73); __m128i DD = _mm_add_epi64(D0, D1); _mm_storeu_si128((__m128i *)(dst + ((x + 0) * 4) + y), AA); _mm_storeu_si128((__m128i *)(dst + ((x + 1) * 4) + y), BB); _mm_storeu_si128((__m128i *)(dst + ((x + 2) * 4) + y), CC); _mm_storeu_si128((__m128i *)(dst + ((x + 3) * 4) + y), DD); } } ~~~ * 輸出畫面： ~~~ tina@tina-X550VB:~/Matrix-multiplication$ gcc main.c -o main tina@tina-X550VB:~/Matrix-multiplication$ ./main 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 naive: 1 us 152 158 164 170 504 526 548 570 856 894 932 970 1208 1262 1316 1370 sse: 1 us 152 158 164 170 504 526 548 570 856 894 932 970 1208 1262 1316 1370 ~~~ * 來試試看 8X8 矩陣，因為 rand() 數字太大，我先用 1~10 驗證： ~~~ tina@tina-X550VB:~/Matrix-multiplication$ gcc main.c -o main tina@tina-X550VB:~/Matrix-multiplication$ ./main scr1 2 1 3 4 4 0 5 5 3 2 6 7 2 8 7 8 0 0 8 7 6 4 0 3 3 6 4 5 1 2 0 1 2 8 7 3 5 8 9 4 9 4 6 3 5 1 7 0 5 0 5 9 4 3 9 0 9 7 4 6 7 9 1 6 scr2 6 9 7 2 0 0 0 0 1 0 5 4 5 3 4 1 9 5 6 2 8 3 0 4 7 9 5 7 6 0 0 2 1 0 5 3 9 8 2 2 0 6 9 9 5 7 2 9 8 4 5 3 7 6 7 1 4 7 2 1 9 4 1 8 naive: 12 us 132 124 112 74 169 94 52 74 213 252 235 182 279 170 85 187 139 148 155 122 207 112 23 118 100 111 125 95 120 56 31 60 197 192 261 189 298 211 125 165 194 172 203 112 185 119 77 60 214 205 202 149 208 122 77 82 178 255 285 207 272 182 73 179 sse: 3 us 132 124 112 74 213 252 235 182 139 148 155 122 100 111 125 95 139 148 155 122 301 218 164 148 492 422 320 369 346 260 178 240 220 167 156 155 298 211 125 165 194 172 203 112 185 119 77 60 214 205 202 149 208 122 77 82 178 255 285 207 272 182 73 179 ~~~ * 延伸到 8X8 的矩陣後，發現 SSE 答案有一些部份有問題，雖然度有明顯變快，但仍然還需除錯，可能有些想法有誤。 ## 作業要求3 * ### 在 GitHub 上 fork [matrix_oo](https://github.com/sysprog21/matrix_oo) 專案，指出原本設計的優缺點，並且提出改進方案 (從物件導向封裝、效能，還有資料驗證等觀點) * 提示: 可執行 `make check`，自動編譯測試程式並驗證 * 在原程式碼上可以看見，一些包裝好的結構，讓我們使用更方便： ~~~ clike= //matrix_naive.c MatrixAlgo NaiveMatrixProvider = { .assign = assign, .equal = equal, .sse = sse, .mul =mul, }; //matrix.h typedef struct { void (*assign)(Matrix *thiz, Mat4x4); bool (*equal)(const Matrix *l, const Matrix *r); bool (*mul)(Matrix *dst, const Matrix *l, const Matrix *r); bool (*sse)(Matrix *dst, const Matrix *l, const Matrix *r); } MatrixAlgo; ~~~ ## 作業要求4 * ### 在 [matrix_oo](https://github.com/sysprog21/matrix_oo) 的基礎上，整合 [Matrix Multiplication using SIMD](https://hackmd.io/s/Hk-llHEyx) 的成果，注意，應該建立新的檔案，如 `matrix_sse.c` 和 `matrix_avx.c` 並且設計新的效能分析工具，需要用到程式碼的 `Stopwatch`，時間刻度用 millisecond (ms) * 在這邊我的 matrix_sse.c 也出現了問題正在想辦法解決： ~~~clike= for (int i = 0; i < 4; i++){ for (int j = 0; j < 4; j++) { __m128d I0 = _mm_loadu_pd((__m128d *)&(PRIV(l)->values[0][i])); __m128d I1 = _mm_loadu_pd((__m128d *)&(PRIV(l)->values[1][i])); __m128d I2 = _mm_loadu_pd((__m128d *)&(PRIV(l)->values[2][i])); __m128d I3 = _mm_loadu_pd((__m128d *)&(PRIV(l)->values[3][i])); __m128d I4 = _mm_loadu_pd((__m128d *)&(PRIV(r)->values[j][0])); __m128d I5 = _mm_loadu_pd((__m128d *)&(PRIV(r)->values[j][1])); __m128d I6 = _mm_loadu_pd((__m128d *)&(PRIV(r)->values[j][2])); __m128d I7 = _mm_loadu_pd((__m128d *)&(PRIV(r)->values[j][3])); __m128d T0 = _mm_mul_pd(I0,I4); __m128d T1 = _mm_mul_pd(I1,I5); __m128d T2 = _mm_mul_pd(I2,I6); __m128d T3 = _mm_mul_pd(I3,I7); __m128d C0 = _mm_add_pd(T0, T1); __m128d C1 = _mm_add_pd(T2, T3); __m128d CC = _mm_add_pd(C0, C1); _mm_storeu_pd((__m128i *)&(PRIV(dst)->values[i][j]),CC); } } ~~~ * 使用 float 印出兩個矩陣的直行橫列： ~~~ 1.000000 5.000000 1.000000 5.000000 1.000000 2.000000 3.000000 4.000000 //確認完 ~~~ * 印出乘法的答案 ~~~ 170146294627689351952634138423208181760.000000 -0.000000 170146294627689351952634138423208181760.000000 -0.000000 nan -0.000000 nan -0.000000 nan -0.000000 nan -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 ~~~ * 我在想或許是因為 _mm_mul_pd 的兩個 register 的問題： ~~~ //返回一个__m128d的寄存器，r0=_A0*_B0, r1=_A1*_B1 extern __m128d _mm_mul_pd(__m128d _A, __m128d _B); ~~~ * 但再 SSE2 裡只有這兩種支援 float 乘法運算，都是有兩個 register 的： ~~~ //返回一个__m128d的寄存器，r0=_A0*_B0, r1=_A1 extern __m128d _mm_mul_sd(__m128d _A, __m128d _B); //返回一个__m128d的寄存器，r0=_A0*_B0, r1=_A1*_B1 extern __m128d _mm_mul_pd(__m128d _A, __m128d _B); ~~~ * ### 克服 [matrix_oo](https://github.com/sysprog21/matrix_oo) 裡頭只支援 4x4 矩陣的限制，考慮 sub-matrix 和 column-major 來提升 cache 使用率 * 驗證過程需要考慮不同大小的矩陣 * 需要考慮到矩陣乘法的有效性，不是任意矩陣都可相互執行乘法 ## 作業要求5 * 需要一併考慮 data alignment 議題 * 學習 crypto API 並設計介面使演算法容易擴充和比較 * 學習的案例: [Linaro Cortex String benchmark for ARMv8](https://wiki.linaro.org/WorkingGroups/Kernel/ARMv8/cortex-strings) * 關於 benchmark 設計可參考: [ssvb/tinymembench](https://github.com/ssvb/tinymembench) 和 [hglm/test-arm-kernel-memcpy](https://github.com/hglm/test-arm-kernel-memcpy) * 允許透過亂數作為資料輸入，驗證矩陣乘法的效能，需要視覺化 (透過 gnuplot 或 R) * 留意記憶體存取，避免 leak，參考 [你所不知道的C語言：技巧篇](https://hackmd.io/s/HyIdoLnjl) 提到的 smart pointer * 指出原有的程式碼是否 [thread-safe](https://en.wikipedia.org/wiki/Thread_safety)，並提出改進方案