cuda编程笔记（11）--学习cuBLAS使用
cuBLAS 是 NVIDIA 提供的 GPU 加速 BLAS 库；使用时需要#include <cublas_v2.h>
如果使用VS，需要添加cublas.lib的链接；如果用命令编译，-l记得加上cublas
cuBLAS 的核心基础概念
cublasStatus_t
-
类型：枚举类型
-
作用：表示 cuBLAS API 的返回状态（错误码）。
-
常用值：
常量 含义 CUBLAS_STATUS_SUCCESS成功 CUBLAS_STATUS_NOT_INITIALIZEDcuBLAS 库未初始化 CUBLAS_STATUS_ALLOC_FAILEDGPU 设备内存分配失败 CUBLAS_STATUS_INVALID_VALUE传入参数无效 CUBLAS_STATUS_ARCH_MISMATCH硬件不支持请求的功能（如 Tensor Core） CUBLAS_STATUS_EXECUTION_FAILED核函数执行失败
返回值检查（典型模式）：
cublasStatus_t status = cublasCreate(&handle);
if (status != CUBLAS_STATUS_SUCCESS) {
printf("cuBLAS initialization failed!\n");
}
cublasHandle_t
-
类型：指向 cuBLAS 库上下文的句柄（类似于会话）。
-
作用：
-
所有 cuBLAS 函数都需要它。
-
cuBLAS 使用 上下文模型，通过
handle管理状态。
-
-
生命周期：
-
通过
cublasCreate()创建。 -
用完后调用
cublasDestroy()销毁。
-
cublasOperation_t
-
类型：枚举类型
-
作用：指定矩阵是否需要 转置。
-
值：
常量 含义 CUBLAS_OP_N不转置（Normal） CUBLAS_OP_T转置（Transpose） CUBLAS_OP_C共轭转置（Conjugate）
使用场景：cublasSgemm() 等矩阵乘法接口。
例如：
cublasOperation_t transA = CUBLAS_OP_N;
cublasOperation_t transB = CUBLAS_OP_T; // B 矩阵转置
cublasCreate()
cublasStatus_t cublasCreate(cublasHandle_t *handle);
-
作用：
-
初始化 cuBLAS 库。
-
创建 cuBLAS 句柄。
-
-
参数：
-
handle：指向cublasHandle_t的指针，返回创建的句柄。
-
-
返回值：
-
CUBLAS_STATUS_SUCCESS表示成功。
-
cublasDestroy()
cublasStatus_t cublasDestroy(cublasHandle_t handle);
-
作用：
-
释放
cublasHandle_t占用的资源。
-
-
参数：
-
handle：之前用cublasCreate()创建的句柄。
-
-
返回值：
-
成功返回
CUBLAS_STATUS_SUCCESS。
-
矩阵乘法（GEMM）
函数原型
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
);
参数说明
-
transa,transb:-
CUBLAS_OP_N: 不转置 -
CUBLAS_OP_T: 转置
-
-
矩阵维度：
-
m × k矩阵 A -
k × n矩阵 B -
m × n矩阵 C
-
-
alpha,beta: 标量，计算公式：
C = alpha * op(A) * op(B) + beta * C
lda, ldb, ldc: leading dimension（矩阵为列主序时填行数）
列主序
这个ld是理解这个接口最关键的难点。
我们在调用cublas系列的矩阵乘法接口的时候，是按照列主序读取我们的矩阵的。
比如我们的M,N,K是2,4,3.那么矩阵A的维度是m*k即2行3列
假如我们输入是行优先的矩阵A:
[1,2,3]
[4,5,6]
由于我们传递的参数的连续的内存，所以按照我们行优先的分布，这个A被展平为
[1,2,3,4,5,6]
到这里都没错，但是接下来lda参数的设置就很关键了，如果这时候你还是填行数m=2，那么此时的A矩阵将会被这样读取
[1,3,5]
[2,4,6]
也就是每次读取两个数，按照列优先的方式去填充。这和我们想要的完全不一样。
那为什么说列主序矩阵填行数就对了？
假如我们输入的是正确的按列主序的矩阵A
[1,4,2,5,3,6]
这时候lda填行数m=2，就会读取成
[1,2,3]
[4,5,6]
这才是我们想要计算的形式。也就是说假如你想要一个矩阵按照行优先的形式进行计算，那你传进来的参数需要按照列优先的内存排布，然后lda填行优先时的行数。
当然这样的输入矩阵内存分布方式，完全不符合人类的数学习惯，所以一般输入还是会按照行优先，那么对应的函数调用就会有所调整，详见下面的讲解。
最小示例：矩阵乘法 C = A × B
#include <cublas_v2.h>
#include <cuda_runtime.h>
#include <iostream>
#include <vector>
#define M 4
#define N 4
#define K 4
int main() {
// 创建 cuBLAS handle
cublasHandle_t handle;
cublasCreate(&handle);
// Host 矩阵
std::vector<float> h_A(M*K), h_B(K*N), h_C(M*N);
for (int i = 0; i < M*K; i++) h_A[i] = 1.0f;
for (int i = 0; i < K*N; i++) h_B[i] = 2.0f;
float *d_A, *d_B, *d_C;
cudaMalloc(&d_A, M*K*sizeof(float));
cudaMalloc(&d_B, K*N*sizeof(float));
cudaMalloc(&d_C, M*N*sizeof(float));
cudaMemcpy(d_A, h_A.data(), M*K*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B.data(), K*N*sizeof(float), cudaMemcpyHostToDevice);
const float alpha = 1.0f;
const float beta = 0.0f;
// 调用 cuBLAS GEMM
cublasSgemm(handle,
CUBLAS_OP_N, CUBLAS_OP_N,
M, N, K,
&alpha,
d_A, K,
d_B, N,
&beta,
d_C, N);
cudaMemcpy(h_C.data(), d_C, M*N*sizeof(float), cudaMemcpyDeviceToHost);
// 打印结果
for (int i = 0; i < M; i++) {
for (int j = 0; j < N; j++) {
std::cout << h_C[i + j*M] << " "; // 注意列主序
}
std::cout << "\n";
}
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
cublasDestroy(handle);
}
如何保持行优先，正确调用cublasgemm？
如果想要保持行优先，且不想动原矩阵，需要这么调整
-
m × k矩阵 A -
k × n矩阵 B -
m × n矩阵 C
我想要得到C=A*B,且A，B，C都是行优先表示的
需要如此调用（注意A，B的顺序换了）
cublasSgemm(handle,
CUBLAS_OP_N, CUBLAS_OP_N, // 不转置
N, M, K,
&alpha,
d_B, N, // lda = N = B 的列数
d_A, K, // ldb = K = A 的列数
&beta,
d_C, N); // ldc = N = C 的列数
原理是什么呢，这样子调用相当于执行了：,所以d_C里存储的是C的转置，但是d_C现在是按列优先解释的，由于存的时候还是存成了一维，所以先存第一列，再存第二列...；我们只要按照行优先的方式去访问d_C,就相当于正常访问了正常结果的
矩阵转置or矩阵加法
cuBLAS 提供了 cublasSgeam：
cublasSgeam 是 cuBLAS 提供的矩阵加法/转置操作接口，可以完成：
其中 op(X) 表示矩阵 X 是否转置。
cublasStatus_t cublasSgeam(
cublasHandle_t handle, // cuBLAS 句柄
cublasOperation_t transa, // A 是否转置 (CUBLAS_OP_N 或 CUBLAS_OP_T)
cublasOperation_t transb, // B 是否转置
int m, // C 的行数
int n, // C 的列数
const float *alpha, // 缩放系数 α
const float *A, // 矩阵 A
int lda, // A 的 leading dimension (步长)
const float *beta, // 缩放系数 β
const float *B, // 矩阵 B
int ldb, // B 的 leading dimension
float *C, // 结果矩阵 C
int ldc // C 的 leading dimension
);
-
handle：cuBLAS 上下文。
-
transa / transb：
-
CUBLAS_OP_N→ 不转置 -
CUBLAS_OP_T→ 转置
-
-
m, n：结果矩阵 C 的维度 (m 行 n 列)【注意这是在cublas列优先存储视角下的,可以看下面的例子】。
-
alpha, beta：系数，通常
alpha = 1.0f，beta = 1.0f。 -
A, lda：矩阵 A 的指针和 leading dimension。
-
B, ldb：矩阵 B 的指针和 leading dimension。
-
C, ldc：结果矩阵 C 的指针和 leading dimension。
⚠ leading dimension (lda, ldb, ldc)：
cuBLAS 默认使用 列主存储（和 Fortran 一致），所以 lda = A 的行数，即 A 每列元素的间隔。
-
可以实现 矩阵转置（单独把 beta 设为 0）。
-
可以实现 矩阵加法。
-
可以实现 带缩放的线性组合。
示例代码（所有矩阵都是行优先）
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <iostream>
int main() {
int m = 2, n = 3;
float alpha = 1.0f, beta = 1.0f;
// Host 数据
float h_A[6] = {1, 2, 3, 4, 5, 6}; // 2x3
float h_B[6] = {10, 20, 30, 40, 50, 60};
float h_C[6];
// Device 指针
float *d_A, *d_B, *d_C;
cudaMalloc((void**)&d_A, 6 * sizeof(float));
cudaMalloc((void**)&d_B, 6 * sizeof(float));
cudaMalloc((void**)&d_C, 6 * sizeof(float));
cudaMemcpy(d_A, h_A, 6 * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, 6 * sizeof(float), cudaMemcpyHostToDevice);
cublasHandle_t handle;
cublasCreate(&handle);
cublasSgeam(handle,
CUBLAS_OP_T, CUBLAS_OP_T, // 对 A、B 都做转置
n, m, // 注意，行列要交换
&alpha,
d_A, m, // lda = m
&beta,
d_B, m,
d_C, n); // ldc = n
cudaMemcpy(h_C, d_C, 6 * sizeof(float), cudaMemcpyDeviceToHost);
std::cout << "Result C:\n";
for (int i = 0; i < m; i++) {
for(int j=0;j<n;j++){
std::cout << h_C[i*m+j] << " ";
}
std::cout<<std::endl;
}
cublasDestroy(handle);
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
return 0;
}
我们可以这么去思考，就是我们函数最终得到的矩阵，得是行优先视角下的A+B=C的转置，然后因为C是列优先存储，我们再通过行优先去访问，结果会是正确的。
明白了这个，我们的A，B参数的设置就可以理解了，就是要调整一下，让结果变成C的转置
转置的注意点
比如代码里，A设置了转置；但是 lda 定义的是 A 在内存中的存储布局，而不是 Aᵀ 的逻辑维度。
对于A来说，我们设置它的一列（cublas默认列优先）有m个元素，转置之后，它的一列有n个元素（注意cublas列优先的原则）；于是参与计算的Aᵀ形状如下
[1,4
2,5
3,6]
同理，Bᵀ形状如下
[10,40
20,50
30,60]
于是结果存到d_C(n行m列)就是
[11,44
22,55
33,66]
它确实是原本C的转置；但是存储到d_C也是按照列优先存储的，所以d_C最终是[11,22,33,44,55,66]。我们按照行优先去访问，就能得到正确的结果：
[11,22,33]
[44,55,66]
当然，也可以不用主动设置转置，因为加法操作的转置有分配律，可以直接倒转维度
cublasSgeam(handle,
CUBLAS_OP_N, CUBLAS_OP_N, //
n, m, // 注意，行列要交换
&alpha,
d_A, n,
&beta,
d_B, n,
d_C, n);
这样计算得到的d_C,按照m行n列行优先访问也是没问题的
但是对于上面介绍的乘法操作，就没有这么随便了，因为乘法没有交换律，也没有转置的分配律，需要仔细考虑维度
向量乘矩阵
-
op(A)表示 A 或 Aᵀ（取决于trans参数） -
x和y是向量 -
alpha、beta是缩放因子
cublasStatus_t cublasSgemv(
cublasHandle_t handle,
cublasOperation_t trans, // 是否转置 A
int m, int n, // 矩阵 A 的维度 (m x n)
const float *alpha, // 缩放因子 alpha
const float *A, int lda, // 矩阵 A
const float *x, int incx, // 向量 x
const float *beta, // 缩放因子 beta
float *y, int incy // 向量 y
);
-
trans：-
CUBLAS_OP_N→ A 不转置，维度 m×n -
CUBLAS_OP_T→ A 转置，维度 n×m
-
-
lda：leading dimension，通常等于m（A 的行数） -
incx、incy：向量步长，通常 = 1
示例代码：
#include <cublas_v2.h>
#include <cuda_runtime.h>
#include <iostream>
int main() {
cublasHandle_t handle;
cublasCreate(&handle);
const int m = 3, n = 2;
float h_A[m*n] = { 1, 2, 3, 4, 5, 6}; // 3x2 矩阵
float h_x[n] = { 1, 1}; // 向量
float h_y[m] = { 0, 0, 0 };
float *d_A, *d_x, *d_y;
cudaMalloc(&d_A, m*n*sizeof(float));
cudaMalloc(&d_x, n*sizeof(float));
cudaMalloc(&d_y, m*sizeof(float));
cudaMemcpy(d_A, h_A, m*n*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_x, h_x, n*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_y, h_y, m*sizeof(float), cudaMemcpyHostToDevice);
float alpha = 1.0f;
float beta = 0.0f;
// y = alpha * A * x + beta * y
cublasSgemv(handle, CUBLAS_OP_T, n, m,//n*m是转置之前的维度，转置后是正常的m*n
&alpha, d_A, n, d_x, 1, &beta, d_y, 1);//这里lda需要填n
cudaMemcpy(h_y, d_y, m*sizeof(float), cudaMemcpyDeviceToHost);
std::cout << "Result y: ";
for (int i = 0; i < m; i++) std::cout << h_y[i] << " ";
std::cout << std::endl;
cudaFree(d_A); cudaFree(d_x); cudaFree(d_y);
cublasDestroy(handle);
return 0;
}
矩阵 × 向量（gemv）其实可以用 矩阵 × 矩阵（gemm）来替代，但有几点需要注意：
-
性能原因：
gemv是 BLAS Level 2 操作，专门针对 矩阵-向量做了优化，尤其在小规模问题上，它的内存访问模式和寄存器利用更高效。 -
内存开销：
用gemm处理 m×n 矩阵 × n×1 向量，cuBLAS 内部仍按矩阵矩阵方式调度，会引入不必要的开销。 -
语义清晰：
gemvAPI 直接说明这是矩阵-向量操作，避免错误。
总结：
-
gemm是通用 API，可以替代gemv。 -
但如果只是单次
矩阵 × 向量，gemv更高效。 -
深度学习场景，通常是批处理（
GEMMdominate），所以主流框架用gemm。
多流并发
cublasHandle_t 绑定 CUDA 流：
-
默认情况下，cuBLAS 的操作运行在 默认流（stream 0） 上，这意味着 cuBLAS 操作会和其他默认流的操作 顺序执行，无法并发。
-
绑定流后（通过
cublasSetStream(handle, stream)）：-
cuBLAS 的运算将提交到指定的 CUDA 流
-
可以和其他流中的 kernel 并发执行，实现任务流水线和异步调度
-
对多 GPU 或 pipeline 加速非常重要（比如计算和数据拷贝重叠）
-
cuBLAS系列函数的本质：
-
cuBLAS API（如
cublasSgemm）是主机端函数，它的作用是：-
通过
handle向 CUDA 驱动提交任务 -
把矩阵乘法封装成 GPU 内核，并异步调度到 GPU 上运行
-
-
计算实际发生在 GPU 上，不是 CPU
-
调用
cublasSgemm后，不会阻塞主机（异步执行），除非你显式调用cudaDeviceSynchronize或访问 GPU 结果
核心区别：
-
cublasSgemm不做计算本身，它只是 发起任务 -
CUDA 流 +
handle控制 任务在哪个流上执行 -
默认流会串行执行，多个流能并发
举个并发例子：
cublasHandle_t handle;
cublasCreate(&handle);
cudaStream_t stream1, stream2;
cudaStreamCreate(&stream1);
cudaStreamCreate(&stream2);
// 第一批数据在 stream1
cublasSetStream(handle, stream1);
cublasSgemm(handle, ... /* batch1 */);
// 第二批数据在 stream2
cublasSetStream(handle, stream2);
cublasSgemm(handle, ... /* batch2 */);
// 两个 GEMM 并行执行
cudaStreamSynchronize(stream1);
cudaStreamSynchronize(stream2);
如果学习过我的Boost.Asio文章的话，cublasHandle_t非常像Boost.Asio里的io_context
北京人形旗下天工造物具身智能开源社区，聚焦具身天工与慧思开物两大平台
更多推荐
11
0- 0
已为社区贡献17条内容
所有评论(0)