2018q3 矩陣乘法 Matrix multiplication

Contributed by < DyslexiaS, siahuat0727, pjchiou >
GitHub
Task : 重現 Matrix Multiplication using SIMD 實驗，並依循 CS:APP 第 6 章指引，分析個別實作的效能

outline

演算法
SSE 指令集與 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 的特性來加速運算，把握兩個重點

Temporal locality : 一旦取用過的變數，重複取用。
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 的概念，先把矩陣拆成較小的矩陣，各自運算後再加起來。

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 →

在這個示意圖中，A、B兩個矩陣都是正方形的，其 sub-matrix 行列皆為原本的一半。這與 naive 版本的功能不符，所以我們將其稍做調整

A、B 矩陣不必為正方形，符合矩陣相乘的限制 A.col==B.row 即可。
sub-matrix 行列大小不必為 N/2 ，但必須為正方形。
sub-matrix 大小亦不必可被 C.row 或 C.col 整除。

time complexity :
$O (n^{3})$
演算法
reference : divide and conquer


























//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$ 階分塊矩陣乘法
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 →

其中每一分塊皆為方陣。傳統的矩陣乘法運算方式，

C_{i j} = A_{i 1} B_{1 j} + A_{i 2} B_{2 j}

，總共使用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 = 2 m$ ，上面所有分塊皆為
$m \times m$ 。
傳統的矩陣乘法 :
- 使用了
  $8 m^{3}$ 個乘法和
  $8 m^{3} - 4 m^{2}$ 個加法。
Strassen :
- Strassen 演算法僅執行至分塊等級，也就是說，
  $m \times m$ 分塊乘法仍採用傳統方法計算
- 則7次分塊乘法共使用
  $7 m^{3}$ 個乘法和
  $7 (m^{3} - m^{2})$ 個加法，再加上 18 次分塊加法使用的
  $18 m^{2}$ 個加法，總計是
  $7 m^{3}$ 個乘法和
  $7 m^{3} + 11 m^{2}$ 個加法
- 若
  $m ≫ 1$ ，不論乘法或加法，Strassen 演算法僅耗費傳統方法
  $7 / 8$ 的計算量。

Reference : 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 。

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 →

我們要做的事情就是從 SIMD 指令集中組合出可以達到需求的指令。以下說明如何得到最上面的 row ，其他以此類推。

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

讀取矩陣值，存入暫存器。宣告型別為 __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 矩陣用後者。

將 V1 分別對 V5~V8 做內積，可得四個內積後的向量 V15, V16, V17, V18

__m128i _mm_mullo_epi32 (__m128i a, __m128i b)。

我們要將這四個向量 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)

將這四個向量相加後，回到 (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 :

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 一文的做法。
所有的指令可以參考 Intel Instrinsics Guide，而我們會用到的指令有以下幾個。(記得先 #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]

。

我們將
$[A]$ 、
$[B]$ 矩陣切成
$8 * 8$ 的 sub-matrix ，然候按照下圖同顏色做乘積。
假設
$[A]_{x y}$ 為
$[A]$ 矩陣中 row x, col y 的元素；
$[B]_{x}$ 為
$[B]$ 矩陣中第 x 個 row 的向量。乘積後我們會有 8 個向量
$V_{1 x} = [A]_{1 x} * [B]_{x}$ ，
$x \in {1, 2, 3, 4, 5, 6, 7, 8}$ 。

$[B]_{x} = ([B]_{x 1}, [B]_{x 2}, [B]_{x 3}, [B]_{x 4}, [B]_{x 5}, [B]_{x 6}, [B]_{x 7}, [B]_{x 8})$
由上圖可看出，加總 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 內對於這兩種指令集讀取連續記憶體的 function 說明中分別寫著

__m128i _mm_load_si128 (__m128i const* mem_addr)

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)

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

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

./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

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

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

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，有一些部分還寫不好），所以比較有概念。
另外，因爲 Strassen 只是 divide-and-conquer 後使得計算量減少，裡面還是會進行矩陣乘法，而這裡的乘法應該要可以自由選擇用 naive 或 simd 之類的實作。但上一組的寫法不能方便的做到這些，所以目前把乘法宣告成：
typedef Matrix (*MatrixMulFunc)(const Matrix, const Matrix, void*)
其中的 void* 保留給需要的 function （如 Strassen 中用來指定用哪種乘法進行運算以及分割的最小大小，或是將來決定是否平行運算之類的）。
對此大概寫了一點框架 matrix-multiplication，不過不確定設計得好不好，有問題可以討論看看～
想說這部分要先達到共識之後才比較方便？
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 指的是？
pjchiou
ctx 指 context。
先說上一組的 sse 和 sse_prefetch，沒錯的話兩者差別不大，一大段的程式碼是複製貼上的，因爲他們是以一個 function 爲單位，不好修改。
所以想說給每種方法保留調整的空間，像是 sse 是否進行 prefetch，那這時需要傳一個 bool 給它。又或是 strassen 要用哪種矩陣乘法還有最小切割單位是多少之類的，這時就可以傳一個 struct 去控制，這些額外的空間就需要 extra_size 去存。
因爲最後要串起來，API 要一致。
siahuat0727

ok 明白了
pjchiou






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

這邊返回 Matrix{0} 的用意是？
pjchiou
這裡是參考你所不知道的C語言：技巧篇提到的 prefer-to-return-a-value-rather-than-modifying-pointers。
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.
// 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) 或類似的處理。
siahuat0727

感覺這個範例跟我們 matrix_free 的使用時機不大一樣，但我一時說不上哪不同。
pjchiou
嗯嗯我也糾結了很久，都可以討論看看。
siahuat0727
在想不一樣的地方是不是因爲我們的 struct 需要 allocate 空間？
覺得比較麻煩的是用 Matrix matmul（const Matrix, const Matrix, void*) 這種寫法每次乘出來的 Matrix 都需要被 free，這有點討厭。
如果用 void matmul（const Matrix, const Matrix, Matrix *dst, void*) 的話雖然要去檢查 dst 空間是不是已經配置過了之類的，但就可以一直重複使用同一個 dst，最後不用時再 free，不知道會不會比較好。
siahuat0727

我認為都可以，隨意選一種就好，我們把重點擺在運算的部份，反正分配記憶體的時間不算在我們要量測的範圍內。我們把其他方法也直接寫在 naive.[hc] 內好嗎?
pjchiou

也是，計時是夾整個 function，第一種方法會把分配的時間算進去，那我們就用第二種然後 function 裡假設 dst 已經被分配了而且大小是正確的吧。合在一起感覺大家一起寫的時候有點麻煩，我先去改成第二種的，如果覺得還是合在一起比較好的話我也都 ok。
siahuat0727