cuBLAS简介在科学计算和数值分析领域，经常需要解决矩阵特征值、线性方程（代数）等问题，形成了如EISPACK、LINPACK、LAPACK数学求解库。这些库解决的是偏上层的数学问题，它们共同依赖了更加底层的、更基础的数学库BLAS（Basic Linear Algebra Subprograms）。BLAS总共分为三个层次：
- 第一层是向量操作，如 𝑦=𝛼𝑥+𝑦 计算问题；
- 第二层向量-矩阵操作，如 𝑦=𝛼𝐴𝑥+𝛽𝑦 计算问题；
- 第三层矩阵-矩阵操作，如 𝐶=𝛼𝐴𝐵+𝛽𝐶 。
𝐶=𝛼𝐴𝐵+𝛽𝐶 就是我们常说的GEMM（GEneral Matrix Multiplication）问题，它也是BLAS第三层次要求解的问题。其中的广义来自于 𝛼 𝛽可以是任意数值，𝐴和𝐵可以是转置的也可以是非转置的，数值精度方面可以是单精度浮点数、双精度浮点数等。BLAS早期实现的时候采用的是数值计算领域的Fortran语言，在Fortran语言中二维数组是列主序的（列优先），所以在BLAS的API中针对AB矩阵的转置情况描述而言，列优先称为Normal（简写为N），行优先则称为Transpose（简写为T），值得注意的是转置与否是针对A B矩阵而言的， 𝐶 𝐷 矩阵必须是列优先的。同时为了区分不同的数值精度，BLAS规定了不同的简写，如单精度浮点为s（缩写自single precision）, 双精度浮点为d （缩写自double precision），单精度复数为c （缩写自complex），双精度浮点复数记为z。所以在附加上数据类型（如单精度浮点）之后函数名字就由gemm变为了sgemm（single precision general matrix multiplation）。同时BLAS约定了矩阵的维度，表示为
D _ {m,n} = \alpha A _ {m,k}B _ {k,n} + \beta C _ {m,n}其中的m，n表示行列数目，k为reduce轴上的数据维度。
深度学习中也有矩阵乘法所表达的网络层（如Pytorch中Linear层、Tensorflow中的Dense层），但是神经网络中该层的公式可以表达为 𝑦=𝑥𝐴+𝑏，其中x为input，维度一般为（N, Cin）y为output维度为（N, Cout），b为bias（一般是boardcast语义，此处我们忽略掉）。我们在忽略掉bias同时采用BLAS的矩阵形式描述，则其公式可以变更为 𝐷=𝐴𝐵 ，我们可以发现深度学习中需要的也是矩阵乘法的语义，和BLAS类似但是又不完全相同，不同点一深度学习中的公式没有BLAS公式中的 𝛼, 不同点二深度学习中公式没有BLAS中的加法右侧项，不同点三输出的D是行优先的而非BLAS中的列优先。除了公式层面的差异外，深度学习中的层为了计算效率会选用比传统的BLAS中的数据类型更低的精度如半精度和整数量化精度。
由于深度学习的输出要求的是行优先的，而BLAS中的约定中输出矩阵D是列优先的，因此当我们借助BLAS类函数完成深度学习的矩阵乘法时，需要对BLAS的公式进行转置变换，即
D^T=\alpha B^TA^T+\beta C^T ，同时设置 𝛼 为1， 𝛽 为0。 ——节选自reed大佬的文章
如上面介绍的，cuBLAS gemm中使用的矩阵是列优先，而C/C++中的矩阵默认为行优先，这导致我们在C/C++中使用cuBLAS时需要用一些取巧的方法完成矩阵计算，我在初次学习使用cuBLAS gemm时掌握了一个常用场景（A矩阵为mxk，B矩阵为kxn）的使用方式后便没继续深究。但最近我发现 A矩阵为mxk，B矩阵为nxk，即B矩阵需要转置后才能参与运算的场景也非常常见，如A \times A^T 一般不会特意存一个A^T，这次打算彻底弄清楚cuBLAS gemm中关于矩阵转置的问题。
1.cuBLAS环境
cuBLAS库中一共包含四个API
- The cuBLAS API, which is simply called cuBLAS API in this document (starting with CUDA 6.0)
- The cuBLASXt API (starting with CUDA 6.0)
- The cuBLASLt API (starting with CUDA 10.1)
- The cuBLASDx API (not shipped with the CUDA Toolkit)
其中前三个都在CUDA Toolkit中，我们所说的cuBLAS gemma在cuBLAS API中，在使用CUDA Toolkit的环境中只需引入#include <cublas_v2.h>头文件并链接上 cublas.so/dll/dylib 动态链接库，即加上-lcublas编译参数
2.cuBLAS gemm参数基本介绍
cuBLAS gemm提供了5种函数，差别在于支持的数据类型
- cublasSgemm 支持float类型（Single precision)
- cublasDgemm 支持double类型
- cublasCgemm 支持复数类型（Complex）
- cublasZgemm 支持DoubleComplex类型
- cublasHgemm 支持半精读（Half precision)
cublasStatus_t cublasSgemm(cublasHandle_t handle,
cublasOperation_t transa, cublasOperation_t transb,
int m, int n, int k,
const float *alpha,
const float *A, int lda,
const float *B, int ldb,
const float *beta,
float *C, int ldc)
cublasStatus_t cublasDgemm(cublasHandle_t handle,
cublasOperation_t transa, cublasOperation_t transb,
int m, int n, int k,
const double *alpha,
const double *A, int lda,
const double *B, int ldb,
const double *beta,
double *C, int ldc)
cublasStatus_t cublasCgemm(cublasHandle_t handle,
cublasOperation_t transa, cublasOperation_t transb,
int m, int n, int k,
const cuComplex *alpha,
const cuComplex *A, int lda,
const cuComplex *B, int ldb,
const cuComplex *beta,
cuComplex *C, int ldc)
cublasStatus_t cublasZgemm(cublasHandle_t handle,
cublasOperation_t transa, cublasOperation_t transb,
int m, int n, int k,
const cuDoubleComplex *alpha,
const cuDoubleComplex *A, int lda,
const cuDoubleComplex *B, int ldb,
const cuDoubleComplex *beta,
cuDoubleComplex *C, int ldc)
cublasStatus_t cublasHgemm(cublasHandle_t handle,
cublasOperation_t transa, cublasOperation_t transb,
int m, int n, int k,
const __half *alpha,
const __half *A, int lda,
const __half *B, int ldb,
const __half *beta,
__half *C, int ldc)上面的gemm接口都可以抽象为 C _ {m,n} = \alpha OP(A _ {m,k}) OP(B _ {k,n}) + \beta C _ {m,n} ，在矩阵乘法场景中一般会设置\alpha=1 ,\beta=0，因此可以简化为C _ {m,n} = OP(A _ {m,k}) OP(B _ {k,n})
参数列表如下
| Param. | Memory | In/out | Meaning |
|---|---|---|---|
| handle | input | handle to the cuBLAS library context. | |
| transa | input | operation op(A)，是否对A矩阵执行转置操作，可选参数为 CUBLAS_OP_N/T/C |
|
| transb | input | operation op(B)，是否对B矩阵执行转置操作，可选参数为 CUBLAS_OP_N/T/C |
|
| m | input | op(A)和C的行数 | |
| n | input | op(A)和C的列数 | |
| k | input | op(A)的列数和op(B)的行数 |
|
| alpha | host or device | input | 即D _ {m,n} = \alpha A _ {m,k}B _ {k,n} + \beta C _ {m,n} 中的\alpha，标量，可以在host也可以在device上 |
| A | device | input | 矩阵A，二维 |
| lda | input | A矩阵的leading dimension | |
| B | device | input | 矩阵B，二维 |
| ldb | input | B矩阵的leading dimension | |
| beta | host or device | input | 即D _ {m,n} = \alpha A _ {m,k}B _ {k,n} + \beta C _ {m,n} 中的\beta，标量，可以在host也可以在device上 |
| C | device | in/out | 矩阵C，二维 |
| ldc | input | C矩阵的leading dimension |
下面拆解几个可能令人困惑的点
-
handle在新版cuBLAS API（cublas_v2.h）中是一个重要概念，用来做上下文管理，如在多GPU场景可以将一个handle与一个设备绑定，每个handle可以互不干扰的执行任务，具体可以参考 New and Legacy cuBLAS API)👈
-
transa/b 表示是否对矩阵做转置操作，由于cuBLAS是列优先CUBLAS_OP_N（Normal)，即为列优先，CUBLAS_OP_T（Transpose）为行优先。
-
lda/b/c即为leading dimension，在官方文档的Data Layout中，我们可以看到对ld的使用方式，
#define IDX2C(i,j,ld) (((j)*(ld))+(i))，其中i,j分别为行列索引，从这也能看出是列主序（行主序应该是 i * ld + j），用来关联一维的物理存储空间和二维的逻辑存储空间，这是本篇的重点。
3.行优先 & 列优先
在深度学习的常用的场景中，我们希望使用 A(m x k) * B (k x n) = C (m x n) 的矩阵乘法。然而由于cuBLAS gemm默认支持列优先，很多人都会推荐这种方案
- 改写成 B.T(n x k) * A.T(k x m) = C.T(n x m)
- cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, n, m, k, &alpha, B, n, A, k, &beta, C, n)
这种方案非常的nice，大部分场景都够用了，但是会让人有许多疑问：
- 为什么我们传进去的参数仍然是A，B 而不是他们的转置？
- 为什么ld需要那样传参？
- 为什么出来的C可以直接在行优先的C++中使用？
这些问题的来源其实都是对行优先和列优先没概念我们假设一个C++中的矩阵A为
A =\begin{bmatrix} 1& 2 & 3\\ 4& 5 & 6\end{bmatrix} ，由于C++是行优先的，所以他的物理存储为 A = \begin{bmatrix} 1 & 2 & 3 & 4 & 5 & 6\end{bmatrix}
而如果A是列优先存储的，他的物理存储为A = \begin{bmatrix} 1 & 4 & 2 & 5 & 3 & 6\end{bmatrix}，这样大概能对行列优先存储有个概念。
接下来我们看看，一个行优先存储的A被列优先读取会怎么样
A =\begin{bmatrix} 1& 2 & 3\\ 4& 5 & 6\end{bmatrix} ，行优先存储为 A = \begin{bmatrix} 1 & 2 & 3 & 4 & 5 & 6\end{bmatrix}
由于矩阵的行数列数这种属性和存储方式是无关的，仅仅改变读取方式为列优先的A_col矩阵也和A具有相同的形状，那么A_col为
A\_col =\begin{bmatrix} 1& 3 & 5\\ 2& 4 & 6\end{bmatrix} , 注意A和A_col使用相同的物理空间，只是逻辑空间使用的映射方式不一样
A[i,j] = A[i * 3 + j] （1）
A_col[i,j] = A[j * 2 + i] （2）
这时我们可以引入leading dimension的概念，其中1式中的3和2式中的2，即为leading dimension，描述连续行（行优先）或列（列优先）之间的元素间距
leading dimension在cublas文档中的定义为：
describes the stride in elements between consecutive rows (for row major layout) or columns (for column major layout)
如果我们把A_col的leading dimension改为3，即
A_col[i,j] = A[j * 3 + i]，那么 A\_col =\begin{bmatrix} 1& 4 \\ 2& 5 \\ 3& 6\end{bmatrix} 刚好为A的转置，
现在就能回答上面的1，2问题了，我们传入B矩阵，并设置ld为k，被列优先读取后会变为B的转置。
leading dimension怎么选才能有这个效果？只需要保证ld的值和行优先时的leading dimension一样（A在行优先的情况下会在连续存储3个行元素后逻辑换行）
关于m,n,k三个参数的设置，需要考虑OP操作，B.T(n x k) * A.T(k x m) = C.T(n x m) 使用cuBLAS gemm可以写为
OP\_N(B)_{n\times k} * OP\_N(A)_{k\times m} = C^T_{n\times m} , m为乘法左项的行数，n为乘法右项的列数，k为左项的列数或右项的行数
所以，m = n, n = m, k = k
得出的矩阵C必然为列优先，C默认会以行数作为leading dimension（C会在连续存储n个数后换列），即ld=n ，而列优先的C^T_{n \times m} 被行优先方式读取且leading dimension为n(即行数为n)，相当于 C_{m\times n} = (C^T_{n \times m})^T_{m \times n}
4.整点复杂的
1.把A,B使用列优先存储后传入cublas gemm | 等同于把列优先存储的A,B传入cublas gemm
A_col = row2col(A)
B_col = row2col(B)相当于列优先存储的A，B被列优先读取，此时A，B都是正常的（未转置），可以完成正常的 AB = C
接口调用可以抽象为 OP\_N(row2col(A_{m\times k}))_{m\times k} \times OP\_N(row2col(B_{k\times n}))_{k\times n} = C_{m\times n}
接口调用为：
cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A_col, m, B_col, k, &beta, C, m);不过可惜的时，这时的C不能直接使用，注意此时C的leading dimension为m（我们需要的C(mxn)的ld应该为n），这意味着不可以直接使用行优先的C(mxn)读取并转置，因为内部元素在逻辑上是“错位的”，我们可以在行优先使用ld=m即形状为C(nxm)的矩阵读取然后转置，或者在物理存储中把列优先转为行优先存储（如使用下面的col2row函数），再用行优先的C(mxn)读取。
// convert row-major matrix to column-major matrix
void row2col(const float *a, int rows, int cols, float *b)
{
int size = rows * cols;
for (int i = 0; i < size; i++)
{
int row = i / cols;
int col = i % cols;
b[col * rows + row] = a[i];
}
}
// convert column-major matrix to row-major matrix
void col2row(const float *a, int rows, int cols, float *b)
{
int size = rows * cols;
for (int i = 0; i < size; i++)
{
int col = i / rows;
int row = i % rows;
b[row * cols + col] = a[i];
}
}2.使用CUBLAS_OP_T按行优先的方式读取A,B
行优先存储的A，B被行优先方式读取，所以A，B都是正常的，可以完成正常的 AB = C
接口可以抽象为 OP\_T(A_{m\times k})_{m\times k} \times OP\_T(B_{k\times n})_{k\times n} = C_{m\times n}
接口调用为：
cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_T, m, n, k, &alpha, A, k, B, n, &beta, C, m);返回结果同1，这时的C不能直接使用，注意此时C的leading dimension为m（我们需要的C(mxn)的ld应该为n），这意味着不可以直接使用行优先的C(mxn)读取并转置，因为内部元素在逻辑上是“错位的”，我们可以在行优先使用ld=m即形状为C(nxm)的矩阵读取然后转置，或者在物理存储中把列优先转为行优先存储，再用行优先的C(mxn)读取。
3.使用未转置的B矩阵，列优先读取B，AxB.T = C
此时的B矩阵为nxk，需要转置为kxn后参与计算，不过我们不需要单独花计算量去做转置操作，可以使用gemm的列优先读取特点在读取过程中实现，具体为，使用OP_N和ld参数来转置B矩阵
可以抽象为： OP\_T(A_{m\times k})_{m\times k} \times OP\_N(B_{n\times k})_{k\times n} = C_{m\times n}
接口调用为：
cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, k, B_nxk, k, &beta, C, m)返回结果同1，2，这时的C不能直接使用，注意此时C的leading dimension为m（我们需要的C(mxn)的ld应该为n），这意味着不可以直接使用行优先的C(mxn)读取并转置，因为内部元素在逻辑上是“错位的”，我们可以在行优先使用ld=m即形状为C(nxm)的矩阵读取然后转置，或者在物理存储中把列优先转为行优先存储，再用行优先的C(mxn)读取。
4.使用未转置的B矩阵，行优先读取B, B x A.T = C.T
可以抽象为： OP\_T(B_{n\times k})_{n\times k} \times OP\_N(A_{m\times k})_{k\times m} = C^T_{n\times m}
接口调用为：
cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, n, m, k, &alpha, B_nxk, k, A, k, &beta, C, n);返回结构可以直接用行数为n的行优先方式读取。
代码地址：cudaLearn/cublas_gemma.cu at main · HeduAiDev/cudaLearn (github.com)
Comments | NOTHING