cublas 矩阵乘法

cublas是cuda用来解决线性代数的问题的一个函数库，而且对于矩阵运算来说，其效率比大部分人自己写核函数高不少，只是cublas不同于C++，是列优先存储，因此参数一不小心设的不对，结果大不相同，所以记录下来；

首先介绍一下什么是列优先存储？

C++中，二维数组也是通过一维数组传进去的，假设我们有1个2*2的数组，设为a=[1,2,3,4],那么C++变成2*2的数组，就为[1,2;3,4],而cublas就变成了[1,3;2,4],这就是列优先存储，先填充列，再填充行。

cublas中提供矩阵乘法的函数有4个：cublasSgemm（单精度实数）、cublasDgemm（双精度实数）、cublasCgemm（单精度复数）、cublasZgemm（双精度复数）。对应于不同的精度，用法基本上是一样的。注意在cuda中，GPU对单精度的效率快了双精度很多，我测试过有两倍不止，因此如果不需要很大的精度的话，尽量使用单精度计算。

cublas用于矩阵乘法大概分为如下几步：创建句柄，分配显存，计算及释放显存和销毁句柄等步骤，详情可以参考博客：

https://www.cnblogs.com/scut-fm/p/3756242.html

首先我们假设有个2*3和3*2的矩阵相乘，分别设为A,B；

A=[4 8 4;7 5 6] ,B=[7 6;6 3;5 6];

需要算C=A*B；

cublasSgemm 计算的计算公式：C=alpha*op(A)op(B)+beta*C

如果要求C=A*B,令alpha=1,beta=0即可；

下面主要对参数设置进行详解：

 如果计算A*B（注意顺序）
 CUBLASAPI cublasStatus_t CUBLASWINAPI cublasSgemm_v2
 (
      cublasHandle_t handle,//句柄
      cublasOperation_t transa,//对A是否转置
      cublasOperation_t transb,//对B是否转置
      int m,//A的行数
      int n, //B的列数
      int k, //A的列数
      const float *alpha,//标量 α 的指针，可以是主机指针或设备指针，只需要计算矩阵乘法时命 α = 1.0f
      const float *A, // 矩阵（数组）A 的指针，必须是设备指针
      int lda,//A的主维，即A的行数
      const float *B,// 矩阵（数组）B 的指针，必须是设备指针
      int ldb,//B的主维
      const float *beta,//标量 β 的指针,可以是主机指针或设备指针，只需要计算矩阵乘法时命 β = 0.0f
      float *C,//矩阵（数组）C 的指针，必须是设备指针
      int ldc //矩阵C的主维
);

因为cublas为列优先存储，当将A拷贝到显存后，如果按正常设置，cublas中矩阵A就变为了[4 4 5;8 7 6],B为

[7 3;6 5;6 6];这显然不是我们所要的结果；

那怎么办呢？下面是一种正确的用法：

cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_T, m,n,k, &alpha, d_A, k, d_B, n, &beta, d_C, m);

最终的结果为结果的转置；

来分析一下具体的原理：

这里主要利用了主维的概念，即矩阵的行数（因为按列存储），负责对blas中矩阵进行变形，可以理解为matlab中resize,但它又与resize有不同之处，它存在着下限且列数始终是固定的：

设A行为m,列为k,当参数 transa 为 CUBLAS_OP_N 时，矩阵 A 在 cuBLAS 中的尺寸实际为 lda × k，此时要求 lda ≥ max(1, m)（否则由该函数直接报错，输出全零的结果）；当参数 transa 为 CUBLAS_OP_T 或 CUBLAS_OP_C 时，矩阵 A 在 cuBLAS 中的尺寸实际为 lda × m，此时要求 lda ≥ max(1, k) 。

可以看出lda的下限至少要与矩阵大小一致，转置也是与转置后的大小一致，如果lda大于下限，矩阵将自动补0；

还是以A=[4 8 4;7 5 6]为例，这是一个2行3列矩阵；

将它拷贝到显存，放入cublas中，注意因为主维设为了3，所以cublas中d_A resize为一个3*2的矩阵（满足下限）,遵循列优先存储，d_A就变为了[4 7;8 5;4 6];然后转置，d_A就变为了[4 8 4;7 5 6];

同理,d_B也就变为了[7 6;6 3;5 6]；

这样算出来的结果就为正常的结果；

还是因为列优先原则，所以当从cublas拷贝出来时,结果为其转置;

当算矩阵乘法时（即主维刚好是下限），其实可以换个思路理解：当转置时，计算时cublas变成了行优先存储，主维就变成了矩阵列数，而C并没有转置，仍为列优先存储，主维为矩阵行数，所以最后拷贝出来时矩阵为结果的转置；

还有一种用法，利用了A*B=(BT*AT)T的概念，

设A=m*k;B=k*n;

cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, n, m, k, &alpha, *B, n, *A, k, &beta, *C, n);

因为cublas按列存储，因此如果我们把blas中d_A矩阵大小设为k*m,那么当从显存中拷贝数据到blas中后，blas中d_A就刚好为A的转置，同理d_B矩阵大小设为n*k,那么就能计算出来BT*AT;当从blas中按列拷贝出来后，结果就正好为A*B;

本文参考链接：https://www.cnblogs.com/cuancuancuanhao/p/7763256.html