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 }