EIE稀疏矩阵乘法硬件模拟

#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]]$，其行号则为$x[p[j],p[j+1]]$，当然论文中采用的是相对索引存储方式，行号的获取略有差异，然后将这些非零元素和激活进行相乘，将结果累加在对应行的累加器上，在上述代码中表现为
$out[x[k]]+=in[j]\ast v[k]k=p[j],p[j]+1,...,p[j+1]-1.$
需要注意的是，第$i$个$PE$存储的是第$i,i+N\_PE,i+2\ast N\_PE,...$等行的权重。