Matrix Multiplication using SIMD

--- tags: sysprog21-hw, 0140454, linachiu, kaizsv, HuangWenChen, aweimeow --- # Matrix Multiplication using SIMD contributed by <`0140454`>, <`linachiu`>, <`kaizsv`>, <`HuangWenChen`>, <`aweimeow`> ## 預期目標 * 參照 [software-pipeline](https://hackmd.io/s/ry7eqDEC)，以 SSE/AVX 實作矩陣乘法 * 針對記憶體操作提出改進方案 (注意 cache 的影響) * 研讀 [Memory Hierarchies and Optimizing Matrix Multiplication](https://people.eecs.berkeley.edu/~demmel/cs267_Spr99/Lectures/Lect_02_1999b.pdf), [Optimizing matrix multiplication](http://web.cs.ucdavis.edu/~bai/ECS231/optmatmul.pdf), [Optimizing Parallel Multiplication Operation for Rectangular and Transposed Matrices](http://hpc.pnl.gov/srumma/srumma-icpads.pdf) * 提出測試矩陣乘法正確性和效能之有效機制 ## 產出 GitHub: [matrix-multiplication](https://github.com/0140454/matrix-multiplication) YouTube Video: [Matrix Multiplication using SIMD](https://www.youtube.com/watch?v=3rU6BX7w8Tk&list=PLKT8ER2pEV3umVSMwd06LY_eSIX-DnU6A); Nov 16, 2016 研究報告: [研究報告](https://hackmd.io/IYBgpmAmBMCsDsBaWBjeA2RAWAzJSiARlgBwiLQCMew8AnAGaxYopA==) ## 待改善項目 - [x] 提出測試矩陣乘法正確性和效能之有效機制 - [x] 學習 Linux crypto API，可以切換不同的矩陣乘法實作 ## 實作 ### Naive version **程式碼** ```clike= void naive_multiply(int *src1, int *src2, int *dst, int src1_w, int src1_h, int src2_w, int src2_h) { for (int i = 0; i < src1_h; ++i) { for (int j = 0; j < src2_w; ++j) { dst[i * src2_w + j] = 0; for (int k = 0; k < src2_h; ++k) { dst[i * src2_w + j] += src1[i * src1_w + k] * src2[k * src2_w + j]; } } } } ``` 從程式碼中可以看出，下面兩張圖分別為 `src1` 與 `src2` 的 access pattern。 ![](https://i.imgur.com/JChvCUm.png) 我們都知道 C 是 row-major 的語言，所以 `src1` 在 cache 的部份是有發揮效益的。但在 `src2` 的部份，則會因為一直以 column order 的方式去讀取，導致 cache 沒辦法有效地被利用。因此，接下來的工作除了要使用 SSE 和 AVX 來實作外，也要注意 cache 相關的問題。 **結果** 執行兩個 1024x1024 矩陣相乘的時間。 ```shell= $ ./naive naive: 17820600 us ``` ### Sub-matrix version 在實作 SIMD 前，先實作另外一個版本。相較於 naive 而言，這個版本會將整個矩陣分為 4x4 的小矩陣進行運算。因為拆成小矩陣運算的關係，`src2` 的 access pattern 變化如下，這樣會有比較大的機會讓 cache 發揮效用。 ![](https://i.imgur.com/KD8uENc.png) 實際執行結果如下，$Speedup = \frac{17820600}{8369778} = 2.13$ ```shell= $ ./submatrix submatrix: 8369778 us ``` ### SSE version **程式碼** 因為太長了可以去github看這裡講想法跟解釋一些用到的function 想法是先把兩個矩陣切成很多4x4的小矩陣，把對應的位置先做好矩陣相乘後再相加 ![](https://i.imgur.com/0XJ6usb.gif) * matrix1 一次橫向讀入128bit = 四個int 存進I0 ，只要給要讀入陣列的起始位置，他會往後抓127bit, I1,I2,I3同理 ``` __m128i I0 = _mm_loadu_si128((__m128i *)(src1 + (x + 0) * src1_w + k)); ``` ![](https://i.imgur.com/qs3XLYI.png) * matrix2 一次縱向讀入128bit = 四個int 存進I4 ，因為陣列縱向的記憶體位置不連續，所以要給他四個int， I5,I6,I7同理 ``` __m128i I4 = _mm_set_epi32 (src2[(k+3) * src2_w + y], src2[(k+2) * src2_w + y], src2[(k+1) * src2_w + y], src2[k * src2_w + y]); ``` ![](https://i.imgur.com/KT5Yig6.png) * 將I0與I4,I5,I6,I7相乘 , I1,I2,I3同理 ``` __m128i T0 = _mm_mullo_epi32(I0, I4); __m128i T1 = _mm_mullo_epi32(I0, I5); __m128i T2 = _mm_mullo_epi32(I0, I6); __m128i T3 = _mm_mullo_epi32(I0, I7); ``` 先講解_mm_mullo_epi32(m128i, m128i) 他是 SSE4.1 的function 所以要 #include \<smmintrin.h\> 原理: ``` _mm_mullo_epi32(m128i, m128i) FOR j := 0 to 1 i := j*64 dst[i+63:i] := a[i+31:i] * b[i+31:i] ENDFOR ``` I0*I4 這時候我們得到T0的值會是 T0(M1 00 * M2 00 , M1 01 * M2 10 , M1 02 * M2 20 , M1 03 * M2 30) 而我們最終要放在 matrix 00位置的值應該要是(M1 00 * M2 00 + M1 01 * M2 10 + M1 02 * M2 20 +...+ M1 0n * M2 n0) * 所以在這裡將他轉置 ``` __m128i T16 = _mm_unpacklo_epi32(T0, T1); __m128i T17 = _mm_unpacklo_epi32(T2, T3); __m128i T18 = _mm_unpackhi_epi32(T0, T1); __m128i T19 = _mm_unpackhi_epi32(T2, T3); __m128i T20 = _mm_unpacklo_epi64(T16, T17); __m128i T21 = _mm_unpackhi_epi64(T16, T17); __m128i T22 = _mm_unpacklo_epi64(T18, T19); __m128i T23 = _mm_unpackhi_epi64(T18, T19); ``` * 就可以把他們相加起來 ``` T20 = _mm_add_epi32(T20,T21); T20 = _mm_add_epi32(T20,T22); T20 = _mm_add_epi32(T20,T23); des0 = _mm_add_epi32(T20,des0); ``` ![](https://i.imgur.com/uw2w9rP.png) * 我們來比較一下效能 ```shell= $ ./sse sse: 4070998 us ``` * 在加上 prefetch 的版本 ```shell= $ ./sse_prefetch sse_prefetch: 3880293 us ``` 因為prefetch可以幫我們預測接下來會用到的cache，並且預先載入，所以效能提升。 >> 改 [prefetch 第二個矩陣](https://github.com/0140454/matrix-multiplication/commit/df7d9e8f7051bfca5f4e57e5428638f59f7c9ea5) 後，效能有所提升。 >> >> ```shell= >> sse: 2864652 us >> >> sse_prefetch: 1989581 us >> ``` >> >> 另外，如果用 AVX 中提到的向量方式，直接運算而不轉置。可以發現其執行時間又會更短一些。 >> >> 程式碼可以到 [branch sse_v2](https://github.com/0140454/matrix-multiplication/tree/sse_v2) 觀看。 >> >> ```shell= >> sse: 2310845 us >> >> sse_prefetch: 1358132 us >> ``` >> [name=吳勃興] ``` prefetch src2: 7755,9769 L1-dcache-load-misses prefetch src1: 1,0313,4778 L1-dcache-load-misses without prefetch: 1,1071,8285 L1-dcache-load-misses ``` ### AVX version 因為 AVX 沒有提供完整的 integer 操作函數，又電腦剛好支援 AVX2，所以這一部份就改用 AVX2 來實作。概念和上面提及的方式差不多，都是將大矩陣分割成小矩陣來運算，因為 AVX 可以同時對 256-bit 的資料做操作，也就是 8 個 integer。因此，在這邊小矩陣的大小為 8x8。另外，小矩陣相乘的時候並不會先進行轉置的動作，而是透過[向量的方式](https://zh.wikipedia.org/wiki/%E7%9F%A9%E9%99%A3%E4%B9%98%E6%B3%95#.EF.BC.8D.E5.90.91.E9.87.8F.E6.96.B9.E6.B3.95.3D)直接運算。主要用到的函數有 * _mm256_setzero_si256 將 `__m256i` 內的元素全部填零 * _mm256_loadu_si256 從指定記憶體位置中讀取 256-bit 的資料到 `__m256i` * _mm256_set1_epi32 將 `__m256i` 中以某個數字填充 * _mm256_mullo_epi32 這是 AVX2 才引入的函數，對 `__m256i` 做乘法 * _mm256_add_epi32 也是 AVX2 才引入的函數，對 `__m256i` 做加法 * _mm256_storeu_si256 將 `__m256i` 中的資料複製到指定的記憶體位置而在效能方面 ```shell= $ ./avx avx: 873441 us $ ./avx_prefetch avx_prefetch: 807615 us ``` ### 各版本初步比較 | 版本 | naive | sub-matrix | sse | sse_prefetch | sse_prefetch_modify | | --- | ----- | ---------- | --- | ------------ | -------------- http://tiny.pl/gp71n <<<<-------------->>>> http://forumloverz.lark.ru http://periscopetv.lark.ru 六、0.3 <= P <= 0.7 EXCELLENT ( 完美 ) ``` * 實際運算: 1. 先求出$X^2$ = 32663 2. 將$X^2$帶入公式計算$Q_{x^2,d} \approx 0.65$ 3. P = 0.35 => 完美 ![](https://i.imgur.com/t5mGluJ.png) * 在亂數範圍為0~32767下 $32550 \leq X^2 \leq 32980$ 為不錯的集中度 ($0.202 \leq Q_{x^2,d} \leq 0.801$ 用[wolfram](http://www.wolframalpha.com)測試的) 2. 亂數獨立度驗證 * RUNS-UP-AND-DOWN TEST RUN：一個序列中分享所有性質的元件群 RUNS UP：一個序列中所有的點都上昇 RUNS DOWN：一個序列中所有的點都下降例如： ``` 0.9501 0.2311 0.6068 0.4860 0.8913 0.7621 0.4565 0.0185 0.8214 0.4447 - + - + - - - + - ``` RUNS UP 的數量：(用「+」來表示 ) = 3 RUNS DOWN 的數量：(用「-」來表示 ) = 4 所有 RUNS 的總數： = 7 RUNS UP 的長度： = 1,1,1 RUNS DOWN 的長度： = 1,1,3,1 APPROXIMATE TEST: $\mu + Z_1\sigma \le R \le \mu + Z_2\sigma$ $\mu = (2n-1)/3$ $\sigma^2 = (16n-29)/90$ R：RUNS 的數量 $Z_1$和 $Z_2$：分別是 -1.96 和 1.96 當 $Z$ < -1.96 時，就表示有太少的 RUNS；當 $Z$ > +1.96 時，就表示有太多的 RUNS。 * 再 1024 x 1024 的矩陣下產生的 runs 次數 ``` runs: min runs = 698204.091697 runs = 698984 max runs = 699896.574969 ``` * 出現次數可以看的出來蠻均勻的 ![](https://i.imgur.com/hssnwNm.png) ## 學習 crypto API 並設計介面使演算法容易擴充和比較在 Linux 的核心中有許多的 cryptography kernel module 可以使用，在開發核心的應用時，如果需要加解密的話，可以使用這些 module 來加速。 >> 可以使用 `ls /lib/modules/$(uname -r)/kernel/crypto` 來觀看有哪些 kernel module。 >> [name=吳勃興] 為了讓 user space 的程式也可以使用相關的 API，所以在 Linux 2.6.38 時引入了 netlink-based crypto API。 Netlink 是一種特殊的 socket，透過他可以讓 process 與 kernel 進行通訊。而 user space 的 crypto API 便是在 kernel 中透過 [af_alg.c](http://lxr.free-electrons.com/source/crypto/af_alg.c) 和 [algif_hash.c](http://lxr.free-electrons.com/source/crypto/algif_hash.c) 來實作他。 >> 可以使用 `cat /proc/crypto` 來觀看有哪些演算法可以使用。 >> [name=吳勃興] >> 要學習的不是 cyrpto API 的用法，而是如何切換 algorithm 的「手法」，對應到本例就是「如何預先定義存取介面，讓許多不同的矩陣乘法演算法得以實作指定的介面，並且讓測試程式不需要事先得知有幾套 algorithm provider」 [name=jserv] >> 了解！這一部分正在規劃 >> [name=吳勃興] 參考`raytracing object.c`[^first]的寫法，在`benchmark.c`開始的時候為每個乘法演算法 append 到串列，這裡使用`__attribute__ ((constructur))`提醒GCC在`main()`開始之前將每個演算法加到一個 linked list 中。如此就可以用簡單的`for each`迴圈走訪每個演算法。 ![](https://i.imgur.com/84EjImI.png) {%gist kaizsv/6c706ce04e4db402dfc2d26c1ea9f660%} ## 閱讀資料 ### [Memory Hierarchies and Optimizing Matrix Multiplication](https://hackmd.io/s/r1Tb8hOel) * Understanding Caches * Optimizing Matrix Multiplication ### [Optimizing Parallel Multiplication Operation for Rectangular and Transposed Matrices](https://hackmd.io/s/rk_mvfree) * paper study * OpenMPI ### [Optimizing matrix multiplication](https://hackmd.io/s/BJ4mw_reg) ## 參考資料 * [Intel Intrinsics Guide](https://software.intel.com/sites/landingpage/IntrinsicsGuide/) * [SSE 引入的標頭檔](http://www.aiuxian.com/article/p-1028293.html) * [Memory part 5: What programmers can do](https://lwn.net/Articles/255364/) * [AVX matrix-multiplication, or something like it](http://www.mindfruit.co.uk/2012/02/avx-matrix-multiplication-or-something.html) * [Aligned vs. unaligned memory access](http://www.alexonlinux.com/aligned-vs-unaligned-memory-access) * [strassen-演算法](https://ccjou.wordpress.com/2013/06/04/%E5%88%86%E6%B2%BB%E7%9F%A9%E9%99%A3%E4%B9%98%E6%B3%95%E2%94%80%E2%94%80strassen-%E6%BC%94%E7%AE%97%E6%B3%95/) * [常見的 LaTeX 語法](http://mohu.org/info/symbols/symbols.htm) * [CheHsuan 共筆](/s/BkbydgVa) * [隨機亂數驗證](http://www.islab.iecs.fcu.edu.tw/Subject/892004.pdf) [^first]:[raytracing_object.c](https://github.com/sysprog21/raytracing/blob/master/objects.c)

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.