cuda 矩阵乘法函数之cublasSgemm

，可以考虑使用, 例如cublasSgeam()(矩阵加法), 进行一次1.0 * AT + 0.0 * B的参数设定, 利用内置的转置功能(注意这里的1和0), 来进行将A转换成AT. 在使用CUDA的cuBLAS库中矩阵乘法函数cublasSgemm时，注意到cuda其中的二维矩阵的储存是“按列储存”，一天都处于蒙蔽状态，查了很多资料，按所得结果情况，总结出如下几条。

一、获得按行存储的结果

由博文：http://blog.csdn.net/xfortius/article/details/9225799收到启发：

比如，我们想求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自带库函数
#inc
4000
lude <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;
}

总结：若想获得按行存储的结果矩阵C=A*B（可复制到host）,可以使用代码：

cublasSgemm(handle,CUBLAS_OP_N,CUBLAS_OP_N,m,n,k,B,ldb,A,lda,C,ldc);

说明：

输入的参数为CUBLAS_OP_N ， CUBLAS_OP_N

相比原参数，A和B的位置调换（C=A*B —>CT=(BT*AT) (其中CT按列存储，实际上就是按行存储的C ！)）

其中参数m为col of B ,n 为row of A , k为A和B相同的(row of B or col of A) ,参数ldb , lda , ldc 分别为col of B,A,C ;

经过验证，结果正确。

二、获得按列存储结果（存储在GPU显存中）

有些矩阵为中间结果，为了避免重复使用上述方法，可以使用下面的方法，其中，GPU显存中的为按列存储，Host内存中的为按行存储。

灵感来源博文：http://www.xuebuyuan.com/1966642.html

由于cublas为了更大的适应Fortan语言，二维数据的存储采用以列优先的方式，这与C/C++中，行优先的存储方式不同。由于本人的研究是数据的来源是C代码得到的，为了加速矩阵的运算效率，利用cublas来完成。本文档提出了一种有效的解决方案。

为了更好的说明，以函数cublasSgemm的实现C= A*B为例。接口cublasSgemm 实现的功能为C = alpha*A*B + beta*C，为了完成C= A*B 的功能，令alpha= 1.0f，beta = 0.0f

利用cublasSgemm的参数 transa（transb） 和 lda（ldb）的设置来共同解决存储方式改变的问题。

第一种：输入矩阵A,B均在CPU上

cublasSgemm(‘t’,’t’,结果矩阵C的行,结果矩阵C的列，左矩阵A的列, alpha，左矩阵,左矩阵的列,右矩阵,右矩阵的列, beta,结果矩阵,结果矩阵的行);

如果求AT*B，只需将A的t改为n（t的转置为t）,并且关于A的参数（行和列）是转置以前的。

第二种：输入矩阵A在显存上以column-major的方式存储，B在CPU上

cublasSgemm(‘n’,’t’,结果矩阵C的行,结果矩阵C的列，左矩阵A的列, alpha，左矩阵A,左矩阵A的行,右矩阵B,右矩阵B的列, beta,结果矩阵C,结果矩阵C的行);

第三种：输入矩阵B在显存上以column-major的方式存储，A在CPU上

cublasSgemm(‘t’,’n’,结果矩阵C的行,结果矩阵C的列，左矩阵A的列, alpha，左矩阵A,左矩阵A的列,右矩阵B,右矩阵B的行, beta,结果矩阵C,结果矩阵C的行);

第四种：输入矩阵A，B均在显存上，以column-major的方式存储

cublasSgemm(‘n’,’n’,结果矩阵C的行,结果矩阵C的列，左矩阵A的列, alpha，左矩阵A,左矩阵A的行,右矩阵B,右矩阵B的行, beta,结果矩阵C,结果矩阵C的行);

得到的结果C都是按列存储的（拷贝出host内存之前转置一下，可以考虑使用, 例如cublasSgeam()(矩阵加法), 进行一次1.0 * AT + 0.0 * B的参数设定, 利用内置的转置功能(注意这里的1和0), 来进行将A转换成AT.）。

经过验证，结果正确。

总结：

如果前边的参数是’t’，那么leading dimesion 就是矩阵的列数，因为此时的矩阵是按照C语言以行优先的方式来存储的；反之如果前边的参数是’n’，那么leading dimesion 就是矩阵的行数，此时的矩阵保持CUBLAS的列优先存储方式。