---
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。
我們都知道 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 發揮效用。
實際執行結果如下，$Speedup = \frac{17820600}{8369778} = 2.13$
```shell=
$ ./submatrix
submatrix: 8369778 us
```
### SSE version
**程式碼**
因為太長了 可以去github看
這裡講想法跟解釋一些用到的function
想法是先把兩個矩陣切成很多4x4的小矩陣，把對應的位置先做好矩陣相乘後再相加
* matrix1 一次橫向讀入128bit = 四個int 存進I0 ，只要給要讀入陣列的起始位置，他會往後抓127bit, I1,I2,I3同理
```
__m128i I0 = _mm_loadu_si128((__m128i *)(src1 + (x + 0) * src1_w + k));
```
* 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]);
```
* 將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);
```
* 我們來比較一下效能
```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 => 完美
* 在亂數範圍為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
```
* 出現次數可以看的出來蠻均勻的
## 學習 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`迴圈走訪每個演算法。
{%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)
Sign in with Wallet
Connect another wallet