# 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
所以如果你的資料是分布在 1~4 那還是會
* 第一次取:0~3 將，0 的資料去掉，留下 1~3
* 第二次取:4~7 將，5~7 的資料去掉，留下 4
* 再將 1~3 4 合起來
由於資料分布不在 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
```
* 最後畫出圖形
有幾比測資誤差滿大的，但大部份都相同，關於誤差很大的問題，到現在還是不懂，只知道這個問題從第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 版本
* sse 版本
## 矩陣乘法 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)
Sign in with Wallet
Connect another wallet