利用cuda的cublas库实现任意矩阵的乘法

使用cublas库的相关函数做矩阵乘法时，应当注意两点：一是cublas的中数据是按列存储的，二是，要注意矩阵是在GPU端产生的还是CPU端产生的。cublas库实现矩阵乘法运算的函数

cublasSgemm和cublasDgemm.

在cublas库中，实现矩阵的乘法有以上两个函数，分别完成单精度和双精度的矩阵乘法运算。完成的是以下运算（具体可以参看SDK文档）

C=alpha*op(A)op(B)+beta*C

特别小心：ld(leading dimension)的意思，对于矩阵在由cublas的随机函数数在GPU端产生的话，其值为矩阵的一列（因为按列存储）所包含的元素的个数（即取值为矩阵的rows），若矩阵为CPU端产生的数据，则ld为矩阵的一行所包含的元素的个数（即取值为矩阵的columns）.

实例：

矩阵A：3x5,B:5x6,C=AB（C为3x6的矩阵）

若A，B的元素均由cublas库的curand产生的随机数构成，调用函数cublasSgemm的话，可以得到C，但是要注意C按列存储的，将C传回GPU的，若要打印结果，注意要转换过来。否则打印出来是C的转置。

若A，B的元素均在CPU端由随机函数产生，再传到GPU端做运算调用函数时C=AB，这个时候其实执行的是A^T*B^T，即A的转置和B的转置相乘，显然不对。

若函数调用时做C=BA的计算，即将函数中A和B调换位置，那么实际是C^T=B^T*A^T(得到6x3的矩阵即C的转置），同时注意函数调用时lda,ldb,ldc对应的是的6,5,6（矩阵的列数）transa==CUBLAS_OP_N.另外的做法就是A，B的位置不要调换，而是参数transa==CUBLAS_OP_T，那么调用函数时执行的是A*B.

下面代码是自己做的一个100x200和200x110的矩阵的乘法，矩阵是从外部导入到CPU的，在传到GPU上，利用cublas库的函数cublasSgemm实现矩阵的乘法

代码不够漂亮仅当初学之用

/*********************************************
该程序的功能是借用cublas来做矩阵并行运算,C=AB
考虑到cublas中矩阵是按列存储的
调用cublasSgemm时是按C=BA调用,实际执行的是B^T x A^T ,由于按列存储，这样乘积得到的矩阵刚好是按行存储的C，打印结果是就很方便了
*********************************************/
#include
   
   
    
    
#include
    
    
     
     
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include
     
     
      
      
#include 
      
      
       
       
#include 
       
       
         #include 
        
          #include 
         
           #include"printMatrix.h" #include"initCUDA.h" #include"load.h" #define A_ROW 100 //随意定义A和B的size #define A_COL 200 #define B_COL 110 #define IDX2C(i,j,leading) (((j)*(leading))+(i)) #define CUDA_SAFE_CALL(call) do { cudaError_t err = call; \ if(err != cudaSuccess) { \ printf("Error at %s:%d\nError: %s\n",__FILE__,__LINE__, cudaGetErrorString(err)); \ exit(EXIT_FAILURE);}} while(0) void gpu_blas_mmul(cublasHandle_t &handle, float *C, const float *A, const float *B, const int m, const int k, const int n); int main() { if(!InitCUDA()) return 0; //分配cpu上的内存 float *host_A=(float*) malloc(sizeof(float)*A_ROW*A_COL); float *host_B=(float*) malloc(sizeof(float)*A_COL*B_COL); float *host_C=(float*) malloc(sizeof(float)*A_ROW*B_COL); //从外部导入数据 load("D:\\E\\myprogramm\\vs2010\\mycudaexercises\\cudaStart\\matrixA.txt",host_A,A_ROW*A_COL); load("D:\\E\\myprogramm\\vs2010\\mycudaexercises\\cudaStart\\matrixB.txt",host_B,A_COL*B_COL); //GPU上内存分配 float *dev_A=0,*dev_B=0,*dev_C=0; clock_t startTime, finishTime; /*size_t A_pitch=0; size_t B_pitch=0; size_t C_pitch=0;*/ cudaMalloc((void**)&dev_A,sizeof(float)* A_ROW* A_COL ); cudaMalloc((void**)&dev_B,sizeof(float)* A_COL* B_COL ); cudaMalloc((void**)&dev_C,sizeof(float)* A_ROW* B_COL ); cublasHandle_t handle;//handle，表示是cublas的上下文句柄,初始化cublas cublasCreate(&handle);//计算前create handle CUDA_SAFE_CALL(cudaMemcpy(dev_A,host_A,sizeof(float)* A_ROW* A_COL,cudaMemcpyHostToDevice)); CUDA_SAFE_CALL(cudaMemcpy(dev_B,host_B,sizeof(float)* A_COL* B_COL,cudaMemcpyHostToDevice)); //cublasSetMatrix(A_ROW, A_COL, sizeof(float),host_A, A_ROW, dev_A,A_ROW); //cublasSetMatrix(A_COL, B_COL, sizeof(float),host_B, A_COL, dev_A, A_COL); startTime = clock();//计算时用time.h库中的clock()函数 gpu_blas_mmul(handle,dev_C,dev_A, dev_B, A_ROW,B_COL, A_COL); cudaDeviceSynchronize(); finishTime=clock(); float time =float((finishTime - startTime))/CLOCKS_PER_SEC; printf("time to compute on gpu is ：%.5f second\n",time);// //传回数据 cudaMemcpy(host_C,dev_C,sizeof(float)* A_ROW* B_COL,cudaMemcpyDeviceToHost); //打印计算结果，注意cublas的函数计算出数据是按列存储的 /*for(int i=0;i 
          
            #include 
           
             //#define IDX2C(i,j,ld) (((j)*(ld))+(i)) void print_matrix(const float *A, const int row_A, const int col_A) { for(int i = 0; i < row_A; i++){ for(int j = 0; j < col_A; j++){ printf("%.6f" ,A[i*col_A+j]); printf(" "); } printf("\n"); } printf("\n"); } #endif //从外部导入数据到vs的load函数 #ifndef LOAD_H_ #define LOAD_H_ #include 
            
              #include 
             
               void load(char* filename,float *x,int size){ std::ifstream indata; indata.open(filename); if( !indata ) { // file couldn't be opened std::cout<<"Error: file could not be opened"< 
              
                > x[i]; indata.close(); } #endif