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

 

 

 

 

发布了16 篇原创文章 · 获赞 5 · 访问量 361
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 大白 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览