【cuda】六、基础库：cuBLAS入门

大家好，欢迎来到IT知识分享网。

cuBLAS 基础

介绍

CUDA Basic Linear Algebra Subprograms（BLAS）提供了高效计算线性代数的方法。

有三级API和cuBLAS 扩展、辅助API：

1. 最基础操作，例如加、减、最大值、复制、转置
2. 矩阵的一般操作，例如特殊类型矩阵的乘法、rank
3. 更复杂一些的例子，例如“使用一般矩阵计算批量的矩阵-矩阵乘积”，‘使用高斯复杂度降低算法计算一般矩阵的矩阵-矩阵乘积’

API介绍：https://docs.nvidia.com/cuda/cublas/index.html

样例代码：https://github.com/NVIDIA/CUDALibrarySamples/tree/master/cuBLAS

功能：

1. 向量和矩阵操作：包括向量加法、向量-标量乘法、向量点积等。
2. 矩阵乘法：支持各种形式的矩阵乘法，包括方阵乘法、矩阵-向量乘法等。
3. 分解和求逆：例如LU分解、Cholesky分解和矩阵求逆等。
4. 求解线性系统：使用不同的方法解决线性方程组。

使用教程

使用前需要在linker中注明cuBLAS.lib

使用 cuBLAS API，应用程序必须在 GPU 内存空间中分配所需的矩阵和向量，用数据填充它们，调用所需的 cuBLAS 函数序列，然后将结果从 GPU 内存空间上传回主机。cuBLAS API 还提供用于从 GPU 写入和检索数据的辅助函数。

需要注意的是，cuBLAS 库使用列存储，同时首项从1开始计算。

为了最大限度地兼容现有的 Fortran 环境，cuBLAS 库使用列主存储和基于 1 的索引。由于 C 和 C++ 使用行优先存储，因此用这些语言编写的应用程序无法使用二维数组的本机数组语义。

应用程序必须通过调用cublasCreate()函数来初始化 cuBLAS 库上下文的句柄。然后，该句柄被显式传递给每个后续的库函数调用。一旦应用程序完成使用库，它必须调用函数cublasDestroy()来释放与 cuBLAS 库上下文关联的资源。

这种方法允许用户在使用多个主机线程和多个 GPU 时显式控制库设置。例如，应用程序可以将cudaSetDevice()不同的设备与不同的主机线程关联起来，并且在每个主机线程中，它可以初始化 cuBLAS 库上下文的唯一句柄，该句柄将使用与该主机线程关联的特定设备。然后，使用不同句柄进行的 cuBLAS 库函数调用将自动将计算分派到不同的设备。

数据类型

cublasHandle_t是指向保存 cuBLAS 库上下文的不透明结构的指针类型。cuBLAS 库上下文必须使用cublasCreate()进行初始化，并且返回的句柄必须传递给所有后续库函数调用。最后应使用cublasDestroy()销毁上下文。

cublasStatus_t该类型用于函数状态返回。所有 cuBLAS 库函数都会返回其状态。具体报错查看：https://docs.nvidia.com/cuda/cublas/index.html#cublasstatus-t

cublasOperation_t类型指示需要对稠密矩阵执行哪种操作。它的值对应于 Fortran 字符‘N’or ‘n’（非转置）、‘T’or ‘t’（转置）和‘C’or ‘c’（共轭转置），这些字符通常用作传统 BLAS 实现的参数。

cublasFillMode_t类型指示密集矩阵的哪一部分（下部或上部）已填充，因此应由函数使用。它的值对应于 Fortran 字符Lor l（下）和Uor u（上），这些字符通常用作旧版 BLAS 实现的参数。

基础元素

使用 cuBLAS 库必备的步骤：

1. 包含必要的头文件

在程序的开头包含 cuBLAS 和 CUDA 运行时库的头文件。

#include <cublas_v2.h> #include <cuda_runtime.h> 

2. 初始化 cuBLAS 和 CUDA

初始化 cuBLAS 库和配置 CUDA 设备。

cublasHandle_t cublas_handle; cudaSetDevice(device_id); // 可选，设置CUDA设备 cublasCreate(&cublas_handle); 

3. 分配内存

在 GPU 上分配必要的内存空间来存储你的数据（比如矩阵或向量）。使用 CUDA 的 cudaMalloc 函数进行分配。

double* d_A; cudaMalloc((void**)&d_A, sizeof(double) * size_of_A); // 其他变量类似 

4. 数据传输

将数据从主机内存复制到 GPU 内存。使用 cudaMemcpy 函数从主机到设备进行数据传输。

cudaMemcpy(d_A, h_A, sizeof(double) * size_of_A, cudaMemcpyHostToDevice); // 其他变量类似 

5. 执行 cuBLAS 操作

调用 cuBLAS 函数执行所需的线性代数运算。例如，执行矩阵乘法或向量加法。

cublasDgemm(cublas_handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, d_A, lda, d_B, ldb, &beta, d_C, ldc) 

6. 从 GPU 获取结果

计算完成后，将结果从 GPU 内存复制回主机内存。

cudaMemcpy(h_C, d_C, sizeof(double) * size_of_C, cudaMemcpyDeviceToHost); 

7. 清理

释放 GPU 上分配的内存，并销毁 cuBLAS 句柄。

cudaFree(d_A); // 释放其他分配的内存 cublasDestroy(cublas_handle); 

8. 关闭 CUDA 设备（可选）

如果需要，可以重置 CUDA 设备以清理所有状态。

cudaDeviceReset(); 

重要提示

• 确保所有 CUDA 和 cuBLAS 调用都成功。可以通过检查每个调用的返回状态来实现。
• 在进行大量的 CUDA/cuBLAS 操作时，考虑使用错误检查宏或函数来简化代码和调试。
• 对于复杂的程序，考虑使用 CUDA 流来管理并行执行和数据传输。

遵循这些步骤可以确保你的 cuBLAS

程序能够正确地执行线性代数运算，同时充分利用 GPU 的计算能力。在实际应用中，你可能需要根据具体的计算任务调整这些步骤，比如选择不同的 cuBLAS 函数或处理不同大小和类型的数据。

复杂案例：计算三角带状矩阵向量乘法

来源于cuBLAS 2 级 API –cublas<t>tbmv ：https://github.com/NVIDIA/CUDALibrarySamples/tree/master/cuBLAS/Level-2/tbmv

#include <cstdio> #include <cstdlib> #include <vector> #include <cublas_v2.h> #include <cuda_runtime.h> #include "cublas_utils.h" using data_type = double; // 定义数据类型为 double int main(int argc, char *argv[]) { 
    cublasHandle_t cublasH = NULL; // 声明一个 cuBLAS 句柄 cudaStream_t stream = NULL; // 声明一个 CUDA 流 const int m = 2; // 定义矩阵A的行数 const int n = 2; // 定义矩阵A的列数 const int k = 1; // 定义超对角线元素的个数（用于三角矩阵的函数） const int lda = m; // 定义矩阵A的领先维度（leading dimension） // 初始化矩阵A和向量x const std::vector<data_type> A = { 
   1.0, 3.0, 2.0, 4.0}; // 矩阵A std::vector<data_type> x = { 
   5.0, 6.0}; // 向量x const int incx = 1; // x的步长 data_type *d_A = nullptr; // 设备端的矩阵A data_type *d_x = nullptr; // 设备端的向量x // cuBLAS相关设置 cublasFillMode_t uplo = CUBLAS_FILL_MODE_UPPER; // 使用上三角形式 cublasOperation_t transa = CUBLAS_OP_N; // 矩阵A不进行转置 cublasDiagType_t diag = CUBLAS_DIAG_NON_UNIT; // 矩阵A的对角线元素不被视为1 printf("A\n"); print_matrix(m, n, A.data(), lda); // 打印矩阵A printf("=====\n"); printf("x\n"); print_vector(x.size(), x.data()); // 打印向量x printf("=====\n"); // 步骤1: 创建 cuBLAS 句柄，绑定一个流 CUBLAS_CHECK(cublasCreate(&cublasH)); CUDA_CHECK(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking)); CUBLAS_CHECK(cublasSetStream(cublasH, stream)); // 步骤2: 将数据拷贝到设备端 CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&d_A), sizeof(data_type) * A.size())); CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&d_x), sizeof(data_type) * x.size())); CUDA_CHECK(cudaMemcpyAsync(d_A, A.data(), sizeof(data_type) * A.size(), cudaMemcpyHostToDevice, stream)); CUDA_CHECK(cudaMemcpyAsync(d_x, x.data(), sizeof(data_type) * x.size(), cudaMemcpyHostToDevice, stream)); // 步骤3: 执行计算 // 使用 cuBLAS 的 tbmv 函数进行三角带状矩阵和向量的乘法 CUBLAS_CHECK(cublasDtbmv(cublasH, uplo, transa, diag, n, k, d_A, lda, d_x, incx)); // 步骤4: 将计算结果从设备拷贝回主机 CUDA_CHECK(cudaMemcpyAsync(x.data(), d_x, sizeof(data_type) * x.size(), cudaMemcpyDeviceToHost, stream)); // 同步 CUDA 流以确保所有操作完成 CUDA_CHECK(cudaStreamSynchronize(stream)); /* * x = | 27.00 24.00 | */ printf("x\n"); print_vector(x.size(), x.data()); // 打印计算后的向量x printf("=====\n"); // 释放资源 CUDA_CHECK(cudaFree(d_A)); // 释放设备上的矩阵A CUDA_CHECK(cudaFree(d_x)); // 释放设备上的向量x CUBLAS_CHECK(cublasDestroy(cublasH)); // 销毁 cuBLAS 句柄 CUDA_CHECK(cudaStreamDestroy(stream)); // 销毁 CUDA 流 CUDA_CHECK(cudaDeviceReset()); // 重置 CUDA 设备 return EXIT_SUCCESS; // 程序正常退出 }