cublasSgemm实现矩阵的相乘

今天想调用cuda的库函数实现矩阵的相乘，但是发现在cublasSgemm中矩阵是按照列主元素进行存储的，也就是一列列的存储的。这和c中一般按照行存储完全相反，后面看了一个帖子http://cudazone.nvidia.cn/forum/forum.php?mod=viewthread&tid=6001&extra=&page=2讲了一个方法，后面理了一下，具体分析如下：

比如，我们想求C=A*B这个矩阵运算，其中A={{1,1},{2,2},{3,3}};B={{1},{1}};C={{4},{5},{6}}，而对于A、B、C进行一维数组表示有A={1,1,2,2,3,3}，B={1,1}，C={4,5,6}；这个在c/c++是和前面表示一样，但是在cublasSgemm中就完全不对了，那么这个一维的A其实表示的是{{1,2,3},{1,2,3}};可以看到两个矩阵其实刚好是转置关系。那么我们要求C=A*B,按照一维数据输入的话结果A表示的是AT，B表示的是BT，所以我们要输入的是AT和BT，这样在公式中得到的才是A*B.假设这是得到矩阵C=A*B，但是C也是按列存储的，我们要的是CT,而CT=BT*AT，而这里的BT其实就是原矩阵B，AT其实就是原矩阵A。可见，我们通过交换AB的顺序就可以得到按行存储的C。

这里还有一点，就是cublasSgemm的参数，把自己绕的有点晕。这里我们输入的是B*A即（2*1）*（3*2），但是真正在函数中执行的是BT*AT=CT(1*2)*(2*3=(1*3)).因此主要不要把参数搞错了。。

而代码如下：

#include <cublas_v2.h> //cuda自带库函数
#include <helper_cuda.h>
#include <stdio.h>
int main(void)
{
float alpha=1.0;
float beta=0.0;
float h_A[6]={1,1,2,2,3,3};
float h_B[2]={1,1};
float h_C[3];
float *d_a,*d_b,*d_c;
checkCudaErrors(cudaMalloc((void**)&d_a,6*sizeof(float)));
checkCudaErrors(cudaMalloc((void**)&d_b,2*sizeof(float)));
checkCudaErrors(cudaMalloc((void**)&d_c,3*sizeof(float)));
checkCudaErrors(cudaMemcpy(d_a,&h_A,6*sizeof(float),cudaMemcpyHostToDevice));
checkCudaErrors(cudaMemcpy(d_b,&h_B,2*sizeof(float),cudaMemcpyHostToDevice));
checkCudaErrors(cudaMemset(d_c,0,3*sizeof(float)));
cublasHandle_t handle;
cublasCreate(&handle);
cublasSgemm(handle,CUBLAS_OP_N,CUBLAS_OP_N,1,3,2,&alpha,d_b,1,d_a,2,&beta,d_c,1);
checkCudaErrors(cudaMemcpy(h_C,d_c,3*sizeof(float),cudaMemcpyDeviceToHost));
for(int i=0;i<3;i++)
{
printf("%f\n",h_C[i]);
}
printf("\n");
return 0;
}