HackMD
Matrix Multiplication using SIMD
contributed by <0140454>, <linachiu>, <kaizsv>, <HuangWenChen>, <aweimeow>
預期目標
- 參照 software-pipeline，以 SSE/AVX 實作矩陣乘法
- 針對記憶體操作提出改進方案 (注意 cache 的影響)
- 提出測試矩陣乘法正確性和效能之有效機制
產出
GitHub: matrix-multiplication
YouTube Video: Matrix Multiplication using SIMD; Nov 16, 2016
研究報告: 研究報告
待改善項目
- 提出測試矩陣乘法正確性和效能之有效機制
- 學習 Linux crypto API，可以切換不同的矩陣乘法實作
實作
Naive version
程式碼
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 矩陣相乘的時間。
$ ./naive
naive: 17820600 us
Sub-matrix version
在實作 SIMD 前，先實作另外一個版本。
相較於 naive 而言，這個版本會將整個矩陣分為 4x4 的小矩陣進行運算。
因為拆成小矩陣運算的關係，src2 的 access pattern 變化如下，這樣會有比較大的機會讓 cache 發揮效用。
實際執行結果如下，
$ ./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);
- 我們來比較一下效能
$ ./sse
sse: 4070998 us
- 在加上 prefetch 的版本
$ ./sse_prefetch
sse_prefetch: 3880293 us
因為prefetch可以幫我們預測接下來會用到的cache，並且預先載入，所以效能提升。
改 prefetch 第二個矩陣 後，效能有所提升。
另外，如果用 AVX 中提到的向量方式，直接運算而不轉置。可以發現其執行時間又會更短一些。
程式碼可以到 branch sse_v2 觀看。
吳勃興
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。
另外，小矩陣相乘的時候並不會先進行轉置的動作，而是透過向量的方式直接運算。
主要用到的函數有
-
_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中的資料複製到指定的記憶體位置
而在效能方面
$ ./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。
吳勃興
為了讓 user space 的程式也可以使用相關的 API，所以在 Linux 2.6.38 時引入了 netlink-based crypto API。
Netlink 是一種特殊的 socket，透過他可以讓 process 與 kernel 進行通訊。而 user space 的 crypto API 便是在 kernel 中透過 af_alg.c 和 algif_hash.c 來實作他。
可以使用
cat /proc/crypto來觀看有哪些演算法可以使用。
吳勃興
要學習的不是 cyrpto API 的用法，而是如何切換 algorithm 的「手法」，對應到本例就是「如何預先定義存取介面，讓許多不同的矩陣乘法演算法得以實作指定的介面，並且讓測試程式不需要事先得知有幾套 algorithm provider」 jserv
了解！這一部分正在規劃
吳勃興
參考raytracing object.c[1]的寫法，在benchmark.c開始的時候為每個乘法演算法 append 到串列，這裡使用__attribute__ ((constructur))提醒GCC在main()開始之前將每個演算法加到一個 linked list 中。如此就可以用簡單的for each迴圈走訪每個演算法。
閱讀資料
Memory Hierarchies and Optimizing Matrix Multiplication
- Understanding Caches
- Optimizing Matrix Multiplication
Optimizing Parallel Multiplication Operation for Rectangular and Transposed Matrices
- paper study
- OpenMPI
Optimizing matrix multiplication
參考資料
- Intel Intrinsics Guide
- SSE 引入的標頭檔
- Memory part 5: What programmers can do
- AVX matrix-multiplication, or something like it
- Aligned vs. unaligned memory access
- strassen-演算法
- 常見的 LaTeX 語法
- CheHsuan 共筆
- 隨機亂數驗證
