# 2018q3 矩陣乘法 Matrix multiplication
Contributed by < [DyslexiaS](https://github.com/DyslexiaS), [siahuat0727](https://github.com/siahuat0727), [pjchiou](https://github.com/pjchiou) >
[GitHub](https://github.com/siahuat0727/matrix-multiplication)
Task : 重現 [Matrix Multiplication using SIMD](https://hackmd.io/s/Hk-llHEyx) 實驗，並依循 CS:APP 第 6 章指引，分析個別實作的效能

---

## outline
- 演算法
- SSE 指令集與 AVX 指令集的對齊問題。

## 演算法
- [x] naive
- [x] cache-friendly
- [x] sub-matrix
- [x] strassen
- [x] SIMD-SSE
- [x] SIMD-AVX

---

### Naive algorithm
最簡單的矩陣相乘演算法
- time complexity : $O(n^3)$
- 演算法

 for (int i = 0; i < a.row; ++i)
        for (int j = 0; j < b.col; ++j) {
            int sum = 0;
            for (int k = 0; k < a.col; ++k)
                sum += a.values[i][k] * b.values[k][j];
            dst->values[i][j] += sum;
        }

---

### cache-friendly

在 naive algorithm 中， `b.values[k][j]` 一直在取用距離為 `b.col` 的變數，這樣對 cache 的使用不利。要利用 cache 的特性來加速運算，把握兩個重點
1. Temporal locality : 一旦取用過的變數，重複取用。
2. Spatial locality : 取用記憶體中相鄰的變數。
- time complexity : $O(n^3)$
- 演算法

for (int i = 0; i < a.row; ++i)
        for (int k = 0; k < a.col; ++k) {
            int r = a.values[i][k];
            for (int j = 0; j < b.col; ++j)
                dst->values[i][j] += r * b.values[k][j];
        }

如此一來，重複取用 `a[i][k] `且每個矩陣都是以步長為 1 取用變數。

---

### Sub-matrix
Divied and conquer 的概念，先把矩陣拆成較小的矩陣，各自運算後再加起來。
在這個示意圖中，`A`、`B`兩個矩陣都是正方形的，其 sub-matrix 行列皆為原本的一半。這與 naive 版本的功能不符，所以我們將其稍做調整
1. `A`、`B` 矩陣不必為正方形，符合矩陣相乘的限制 `A.col==B.row` 即可。
2. sub-matrix 行列大小不必為 N/2 ，但必須為正方形。
3. sub-matrix 大小亦不必可被 C.row 或 C.col 整除。

- time complexity : $O(n^3)$
- 演算法
- reference : [divide and conquer](https://www.geeksforgeeks.org/strassens-matrix-multiplication/)

//Calculate sub-matrix at (c_row,c_col) and the size is stride.
void matmul_stride(const Matrix l,
                   const Matrix r,
                   const Matrix dst,
                   int c_row,
                   int c_col,
                   int stride)
{
    for (int k = 0; k < l.col; k += stride)
        for (int i = 0; i < stride && i + c_row < l.row; i++)
            for (int j = 0; j < stride && j + c_col < r.col; j++)
                for (int m = k; m < k + stride && m < l.col; m++)
                    dst.values[i + c_row][j + c_col] +=
                        l.values[i + c_row][m] * r.values[m][j + c_col];
}

void sub_matmul(const Matrix l, const Matrix r, const Matrix dst, void *ctx)
{
    CHECK(l, r, dst);

    assert(ctx != NULL);
    int stride = ((SubMatrixInfo*)ctx)->stride;
    for (int i = 0; i < l.row; i += stride)
        for (int j = 0; j < r.col; j += stride)
            matmul_stride(l, r, dst, i, j, stride);
}

---

### Strassen
兩個矩陣相乘所要用到的乘法及加法次數：
因為 $n \times n$ 矩陣 $C$ 有 $n^2$ 個元素，計算每一元素需要 n 個乘法和 (n-1) 個加法，故知 $n\times n$ 階矩陣乘積共使用了 $n^3$ 個乘法和 $n^3-n^2$ 個加法。
這個迷思直到更快速的 divide-and-conquer 矩陣乘法 strassen 被提出才打破。

- 考慮 $2\times 2$ 階分塊矩陣乘法

其中每一分塊皆為方陣。傳統的矩陣乘法運算方式，$C_{ij}=A_{i1}B_{1j}+A_{i2}B_{2j}$，總共使用8個分塊乘法和4個分塊加法。
Strassen 演算法使用7個分塊乘法和18個分塊加法，如下：

\begin{aligned} P_1&=(A_{11}+A_{22})(B_{11}+B_{22})\\ P_2&=(A_{21}+A_{22})B_{11}\\ P_3&=A_{11}(B_{12}-B_{22})\\ P_4&=A_{22}(B_{21}-B_{11})\\ P_5&=(A_{11}+A_{12})B_{22}\\ P_6&=(A_{21}-A_{11})(B_{11}+B_{12})\\ P_7&=(A_{12}-A_{22})(B_{21}+B_{22})\\ C_{11}&=P_1+P_4-P_5+P_7\\ C_{12}&=P_3+P_5\\ C_{21}&=P_2+P_4\\ C_{22}&=P_1+P_3-P_2+P_6. \end{aligned}

- 直接代入化簡即可驗證上述公式的正確性
假設 $n=2m$，上面所有分塊皆為 $m\times m$ 。
- 傳統的矩陣乘法 :
- 使用了 $8m^3$ 個乘法和 $8m^3-4m^2$ 個加法。
- Strassen :
- Strassen 演算法僅執行至分塊等級，也就是說，$m\times m$ 分塊乘法仍採用傳統方法計算
- 則7次分塊乘法共使用 $7m^3$ 個乘法和 $7(m^3-m^2)$ 個加法，再加上 18 次分塊加法使用的 $18m^2$ 個加法，總計是 $7m^3$ 個乘法和 $7m^3+11m^2$ 個加法
- 若 $m\gg 1$，不論乘法或加法，Strassen 演算法僅耗費傳統方法 $7/8$ 的計算量。

Reference : [strassen](https://ccjou.wordpress.com/2013/06/04/分治矩陣乘法──strassen-演算法/)

---

### SIMD-SSE

- SIMD 為 ==Single Instruction Multiple Data== 的縮寫，以 Intel 的處理器為例，其發展的 MMX SIMD 指令集，能以一道指令處理多個資料，其目的就是為了提高多媒體資料的處理能力。其後，Intel 又發展出 SSE SIMD 指令集，為 MMX 的擴充指令集。
- AMD 的處理器亦有其 SIMD 指令集，但我用的是 Intel 的處理器，所以我以 Intel 的 `SSE SIMD` 指令集做為例子。
- SSE 有 8 個 128 bit 的暫存器，每個暫存器可以用來存
- 4 個 32 bit 的 float
- 2 個 64 bit 的 double
- 4 個 32 bit 的 integer
- 8 個 16 bit 的 short
- 16 個 8 bit 的 char

利用 SIMD 指令集做矩陣乘法時，其流程與 sub-matrix 無異，但有一個地方要注意的是，矩陣乘法的本質就是向量的內積，因為我們的暫存器大小是 8 個 128 bit ，剛好可以用來==存 8 個有 4 個 integer 的向量==，所以 sub-matrix 的大小就是 $4*4$ ，這樣可以充分利用暫存器。

我們分別用以下 4 個 vector 來存 A、B 兩個矩陣中的 sub-matrix 。

我們要做的事情就是從 [SIMD 指令集](https://software.intel.com/sites/landingpage/IntrinsicsGuide/)中組合出可以達到需求的指令。以下說明如何得到最上面的 row ，其他以此類推。

==使用前注意每個 function 必須 include 哪個 header file==

1. 讀取矩陣值，存入暫存器。宣告型別為 `__128i` 的變數來代表這個 128 bit 長可存四個 32 bit integer 的變數。用以下兩個 function 來讀取矩陣值。
- __m128i _mm_load_si128 (__m128i const* mem_addr)
- __m128i _mm_set_epi32 (int e3, int e2, int e1, int e0)

因為 A 矩陣中的向量值是存在連續的記憶體內，用前者； B 矩陣用後者。

2. 將 V1 分別對 V5~V8 做內積，可得四個內積後的向量 `V15`, `V16`, `V17`, `V18`
- __m128i _mm_mullo_epi32 (__m128i a, __m128i b)。
3. 我們要將這四個向量 V15, V16, V17, V18 組成的矩陣做轉置，但這四個向量是分別獨立的四個 __m128i 變數，只是在我們腦海裡是一個矩陣，我們要利用 SIMD 的四個指令達到等價效果。
- __m128i _mm_unpackhi_epi32 (__m128i a, __m128i b)

將兩個 m128i 變數的前兩筆 32 bits 資料交錯擺放。
r[0] = a[0] ; r[1] = b[0];
r[2] = a[1] ; r[3] = b[1];

- __m128i _mm_unpackhi_epi64 (__m128i a, __m128i b)

將兩個 m128i 變數的前兩筆 64 bits 資料交錯擺放。
r[0] = a[0]; r[1] = a[1]; 
r[2] = b[0]; r[3] = b[1];

- __m128i _mm_unpacklo_epi32 (__m128i a, __m128i b)

將兩個 m128i 變數的後兩筆 32 bits 資料交錯擺放。
r[0] = a[2] ; r[1] = b[2];
r[2] = a[3] ; r[3] = b[3];

- __m128i _mm_unpacklo_epi64 (__m128i a, __m128i b)

將兩個 m128i 變數的後兩筆 64 bits 資料交錯擺放。
r[0] = a[2]; r[1] = a[3]; 
r[2] = b[2]; r[3] = b[3];

先將這四個向量分為兩個一組後，依序做以下四個函式，可以得到另外四個向量。
- S1 = _mm_unpacklo_epi32(V15,V16)
- S2 = _mm_unpackhi_epi32(V15,V16)
- S3 = _mm_unpacklo_epi32(V17,V18)
- S4 = _mm_unpackhi_epi32(V17,V18)
接著再做四個函式，會得到等價於轉置的結果
- D1 = _mm_unpacklo_epi64(S1,S3)
- D2 = _mm_unpackhi_epi64(S1,S3)
- D3 = _mm_unpacklo_epi64(S2,S4)
- D4 = _mm_unpackhi_epi64(S2,S4)

4. 將這四個向量相加後，回到 (2.) 改以 V2 與 V5~V8 做內積，以此類推，最終就能得到矩陣乘積的結果。

此處完整的程式碼如下:

void SIMD_matmul(const Matrix l, const Matrix r, const Matrix dst, void *ctx)
{
    for (int i = 0; i < l.row; i += 4)
        for (int j = 0; j < r.col; j += 4)
            SIMD_matmul4(l, r, dst, i, j);
}

void SIMD_matmul4(const Matrix l,
                  const Matrix r,
                  const Matrix dst,
                  int c_row,
                  int c_col)
{
    __m128i I[4], R[4], D[4], T[4], S[4], Sum[4];

    for (int i = 0; i < 4; i++)
        Sum[i] = _mm_setzero_si128();

    for (int k = 0; k < l.col; k += 4) {
        // Read matrix A
        for (int i = 0; i < 4; i++)
            I[i] = _mm_load_si128((__m128i *) (&l.values[i + c_row][k]));

        // Read matrix B
        for (int i = 0; i < 4; i++)
            R[i] = _mm_set_epi32(
                r.values[k][i + c_col], r.values[k + 1][i + c_col],
                r.values[k + 2][i + c_col], r.values[k + 3][i + c_col]);


        for (int i = 0; i < 4; i++) {
            // Inner product of vector from matrix A and B
            for (int j = 0; j < 4; j++)
                T[j] = _mm_mullo_epi32(I[i], R[j]);

            for (int j = 0; j < 2; j++) {
                S[j] = _mm_unpacklo_epi32(T[j * 2], T[j * 2 + 1]);
                S[j + 2] = _mm_unpackhi_epi32(T[j * 2], T[j * 2 + 1]);
            }

            for (int j = 0; j < 2; j++) {
                D[j] = _mm_unpacklo_epi64(S[2 * j], S[2 * j + 1]);
                D[j + 2] = _mm_unpackhi_epi64(S[2 * j], S[2 * j + 1]);
            }

            for (int j = 0; j < 4; j++)
                Sum[i] = _mm_add_epi32(Sum[i], D[j]);
        }
    }

    for (int i = 0; i < 4; i++)
        _mm_store_si128((__m128i *) (&dst.values[c_row + i][c_col]), Sum[i]);
}

Reference :
- [Intel Instrinsics Guide](https://software.intel.com/sites/landingpage/IntrinsicsGuide/)
- [Wiki-SSE](https://zh.wikipedia.org/wiki/SSE)

### SIMD-AVX
- AVX(Advanced Vector Extensions) 指令集為 SSE 指令集的擴充。與 SSE 不同的是其有 16 個 256 bit 的暫存器。
- 既然暫存器的大小是 256 bit，我們可以一次計算 $8*8$ 大小的矩陣。
- 這裡用的方法與 cache-friendly 一樣，==與 sub-matrix 不同==。可以參考 [AVX matrix-multiplication, or something like it](http://www.mindfruit.co.uk/2012/02/avx-matrix-multiplication-or-something.html) 一文的做法。
- 所有的指令可以參考 [Intel Instrinsics Guide](https://software.intel.com/sites/landingpage/IntrinsicsGuide/)，而我們會用到的指令有以下幾個。(記得先 `#include <immintrin.h>`)
- __m256i _mm256_setzero_si256 (void)
將 256 bit 全部設為 0。
- __m256i _mm256_loadu_si256 (__m256i const * mem_addr)
從 mem_addr 讀取 256 bit 的整數資料，不用考慮對齊。
- __m256i _mm256_set1_epi32 (int a)
將 256 bit 內存的 8 個 integer 都設成 `a` 。
- __m256i _mm256_mullo_epi32 (__m256i a, __m256i b)
向量 `a` 與 `b` 內積。
- __m256i _mm256_add_epi32 (__m256i a, __m256i b)
向量 `a` 與 `b` 相加。
- void _mm256_storeu_si256 (__m256i * mem_addr, __m256i a)
將 `a` 的內容存入 MEM[mem_addr] ，不用考慮對齊。

假設我們要計算 $[C]=[A]*[B]$ ，以下說明如何計算 $[C]$ 最上方 row 的結果，類推即可得完整的矩陣$[C]$。

1. 我們將 $[A]$ 、 $[B]$ 矩陣切成 $8*8$ 的 sub-matrix ，然候按照下圖同顏色做乘積。

2. 假設 $[A]_{xy}$ 為 $[A]$ 矩陣中 row x, col y 的元素；$[B]_x$ 為 $[B]$ 矩陣中第 x 個 row 的向量。乘積後我們會有 8 個向量 $V_{1x}=[A]_{1x}*[B]_{x}$ ， $x\in\{1,2,3,4,5,6,7,8\}$。
==$[B]_x=([B]_{x1},[B]_{x2},[B]_{x3},[B]_{x4},[B]_{x5},[B]_{x6},[B]_{x7},[B]_{x8})$==

3. 由上圖可看出，加總 8 個向量即可得 $[C]$ 矩陣第一個 row 的結果，剩下的 row 依序將 $A_{12}, A_{13}...,A_{18}$ 代入即可得。

完整程式碼如下:

void SIMD_AVX_matmul8(const Matrix l,
                      const Matrix r,
                      const Matrix dst,
                      int c_row,
                      int c_col)
{
    __m256i I[8], R[8], S[8], Sum[8];

    for (int i = 0; i < 8; i++)
        Sum[i] = _mm256_setzero_si256();

    for (int k = 0; k < l.col; k += 8) {
        for (int i = 0; i < 8; i++)
            R[i] = _mm256_loadu_si256((__m256i *) (&r.values[k + i][c_col]));


        for (int i = 0; i < 8; i++) {
            for (int j = 0; j < 8; j++) {
                I[j] = _mm256_set1_epi32(l.values[c_row + i][k + j]);
                S[j] = _mm256_mullo_epi32(R[j], I[j]);
            }

            for (int j = 0; j < 8; j++)
                Sum[i] = _mm256_add_epi32(Sum[i], S[j]);
        }
    }

    for (int i = 0; i < 8; i++)
        _mm256_storeu_si256((__m256i *) (&dst.values[c_row + i][c_col]),
                            Sum[i]);
}

void SIMD_AVX_matmul(const Matrix l,
                     const Matrix r,
                     const Matrix dst,
                     void *ctx)
{
    CHECK(l, r, dst);
    for (int i = 0; i < l.row; i += 8)
        for (int j = 0; j < r.col; j += 8)
            SIMD_AVX_matmul8(l, r, dst, i, j);
}

## SSE 指令集與 AVX 指令集的對齊問題
從 [Intel Instrinsics Guide](https://software.intel.com/sites/landingpage/IntrinsicsGuide/) 內對於這兩種指令集讀取連續記憶體的 function 說明中分別寫著
- __m128i _mm_load_si128 (__m128i const* mem_addr)
:::info
Load 128-bits of integer data from memory into dst. mem_addr must be aligned on a ==16-byte== boundary or a general-protection exception may be generated.
:::
- __m256i _mm256_load_si256 (__m256i const * mem_addr)
:::info
Load 256-bits of integer data from memory into dst. mem_addr must be aligned on a ==32-byte== boundary or a general-protection exception may be generated.
:::

兩者對於對齊的要求不同，在我的電腦(x86_64, 64 bit)下，每一次成功的 malloc 都會得到一個對齊 16-byte 的記憶體位置，但不見得會對齊 32-byte 。

從我的電腦連續 malloc 10000 次，每次都向系統要求 1~10 個 `sizeof(int)`，實驗的結果如下圖，每一次都可以得到對齊 16-byte 的記憶體位置，但只有大約一半次數會對齊 32-byte 。

實驗的程式碼如下

#include <stdint.h>                                                                                                                                                                                         
#include <stdio.h>
#include <stdlib.h>

#define VSIZE 10000

struct AlignTester {
    int *ptr;
    int iSize;
};

typedef struct AlignTester Tester;

int main()
{
    int iAlign16 = 0, iAlign32 = 0;
    Tester box[VSIZE];

    for (int i = 0; i < VSIZE; i++) {
        box[i].iSize = (rand() % 10)+1;
        box[i].ptr = (int *) malloc(sizeof(int) * box[i].iSize);
        if (((intptr_t) box[i].ptr) % 16 == 0)
            iAlign16++;
        if (((intptr_t) box[i].ptr) % 32 == 0)
            iAlign32++;
    }

    printf("1\t%.2f\t2\t%.2f\t0.5\n", (double) iAlign16 / VSIZE * 100.0,
           (double) iAlign32 / VSIZE * 100.0);

    for (int i = 0; i < VSIZE; i++)
        free(box[i].ptr);

    return(0);
}

對於這種情況，SSE 指令集不會有問題，但 AVX 指令集會有一半的機率遇到 segmentation fault，解決這個問題可以有兩種做法。
- 改用不需要對齊的 function ，如用上一節使用的 `_mm256_loadu_mi256`。
- 把 `malloc` 換成 `posix_memalign` 函式，可以使取得的空間對齊 32-byte 。但是如果矩陣的大小並不是 8 的倍數，也就是會有 sub-matrix 不是完整的 $8*8$ ，還是會出現 segmentation fault 。

## 當使用 AVX 指令集時，矩陣大小並非 8 的倍數
有兩種情況
- 右方接近邊界，無法填滿 256 個 bit 。
- 下方接近邊界，無法產生 8 個向量加總。

後者較為單純，邏輯不變，只要修改迴圈上限即可；前者需將右方不足的部分填 0 ，將其視為 $8*8$ 的矩陣計算，為了避免讀取超出矩陣範圍造成 segmentation fault ，額外使用以下幾個指令做處理。
- __m256i _mm256_maskload_epi32 (int const* mem_addr, __m256i mask)
- __m256i _mm256_setr_epi32 (int e7, int e6, int e5, int e4, int e3, int e2, int e1, int e0)
- void _mm256_maskstore_epi32 (int* mem_addr, __m256i mask, __m256i a)

可接受矩陣大小非 8 的倍數版本如下:

void SIMD_AVX_matmul8(const Matrix l,
                      const Matrix r,
                      const Matrix dst,
                      int c_row,
                      int c_col)
{
    __m256i I[8], R[8], S[8], Sum[8];
    __m256i mask;

    for (int i = 0; i < 8; i++)
        Sum[i] = _mm256_setzero_si256();

    for (int k = 0; k < l.col; k += 8) {
        for (int i = 0; i < 8 && k + i < r.row; i++) {
            if (c_col + 8 <= r.col)
                R[i] =
                    _mm256_loadu_si256((__m256i *) (&r.values[k + i][c_col]));
            else {
                int iMask[8];

                for (int j = 0; j < 8; j++)
                    iMask[j] = c_col + j < r.col ? 1 << 31 : 0;
                mask =
                    _mm256_setr_epi32(iMask[0], iMask[1], iMask[2], iMask[3],
                                      iMask[4], iMask[5], iMask[6], iMask[7]);
                R[i] = _mm256_maskload_epi32((int *) (&r.values[k + i][c_col]),
                                             mask);
            }
        }


        for (int i = 0; i < 8 && c_row + i < l.row; i++) {
            for (int j = 0; j < 8 && k + j < l.col; j++) {
                I[j] = _mm256_set1_epi32(l.values[c_row + i][k + j]);
                S[j] = _mm256_mullo_epi32(R[j], I[j]);
            }

            for (int j = 0; j < 8 && k + j < l.col; j++)
                Sum[i] = _mm256_add_epi32(Sum[i], S[j]);
        }
    }

    for (int i = 0; i < 8 && c_row + i < l.row; i++) {
        if (c_col + 8 <= r.col)
            _mm256_storeu_si256((__m256i *) (&dst.values[c_row + i][c_col]),
                                Sum[i]);
        else
            _mm256_maskstore_epi32(&dst.values[c_row + i][c_col], mask, Sum[i]);
    }
}

void SIMD_AVX_matmul(const Matrix l,
                     const Matrix r,
                     const Matrix dst,
                     void *ctx)
{
    CHECK(l, r, dst);
    for (int i = 0; i < l.row; i += 8)
        for (int j = 0; j < r.col; j += 8)
            SIMD_AVX_matmul8(l, r, dst, i, j);
}

---

## 實作成果

實驗環境

$ lscpu
Architecture:        x86_64
CPU op-mode(s):      32-bit, 64-bit
Byte Order:          Little Endian
CPU(s):              8
On-line CPU(s) list: 0-7
Thread(s) per core:  2
Core(s) per socket:  4
Socket(s):           1
NUMA node(s):        1
Vendor ID:           GenuineIntel
CPU family:          6
Model:               158
Model name:          Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz
Stepping:            9
CPU MHz:             826.723
CPU max MHz:         3800.0000
CPU min MHz:         800.0000
BogoMIPS:            5616.00
Virtualization:      VT-x
L1d cache:           32K
L1i cache:           32K
L2 cache:            256K
L3 cache:            6144K
NUMA node0 CPU(s):   0-7

以下下實驗矩陣大小均爲 1024

1. make

$ make 

gcc --std=gnu99 -O2 -DNDEBUG -msse4.1 -Wall   -c -o main.o main.c
gcc --std=gnu99 -O2 -DNDEBUG -msse4.1 -Wall -mavx -march=native -c matmul.c
gcc --std=gnu99 -O2 -DNDEBUG -msse4.1 -Wall   -c -o matrix.o matrix.c
gcc --std=gnu99 -O2 -DNDEBUG -msse4.1 -Wall   -c -o strassen.o strassen.c
gcc --std=gnu99 -O2 -DNDEBUG -msse4.1 -Wall -o main main.o matmul.o matrix.o strassen.o

2. `./main`
- 由使用者自行選擇要使用的矩陣乘法

$ ./main 
------------------------------
Print matrix?
0: No
1: Yes
Please input value from 0 to 1: 0
------------------------------
Square matrix?
0: No
1: Yes
Please input value from 0 to 1: 1
------------------------------
Matrix size?
Please input value from 1 to 4096: 1024
Matrix multiplication: (1024 x 1024) x (1024 x 1024)

------------------------------
Choose a matrix multiplication method
0: quit
1: naive
2: cache friendly naive
3: submatrix
4: simd
5: simd_avx
6: strassen
Please input value from 0 to 6: 3
------------------------------
Stride?
Please input value from 2 to 4096: 4
Submatrix(stride=4)
correct!
CPU time: 3.106043 seconds

3. make test
- 提供使用者隨意輸入矩陣大小
- 由程式亂數創出該大小的矩陣
- 得知每個相乘方法的時間長度

$ make test

./main --test < test_script/test.txt
Matrix multiplication: (1024 x 1024) x (1024 x 1024)

Naive
correct!
CPU time: 2.213579 seconds

-----------------------------------------------
Cache_friendly_naive
correct!
CPU time: 0.595068 seconds

-----------------------------------------------
Submatrix(stride=4)
correct!
CPU time: 1.246918 seconds

-----------------------------------------------
Submatrix(stride=6)
correct!
CPU time: 1.169062 seconds

-----------------------------------------------
SIMD
correct!
CPU time: 0.654296 seconds

-----------------------------------------------
SIMD_AVX
correct!
CPU time: 0.238938 seconds

-----------------------------------------------
Strassen(threshold=128)+Submatrix(stride=6)
correct!
CPU time: 0.543757 seconds

-----------------------------------------------
Strassen(threshold=128)+SIMD_AVX
correct!
CPU time: 0.195370 seconds

-----------------------------------------------

4. make perf
用 fork() 的方式，達到測試每個方法的 cache miss
- child 去做 perf 的指令
- 做完之後 parent 就 kill 掉他

**Code:**
- 利用 makefile 在編譯時加上 -DPERF

#if defined(PERF)
        int pid = getpid();
        int cpid = fork();
        if (cpid == 0) {
            // child process .  Run perf stat
            char buf[128];
            sprintf(buf, "perf stat -e cache-misses,cache-references,"
                    "instructions,cycles,branches,branch-misses -p %d", pid);
            execl("/bin/sh", "sh", "-c", buf, NULL);
        }
        setpgid(cpid, 0);
#endif
        clock_t tic = clock();
        matmuls[mm](m, n, o, &matmul_ctx);
        clock_t toc = clock();
#if defined(PERF)
        kill(-cpid, SIGINT);
        wait(&(int){0});
#endif

**Output:**

Matrix multiplication: (1024 x 1024) x (1024 x 1024)


 Performance counter stats for process id '4064':

        61,190,796      cache-misses:u            #    3.626 % of all cache refs    
     1,687,460,171      cache-references:u                                          
     8,428,134,785      instructions:u            #    0.83  insn per cycle         
    10,161,492,593      cycles:u                                                    
     1,055,317,279      branches:u                                                  
         1,028,966      branch-misses:u           #    0.10% of all branches        

       2.735073074 seconds time elapsed

Naive
correct!
CPU time: 2.797990 seconds

###############################################

 Performance counter stats for process id '4064':

         5,209,865      cache-misses:u            #    4.085 % of all cache refs    
       127,541,989      cache-references:u                                          
     7,209,871,205      instructions:u            #    3.39  insn per cycle         
     2,128,038,207      cycles:u                                                    
     1,032,419,939      branches:u                                                  
         1,005,070      branch-misses:u           #    0.10% of all branches        

       0.576083717 seconds time elapsed

Cache-friendly-naive
correct!
CPU time: 0.604449 seconds

###############################################

 Performance counter stats for process id '4064':

        18,270,261      cache-misses:u            #    4.643 % of all cache refs    
       393,530,544      cache-references:u                                          
    12,839,811,376      instructions:u            #    2.56  insn per cycle         
     5,018,055,509      cycles:u                                                    
     2,010,377,946      branches:u                                                  
            77,591      branch-misses:u           #    0.00% of all branches        

       1.354679373 seconds time elapsed

Submatrix(stride=4)
correct!
CPU time: 1.382032 seconds

###############################################

 Performance counter stats for process id '4064':

        25,508,482      cache-misses:u            #    8.113 % of all cache refs    
       314,407,369      cache-references:u                                          
    11,108,525,213      instructions:u            #    2.36  insn per cycle         
     4,704,450,795      cycles:u                                                    
     1,647,827,769      branches:u                                                  
         4,975,617      branch-misses:u           #    0.30% of all branches        

       1.287659037 seconds time elapsed

Submatrix(stride=6)
correct!
CPU time: 1.314957 seconds

###############################################

 Performance counter stats for process id '4064':

        15,755,922      cache-misses:u            #    5.458 % of all cache refs    
       288,694,157      cache-references:u                                          
     4,489,265,404      instructions:u            #    1.55  insn per cycle         
     2,902,179,862      cycles:u                                                    
       522,381,422      branches:u                                                  
            81,643      branch-misses:u           #    0.02% of all branches        

       0.795661427 seconds time elapsed

SIMD
correct!
CPU time: 0.828697 seconds

###############################################

 Performance counter stats for process id '4064':

         6,045,090      cache-misses:u            #    9.784 % of all cache refs    
        61,783,562      cache-references:u                                          
     1,931,992,621      instructions:u            #    1.90  insn per cycle         
     1,017,711,196      cycles:u                                                    
       279,494,672      branches:u                                                  
         1,976,689      branch-misses:u           #    0.71% of all branches        

       0.281439317 seconds time elapsed

SIMD-AVX
correct!
CPU time: 0.312442 seconds

###############################################

## 實驗結果視覺化

### 不同矩陣乘法的運算時間


### Strassen 演算法加速

Strassen 演算法在減少矩陣乘法計算量的同時會帶來一些 overhead，以下以結合 cache-friendly 方法爲例，實驗 Strassen 演算法中 threshold（不再進行切割的矩陣大小） 的選擇對運算時間的影響：


Strassen 結合不同演算法帶來的加速效果：


可以看到 naive 和 cache-friendly 的方法在結合 Strassen 之後時間相差不大，表示在 64*64 的矩陣乘法中，使用 naive 和 cache-friendly 的方法在運算時間上差異不大。

### Submatrix 切不同大小對運算時間的影響



## Matrix Multiplication using SIMD 討論區

> 覺得可以直接先看上一組的影片，講得很清楚。
剛好之前有做過類似的（[matrix-multiplication-parallel](https://github.com/siahuat0727/matrix-multiplication-parallel)，有一些部分還寫不好），所以比較有概念。
另外，因爲 Strassen 只是 divide-and-conquer 後使得計算量減少，裡面還是會進行矩陣乘法，而這裡的乘法應該要可以自由選擇用 naive 或 simd 之類的實作。但上一組的寫法不能方便的做到這些，所以目前把乘法宣告成：
`typedef Matrix (*MatrixMulFunc)(const Matrix, const Matrix, void*)`
其中的 void* 保留給需要的 function （如 Strassen 中用來指定用哪種乘法進行運算以及分割的最小大小，或是將來決定是否平行運算之類的）。
對此大概寫了一點框架 [matrix-multiplication](https://github.com/siahuat0727/matrix-multiplication/tree/a866e5808e7b9a19fff01d62308150765825610f)，不過不確定設計得好不好，有問題可以討論看看～
想說這部分要先達到共識之後才比較方便？
> [name=siahuat0727]

---

MatrixMulFuncEle *matmul_listadd(MatrixMulFuncEle *list, MatrixMulFunc func,
        char *name, void *ctx, size_t extra_size)
{ 
    MatrixMulFuncEle *node = malloc(sizeof(MatrixMulFuncEle) + extra_size);
    return_val_if_fail(node != NULL, 0);
    node->matmul = func;
    node->name = name;
    memcpy(node->ctx, ctx, extra_size);
    return node;
}

> 這裡 extra_size 打算用來做什麼？
> ctx 指的是？
> [name=pjchiou]
> ctx 指 context。
> 先說上一組的 sse 和 sse_prefetch，沒錯的話兩者差別不大，一大段的程式碼是複製貼上的，因爲他們是以一個 function 爲單位，不好修改。
> 所以想說給每種方法保留調整的空間，像是 sse 是否進行 prefetch，那這時需要傳一個 bool 給它。又或是 strassen 要用哪種矩陣乘法還有最小切割單位是多少之類的，這時就可以傳一個 struct 去控制，這些額外的空間就需要 extra_size 去存。
> 因爲最後要串起來，API 要一致。
[name=siahuat0727]
>> ok 明白了
>> [name=pjchiou]

---

Matrix matrix_free(const Matrix m)
{
    free(m.values[0]);
    free(m.values);
    return (Matrix){0};
}

> 這邊返回 `Matrix{0}` 的用意是？
> [name=pjchiou]
> 這裡是參考 [你所不知道的C語言：技巧篇](https://hackmd.io/s/HyIdoLnjl#%E5%BE%9E%E7%9F%A9%E9%99%A3%E6%93%8D%E4%BD%9C%E8%AB%87%E8%B5%B7) 提到的 [prefer-to-return-a-value-rather-than-modifying-pointers](https://github.com/mcinglis/c-style#prefer-to-return-a-value-rather-than-modifying-pointers)。
>:::info
>**Prefer to return a value rather than modifying pointers**
>This encourages immutability, cultivates pure functions, and makes things simpler and easier to understand. It also improves safety by eliminating the possibility of a NULL argument.
>```clike
>// Bad: unnecessary mutation (probably), and unsafe
>void drink_mix( Drink * const drink, Ingredient const ingr ) {
> assert( drink != NULL );
> color_blend( &( drink->color ), ingr.color );
> drink->alcohol += ingr.alcohol;
>}
>
>// Good: immutability rocks, pure and safe functions everywhere
>Drink drink_mix( Drink const drink, Ingredient const ingr ) {
> return ( Drink ){
> .color = color_blend( drink.color, ingr.color ),
> .alcohol = drink.alcohol + ingr.alcohol
> };
>}
>```
>This isn't always the best way to go, but it's something you should always consider.
>:::
>因爲應該要把 matrix 的 row 和 column 也設成 0，所以用的時候會是 `m = matrix_free(m)` 而不是 `matrix_free(&m)`
>其實我也是第一次這麼寫，想說試試看？
>不過我這裡還是少了 `assert(m.values!=NULL)` 或類似的處理。
[name=siahuat0727]
>> 感覺這個範例跟我們 `matrix_free` 的使用時機不大一樣，但我一時說不上哪不同。
>> [name=pjchiou]
>> 嗯嗯我也糾結了很久，都可以討論看看。
>> [name=siahuat0727]
>> 在想不一樣的地方是不是因爲我們的 struct 需要 allocate 空間？
>> 覺得比較麻煩的是用 `Matrix matmul（const Matrix, const Matrix, void*)` 這種寫法每次乘出來的 Matrix 都需要被 free，這有點討厭。
>> 如果用 `void matmul（const Matrix, const Matrix, Matrix *dst, void*)` 的話雖然要去檢查 dst 空間是不是已經配置過了之類的，但就可以一直重複使用同一個 dst，最後不用時再 free，不知道會不會比較好。
>> [name=siahuat0727]
>>> 我認為都可以，隨意選一種就好，我們把重點擺在運算的部份，反正分配記憶體的時間不算在我們要量測的範圍內。我們把其他方法也直接寫在 `naive.[hc]` 內好嗎?
>>> [name=pjchiou]
>>>> 也是，計時是夾整個 function，第一種方法會把分配的時間算進去，那我們就用第二種然後 function 裡假設 dst 已經被分配了而且大小是正確的吧。合在一起感覺大家一起寫的時候有點麻煩，我先去改成第二種的，如果覺得還是合在一起比較好的話我也都 ok。
>>>> [name=siahuat0727]
---