SIMD向量化运算
隨著機器學習等人工智能技術的飛速發展，矩陣乘法的應用越來越多，intel芯片先后提供了不同系列的向量指令，包括mmx、sse、avx等，支持simd操作。后來為了更好地支持矩陣乘法，又增加了fma（Fused Multiply-Add）指令。fma指令需要三個向量參數va,vb,vcva,vb,vc，其效果等價于表達式(va?vb)+vc(va?vb)+vc，其中的乘法和加法都是面向向量中的元素的，也就是fma指令的結果是一個同樣長度的向量。fma指令的出現為矩陣乘法提供了方便，但是其效果同樣可以用avx指令系列中的乘法和加法的組合來實現，本文使用例子來分析不同向量指令在矩陣乘中的性能和精度。?
例子主要計算了一個矩陣WW和向量xx的乘積，WW的列數等于xx的長度，結果仍然是一個向量，長度等于WW的行數。代碼的實現如下。
#include <stdio.h>
#include <time.h>
#include <x86intrin.h>
int main() {
? const int col = 1024, row = 64, num_trails = 1000000;
? float w[row][col];
? float x[col];
? float y[row];
? float scratchpad[8];
? for (int i=0; i<row; i++) {
? ? for (int j=0; j<col; j++) {
? ? ? w[i][j]=(float)(rand()%1000)/800.0f;
? ? } ??
? }
? for (int j=0; j<col; j++) {
? ? x[j]=(float)(rand()%1000)/800.0f;
? }
? clock_t t1, t2;?
// The original matrix multiplication version
? t1 = clock();
? for (int r = 0; r < num_trails; r++)
? ? for(int j = 0; j < row; j++)
? ? { ??
? ? ? float sum = 0;
? ? ? float *wj = w[j];
? ? ? for(int i = 0; i < col; i++)
? ? ? ? sum += wj[i] * x[i];
? ? ? y[j] = sum;
? ? } ??
? t2 = clock();
? float diff = ((float)t2 - (float)t1) / CLOCKS_PER_SEC;
? printf("\nTime taken: %.2f second.\n", diff);
? for (int i=0; i<row; i++) {
? ? printf("%.4f, ", y[i]);
? }
? printf("\n");
// The avx matrix multiplication version.
? const int col_reduced_8 = col - col % 8;
? __m256 op0, op1, tgt, tmp_vec;
? t1 = clock();
? for (int r = 0; r < num_trails; r++)
? ? for (int i=0; i<row; i++) {
? ? ? float res = 0;
? ? ? tgt = _mm256_setzero_ps();
? ? ? for (int j = 0; j < col_reduced_8; j += 8) {
? ? ? ? op0 = __builtin_ia32_loadups256(&x[j]);
? ? ? ? op1 = __builtin_ia32_loadups256(&w[i][j]);
? ? ? ? tmp_vec = __builtin_ia32_mulps256(op0, op1);
? ? ? ? tgt = __builtin_ia32_addps256(tmp_vec, tgt);
? ? ? }
? ? ? __builtin_ia32_storeups256(scratchpad, tgt);
? ? ? for (int k=0; k<8; k++)
? ? ? ? res += scratchpad[k];
? ? ? for (int l=col_reduced_8; l<col; l++) {
? ? ? ? res += w[i][l] * x[l];
? ? ? }
? ? ? y[i] = res;
? ? }
? t2 = clock();
? diff = ((float)t2 - (float)t1) / CLOCKS_PER_SEC;
? printf("\nTime taken: %.2f second.\n", diff);
? for (int i=0; i<row; i++) {
? ? printf("%.4f, ", y[i]);
? }
? printf("\n");
?// The fma matrix multiplication version.
? t1 = clock();
? for(int r = 0; r < num_trails; r++)
? ? for(int i = 0; i < row; i++)
? ? {
? ? ? float rlt = 0;
? ? ? tgt = _mm256_setzero_ps();
? ? ? for(int j = 0; j < col_reduced_8; j += 8)
? ? ? {
? ? ? ? op0 = __builtin_ia32_loadups256(&x[j]);
? ? ? ? op1 = __builtin_ia32_loadups256(&w[i][j]);
? ? ? ? tgt = _mm256_fmadd_ps(op0, op1, tgt);
? ? ? }
? ? ? __builtin_ia32_storeups256(scratchpad, tgt);
? ? ? for(int k = 0; k < 8; k++)
? ? ? {
? ? ? ? rlt += scratchpad[k];
? ? ? }
? ? ? for(int l = col_reduced_8; l < col; l++)
? ? ? {
? ? ? ? rlt += w[i][l] * x[l];
? ? ? }
? ? ? y[i] = rlt;
? ? }
? t2 = clock();
? diff = ((float)t2 - (float)t1) / CLOCKS_PER_SEC ;
? printf("\nTime taken: %.2f second.\n", diff);
? for(int i=0; i<row; i++)
? {
? ? printf("%.4f, ", y[i]);
? }
? printf("\n");?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
在ubuntu系統中，程序的編譯命令是:?
gcc -O2 -mfma test.c -o test?
需要注意的是，只有在支持fma的芯片結構下，程序才能夠執行。可以通過命令：?
cat /proc/cpuinfo | grep fma?
來判斷芯片是否支持fma。?
其執行結果為：?
Time taken: 93.56 second.?
409.8341, 413.4546, 398.7332, 399.8303, 404.1195, 402.3861, 394.6979, 412.6429, 409.0014, 390.9019, 400.3911, 392.7900, 400.5019, 418.6781, 399.3336, 404.0719, 414.9839, 411.6887, 396.0086, 406.6972, 384.5781, 399.3724, 400.0473, 391.6383, 401.3511, 400.8543, 418.4066, 406.6425, 405.5102, 408.4534, 403.0285, 406.3510, 410.2005, 414.9617, 417.3602, 406.4511, 397.1705, 406.1265, 393.3314, 407.1777, 389.9053, 397.3145, 401.7866, 413.3134, 415.7482, 414.2341, 403.3439, 405.4922, 395.4076, 399.6389, 409.6675, 419.8184, 412.3336, 399.8252, 403.3434, 387.4861, 402.2747, 399.8241, 414.1568, 405.4861, 406.6151, 410.4040, 408.9755, 398.9610,
Time taken: 10.94 second.?
409.8341, 413.4549, 398.7335, 399.8304, 404.1191, 402.3860, 394.6979, 412.6424, 409.0016, 390.9022, 400.3909, 392.7900, 400.5020, 418.6781, 399.3336, 404.0718, 414.9842, 411.6884, 396.0087, 406.6971, 384.5780, 399.3723, 400.0472, 391.6382, 401.3510, 400.8541, 418.4067, 406.6424, 405.5103, 408.4536, 403.0287, 406.3513, 410.2007, 414.9618, 417.3603, 406.4513, 397.1708, 406.1266, 393.3315, 407.1776, 389.9049, 397.3150, 401.7864, 413.3134, 415.7483, 414.2341, 403.3439, 405.4922, 395.4075, 399.6392, 409.6674, 419.8183, 412.3336, 399.8253, 403.3433, 387.4865, 402.2746, 399.8239, 414.1567, 405.4861, 406.6153, 410.4034, 408.9752, 398.9612,
Time taken: 12.08 second.?
409.8341, 413.4549, 398.7335, 399.8304, 404.1191, 402.3860, 394.6979, 412.6424, 409.0016, 390.9022, 400.3909, 392.7900, 400.5021, 418.6781, 399.3336, 404.0718, 414.9842, 411.6884, 396.0087, 406.6971, 384.5780, 399.3722, 400.0472, 391.6382, 401.3510, 400.8541, 418.4067, 406.6424, 405.5102, 408.4536, 403.0287, 406.3513, 410.2007, 414.9618, 417.3603, 406.4513, 397.1708, 406.1266, 393.3315, 407.1776, 389.9050, 397.3150, 401.7864, 413.3134, 415.7483, 414.2341, 403.3439, 405.4922, 395.4075, 399.6392, 409.6674, 419.8183, 412.3336, 399.8253, 403.3433, 387.4865, 402.2746, 399.8239, 414.1568, 405.4861, 406.6153, 410.4034, 408.9752, 398.9612,
可見，avx對乘加的組合實現性能還略高于fma指令。而精度兩者相似，略低于原始的運算。
---------------------?
作者：softee?
來源：CSDN?
原文：https://blog.csdn.net/softee/article/details/55057476?
版權聲明：本文為博主原創文章，轉載請附上博文鏈接！
總結
- 上一篇： 哪些亚比锁属性的啊
- 下一篇： makefile 中 =, :=, ?=