cuBLAS矩阵运算库实战：从CUDA入门到性能调优

如果你写过一阵子CUDA，大概率会有这种体会：核函数写得越多，越觉得矩阵运算这块又关键又烦人。自己实现一个矩阵乘法，要处理线程块划分、共享内存、bank conflict、访存合并，性能还不一定拼得过库。后来我发现，其实N卡生态里早就把这事儿干完了——cuBLAS就是NVIDIA官方提供的高性能BLAS库，专门在CUDA平台上做线性代数运算。它的定位非常简单粗暴：你们别自己造轮子了，直接调我，我针对不同GPU架构做了深度优化，性能基本是当前硬件条件的上限。

这篇文章就围绕cuBLAS展开，从它是什么、能解决什么问题，到怎么在自己的C++工程里集成和调用，再到性能调优时需要注意的细节，一次性理清楚。如果你正准备在N卡GPU上做科学计算、机器学习推理、图像处理或者任何涉及大规模矩阵运算的项目，这篇内容应该能帮你少走不少弯路。

1. cuBLAS是什么，它到底解决了什么问题

1.1 为什么矩阵运算会成为性能瓶颈

几乎所有数值计算任务，最后都能归结到矩阵运算上——神经网络的前向传播就是一连串矩阵乘法加激活，有限元求解的核心是稀疏矩阵迭代，推荐系统的协同过滤要算相似度矩阵。这些计算有个共同特点：数据规模大、计算模式规整、重复性强。

在CPU上做这件事，受限于核心数量和内存带宽，规模一大就顶不住了。用GPU当然是对的路径，但真正的问题在于：GPU编程门槛不低。你不仅要理解线程层次、内存层次、同步机制，还得针对具体的矩阵尺寸、数据类型、GPU架构去调参。而且写出来的代码性能往往不稳定，换个GPU型号可能就废了。

1.2 cuBLAS在CUDA生态里的定位

cuBLAS是NVIDIA官方发布的BLAS（Basic Linear Algebra Subprograms）库的GPU版本。BLAS是线性代数领域的事实标准接口，分三级：一级做向量-向量运算，二级做矩阵-向量运算，三级做矩阵-矩阵运算。cuBLAS实现了这套标准接口，同时针对CUDA的硬件特性做了极致的优化。

说白了，cuBLAS就是NVIDIA替你把最苦最累的那部分活儿干完了。它的核心价值有几个方面：

• 性能极佳。cuBLAS的GEMM（通用矩阵乘法）内核在整个CUDA生态里是性能标杆，从Kepler架构到Hopper架构，每一代都有对应的优化版本。
• 接口标准。它遵循BLAS命名规范，学过数值计算的人几乎零成本上手。
• 持续更新。CUDA Toolkit里换了版本，cuBLAS会跟着适配新架构的新特性。

我用一个具体例子来说明这个差距。我自己写过一个简单的SGEMM内核，用共享内存和分块做了优化，在同样规模的1024x1024矩阵乘法上，测出来的性能大概只有cuBLAS的40%到60%。而且cuBLAS在新卡上的表现更夸张，因为NVIDIA会为最新架构单独调内核，靠自己的代码很难追得上这种差距。

1.3 什么人需要关注cuBLAS

如果你属于下面这几类人，cuBLAS就是绕不开的东西：

第一类是科学计算和工程仿真方向的开发者，需要做大规模稠密矩阵运算，比如求解线性方程组、特征值分解。第二类是机器学习框架的开发者和使用者，很多框架的底层卷积和全连接层，最终都是通过cuBLAS或它的变体来完成的。第三类是自己写CUDA程序但又不想重复造轮子的同学，想聚焦在自己的核心算法上，矩阵运算直接调库。

即使你打算自己写卷积或者特殊的矩阵运算，也应该先看一下cuBLAS的实现思路，比如它的分块策略、双缓冲、向量化优化技巧，对你的自研内核设计会有很大启发。

2. 工具准备与环境配置：让cuBLAS跑起来

2.1 安装CUDA Toolkit时需要注意的细节

首先把CUDA Toolkit装上。去官网下载对应你系统版本的安装包，Windows和Linux都有。安装过程本身不复杂，但在Linux上有几个容易踩的坑。

第一，驱动版本和CUDA版本要匹配。如果你用发行版自带的包管理器安装驱动，版本可能偏旧。我个人的习惯是先把驱动升级到NVIDIA官方最新版本，再装CUDA Toolkit，这样兼容性问题最少。

第二，环境变量要配好。以Linux为例，需要把CUDA的bin目录和lib目录加到PATH和LD_LIBRARY_PATH里。我见过不少同学装好了CUDA，但是一编译就说找不到库，十有八九是环境变量没配。

export PATH=/usr/local/cuda/bin:$PATH export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH

第三，检查安装是否成功。装完之后，跑一下nvcc --version确认编译器版本，再用nvidia-smi看一眼驱动状态。建议再跑一下CUDA自带的deviceQuery示例程序，确保显卡能被识别、CUDA上下文能创建出来。

2.2 确认你的GPU架构与算力

cuBLAS虽然能跑在各种N卡上，但不同架构的支持程度不一样。算力（Compute Capability）是一个关键指标，它决定了你能用哪些cuBLAS API和底层特性。

比如旧卡算力在3.x、5.x，只能用比较基础的API。如果你想用cuBLAS的batched GEMM或者使用Tensor Core的算法变体，就得确认你的卡支持对应的特性。可以在代码里调用cudaGetDeviceProperties查询，也可以直接看nvidia-smi的输出。

这里有个点要提醒：如果项目要部署到多台不同配置的机器上，尽量在代码里做动态能力检查，不要假设所有机器都能用最激进的优化选项。

2.3 在CMake工程中链接cuBLAS

工程的构建配置是很多新手卡壳的地方。我自己比较习惯用CMake管理CUDA项目，找包和链接都干净利落。

cmake_minimum_required(VERSION 3.18) project(cublas_demo LANGUAGES CXX CUDA) find_package(CUDAToolkit REQUIRED) add_executable(gemm_demo main.cu) target_link_libraries(gemm_demo PRIVATE CUDA::cublas)

CUDA::cublas是CMake提供的导入目标，会自动设置头文件路径、库路径和编译选项。这里补充一点：如果CMake版本比较老，找不到CUDAToolkit模块，可以直接手动指定CUDA_INCLUDE_DIRS和CUDA_LIBRARIES，但最好还是升级CMake以省心。

Windows下用Visual Studio的话也很方便，在项目属性里配置CUDA Toolkit的包含目录和库目录，然后把cublas.lib加进链接器输入就可以。不过相比命令行构建，VS里配置确实更加繁琐。

2.4 验证一个最小的cuBLAS示例

配置完环境后，我建议先跑一个最简单的示例验证链路通不通——直接调用cublasSgemm计算两个小矩阵的乘积。我会在下一节详细拆这个示例，这里先把编译命令给出。

nvcc -o gemm_demo gemm_demo.cu -lcublas

如果编译通过并且运行结果正确，说明你的工具链是通的，后面可以放心写更复杂的代码。如果链接报错，优先检查-lcublas是否加上了，以及LD_LIBRARY_PATH在运行时有没有配好。

3. 第一个cuBLAS程序：从SGEMM到完整流程

3.1 cublasHandle_t与库的初始化机制

cuBLAS的API设计跟很多C库不太一样，它的大部分操作都需要一个handle参数。这个cublasHandle_t可以理解为cuBLAS库的一个上下文对象，内部保存了算法配置、流绑定、工作空间等信息。

为什么要设计handle？因为创建一个cuBLAS上下文是有成本的，包括分配内部数据结构、初始化算法参数等。如果每个小操作都创建销毁，开销会非常大。正确的做法是：在程序开始时创建一个handle，中间反复使用，程序结束时销毁。

cublasHandle_t handle; cublasStatus_t status = cublasCreate(&handle); if (status != CUBLAS_STATUS_SUCCESS) { // 处理错误 } // ... 使用handle执行各种计算 ... cublasDestroy(handle);

这里还有一点，cuBLAS默认是在默认CUDA流上执行的。如果你想在多个流上并行执行不同的cuBLAS操作，可以调用cublasSetStream把handle绑定到特定流上。这个技巧在后面讲并发优化的时候会重点说，现在先记住handle和流是可以解耦的就够了。

3.2 主机端与设备端的内存管理对比

cuBLAS操作的是GPU显存，不是主机内存。这就意味着，你要先把数据从CPU拷到GPU，计算完再把结果拷回来。内存的分配和拷贝是GPU编程的基础，但常见的问题恰恰出在这里。

分配显存用cudaMalloc，拷贝用cudaMemcpy，释放用cudaFree。这里我特别强调一个容易出错的地方：cuBLAS的很多API支持非默认存储顺序，特别是列主序（column-major）和行主序（row-major）的区分。

cuBLAS遵循Fortran的惯例，默认按列主序解释矩阵。你在C++里习惯写行主序，直接传进cuBLAS就会得到错误的结果。解决的办法有三种：

• 在传入cuBLAS之前，把行主序矩阵转置成列主序。
• 调用cuBLAS的cublas<t>geam来对矩阵进行转置操作。
• 使用cuBLAS提供的cublasLt接口或者设置对应的transpose参数，在计算时指定转置。

个人最推荐第三种方案。因为显式的转置不仅浪费内存带宽，还增加了一次额外的拷贝。cuBLAS本身就有转置参数，完全可以在计算时把这个问题解决掉。

3.3 一个可运行的SGEMM完整示例

下面是一个完整的SGEMM示例，计算C = A * B，其中A、B、C都是1024x1024的float矩阵。

#include <cstdio> #include <cstdlib> #include <cuda_runtime.h> #include <cublas_v2.h> #define CHECK_CUDA(call) do { \ cudaError_t err = (call); \ if (err != cudaSuccess) { \ fprintf(stderr, "CUDA error at %s:%d: %s\n", \ __FILE__, __LINE__, cudaGetErrorString(err)); \ exit(EXIT_FAILURE); \ } \ } while (0) #define CHECK_CUBLAS(call) do { \ cublasStatus_t status = (call); \ if (status != CUBLAS_STATUS_SUCCESS) { \ fprintf(stderr, "cuBLAS error at %s:%d: %d\n", \ __FILE__, __LINE__, status); \ exit(EXIT_FAILURE); \ } \ } while (0) int main() { const int M = 1024; const int N = 1024; const int K = 1024; // 分配主机内存并初始化 float *h_A = (float*)malloc(M * K * sizeof(float)); float *h_B = (float*)malloc(K * N * sizeof(float)); float *h_C = (float*)malloc(M * N * sizeof(float)); 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; CHECK_CUDA(cudaMalloc((void**)&d_A, M * K * sizeof(float))); CHECK_CUDA(cudaMalloc((void**)&d_B, K * N * sizeof(float))); CHECK_CUDA(cudaMalloc((void**)&d_C, M * N * sizeof(float))); // 数据拷贝到显存 CHECK_CUDA(cudaMemcpy(d_A, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice)); CHECK_CUDA(cudaMemcpy(d_B, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice)); // 创建cuBLAS handle cublasHandle_t handle; CHECK_CUBLAS(cublasCreate(&handle)); // 执行 C = A * B float alpha = 1.0f; float beta = 0.0f; // 注意：cuBLAS默认列主序，这里用CUBLAS_OP_N表示不转置，按列主序解释数据 CHECK_CUBLAS(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, M, N, K, &alpha, d_A, M, d_B, K, &beta, d_C, M)); CHECK_CUDA(cudaDeviceSynchronize()); // 结果拷回主机 CHECK_CUDA(cudaMemcpy(h_C, d_C, M * N * sizeof(float), cudaMemcpyDeviceToHost)); // 验证：结果应该约等于 1.0 * 2.0 * K = 2048 bool ok = true; for (int i = 0; i < 16; i++) { if (fabs(h_C[i] - 2048.0f) > 1e-3) { ok = false; break; } } printf("Result: %s\n", ok ? "PASS" : "FAIL"); // 清理 cublasDestroy(handle); cudaFree(d_A); cudaFree(d_B); cudaFree(d_C); free(h_A); free(h_B); free(h_C); return 0; }

这段代码别看简单，它把cuBLAS的所有关键点都覆盖到了：handle的创建与销毁、显存的分配与回收、错误检查的宏定义、SGEMM参数的完整调用。你自己写项目的时候，直接把这个骨架扩展就行。

值得一提的是，cublasSgemm调用本身是异步的。它把计算任务提交到GPU上就返回了，真正的计算需要等GPU执行完。所以如果你紧接着就要用h_C的结果，必须先调用cudaDeviceSynchronize或者cudaMemcpy。代码里我在计算后放了cudaDeviceSynchronize，实际上cudaMemcpy本身也会隐式同步，但显式写出来更清晰，不容易出错。

3.4 错误处理：为什么你必须检查每个返回值

很多CUDA初学者不检查错误，程序跑挂了也不知道是哪里出的问题。CUDA和cuBLAS的API都有一个共同的特点：错误是在调用时返回的，但有的错误会污染后续的所有调用。

具体来说，如果一次cudaMalloc失败了，但你没有检查，后续的cudaMemcpy就会去访问非法指针，而错误状态会一直累积。这时候程序的表现非常诡异，可能崩溃，也可能返回一个看似正常的错误码。

我在自己的项目里会强制使用宏来做错误检查。上面示例里的CHECK_CUDA和CHECK_CUBLAS我很推荐引入到你的工程里。加上__FILE__和__LINE__，出问题的时候一眼就能定位到代码位置，排查效率翻倍。

4. 常用API与核心场景：不止矩阵乘法

4.1 cuBLAS的三级API对应关系

cuBLAS覆盖了BLAS的全部三级接口。这里把它们对应起来，方便你查用：

一级（向量-向量）：

• cublasSdot：点积
• cublasSaxpy：y = alpha * x + y
• cublasSnrm2：二范数

二级（矩阵-向量）：

• cublasSgemv：矩阵乘向量
• cublasSger：外积更新

三级（矩阵-矩阵）：

• cublasSgemm：通用矩阵乘法
• cublasSsymm：对称矩阵乘
• cublasStrmm：三角矩阵乘

每个函数根据数据类型不同还有前缀变化，比如cublasDgemm是双精度版本，cublasHgemm是半精度版本。命名规则就是BLAS前面加上cublas，数据类型缩写夹在中间，非常好记。

如果你要处理复数，cuBLAS也支持，cublasCgemm对应单精度复数，cublasZgemm对应双精度复数。做信号处理或者量子力学模拟的同学，这些接口直接用就行。

4.2 矩阵乘法的最优参数选择策略

一个很常见的问题是：cublasSgemm的参数到底该怎么填？尤其是lda、ldb、ldc（leading dimension，主维度跨度）这三个参数，最容易搞混。

在cuBLAS中，如果你的矩阵是列主序的，lda就是矩阵的行数M（不是列数）。为什么？因为列主序中相邻列之间间隔的是行数。这里我要画个重点：方案A是把数据转置后再调用，这是很多初学者会用但性能并不好的方式。方案B是直接利用API本身的转置参数，这才是最优实践。

比如C++内存里有一个行主序的MxK矩阵，你想计算C = A^T * B，在cuBLAS中你只需要把第一个转置参数设成CUBLAS_OP_T即可，不需要在代码里事先做任何转置操作。这样既少了一次额外拷贝，又避免了带宽浪费。

我还想提醒一点：cublasSgemm内部会为不同尺寸的矩阵选择不同的内核。对于特别小的矩阵（比如几行几列），建议走cublasSgemmStridedBatched批量处理，把大量小矩阵合批计算，吞吐会高很多。

4.3 批量矩阵运算的隐藏利器

机器学习里经常会遇到这种情况：一个batch里有几百张图，每张图都要做一轮矩阵运算。如果一张一张地调用cublasSgemm，每次调用都有启动开销，几百次下来性能损耗非常明显。

cuBLAS提供的batched API就是为这种场景设计的。cublasSgemmBatched接收一个指针数组，每个指针指向一个独立的矩阵；cublasSgemmStridedBatched更进一步，通过一个固定的stride来定位不同batch的数据，省去了指针数组的管理成本。

我在实现某个跨平台系统的推荐模块时，用户特征矩阵就是典型的batch场景。一开始用循环逐个调用，效率不理想；换成cublasSgemmStridedBatched之后，单位时间内处理的请求数提了几倍。这个改造其实不难，核心就是计算好stride，让数据在显存里连续排布。

4.4 cuBLAS的扩展生态：cuBLASLt与cuSOLVER

cuBLASLt是cuBLAS的轻量级高性能版本，它的特色是提供了更灵活的GEMM配置，可以精细控制矩阵布局、数据类型和算法的组合。如果你对性能有极致的追求，或者要处理混合精度的场景，值得深入研究这个接口。

cuSOLVER则是和cuBLAS配套的求解器库，提供LU分解、QR分解、特征值求解等高级功能。很多线性代数算法最终会落到这些分解上。我做有限元分析相关项目时，装配刚度矩阵之后就交给cuSOLVER来求解，比自己重写一套求解器可靠得多。

如果想把cuBLAS和cuSOLVER串起来，一个典型的路径是：先用cuBLAS做矩阵乘法得到中间结果，然后调用cuSOLVER的cusolverDnSgesv求解线性方程组。两者通过cudaStream_t可以共享一个流，避免额外的同步开销。

5. 性能优化与调优技巧：发挥GPU的每一分算力

5.1 避免主机与设备之间的频繁拷贝

GPU性能杀手排行榜上，主机和显存之间的PCIe传输绝对能排进前三，它是目前整个系统的最大短板之一，我特别强调这个问题是因为它太常见了。

一个典型反面案例是：循环计算里面，每个迭代都分配一次显存、拷贝一次数据、计算一次、再拷贝回主机。这样计算本身可能只要几毫秒，但数据拷贝占了大部分时间。

正确的思路是流水线化。把需要反复计算的数据常驻显存，尽量让数据在GPU侧流转。如果必须从主机侧读入数据，也应该用cudaMemcpyAsync配合内存映射，把拷贝时间和计算时间重叠起来。这个优化的核心原则就是：不要让GPU空等数据传输。

5.2 用多个CUDA流并行执行cuBLAS计算

一个handle只能绑定一个流，但你可以创建多个handle，每个绑定不同的流，让不同的操作并发执行。这个方法在同时处理多个独立任务的时候极其有效。

举个例子，某图像处理Demo里要同时对4张图像做矩阵变换，四张图之间没有依赖关系。我为它们创建了4个流、4个cuBLAS handle，把4个操作分别提交。在计算能力足够的GPU上，这些操作会并行执行，总耗时接近单张图像的处理时间。

需要注意两点：第一，流的创建不是没有代价的，流太多反而可能导致调度开销增加，建议跟GPU的实际SM数量匹配。第二，多流并发时，同一个handle不能被多个流同时使用，否则会产生数据竞争。要么给每个流单独创建handle，要么用锁保证同一个handle在同一时间只被一个流使用。

5.3 混合精度计算：Tensor Core的正确打开方式

NVIDIA从Volta架构开始引入Tensor Core，这是专门为矩阵乘加运算设计的硬件单元。它的特色是大幅提升吞吐量。Tensor Core原生支持的是混合精度：输入是FP16，累加用FP32。也就是说，你可以把A和B用半精度存，但累加过程用全精度，结果精度损失很小但速度飞快。

cuBLAS对混合精度的支持很有意思，cublasGemmEx这个API可以单独指定A、B、C的精类型、计算类型和算法选择。比如：

cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, M, N, K, &alpha, d_A, CUDA_R_16F, M, d_B, CUDA_R_16F, K, &beta, d_C, CUDA_R_16F, M, CUBLAS_COMPUTE_32F, CUBLAS_GEMM_DEFAULT_TENSOR_OP);

这段代码的含义是：输入输出都是FP16，但中间累加用FP32。对很多深度学习推理场景来说，这种配置是精度和性能的最佳平衡点。

需要提醒的是，如果你的卡不支持Tensor Core，调用CUBLAS_GEMM_DEFAULT_TENSOR_OP会失败或者退回普通内核。所以在代码里最好判断一下设备算力，再决定是否走Tensor Core路径。

5.4 规模小到一定程度，就该告别cuBLAS了

虽然cuBLAS已经很快了，但它毕竟有固定的调用开销。当矩阵尺寸小到一定程度（比如低于50x50），调用cuBLAS可能比自己写一个简单核函数还要慢。因为调用开销占比太高了。

这种场景下，你有两个选择：一是把多个小矩阵打包成一个大矩阵做批量计算，用cublasSgemmStridedBatched；二是直接写一个针对小矩阵的专用CUDA核函数，用常量内存或者寄存器缓存数据。

我自己在处理某跨平台系统里的小规模矩阵运算时，实测下来，对于16x16的矩阵，自己写的核函数比调用cuBLAS快了大约2到3倍。因为核函数专门为了这个尺寸做了展开优化，没有额外的API调度开销。建议你也试试，把cuBLAS和大规模计算绑定，把小矩阵场景交给定制内核。

5.5 性能分析的黄金工具：Nsight Compute与Nsight Systems

优化这件事，光靠猜肯定不行。NVIDIA官方的性能分析工具是必备利器。

Nsight Systems适合看全局：你可以看到kernel的启动时间、数据传输时间、空闲间隙，从宏观上找出瓶颈在哪里。Nsight Compute则深入到单个kernel内部，给出SM占用率、访存效率、指令流水线耗时等指标，能直接定位到是哪段代码拖慢了速度。

我习惯的调优流程是这样：先用Nsight Systems跑一遍完整程序，看看GPU利用率高不高。如果利用率低，优先处理数据搬运和kernel依赖问题。如果利用率已经很高了，再用Nsight Compute针对最耗时的kernel做微观分析，看它的计算吞吐是不是接近理论峰值。

记住一个经验：做性能分析之前，先把非必要的日志输出全部关掉，否则IO操作会严重干扰分析结果。

6. 常见问题与排查技巧实录

6.1 我最常遇到的几个cuBLAS错误

cuBLAS的错误码一开始看的人一头雾水，但踩过几次坑之后，其实非常规律。给你整理一个速查表：

其中CUBLAS_STATUS_INVALID_VALUE是最常见的，浪费掉的时间也最多。我教一个排查技巧：从调用链最内层开始检查，先看ld参数是否大于等于对应的维度，再看M、N、K是否为正数，最后检查d_A、d_B、d_C有没有被错误释放。

6.2 为什么结果不对：列主序与行主序的陷阱

如果你把行主序矩阵直接传进cuBLAS，通常会得到一份很奇怪的结果。这不是cuBLAS算错了，而是你把矩阵的行列顺序搞反了，这个问题的出现频率远超你想象。

解决方式上面提过，最优雅的是用transpose参数，但这里补充一个更隐蔽的坑：当你用CUBLAS_OP_T时，对应的ld含义也会变化。比如原来矩阵A是MxK行主序，你传CUBLAS_OP_T表示按转置后的KxM来读，那么这时lda实际上应该传K（原始矩阵的行数）而不是M。

还有个更隐蔽的问题：非方阵转置后，行数变成了原来的列数。很多人参数填得稀里糊涂，导致程序时对时错。我的建议是：写一个小规模的测试用例，比如3x5的矩阵，手动计算期望结果，跟cuBLAS输出做对比。跑通之后再去改大尺寸参数，基本就不会再出错。

6.3 性能忽高忽低，首先查流的同步问题

有一种现象很让人费解：函数跑起来性能不稳定，第一次调用的耗时特别长，后面几次恢复正常。这通常是因为第一次调用时，CUDA上下文和cuBLAS内部的工作空间还没有初始化，加上GPU频率还没拉升，所以显得特别慢。

如果性能持续忽高忽低，就要重点看有没有隐式同步。最常见的情况是：每个cudaMemcpy调用都同步了全局执行，把流水线打断了。解决办法是改用cudaMemcpyAsync，并且确保拷贝的源、目标内存是页锁定内存（pinned memory），也就是用cudaMallocHost分配的主机端内存。

页锁定内存这个概念也值得多说一句：普通的malloc分配的内存，操作系统可能会换页到磁盘上；而cudaMallocHost分配的内存锁在物理内存里，不会被换出。这样DMA控制器就可以直接访问它，不需要CPU介入。这个细节对传输性能的影响非常显著。

6.4 显存不够用怎么办

大矩阵运算极其消耗显存。比如三个4096x4096的float矩阵，每个就有64MB，三个就是192MB，看起来还能接受。如果尝试65536x65536，一个矩阵就是16GB，显存直接爆掉。

解决办法通常有这几条路：

• 使用FP16或BF16等低精度类型存储数据，显存占用直接减半。
• 分块计算，但要注意分块策略和汇总方式。
• 考虑用统一虚拟内存或cudaMallocManaged，但这会牺牲性能，因为数据可能被换到主机内存。

我个人经验是，能降精度就优先降精度。在科学计算场景里，很多时候FP32和FP16的结果差距在可接受范围内，但速度和显存占用差距是巨大的。

6.5 多线程环境的坑：handle不是线程安全的

如果你的程序用多线程分别往GPU提交任务，每个线程都创建一个cuBLAS handle，这个模式本身没问题。但如果你共用一个handle，问题就大了——cuBLAS的handle在内部有状态，多个线程同时调用同一个handle，轻则结果错乱，重则直接崩溃。

正确做法是：每个线程独立创建自己的handle，线程结束时销毁。代价是创建handle的开销会略微增加，但这在多线程场景几乎可以忽略。

另一个细节是，多个线程共享同一个CUDA context时，流的创建也需要做同步管理。如果两个线程同时创建流，可能触发CUDA内部锁竞争。建议预先创建好所有需要的流，之后每个线程只使用自己拿到的流，避免运行时再创建。

7. cuBLAS之外的思考：什么情况下自己编写CUDA核函数

7.1 常规矩阵运算交给cuBLAS，特殊需求交给自己

cuBLAS不是万能的。当你需要的不是标准的C = A * B，而是带特殊约简规则的运算，比如每行取最大值再做归一化，或者矩阵元素之间有自定义的依赖关系，这时候就需要自己写CUDA核函数。

判断标准很简单：如果是标准BLAS操作，哪怕它看起来很复杂，都先查一下cuBLAS是否已经提供了。如果没有提供，再考虑自己写。这是性价比最高的路径，避免重复造轮子，也避免被cuBLAS的性能差距碾压。

7.2 自研核函数时的性能基线参考

写自己的核函数，首先要心里有数：性能目标定在哪里。我通常以cuBLAS同类型操作作为基线。比如自己实现SGEMM，性能应该是cuBLAS的多少。如果低于60%到70%，说明优化不到位。

优化方向按照优先级排列：先确保访存合并，每个线程访问的连续地址最好落在同一缓存行；然后用共享内存做数据复用；再考虑向量化加载float4；最后考虑指令级优化和循环展开。

这里有一个从cuBLAS源码和NVIDIA博客里学到的通用技巧：分块大小通常选8x8或16x16，配合共享内存双缓冲，可以让计算和访存重叠。我按照这个思路优化自己的内核后，性能从cuBLAS的30%提升到了65%左右，虽然跟cuBLAS还是有差距，但在自定义算法里已经算很能打了。

7.3 混合方案：自己写预/后处理，矩阵计算交给cuBLAS

更成熟的工程架构，其实是混合方案。自定义的预处理核函数负责数据标准化、填充、格式转换；中间的矩阵乘法、批量变换交给cuBLAS；后处理核函数再负责激活、池化、归一化。

这样分工的优点很明显：既有定制逻辑的灵活性，又有关键运算的高性能。我参与研发的某图像处理Demo就是这么设计卷积层的，预处理和后处理用自研核函数，中间的主矩阵运算是cuBLAS，最终性能和开发效率都让人比较满意。

8. 最后的个人经验与建议

cuBLAS的学习曲线并不陡峭，但有一个习惯我特别建议你培养：每个新函数第一次用的时候，都写一个最小的验证程序，跟CPU参考实现比对结果。不要嫌麻烦，这个习惯能省掉后面大量调试时间。矩阵运算这东西，错一个符号就是完全不同的结果，而且很难肉眼发现。

在实际项目中，我还习惯于把cuBLAS的调用封装成自己的模块，对外提供更简洁的接口，比如MatrixMultiply(A, B, C, transA, transB)这种形式。这样做的好处是，以后想换底层实现（比如换成cuBLASLt或者自研内核），上层业务代码完全不用动。

还有一个容易被忽略的点：编译选项。nvcc编译CUDA代码时，最好加上-O3和对应的架构优化参数，比如-arch=sm_80。不同架构的参数差别不小，如果用了老架构参数，新卡也没法发挥性能。建议用-arch=native或者提前明确目标设备的算力。

回到开头那句话：矩阵运算又关键又烦人。有了cuBLAS，烦人的部分基本被解决了，你能把精力集中在真正有业务价值的地方。而从cuBLAS的源码、文档和性能特征里学到的东西，也会反向帮助你写出更好的自研内核。如果你现在正要开始一个涉及GPU矩阵运算的项目，我的建议很简单：先别急着写核函数，把cuBLAS请出来，让它先跑起来再说。至于性能优化，等遇到瓶颈了再回头折腾也不迟。