cublas矩阵乘

神经网络中有大量的矩阵乘法运算，使用cuda来进行矩阵的乘法运算，可以大大提高神经网络的训练速度，于是学着使用cuda，由于NVIDIA已经提供了非常好的矩阵运算库cublas，所以应该是学着使用cublas，在使用中遇到了一些问题，记录一下，方便以后的查询。

cublas中执行矩阵乘法运算的函数主要是：

cublasSgemm /*用来处理单精度矩阵，也就是float型的*/

cublasDgemm /*用来处理双精度矩阵，也就是double型的*/

首先要注意的是cublas使用的是以列为主的存储方式，和c/c++中的以行为主的方式是不一样的，例如对于一个c/c++中的5*10的矩阵，将其转化为一个一维数组以后，[0,1]个数据 在新的一维数组中的位置是[1] /*从0开始计数*/，但将其转化为以列为主是它在一位数组中的位置是[11]，于是在对数组进行转化是需要进行一个小小的转化，在nvidia提供的文档中已经提供了响应的宏

#define  IDX2C(i,j,leading) (((j)*(leading))+(i))

　　使用这个宏就可以很方便的将我们习惯的行为主的数据转化为列为主的数据了

现在来看看来进行矩阵乘的函数，

cublasSgemm(handle,
                      CUBLAS_OP_N,
                        CUBLAS_OP_N,
                        mat1->_rows,
                        mat2->_cols,
			mat2->_rows,
                        &alpha,
                        d_a,
                        mat1->_rows,
                        d_b,
                        mat2->_rows,
                        &beta,
                        d_c,
                         mat1->_rows);        

　　根据nvidia的文档，该函数主要进行的运算是

C=alpha*F(a)*F(b)+beta*C

第2,3个参数表示是否对矩阵进行转置，第三个参数是矩阵1的行数，第四个参数是矩阵2的列数，第五个参数是矩阵2的行数，第八个参数是矩阵1的leading dimension（我不知道这个怎么翻译，主维数吗），第十个参数是矩阵2的leading dimension，第十三个参数是输出矩阵也就是矩阵C3的主维数，其实也就是矩阵1的主维数。cublasDgemm和cublasSgemm的参数是一样，使用方式是一样的。

  1 #include <algorithm>
  2 #include <iostream>
  3 #include <time.h>
  4 #include <cublas.h>
  5 #include <cublas_v2.h>
  6 #include <assert.h>
  7 using namespace std;
  8 
  9 
 10 #define  IDX2C(i,j,leading) (((j)*(leading))+(i))
 11 
 12 
 13 typedef struct _data *PDATA;
 14 typedef struct _data
 15 {
 16     int _rows;
 17     int _cols;
 18     float *data;
 19 } Data;
 20 
 21 void free_mat(PDATA mat)
 22 {
 23     free(mat->data);
 24     free(mat);
 25 }
 26 
 27 PDATA mat_product(PDATA mat1,PDATA mat2)
 28 {
 29     if(mat1->_cols!=mat2->_rows)
 30     {
 31         printf("this is not right\n");
 32             return NULL;
 33     }
 34     PDATA mat3=new Data;
 35     mat3->data=(float *)malloc(sizeof(float)*(mat1->_rows)*(mat2->_cols));
 36     mat3->_rows=mat1->_rows;
 37     mat3->_cols=mat2->_cols;
 38     /*
 39      *INIT the matrix we want calculate 
 40      * col primary
 41      */
 42     {
 43         float *d_a,*d_b,*d_c;
 44         cublasInit();
 45         cublasAlloc((mat1->_cols)*(mat1->_rows),sizeof(float),(void **)&d_a);
 46         cublasAlloc((mat2->_cols)*(mat2->_rows),sizeof(float),(void **)&d_b);
 47         cublasAlloc((mat3->_rows)*(mat3->_cols),sizeof(float),(void **)&d_c);
 48         cudaMemcpy(d_a,mat1->data,sizeof(float)*(mat1->_cols)*(mat1->_rows),cudaMemcpyHostToDevice);
 49         cudaMemcpy(d_b,mat2->data,sizeof(float)*(mat2->_rows)*(mat2->_cols),cudaMemcpyHostToDevice);
 50         cublasHandle_t handle;
 51         cublasCreate(&handle);
 52         float alpha=1.0;
 53         float beta=0.0;
 54         cublasSgemm(handle,CUBLAS_OP_N,CUBLAS_OP_N,mat1->_rows,mat2->_cols,
 55             mat2->_rows,&alpha,d_a,mat1->_rows,d_b,mat2->_rows,&beta,d_c,mat1->_rows);
 56         cudaMemcpy(mat3->data,d_c,sizeof(float)*(mat3->_rows)*(mat3->_cols),cudaMemcpyDeviceToHost);
 57         cublasShutdown();
 58     }
 59     /* need to trans the mat3*/
 60     return mat3;
 61 }
 62 
 63 void ele_mat_show(PDATA mat)
 64 {
 65     for (int i=0;i<mat->_rows;i++){
 66         for (int j=0;j<mat->_cols;j++){
 67             cout<<mat->data[IDX2C(i,j,mat->_rows)]<<"\t";
 68         }
 69         cout<<endl;
 70     }
 71 }
 72 float myrand()
 73 {
 74     return rand()%10;
 75 }
 76 int main()
 77 {
 78     clock_t start,end;
 79 
 80 #if 0
 81     for (int i=0;i<M*N;i++)
 82     {
 83         cout<<c[i]<<"\t";
 84     }
 85     cout<<endl;
 86 #endif
 87 
 88     PDATA mat1,mat2,mat3;
 89     /* remember to initialize the point*/
 90     mat1=(PDATA)malloc(sizeof(Data));
 91     mat2=(PDATA)malloc(sizeof(Data));
 92     mat3=(PDATA)malloc(sizeof(Data));
 93     mat1->_rows=5;
 94     mat1->_cols=10;
 95     mat1->data=(float *)malloc(sizeof(float)*mat1->_rows*mat1->_cols);
 96     for (int i=0;i<5;i++)
 97         for (int j=0;j<10;j++)
 98             mat1->data[IDX2C(i,j,mat1->_rows)]=i*j;
 99     //generate(mat1->data,mat1->data+(mat1->_cols)*(mat1->_rows),myrand);
100     ele_mat_show(mat1);
101     mat2->_rows=10;
102     mat2->_cols=2;
103     mat2->data=(float *)malloc(sizeof(float)*mat2->_rows*mat2->_cols);
104     for (int i=0;i<10;i++)
105         for (int j=0;j<2;j++)
106             mat2->data[IDX2C(i,j,mat2->_rows)]=i*j;
107     //generate(mat2->data,mat2->data+(mat2->_cols)*(mat2->_rows),myrand);
108     ele_mat_show(mat2);
109     mat3=mat_product(mat1,mat2);
110     for (int i=0;i<mat3->_rows;i++)
111     {
112         for (int j=0;j<mat3->_cols;j++)
113         {
114             cout<<mat3->data[i+j*mat3->_rows]<<'\t';
115         }
116         cout<<endl;
117     }
118 
119     return 0;
120 }