2017q1 team6 matrix

# 2017q1 team6 matrix constributed by <`henry0929016816` `vtim9907`> ## Column order problem 這是一個簡單的 sub-matrix 版本，把巨大的 array 切成小的 4 * 4 矩陣後再做運算。 ```clike= for (int i = 0; i < dst->row; i += 4) { for (int j = 0; j < dst->col; j += 4) { for (int k = 0; k < N; k += 4) { for (int i2 = 0; i2 < 4; ++i2) { for (int j2 = 0; j2 < 4; ++j2) { for (int k2 = 0; k2 < 4; ++k2) { PRIV(dst)->values[i + i2][j + j2] += PRIV(l)->values[i + i2][k + k2] * PRIV(r)->values[k + k2][j + j2]; } } } } } } ``` 在每次跑最內圈的 for 時，取左陣列和右陣列的元素出來看 memory address 到底長怎樣，並把左右陣列 print 出的結果分開來看： ``` address PRIV(l)->values[i + i2][k + k2] = 0x7e2470 address PRIV(l)->values[i + i2][k + k2] = 0x7e2474 address PRIV(l)->values[i + i2][k + k2] = 0x7e2478 address PRIV(l)->values[i + i2][k + k2] = 0x7e247c address PRIV(r)->values[k + k2][j + j2] = 0x803ffc address PRIV(r)->values[k + k2][j + j2] = 0x80443c address PRIV(r)->values[k + k2][j + j2] = 0x80487c address PRIV(r)->values[k + k2][j + j2] = 0x804cbc ``` 很明顯，右陣列因為是以 column order 的方式讀取，所以 memory address 不連續，沒辦法有效利用 cache。不過因為這個版本並沒有特別用到硬體支援的技術，譬如 SEE、 AVX，所以 Data alignment 與這個版本比較沒有關係。 > 我再確定一下 @@ ## data alignment * 在講到 data alignment 之前，我們要先提到在程式語言裡，一個 data object 所擁有的特性是甚麼。一個 data object 具有兩個特性 * value * storage location (address) * 而 data alignment 的意思就是， data 的 address 可以公平的被 1, 2, 4, 8,這些數字整除，從這些數字可以發現他們都是 2 的 N 次方 (N為從 0 開始的整數)。換句話說，這些 data object 可以用 1 ,2 ,4 ,8 byte 去 alignment。 * 懂得 data alignment 是甚麼之後，就要開始講電腦的 cpu 如何抓取資料，cpu 不會一次只抓取 1 byte 的資料，因為這樣太慢了，如果有個資料型態是 int 的資料，如果你只抓取 1 byte ,就必須要抓 4 次(int 為 4 byte)，有夠慢。所以 cpu 通常一次會取 4 byte(要看電腦的規格 32 位元的 cpu 一次可以讀取 32 bit 的資料，64 位元一次可以讀取 64 bit)，並且是按照順序取的，例如: * 第一次取:0~3 * 第二次取:4~7 ![](https://i.imgur.com/aDCYyWc.png) 所以如果你的資料是分布在 1~4 那還是會 * 第一次取:0~3 將，0 的資料去掉，留下 1~3 * 第二次取:4~7 將，5~7 的資料去掉，留下 4 * 再將 1~3 4 合起來 ![](https://i.imgur.com/wIfEVy9.png) 由於資料分布不在 4 的倍數，導致了存取速度降低，~~為了避免這個情況，我們一開始就要先做好資料的處理，連續資料的 address 都能是 4 的倍數~~而 compiler 一開始再分配記憶體時，就會按照宣告的型態去做 alignment ，例如 int 就是 4 byte alignment。 * struct 會自動做 alignment，假設創了一個 struct，如下面 code 所示 ```c= struct s1 { char c; int a; } ``` 照我們學過的 char 為 1 byte，而 int 為 4 byte ，兩個加起來應該為 5 byte，然而實際上為 8 byte，由於 int 為 4 byte ，所以 char 為了要能 alignment 4 byte 就多加了 3 byte 進去，使得 cpu 存取速度不會因 address 不是在 4 的倍數上，存取變慢。 * 原本想說 struct 會自動 alignment 那我們幹嘛討論呢? 後來實驗了一下，發現在 struct 裡放陣列他並不會自動 alignment ```c= #include<stdio.h> #include<stdlib.h> #include<string.h> typedef struct _s1 { char a[5]; }s1; int main() { int i; s1 p[10]; printf("struct s1 size: %ld byte\n",sizeof(s1)); for(i=0;i<10;i++) { printf("the struct p[%d] address =%p\n",i,p+i); } } ``` 得到執行結果為 struct s1 size: 5 byte the struct p[0] address =0x7ffe72bfd270 the struct p[1] address =0x7ffe72bfd275 the struct p[2] address =0x7ffe72bfd27a the struct p[3] address =0x7ffe72bfd27f the struct p[4] address =0x7ffe72bfd284 the struct p[5] address =0x7ffe72bfd289 the struct p[6] address =0x7ffe72bfd28e the struct p[7] address =0x7ffe72bfd293 the struct p[8] address =0x7ffe72bfd298 the struct p[9] address =0x7ffe72bfd29d ~~所以如果我們現在探討的議題，matrix 不是 4x4 的大小，可能是 3x3，那我們就要考慮如何 aligment 才能使資料的位址符合 cpu 抓取的的方式。~~ :::warning 注意 :zap: 以上是我搞錯了，我以為只要不是對齊 4 byte 就沒有自動 alignment ，然而其實還是有 alignment 的，由於 char type 的 data 大小只佔 1 byte 所以只要 1 byte alignment ，也就是不用使用 padding 讓其變成 4 的倍數。而我們原本的 matrix 的實作也可以不用考慮到 data alignment 由於 compiler 會自動幫我們以 data 的大小做 alignment ，假設有 int type（4 byte ) 他再創建的時候就會是已 4 byte alignment。 ::: * 現代的處理器再處理一般的讀取跟寫入資料，不用考慮 data alignment 的議題，例如 intel x86 系統是允許 data unalignment 的以下用我的電腦實驗 data alignment 跟 data unalignment 的差異。 * 首先我創立兩個 struct 兩個放的東西是相同的，唯一不同的是 t1 有加 pack 這條指令告訴 compiler 說 test1 裡的 data 只要 1 byte alignment 就好，t2 則是會按照宣告的 type 作 alignment 所以 t2 裡會有 padding。 :::warning 注意 :zap: 以下實驗為不好的實驗，struct 裡的資料量太小，可能看不出兩個 struct 執行時間的差異，另外沒有去清除 cache 導致讀取的資料可能會再 cache 裡而使得處理速度一樣快的錯覺 ```c= ... ... #pragma pack(push) #pragma pack(1) typedef struct _test1 { int t; char a; int b; }test1; #pragma pack(pop) typedef struct _test2 { int t; char a; int b; }test2; ``` 分別讓兩個 struct 執行對 b 的資料操作 ``` ... for(int i=0;i<4000;i++) { t1.b=0; t1.b+=1; } ... ... for(int i=0;i<4000;i++) { t2.b=0; t2.b+=1; } ... ``` 發現執行的時間沒有差異 * 執行結果 struct test1 size:9 struct test2 size:12 alignment excution time : 0.000011(s) unalignment excution time : 0.000012(s) 證明了現代處理器其實不用考慮 data alignment 議題，然而老舊的處理器可能就不是這樣了（手邊沒有可以實驗的） ::: :::info 重新設計實驗 * 方式一：將 struct 資料結構改大 ```c= #pragma pack(1) typedef struct _test1 { char c[3]; int num[256]; }test1; #pragma pack(pop) typedef struct _test2 { char c[3]; int num[256]; }test2; ``` * 方式二：去掉 cache 的影響在 Makefile 裡加入這一條指令 `echo 3 > /proc/sys/vm/drop_caches` 並且讓這段 code 跑個 200 次看是否平均起來的執行速度差不多 ```c= for i in `seq 0 1 200`; do \ echo 3 > /proc/sys/vm/drop_caches ;\ printf "%d," $$i;\ ./pointer; \ done > clock_gettime.txt ``` * 最後畫出圖形 ![](https://i.imgur.com/yUS7zcw.png) 有幾比測資誤差滿大的，但大部份都相同，關於誤差很大的問題，到現在還是不懂，只知道這個問題從第1週就一直困惑我滿久的，知道跟 cache 有關但不知道為什麼會這個樣子。 ::: * 要考慮 data alignment 議題的只有當我們要實作 avx 板的矩陣乘法，由於 avx 操作的資料量比較大（32 byte)，所以我們要幫資料對齊 32 byte，當然不對齊也是可以操作，只是會比較慢，以下我會實驗看看 32 byte alignment 跟 unalignment 的差異。 >[color=blue]今天剛好有看到 c11 有提供一個新的函式叫 void *aligned_alloc(size_t alignment, size_t size); 可以自訂要 align 的 size，等晚點再來看的詳細說明） ## 矩陣乘法 SSE 版因為 SSE 的暫存器 (XMM0 ~ XMM7) 有 128 bits 的長度，也就是 4 個 integer 的大小，根據這個特性，在做矩陣乘法時，選擇將大矩陣切割成多個 4 * 4 的小矩陣來做會比較有效率。 ### data alignment 因為原本要做乘法的兩個矩陣都是用 malloc 配置出來的，在 [Matrix Multiplication using SIMD](https://hackmd.io/s/Hk-llHEyx#對齊下的-sse) 有看到 malloc 本來配置出來的記憶體位置就有做 alignment，根據 malloc 的 man page 裡提到 : ``` The malloc() and calloc() functions return a pointer to the allocated memory, which is suitably aligned for any built-in type. ``` 但單看這句英文，其實我意會不太過來實際上到底 malloc 做了怎樣的 data alignment，所以又查了 [The GNU C Library - Malloc Example](ftp://ftp.gnu.org/old-gnu/Manuals/glibc-2.2.3/html_node/libc_25.html)，裡面有特別提到 : ``` In the GNU system, the address is always a multiple of eight on most systems, and a multiple of 16 on 64-bit systems. ``` 這句我就有看懂，所以說在 "大多數" 系統下，malloc 會以 8 bytes 進行對齊；而在 64-bit 的系統下，malloc 則會以 16 bytes 進行對齊。做了小實驗，malloc 在我的系統(64-bit)下是以 16 bytes 對齊： ```clike= for(int i=0;i<10000;++i){ char *z; z=malloc(sizeof(char)); } ``` 結果用 gdb 測過之後，發現位址的結尾的確都是 0 結尾，表示真的是以 16-byte 做對齊。 >[color=red]其實這邊最大的疑惑是在 [alignment 系列 function 的 man page](http://man7.org/linux/man-pages/man3/posix_memalign.3.html) 裡，最下面的 NOTES 提到 " The glibc malloc(3) always returns 8-byte aligned memory addresses "，所以才不太確定到底是只會對齊 8-byte，才去查其他文獻和做實驗，證明 malloc 在 64-bit 系統下會以 16-byte 對齊，所以也許上面貼的 man page 是單指 32-bit 的系統吧!? 確認了 malloc 的機制後，可以確定矩陣不需要額外用其他 function 做 data alignment，用原本 malloc 的配置方式，就可以直接按原計畫做 SSE 版本的矩陣乘法，就直接實作。 ### 先做矩陣轉置 ? 原本有在思考要不要先把第二個矩陣先做轉置，這樣就不會有 column major 的問題，不過想想，這樣有點太大費周章...，把轉置的時間也算進去的話，效能大概也會變慢，還是算了@@ ### sse_matrix 實做的部份我把 4 * 4 的小矩陣乘法做了 loop unrolling，所以 code 有一點點長，不過因為想到一種比較特別的矩陣乘法實做方式，滿精簡的，所以實際用到的指令其實比我預想的少。 code : ```clike= for (int i = 0; i < dst->row; i += 4) { for (int j = 0; j < dst->col; j += 4) { __m128i dst0 = _mm_setzero_si128(); __m128i dst1 = _mm_setzero_si128(); __m128i dst2 = _mm_setzero_si128(); __m128i dst3 = _mm_setzero_si128(); for (int k = 0; k < l->col; k += 4) { __m128i r1 = _mm_load_si128((__m128i *)(&PRIV(r)->values[k][j])); __m128i r2 = _mm_load_si128((__m128i *)(&PRIV(r)->values[k+1][j])); __m128i r3 = _mm_load_si128((__m128i *)(&PRIV(r)->values[k+2][j])); __m128i r4 = _mm_load_si128((__m128i *)(&PRIV(r)->values[k+3][j])); __m128i b1 = _mm_set1_epi32(PRIV(l)->values[i][k]); __m128i b2 = _mm_set1_epi32(PRIV(l)->values[i][k+1]); __m128i b3 = _mm_set1_epi32(PRIV(l)->values[i][k+2]); __m128i b4 = _mm_set1_epi32(PRIV(l)->values[i][k+3]); __m128i b5 = _mm_set1_epi32(PRIV(l)->values[i+1][k]); __m128i b6 = _mm_set1_epi32(PRIV(l)->values[i+1][k+1]); __m128i b7 = _mm_set1_epi32(PRIV(l)->values[i+1][k+2]); __m128i b8 = _mm_set1_epi32(PRIV(l)->values[i+1][k+3]); __m128i b9 = _mm_set1_epi32(PRIV(l)->values[i+2][k]); __m128i b10 = _mm_set1_epi32(PRIV(l)->values[i+2][k+1]); __m128i b11 = _mm_set1_epi32(PRIV(l)->values[i+2][k+2]); __m128i b12 = _mm_set1_epi32(PRIV(l)->values[i+2][k+3]); __m128i b13 = _mm_set1_epi32(PRIV(l)->values[i+3][k]); __m128i b14 = _mm_set1_epi32(PRIV(l)->values[i+3][k+1]); __m128i b15 = _mm_set1_epi32(PRIV(l)->values[i+3][k+2]); __m128i b16 = _mm_set1_epi32(PRIV(l)->values[i+3][k+3]); __m128i t0 = _mm_add_epi32( _mm_add_epi32( _mm_mullo_epi32(b1, r1), _mm_mullo_epi32(b2, r2)), _mm_add_epi32( _mm_mullo_epi32(b3, r3), _mm_mullo_epi32(b4, r4))); __m128i t1 = _mm_add_epi32( _mm_add_epi32( _mm_mullo_epi32(b5, r1), _mm_mullo_epi32(b6, r2)), _mm_add_epi32( _mm_mullo_epi32(b7, r3), _mm_mullo_epi32(b8, r4))); __m128i t2 = _mm_add_epi32( _mm_add_epi32( _mm_mullo_epi32(b9, r1), _mm_mullo_epi32(b10, r2)), _mm_add_epi32( _mm_mullo_epi32(b11, r3), _mm_mullo_epi32(b12, r4))); __m128i t3 = _mm_add_epi32( _mm_add_epi32( _mm_mullo_epi32(b13, r1), _mm_mullo_epi32(b14, r2)), _mm_add_epi32( _mm_mullo_epi32(b15, r3), _mm_mullo_epi32(b16, r4))); dst0 = _mm_add_epi32(t0, dst0); dst1 = _mm_add_epi32(t1, dst1); dst2 = _mm_add_epi32(t2, dst2); dst3 = _mm_add_epi32(t3, dst3); } _mm_store_si128((__m128i *)(&PRIV(dst)->values[i][j]), dst0); _mm_store_si128((__m128i *)(&PRIV(dst)->values[i+1][j]), dst1); _mm_store_si128((__m128i *)(&PRIV(dst)->values[i+2][j]), dst2); _mm_store_si128((__m128i *)(&PRIV(dst)->values[i+3][j]), dst3); } } ``` 過程簡單講大概如下：假設矩陣 A x B = C $\begin{bmatrix} A_{1} & A_{2} & A_{3} & A_{4} \\ A_{5} & A_{6} & A_{7} & A_{8} \\ A_{9} & A_{10} & A_{11} & A_{12} \\ A_{13} & A_{14} & A_{15} & A_{16} \\ \end{bmatrix}$ X $\begin{bmatrix} B_{1} & B_{2} & B_{3} & B_{4} \\ B_{5} & B_{6} & B_{7} & B_{8} \\ B_{9} & B_{10} & B_{11} & B_{12} \\ B_{13} & B_{14} & B_{15} & B_{16} \\ \end{bmatrix}$ = $\begin{bmatrix} C_{1} & C_{2} & C_{3} & C_{4} \\ C_{5} & C_{6} & C_{7} & C_{8} \\ C_{9} & C_{10} & C_{11} & C_{12} \\ C_{13} & C_{14} & C_{15} & C_{16} \\ \end{bmatrix}$ 分割步驟來看，先算 $C_1$ ~ $C_4$ ： * 使用 _mm_set1_epi32 可以塞四個相同的 int 到一個 __m128i 裡，如下 $D_1$ = $\begin{bmatrix} A_{1} & A_{1} & A_{1} & A_{1} \end{bmatrix}$ $D_2$ = $\begin{bmatrix} A_{2} & A_{2} & A_{2} & A_{2} \end{bmatrix}$ $D_3$ = $\begin{bmatrix} A_{3} & A_{3} & A_{3} & A_{3} \end{bmatrix}$ $D_4$ = $\begin{bmatrix} A_{4} & A_{4} & A_{4} & A_{4} \end{bmatrix}$ * 而 $R$ 就先存 B 的列 $R_1$ = $\begin{bmatrix} B_{1} & B_{2} & B_{3} & B_{4} \end{bmatrix}$ $R_2$ = $\begin{bmatrix} B_{5} & B_{6} & B_{7} & B_{8} \end{bmatrix}$ $R_3$ = $\begin{bmatrix} B_{9} & B_{10} & B_{11} & B_{12} \end{bmatrix}$ $R_4$ = $\begin{bmatrix} B_{13} & B_{14} & B_{15} & B_{16} \end{bmatrix}$ * 然後就可以用 _mm_mullo_epi32 做相乘，然後全部加起來，就會得到 $C_1$ ~ $C_4$ _mm_mullo_epi32（$R_1$, $D_1$）+ _mm_mullo_epi32（$R_2$, $D_2$）+ _mm_mullo_epi32（$R_3$, $D_3$）+ _mm_mullo_epi32（$R_4$, $D_4$） = $C_1$ ~ $C_4$ 用同樣概念做四次，就完成了一次 4 * 4 的小矩陣相乘！ * 效能 * submatrix 版本 ![](https://i.imgur.com/CweT6gq.png) * sse 版本 ![](https://i.imgur.com/9BnvZMK.png) ## 矩陣乘法 AVX 版 ### avx set 版乘法 ```c= matrix avx_multiple(matrix a,matrix b) { matrix result; __m256 r[2][4]; for(int i=0;i<4;i++) { __m256 a_row_01=_mm256_set_ps(a.value[0][i],a.value[0][i],a.value[0][i],a.value[0][i], a.value[1][i],a.value[1][i],a.value[1][i],a.value[1][i]); __m256 a_row_23=_mm256_set_ps(a.value[2][i],a.value[2][i],a.value[2][i],a.value[2][i], a.value[3][i],a.value[3][i],a.value[3][i],a.value[3][i]); __m256 b_col=_mm256_set_ps(b.value[i][0],b.value[i][1],b.value[i][2],b.value[i][3], b.value[i][0],b.value[i][1],b.value[i][2],b.value[i][3]); r[0][i]=_mm256_mul_ps (a_row_01,b_col); r[1][i]=_mm256_mul_ps (a_row_23,b_col); } __m256 sum0= _mm256_add_ps (r[0][0],r[0][1]); sum0 = _mm256_add_ps (sum0,r[0][2]); sum0 = _mm256_add_ps (sum0,r[0][3]); __m256 sum1= _mm256_add_ps (r[1][0],r[1][1]); sum1 = _mm256_add_ps (sum1,r[1][2]); sum1 = _mm256_add_ps (sum1,r[1][3]); result.value[0][0]=sum0[7]; result.value[0][1]=sum0[6]; result.value[0][2]=sum0[5]; result.value[0][3]=sum0[4]; result.value[1][0]=sum0[3]; result.value[1][1]=sum0[2]; result.value[1][2]=sum0[1]; result.value[1][3]=sum0[0]; result.value[2][0]=sum1[7]; result.value[2][1]=sum1[6]; result.value[2][2]=sum1[5]; result.value[2][3]=sum1[4]; result.value[3][0]=sum1[3]; result.value[3][1]=sum1[2]; result.value[3][2]=sum1[1]; result.value[3][3]=sum1[0]; return result; } ``` * 執行時間跟 naive 板相比 avx execution time 0.000579(s) naive execution time 0.001123(s) 可以發現有比 naive 版本的乘法還要快 * 首先，先來介紹 code 的概念，我們都知道 matrix 乘法是行列乘法 $\begin{bmatrix} A_{1} & A_{2} \\ A_{3} & A_{4} \end{bmatrix}$X$\begin{bmatrix} B_{1} & B_{2} \\ B_{3} & B_{4} \end{bmatrix}$=$\begin{bmatrix} A_{1}B_{1}+A_{2}B_{3} & A_{1}B_{2}+A_{2}B_{4} \\ A_{3}B_{1}+A_{4}B_{3} & A_{3}B_{2}+A_{4}B_{4} \end{bmatrix}$ 而我的寫法就是分別拆開相乘得到的結果，再相加得到答案 $\begin{bmatrix} A_{1}B_{1} & A_{1}B_{2} \\ A_{3}B_{1} & A_{3}B_{2} \end{bmatrix}+\begin{bmatrix} A_{2}B_{3} & A_{2}B_{4} \\ A_{4}B_{3} & A_{4}B_{4} \end{bmatrix}=\begin{bmatrix} A_{1}B_{1}+A_{2}B_{3} & A_{1}B_{2}+A_{2}B_{4} \\ A_{3}B_{1}+A_{4}B_{3} & A_{3}B_{2}+A_{4}B_{4} \end{bmatrix}$ * `_mm256_set_ps (float e7, float e6, float e5, float e4, float e3, float e2, float e1, float e0)` 是設定 avx 的 register 裡的值要為多少的函式，由於我們是 float type 的 matrix 所以要將 8 個值傳進去，傳值的方式還是一個一個傳，而這種作法很慢，如果能一次就讀取 8 個值想必會快的多，而這樣就要使用 _mm256_load_ps 這個函式。 * `_mm256_load_ps (float const * mem_addr)` 由參數的 address 開始往後讀取 31 byte 的資料，如果傳入的 address 不是 32 byte alignment 會出錯 (core dumped)，由於讀取到錯誤的記憶體，所以這裡就要考慮到 data alignment 的議題了（終於啊！） >[name=henry0929016816]關於這裡我還是有點不懂為何會讀取到錯誤的記憶體 ## memory leak ## 參考資料 [Data Alignment](http://www.songho.ca/misc/alignment/dataalign.html) [Alignment in C](https://wr.informatik.uni-hamburg.de/_media/teaching/wintersemester_2013_2014/epc-14-haase-svenhendrik-alignmentinc-paper.pdf) [pragma](http://topalan.pixnet.net/blog/post/22334686)