矩阵乘法优化示例(访存优化和SIMD的使用)
这个矩阵乘法优化的整个代码都是根据张先轶老师一个关于openblas的讲座写出来的。https://www.leiphone.com/news/201704/Puevv3ZWxn0heoEv.html
不多说了直接上代码，一个矩阵类Matrix，按列存储
#include <iostream>
#include <cstdlib>
#include <immintrin.h>
using namespace std;
class Matrix{
private:
double **_Matrix;
size_t _Row,_Column;
public:
Matrix():_Matrix(nullptr),_Row(0),_Column(0){}//默认构造
Matrix(size_t r,size_t c):_Row(r),_Column(c){//构造r行、c列的矩阵
if(!_Column||!_Row) return;
_Matrix=(double**)malloc(_Column*sizeof(double*));
double **p=_Matrix,**end=_Matrix+_Column;
do{
*(p++)=(double*)malloc(_Row*sizeof(double));
}while(p!=end);
}
Matrix(size_t r,size_t c,const double init):_Row(r),_Column(c){//构造r行、c列的矩阵并用init初始化
if(!_Column||!_Row) return;
_Matrix=(double**)malloc(_Column*sizeof(double*));
double **pr=_Matrix,**endr=_Matrix+_Column,*p,*end;
do{
p=*(pr++)=(double*)malloc(_Row*sizeof(double));
end=p+_Row;
do{
*(p++)=init;
}while(p!=end);
}while(pr!=endr);
}
Matrix(const Matrix& B){//拷贝构造
_Row=B._Row;
_Column=B._Column;
_Matrix=(double**)malloc(_Column*sizeof(double*));
double **pbr=B._Matrix,**endbr=B._Matrix+_Column,*pb,*endb,
**par=_Matrix,**endar=_Matrix+_Column,*pa,*enda;
do{
pa=*(par++)=(double*)malloc(_Row*sizeof(double));
enda=pa+_Row;
pb=*(pbr++);
endb=pb+_Row;
do{
*(pa++)=*(pb++);
}while(pa!=enda);
}while(par!=endar);
}
~Matrix(){//析构
if(!_Matrix) return;
double **p=_Matrix,**end=_Matrix+_Column;
do{
free(*(p++));
}while(p!=end);
_Column=_Row=0;
free(_Matrix);
}
double& operator()(size_t i,size_t j){return _Matrix[j][i];}//访问第i行、第j列的元素
const double operator()(size_t i,size_t j)const{return _Matrix[j][i];}//访问第i行、第j列的元素
Matrix& operator=(Matrix&& B){//移动赋值
if(_Matrix){
double **p=_Matrix,**end=_Matrix+_Row;
do{
free(*(p++));
}while(p!=end);
free(_Matrix);
}
_Row=B._Row;
_Column=B._Column;
_Matrix=B._Matrix;
B._Matrix=nullptr;
return *this;
}
};
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
接下来就是实现乘法运算了，习惯上讲我们都是按行计算乘法结果的，即计算结果时外层循环是行号，内层循环是列号，乘积是一行一行计算的。
for(int i=0;i<M;i++){
for(int j=0;j<N;j++){
for(int k=0;k<C;k++) c(i,j)+=a(i,k)*b(k,j);
}
}
1
2
3
4
5
但是这里为了展示访存对计算效率的影响，我们先按列计算。
//这个函数以及后面的函数都是Matrix的成员函数，应当写在Matrix类里。这里把它单独拆出来了，如果要编译测试，这段代码需要放回Matrix类中才能通过。
Matrix multi1(const Matrix &B){
if(_Column!=B._Row) return *this;
Matrix tmp(_Row,B._Column,0);
int i,j(0),k;
do{
i=0;
do{
k=0;
do{
tmp(i,j)+=(*this)(i,k)*B(k,j);
k++;
}while(k<_Column);
i++;
}while(i<_Row);
j++;
}while(j<B._Column);
return tmp;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
好了测试一下效率，这里所有测试全部使用O3优化
#define N 960
int main(int argc,char* argv[]){
Matrix A(N,N,1),B(N,N,2);
double dt;
clock_t start=clock();
A.multi1(B);
dt=static_cast<double>(clock()-start)/CLOCKS_PER_SEC;
cout<<"用时：\n"<<dt<<"s"<<endl;
return 0;
}
1
2
3
4
5
6
7
8
9
10
运行结果是
用时：
3.38886s
1
2
接下来我们换一个计算方式，我们按行计算
Matrix multi2(const Matrix &B){
if(_Column!=B._Row) return *this;
Matrix tmp(_Row,B._Column,0);
int i(0),j,k;
do{
j=0;
do{
k=0;
do{
tmp(i,j)+=(*this)(i,k)*B(k,j);
k++;
}while(k<_Column);
j++;
}while(j<B._Column);
i++;
}while(i<_Row);
return tmp;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
运行结果是
用时：
2.24073s
1
2
其实我们什么多余的工作都没有做，只是把按列求解改成了按行求解，但是用时就从3.38886s减少到了2.24073s。按列求解时列号在外层循环，每一次内层循环中行号先改变，也就是说每求完一个乘积都会重新读取一行。而按列存储的矩阵的每一行在内存中都是不连续的，所以行访问会很慢，这样反复多次进行行访问就会大大拖慢运行速度。按行求解时行号在外层循环，每求完一次乘积会重读一列，但之前读出的一行是不变的。这时需要的那一行可能还在Cache里，所以可以不需要再从那段不连续的内存中读取了。
接下来我们再做一个简单优化，按行求解时我们每一列都是通过下标访问的，这样都会调用()重载，但实际上每一列都是一个连续的一维数组，因此这里我们用一个指针来代替()重载。
Matrix multi3(const Matrix &B){
if(_Column!=B._Row) return *this;
Matrix tmp(_Row,B._Column,0);
int i(0),j(0),k;
double **pr,*p;
do{
j=0;
pr=B._Matrix;
do{
k=0;p=*(pr++);
do{
tmp(i,j)+=(*this)(i,k)**(p++);
k++;
}while(k<_Column);
j++;
}while(j<B._Column);
i++;
}while(i<_Row);
return tmp;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
运行结果是
用时：
2.24051s
1
2
很遗憾，并没有很大的提升，可能编译器已经做过这部分优化了，所以我们这么写并没有什么效果。
我们知道X86-64的CPU都有16个寄存器，而矩阵乘法本来就是个叠加的过程，所以我们可以申请16个寄存器变量，每次计算都直接叠加到上面，从而将乘积的结果16个16个的算出来。这样我们读取的行列都可以复用很多次，从而减少访存的时间。为了使每次读取的行和列都可以多次复用，我们将乘积按照4×4的小块来计算。这样每次从A中取出四行，B中取出四列，然后每行每列都会复用四次，这个过程中所有的数据复用都是在Cache中提取的，因此可以有效减少访存的次数。由前面的例子可以知道乘积在按行计算时效率更高，因此在子块的计算也按行来做，也就是说每算完一个4×4的子块，就接着计算右侧紧挨着的4×4的子块。这时我们的乘法就需要一个kernel函数来协助完成了。
另外由于本身行访问是不连续的，但是每一行又是一直在复用的，直到这一行在计算中不再使用。那么我们还可以使用一种packing的打包优化。我们每计算四行都把对应的数据打包到四个一维数组中，当乘积中关于这四行的所有项计算完之后我们再打包下面四行，如此反复，就可以让每一行的读取都变成连续的了，这又可以进一步提高性能。
void multi4kernel(double **c,double **a,double **b,int row,int col){
register double t0(0),t1(0),t2(0),t3(0),t4(0),t5(0),t6(0),t7(0),
t8(0),t9(0),t10(0),t11(0),t12(0),t13(0),t14(0),t15(0);
double *a0(a[0]),*a1(a[1]),*a2(a[2]),*a3(a[3]),
*b0(b[col]),*b1(b[col+1]),*b2(b[col+2]),*b3(b[col+3]),*end=b0+_Row;
do{
t0+=*(a0)**(b0);
t1+=*(a0)**(b1);
t2+=*(a0)**(b2);
t3+=*(a0++)**(b3);
t4+=*(a1)**(b0);
t5+=*(a1)**(b1);
t6+=*(a1)**(b2);
t7+=*(a1++)**(b3);
t8+=*(a2)**(b0);
t9+=*(a2)**(b1);
t10+=*(a2)**(b2);
t11+=*(a2++)**(b3);
t12+=*(a3)**(b0++);
t13+=*(a3)**(b1++);
t14+=*(a3)**(b2++);
t15+=*(a3++)**(b3++);
}while(b0!=end);
c[col][row]=t0;
c[col+1][row]=t1;
c[col+2][row]=t2;
c[col+3][row]=t3;
c[col][row+1]=t4;
c[col+1][row+1]=t5;
c[col+2][row+1]=t6;
c[col+3][row+1]=t7;
c[col][row+2]=t8;
c[col+1][row+2]=t9;
c[col+2][row+2]=t10;
c[col+3][row+2]=t11;
c[col][row+3]=t12;
c[col+1][row+3]=t13;
c[col+2][row+3]=t14;
c[col+3][row+3]=t15;
}
Matrix multi4(const Matrix &B){
if(_Column!=B._Row) return *this;
Matrix tmp(_Row,B._Column,0);
double *tr[4];
int i(0),j(0);
do{
tr[i++]=(double*)malloc(sizeof(double)*_Column);
}while(i<4);
do{
i=0;
do{
tr[0][i]=_Matrix[i][j];//packing过程，把行数据打包到连续空间
tr[1][i]=_Matrix[i][j+1];
tr[2][i]=_Matrix[i][j+2];
tr[3][i]=_Matrix[i][j+3];
}while((++i)<_Column);
i=0;
do{
multi4kernel(tmp._Matrix,tr,B._Matrix,j,i);
i+=4;
}while(i<B._Column);
j+=4;
}while(j<_Row);
return tmp;
}
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
运行结果是
用时：
0.348799s
1
2
从3.38886s到0.348799s性能提高了10倍。
由于这套代码是优化比较好的，所以我用的这套代码演示。实际上我之前还写过另一套代码，那里有一个只使用4×4的kernel函数计算的一个乘法函数，和一个在kernel函数基础上在进行行packing的乘法函数，用来单独探究packing和kernel的优化效果，不过那个代码有些问题这里就没有贴，从那套代码的测试结果来看，packing对乘法效率的提高更加显著一些。
到这里访存上已经没有什么可以优化的，要说再优化就是用一种block的方式，根据Cache的大小把A、B进行分块计算。因为矩阵很大的时候四行、四列未必能完全存在Cache中，因此需要把A、B分区，得到合适大小的子块，然后运算。基本原理就是把A切分成[A1,A2,A3…]这种行向量的形式，而B则分成[B1;B2;B3…]这种列向量的形式，这时乘法就可以分区进行了，从而减小了对Cache容量的需求。这部分代码我没有仔细想过具体实现，所以也没写。
既然访存不能优化了那就只能从计算上来优化了，这里就要用到前面说到的SIMD技术了，我们利用SSE指令集来加速乘法的运算。
void multi5kernel(double **c,double **a,double **b,int row,int col){
__m128d t01_0,t01_1,t01_2,t01_3,t23_0,t23_1,t23_2,t23_3,
a0,a1,b0,b1,b2,b3;
t01_0=t01_1=t01_2=t01_3=t23_0=t23_1=t23_2=t23_3=_mm_set1_pd(0);
double *pb0(b[col]),*pb1(b[col+1]),*pb2(b[col+2]),*pb3(b[col+3]),*pa0(a[0]),*pa1(a[1]),*endb0=pb0+_Column;
do{
a0=_mm_load_pd(pa0);
a1=_mm_load_pd(pa1);
b0=_mm_set1_pd(*(pb0++));
b1=_mm_set1_pd(*(pb1++));
b2=_mm_set1_pd(*(pb2++));
b3=_mm_set1_pd(*(pb3++));
t01_0+=a0*b0;
t01_1+=a0*b1;
t01_2+=a0*b2;
t01_3+=a0*b3;
t23_0+=a1*b0;
t23_1+=a1*b1;
t23_2+=a1*b2;
t23_3+=a1*b3;
pa0+=2;
pa1+=2;
}while(pb0!=endb0);
_mm_store_pd(&c[col][row],t01_0);
_mm_store_pd(&c[col+1][row],t01_1);
_mm_store_pd(&c[col+2][row],t01_2);
_mm_store_pd(&c[col+3][row],t01_3);
_mm_store_pd(&c[col][row+2],t23_0);
_mm_store_pd(&c[col+1][row+2],t23_1);
_mm_store_pd(&c[col+2][row+2],t23_2);
_mm_store_pd(&c[col+3][row+2],t23_3);
}
Matrix multi5(const Matrix &B){
if(_Column!=B._Row) return *this;
Matrix tmp(_Row,B._Column,0);
double *ta[2];
ta[0]=(double*)malloc(sizeof(double)*2*_Column);
ta[1]=(double*)malloc(sizeof(double)*2*_Column);
tb=(double*)malloc(sizeof(double)*4*B._Row);
end=tb+4*B._Row;
int i(0),j(0),k,t;
do{
k=0;i=0;
do{
ta[0][k]=_Matrix[i][j];
ta[1][k++]=_Matrix[i][j+2];
ta[0][k]=_Matrix[i][j+1];
ta[1][k++]=_Matrix[i++][j+3];
}while(i<_Column);
i=0;
do{
multi5kernel(tmp._Matrix,ta,B._Matrix,j,i);
i+=4;
}while(i<B._Column);
j+=4;
}while(j<_Row);
free(tb);
free(ta[0]);
free(ta[1]);
return tmp;
}
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
这里就在前面的优化基础上将kernel中的计算通过SIMD来执行，从而提高计算速度。
运行结果是
用时：
0.233691s
1
2
我们看到相比于前面的优化SIMD又提高了33%的效率，这还是128bit的情况。
使用AVX2.0可以一次处理256bit的数据，即4个双精度浮点数
void multi6kernel(double **c,double **a,double **b,int row,int col){
__m256d t04_0,t04_1,t04_2,t04_3,t58_0,t58_1,t58_2,t58_3,
a0,a1,b0,b1,b2,b3;
t04_0=t04_1=t04_2=t04_3=t58_0=t58_1=t58_2=t58_3=_mm256_set1_pd(0);
double *pb0(b[col]),*pb1(b[col+1]),*pb2(b[col+2]),*pb3(b[col+3]),*pa0(a[0]),*pa1(a[1]),*endb0=pb0+_Column;
do{
a0=_mm256_loadu_pd(pa0);
a1=_mm256_loadu_pd(pa1);
b0=_mm256_set1_pd(*(pb0++));
b1=_mm256_set1_pd(*(pb1++));
b2=_mm256_set1_pd(*(pb2++));
b3=_mm256_set1_pd(*(pb3++));
t04_0+=a0*b0;
t04_1+=a0*b1;
t04_2+=a0*b2;
t04_3+=a0*b3;
t58_0+=a1*b0;
t58_1+=a1*b1;
t58_2+=a1*b2;
t58_3+=a1*b3;
pa0+=4;
pa1+=4;
}while(pb0!=endb0);
_mm256_storeu_pd(&c[col][row],t04_0);
_mm256_storeu_pd(&c[col+1][row],t04_1);
_mm256_storeu_pd(&c[col+2][row],t04_2);
_mm256_storeu_pd(&c[col+3][row],t04_3);
_mm256_storeu_pd(&c[col][row+4],t58_0);
_mm256_storeu_pd(&c[col+1][row+4],t58_1);
_mm256_storeu_pd(&c[col+2][row+4],t58_2);
_mm256_storeu_pd(&c[col+3][row+4],t58_3);
}
Matrix multi6(const Matrix &B){
if(_Column!=B._Row) return *this;
Matrix tmp(_Row,B._Column,0);
double *ta[2];
ta[0]=(double*)malloc(sizeof(double)*4*_Column);
ta[1]=(double*)malloc(sizeof(double)*4*_Column);
int i(0),j(0),k,t;
do{
k=0;i=0;
do{
ta[0][k]=_Matrix[i][j];
ta[1][k++]=_Matrix[i][j+4];
ta[0][k]=_Matrix[i][j+1];
ta[1][k++]=_Matrix[i][j+5];
ta[0][k]=_Matrix[i][j+2];
ta[1][k++]=_Matrix[i][j+6];
ta[0][k]=_Matrix[i][j+3];
ta[1][k++]=_Matrix[i++][j+7];
}while(i<_Column);
i=0;
do{
multi6kernel(tmp._Matrix,ta,B._Matrix,j,i);
i+=4;
}while(i<B._Column);
j+=8;
}while(j<_Row);
free(ta[0]);
free(ta[1]);
return tmp;
}
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
测试一下
用时：
0.105505s
1
2
可以看到使用AVX2.0之后效率有提高了将近一倍。如果使用支持AVX-512的CPU，并采用对应的指令集，可能效率提高的更多。
我们再把最开始那个最慢的结果拿过来比较一下，从3.38886s到0.105505s，效率提高了30多倍，基本上从C++这个层面上的矩阵乘法优化已经再没有什么可以做的了。剩下的优化就是在汇编层面上进行一定的汇编代码的执行顺序优化，根据张先轶老师的说法，这一步可以进一步提高10%左右的效率，但是做起来就比较麻烦了。如果只是自己写着玩也没必要做这么细了，如果真的对这种线性代数运算的效率要求十分苛刻，那么直接使用openblas或者intel的数学核心库可能是更好的选择。
上一篇：数值计算优化方法C/C++(三)——SIMD
下一篇：数值计算优化方法C/C++(五)——矩阵转置优化示例(访存优化和SIMD的使用)
————————————————
版权声明：本文为CSDN博主「遗世独立的理想乡_」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/artorias123/article/details/86527456
763
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
