#include<iostream>
#include<cstdlib>
#include<ctime>
#include <iomanip>
using namespace std;
typedef int dtype;
typedef struct CSC{
int a; //存储非零元素个数
int n; //存储列数
dtype* v; //存储非零值
int* x; //存储每一列中非零值与上一个非零值之间0的个数
int* p; //存储每一列之前非零元素的个数
}csc_t;
void gen_sparse_matrix(int n,dtype** & mat,double s){
mat=new int*[n];
for(int i=0;i<n;i++)
mat[i]=new int[n];
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
{
int x=rand()%1000;
if(x>1000*s)
mat[i][j]=0;
else
mat[i][j]=x+1;
}
return;
}
void print_matrix(int n,dtype** mat){
int i,j;
for(i=0;i<n;i++){
for(j=0;j<n;j++)
cout<<setw(8)<<setfill(' ')<<mat[i][j]<<",";
cout<<endl;
}
}
void print_csc(csc_t* S){
cout<<"NNZ="<<S->a<<endl;
cout<<"Column Num is "<<S->n<<endl;
cout<<"value of non-zero element\n";
for(int i=0;i<S->a;i++)
cout<<S->v[i]<<",";
cout<<endl;
cout<<"rowindex of value\n";
for(int i=0;i<S->a;i++)
cout<<S->x[i]<<",";
cout<<endl;
cout<<"column pointer\n";
for(int i=0;i<S->n+1;i++)
cout<<S->p[i]<<",";
cout<<endl;
return;
}
void storage_matrix(int N_PE,int n,dtype** mat,csc_t* & S){
S=new csc_t[N_PE];
for(int i=0;i<N_PE;i++){
//为每个PE存储相应的值,包括矩阵中的i,i+N_PE,i+2N_PE,...等行。
int cnt=0;
int row=i;
while(row<n){
for(int j=0;j<n;j++)
if(mat[row][j]!=0)
cnt++;
row+=N_PE;
}
S[i].a=cnt;
S[i].n=n;
S[i].v=new dtype[cnt];
S[i].x=new int[cnt];
S[i].p=new int[n+1];
int k=0; //记录当前已经访问的非零元素个数
int l=0; //记录当前已经访问的列的个数
for(int j=0;j<n;j++){
S[i].p[l++]=k;
for(int r=i;r<n;r+=N_PE)
if(mat[r][j]!=0){
S[i].v[k]=mat[r][j];
S[i].x[k]=r;
k++;
}
}
S[i].p[l]=S[i].a;
}
return;
}
void print_PE(int n,dtype** mat,int i,int N_PE){
for(int r=i;r<n;r+=N_PE){
for(int j=0;j<n;j++)
cout<<setw(8)<<setfill(' ')<<mat[r][j]<<",";
cout<<endl;
}
return;
}
void spmv(int N_PE,csc_t* S,dtype* in,dtype* out){
int n=S[0].n;
for(int i=0;i<n;i++)
out[i]=0;
for(int j=0;j<n;j++)
{
if(in[j]!=0)
{
for(int i=0;i<N_PE;i++)
{
//寻找第j列的非零元素,即v[p[j]:p[j+1]],行号分别是x[p[j]:p[j+1]]
for(int k=S[i].p[j];k<S[i].p[j+1];k++)
out[S[i].x[k]]+=S[i].v[k]*in[j];
}
}
}
return;
}
void gpmv(int n,dtype** mat,dtype* in,dtype* out){
for(int i=0;i<n;i++)
{
out[i]=0;
for(int j=0;j<n;j++)
out[i]+=mat[i][j]*in[j];
}
return;
}
int main(){
int n;
int N_PE;
double s;
cout<<"input n:";
cin>>n;
cout<<"input s:";
cin>>s;
cout<<"input N_PE:";
cin>>N_PE;
dtype** M=NULL;
csc_t* S=NULL;
gen_sparse_matrix(n,M,s);
print_matrix(n,M);
storage_matrix(N_PE,n,M,S);
for(int i=0;i<N_PE;i++){
cout<<"*************************************"<<endl<<"N_PE="<<i<<endl;
print_PE(n,M,i,N_PE);
print_csc(&S[i]);
}
//定义一个输入向量和输出向量
dtype* in=new dtype[n];
dtype* out1=new dtype[n];
dtype* out2=new dtype[n];
for(int i=0;i<n;i++){
in[i]=rand()%10-5;
out1[i]=0;
out2[i]=0;
}
//普通矩阵向量乘法
gpmv(n,M,in,out1);
spmv(N_PE,S,in,out2);
for(int i=0;i<n;i++)
cout<<out1[i]<<" , "<<out2[i]<<endl;
return 0;
}
该代码是对论文EIE: Efficient Inference Engine on Compressed Deep Neural Network所述的初步模拟,还有待完善。该算法使用CSC存储稀疏矩阵,并采用多个PE单元进行并行计算,具体的方法为:
读入一个激活向量,从该向量中扫描,若扫描到非零元素,则将其值和列号j同时送入所有PE单元,PE单元则查找自己存储的属于第j列的非零元素,即
v
[
p
[
j
]
:
p
[
j
+
1
]
]
v[p[j]:p[j+1]]
v[p[j]:p[j+1]],其行号则为
x
[
p
[
j
]
,
p
[
j
+
1
]
]
x[p[j],p[j+1]]
x[p[j],p[j+1]],当然论文中采用的是相对索引存储方式,行号的获取略有差异,然后将这些非零元素和激活进行相乘,将结果累加在对应行的累加器上,在上述代码中表现为
o
u
t
[
x
[
k
]
]
+
=
i
n
[
j
]
∗
v
[
k
]
k
=
p
[
j
]
,
p
[
j
]
+
1
,
.
.
.
,
p
[
j
+
1
]
−
1.
out[x[k]]+=in[j]*v[k]\,\,\,\,\,\,\,\,k=p[j],p[j]+1,...,p[j+1]-1.
out[x[k]]+=in[j]∗v[k]k=p[j],p[j]+1,...,p[j+1]−1.
需要注意的是,第
i
i
i个
P
E
PE
PE存储的是第
i
,
i
+
N
_
P
E
,
i
+
2
∗
N
_
P
E
,
.
.
.
i,i+N\_PE,i+2*N\_PE,...
i,i+N_PE,i+2∗N_PE,...等行的权重。