利用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