2017q1 team6 matrix

constributed by <henry0929016816 vtim9907>

Column order problem

這是一個簡單的 sub-matrix 版本，把巨大的 array 切成小的 4 * 4 矩陣後再做運算。














    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
  Image Not Showing Possible Reasons
  - The image file may be corrupted
  - The server hosting the image is unavailable
  - The image path is incorrect
  - The image format is not supported
  Learn More →

所以如果你的資料是分布在 1~4 那還是會

第一次取:0~3 將，0 的資料去掉，留下 1~3
第二次取:4~7 將，5~7 的資料去掉，留下 4
再將 1~3 4 合起來
Image Not Showing Possible Reasons
- The image file may be corrupted
- The server hosting the image is unavailable
- The image path is incorrect
- The image format is not supported
Learn More →
由於資料分布不在 4 的倍數，導致了存取速度降低，~~為了避免這個情況，我們一開始就要先做好資料的處理，連續資料的 address 都能是 4 的倍數~~而 compiler 一開始再分配記憶體時，就會按照宣告的型態去做 alignment ，例如 int 就是 4 byte alignment。
struct 會自動做 alignment，假設創了一個 struct，如下面 code 所示




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


















#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 抓取的的方式。~~

注意

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

以上是我搞錯了，我以為只要不是對齊 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。

注意

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

以下實驗為不好的實驗，struct 裡的資料量太小，可能看不出兩個 struct 執行時間的差異，另外沒有去清除 cache 導致讀取的資料可能會再 cache 裡而使得處理速度一樣快的錯覺


















...
...
#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 議題，然而老舊的處理器可能就不是這樣了（手邊沒有可以實驗的）

重新設計實驗

方式一：將 struct 資料結構改大













#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 次看是否平均起來的執行速度差不多





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 的差異。

今天剛好有看到 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 有看到 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，裡面有特別提到 :

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 對齊：




    for(int i=0;i<10000;++i){
        char *z;
        z=malloc(sizeof(char));
    }

結果用 gdb 測過之後，發現位址的結尾的確都是 0 結尾，表示真的是以 16-byte 做對齊。

其實這邊最大的疑惑是在 alignment 系列 function 的 man page 裡，最下面的 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 :





































































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{matrix} 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{matrix}]

[\begin{matrix} 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{matrix}]

[\begin{matrix} 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{matrix}]

分割步驟來看，先算

C_{1}

C_{4}

：

使用 _mm_set1_epi32 可以塞四個相同的 int 到一個 __m128i 裡，如下

$D_{1}$ =
$[\begin{matrix} A_{1} & A_{1} & A_{1} & A_{1} \end{matrix}]$

$D_{2}$ =
$[\begin{matrix} A_{2} & A_{2} & A_{2} & A_{2} \end{matrix}]$

$D_{3}$ =
$[\begin{matrix} A_{3} & A_{3} & A_{3} & A_{3} \end{matrix}]$

$D_{4}$ =
$[\begin{matrix} A_{4} & A_{4} & A_{4} & A_{4} \end{matrix}]$
而
$R$ 就先存 B 的列

$R_{1}$ =
$[\begin{matrix} B_{1} & B_{2} & B_{3} & B_{4} \end{matrix}]$

$R_{2}$ =
$[\begin{matrix} B_{5} & B_{6} & B_{7} & B_{8} \end{matrix}]$

$R_{3}$ =
$[\begin{matrix} B_{9} & B_{10} & B_{11} & B_{12} \end{matrix}]$

$R_{4}$ =
$[\begin{matrix} B_{13} & B_{14} & B_{15} & B_{16} \end{matrix}]$
然後就可以用 _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 版乘法

















































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{matrix} A_{1} & A_{2} \\ A_{3} & A_{4} \end{matrix}]$ X
$[\begin{matrix} B_{1} & B_{2} \\ B_{3} & B_{4} \end{matrix}]$ =
$[\begin{matrix} 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{matrix}]$

而我的寫法就是分別拆開相乘得到的結果，再相加得到答案

$[\begin{matrix} A_{1} B_{1} & A_{1} B_{2} \\ A_{3} B_{1} & A_{3} B_{2} \end{matrix}] + [\begin{matrix} A_{2} B_{3} & A_{2} B_{4} \\ A_{4} B_{3} & A_{4} B_{4} \end{matrix}] = [\begin{matrix} 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{matrix}]$
_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 的議題了（終於啊！）

henry0929016816關於這裡我還是有點不懂為何會讀取到錯誤的記憶體

memory leak

參考資料

Data Alignment
Alignment in C
pragma