人工智能 cuda

Cublas相关

cublas学习笔记

Posted by donpromax on September 28, 2025

https://docs.nvidia.com/cuda/cublas/

Data Layout

cublas遵循 Fortran 风格的列优先（column-major order）内存布局。因此使用下面的宏可以计算对应的C/C++风格的0-based一维数组索引：

// i: 行索引, j: 列索引, ld: 矩阵的 leading dimension
#define IDX2C(i,j,ld) (((j)*(ld))+(i))

simple sgemm case

cublasSgemm 是 cuBLAS 库 中最常用的 BLAS Level 3 函数之一，用于执行 单精度浮点数矩阵乘法（Single-precision General Matrix-Matrix Multiplication）。

cublasSgemm

函数原型（API 定义）

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返回CUBLAS_STATUS_SUCCESS时表示成功。详细定义见：https://docs.nvidia.com/cuda/cublas/#cublasstatus-t

返回值 含义
CUBLAS_STATUS_SUCCESS 成功
CUBLAS_STATUS_NOT_INITIALIZED cuBLAS handle 未创建
CUBLAS_STATUS_INVALID_VALUE m < 0 或 n < 0 或 k < 0

参数详解

| 参数 | 类型 | 说明 |
| — | — | — |
| handle | cublasHandle_t | cuBLAS 上下文句柄，必须通过 cublasCreate() 创建 |
| transa | cublasOperation_t | 矩阵 A 是否转置： CUBLAS_OP_N: 不转置 CUBLAS_OP_T: 转置 CUBLAS_OP_C: 共轭转置（对实数等同于 T） |
| transb | cublasOperation_t | 矩阵 B 是否转置（同上） |
| m | int | 矩阵 C 的行数（也是 A 的行数，若 A 不转置） |
| n | int | 矩阵 C 的列数（也是 B 的列数，若 B 不转置） |
| k | int | 内维度（A 的列数，若 A 不转置；B 的行数，若 B 不转置） |
| alpha | const float* | 指向标量 $ \alpha $ 的指针（可位于主机或设备内存） |
| A | const float* | 指向矩阵 A 的设备指针 |
| lda | int | 矩阵 A 的 leading dimension（主维度），即按列主序存储时 A 的行数（或分配的行跨度） |
| B | const float* | 指向矩阵 B 的设备指针 |
| ldb | int | 矩阵 B 的 leading dimension |
| beta | const float* | 指向标量 $ \beta $ 的指针 |
| C | float* | 指向矩阵 C 的设备指针（结果矩阵） |
| ldc | int | 矩阵 C 的 leading dimension |

功能说明

cublasSgemm 执行如下操作：

$ C = \alpha \cdot \text{op}(A) \cdot \text{op}(B) + \beta \cdot C $

其中：

  • $ \text{op}(A) $ 是 $ A $ 或 $ A^T $（或 $ A^H $）
  • $ \text{op}(B) $ 是 $ B $ 或 $ B^T $
  • $ \alpha $ 和 $ \beta $ 是标量
  • 所有矩阵都为 单精度浮点型（float）
  • 所有矩阵在内存中按 列主序（column-major order） 存储（Fortran 风格）

使用示例

示例：计算 $ C = 1.0 \cdot A \times B + 0.0 \cdot C $

cublasHandle_t handle;
float alpha = 1.0f, beta = 0.0f;
float *d_A, *d_B, *d_C;  // 设备内存中的矩阵 A, B, C
int m = 6, n = 4, k = 5;

// 假设 d_A: m×k, d_B: k×n, d_C: m×n
// lda = m, ldb = k, ldc = m

cublasCreate(&handle);

cublasSgemm(
    handle,
    CUBLAS_OP_N,        // A 不转置
    CUBLAS_OP_N,        // B 不转置
    m, n, k,            // C: m×n, A: m×k, B: k×n
    &alpha,             // α = 1.0
    d_A, m,             // A, leading dimension = m
    d_B, k,             // B, leading dimension = k
    &beta,              // β = 0.0
    d_C, m              // C, leading dimension = m
);

cublasDestroy(handle);

等价于：

$ C = A \times B $

支持的数据类型变体

Sgemm 中的 S 表示 Single precision。其他变体包括：

函数 数据类型 描述
cublasSgemm float 单精度实数矩阵乘法
cublasDgemm double 双精度实数矩阵乘法
cublasCgemm cuComplex 单精度复数矩阵乘法
cublasZgemm cuDoubleComplex 双精度复数矩阵乘法

示例：

cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha_d, d_A, lda, d_B, ldb, &beta_d, d_C, ldc);

leading dimension（lda, ldb, ldc）详解

  • 在列主序中，矩阵按列存储。
  • lda 是分配给矩阵 A 的行数跨度（可能大于实际行数，用于对齐）。
  • 例如：A 是 $ m \times k $，按列存储，每列有 lda 个元素（lda >= m）。

即使 m=6，如果分配时对齐为 8，则 lda=8。

example full code

SGEMM的简单示例：

$ \begin{bmatrix}1 & 3 \ 2 & 4\end{bmatrix}
\times
\begin{bmatrix}5 & 7 \ 6 & 8\end{bmatrix}
=
\begin{bmatrix}23 & 31 \ 34 & 46\end{bmatrix} $

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

// 定义列优先0-based索引宏
#define IDX2C(i, j, ld) (((j) * (ld)) + (i))

int main() {
    const int m = 2, n = 2, k = 2;

    float A[m * k];
    float B[k * n];
    float C[m * n] = {0};

    // 初始化 A (2x2):   [1  3]
    //                  [2  4]
    A[IDX2C(0,0,m)] = 1.0f;
    A[IDX2C(1,0,m)] = 2.0f;
    A[IDX2C(0,1,m)] = 3.0f;
    A[IDX2C(1,1,m)] = 4.0f;

    // 初始化 B (2x2):   [5  7]
    //                  [6  8]
    B[IDX2C(0,0,k)] = 5.0f;
    B[IDX2C(1,0,k)] = 6.0f;
    B[IDX2C(0,1,k)] = 7.0f;
    B[IDX2C(1,1,k)] = 8.0f;

    // 设备指针
    float *d_A, *d_B, *d_C;
    cublasHandle_t handle;

    cudaMalloc(&d_A, m * k * sizeof(float));
    cudaMalloc(&d_B, k * n * sizeof(float));
    cudaMalloc(&d_C, m * n * sizeof(float));

    // Host → Device 
    cudaMemcpy(d_A, A, m * k * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, B, k * n * sizeof(float), cudaMemcpyHostToDevice);

    // 创建 cuBLAS 句柄
    cublasCreate(&handle);

    // C = alpha * A * B + beta * C
    const float alpha = 1.0f, beta = 0.0f;
    cublasSgemm(handle,
                CUBLAS_OP_N, CUBLAS_OP_N,  // 不转置
                m, n, k,                   // 维度
                &alpha,
                d_A, m,                    // A, leading dimension = m
                d_B, k,                    // B, leading dimension = k
                &beta,
                d_C, m);                   // C, leading dimension = m

    // Device → Host
    cudaMemcpy(C, d_C, m * n * sizeof(float), cudaMemcpyDeviceToHost);

    // 使用 IDX2C 宏打印结果矩阵 C
    printf("C = \n");
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            printf("%.2f ", C[IDX2C(i, j, m)]);
        }
        printf("\n");
    }

    cublasDestroy(handle);
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);

    return 0;
}

编译：

nvcc -o sgemm sgemm.cu -lcublas

output:

C = 
23.00 31.00 
34.00 46.00

official scal example

https://docs.nvidia.com/cuda/cublas/#example-code

cublasSscal

cublasSscal 是 cuBLAS 库中的一个 BLAS Level 1 函数，用于对 单精度浮点数向量 执行 缩放操作（scaling），即$ x \leftarrow \alpha \cdot x $

函数原型（API 定义）

cublasStatus_t cublasSscal(
    cublasHandle_t handle,
    int n,
    const float* alpha,
    float* x,
    int incx
);

cublasStatus_t返回CUBLAS_STATUS_SUCCESS时表示成功。详细定义见：https://docs.nvidia.com/cuda/cublas/#cublasstatus-t

返回值 含义
CUBLAS_STATUS_SUCCESS 成功
CUBLAS_STATUS_NOT_INITIALIZED cuBLAS handle 未创建
CUBLAS_STATUS_INVALID_VALUE n < 0 或 incx < 0

参数详解

| 参数 | 类型 | 说明 |
| — | — | — |
| handle | cublasHandle_t | cuBLAS 库的句柄，必须先通过 cublasCreate() 创建 |
| n | int | 向量 x 中参与运算的元素个数 |
| alpha | const float* | 指向标量 alpha 的设备或主机指针，用于缩放：$ x := \alpha \cdot x $ |
| x | float* | 指向单精度实数向量 x 的设备指针 |
| incx | int | 向量 x 的增量（stride），即访问元素的间隔步长 |

功能说明

该函数执行如下操作：

$ x_i \leftarrow \alpha \cdot x_i \quad \text{for } i = 0, \text{incx}, 2\cdot\text{incx}, \dots, (n-1)\cdot\text{incx} $

即：从 x[0] 开始，每隔 incx 个元素取一个，共取 n 个元素，全部乘以标量 alpha。

使用示例

示例 1：缩放整个向量（incx = 1）

float alpha = 2.0f;
float *d_x;  // device vector of size 5
int n = 5;
int incx = 1;

cublasSscal(handle, n, &alpha, d_x, incx);
// 结果：d_x[i] *= 2.0f, for i = 0..4

示例 2：缩放偶数索引元素（incx = 2）

cublasSscal(handle, 3, &alpha, d_x, 2);
// 操作：d_x[0], d_x[2], d_x[4] 乘以 alpha（共 3 个元素）

内存位置说明

  • alpha 可以位于 主机或设备内存 中（cuBLAS 会自动检测）。
  • x 必须位于 GPU 设备内存 中。

支持的数据类型变体

Sscal 中的 S 表示 Single precision（单精度浮点）。cuBLAS 提供了多种数据类型的版本：

函数 数据类型 描述
cublasSscal float 单精度实数向量缩放
cublasDscal double 双精度实数向量缩放
cublasCscal cuComplex 单精度复数向量缩放alpha为复数
cublasZscal cuDoubleComplex 双精度复数向量缩放alpha为复数
cublasCsscal cuComplex 单精度复数向量缩放alpha 为实数
cublasZdscal cuDoubleComplex 双精度复数向量缩放alpha 为实数

示例：

cublasDscal(handle, n, &alpha_double, d_x_double, incx);

example full code

// Example 2 Modified: Use cudaMemcpy2D instead of cublasSetMatrix/cublasGetMatrix
//-----------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda_runtime.h>
#include "cublas_v2.h"

#define M 6
#define N 5
#define IDX2C(i,j,ld) (((j)*(ld))+(i))  // 0-based row-major indexing in column-major storage

static __inline__ void modify(cublasHandle_t handle, float *m, int ldm, int n, int p, int q, float alpha, float beta) {
    // Scale row p from column q to n-1 (horizontal)
    cublasSscal(handle, n - q, &alpha, &m[IDX2C(p, q, ldm)], ldm); // stride is ldm for row
    // Scale column q from row p to ldm-1 (vertical)
    cublasSscal(handle, ldm - p, &beta, &m[IDX2C(p, q, ldm)], 1); // stride is 1 for column
}

int main(void) {
    cudaError_t cudaStat;
    cublasStatus_t stat;
    cublasHandle_t handle;
    int i, j;
    float* devPtrA;
    float* a = NULL;

    // Host memory allocation
    a = (float*)malloc(M * N * sizeof(*a));
    if (!a) {
        printf("Host memory allocation failed\n");
        return EXIT_FAILURE;
    }

    // Initialize matrix with 0-based indexing, column-major layout
    for (i = 0; i < M; i++) {
        for (j = 0; j < N; j++) {
            a[IDX2C(i, j, M)] = (float)(i * N + j + 1);
        }
    }

    // Device memory allocation
    cudaStat = cudaMalloc((void**)&devPtrA, M * N * sizeof(*a));
    if (cudaStat != cudaSuccess) {
        printf("Device memory allocation failed: %s\n", cudaGetErrorString(cudaStat));
        free(a);
        return EXIT_FAILURE;
    }

    // Create cuBLAS handle
    stat = cublasCreate(&handle);
    if (stat != CUBLAS_STATUS_SUCCESS) {
        printf("CUBLAS initialization failed\n");
        free(a);
        cudaFree(devPtrA);
        return EXIT_FAILURE;
    }

    // --- 替换 cublasSetMatrix: 使用 cudaMemcpy2D 上传数据 ---
    cudaStat = cudaMemcpy2D(
        devPtrA,              // dst: device memory pointer
        M * sizeof(float),    // dpitch: leading dimension on device (in bytes)
        a,                    // src: host matrix
        M * sizeof(float),    // spitch: leading dimension on host (in bytes)
        M * sizeof(float),    // width  of matrix in bytes (number of rows in float bytes)
        N,                    // height of matrix (number of columns)
        cudaMemcpyHostToDevice
    );
    if (cudaStat != cudaSuccess) {
        printf("Data upload to device failed: %s\n", cudaGetErrorString(cudaStat));
        free(a);
        cudaFree(devPtrA);
        cublasDestroy(handle);
        return EXIT_FAILURE;
    }

    // Perform modification on GPU
    modify(handle, devPtrA, M, N, 1, 2, 16.0f, 12.0f);

    // --- 替换 cublasGetMatrix: 使用 cudaMemcpy2D 下载数据 ---
    cudaStat = cudaMemcpy2D(
        a,                    // dst: host memory pointer
        M * sizeof(float),    // dpitch: host leading dimension (in bytes)
        devPtrA,              // src: device matrix
        M * sizeof(float),    // spitch: device leading dimension (in bytes)
        M * sizeof(float),    // width in bytes (rows)
        N,                    // height (columns)
        cudaMemcpyDeviceToHost
    );
    if (cudaStat != cudaSuccess) {
        printf("Data download from device failed: %s\n", cudaGetErrorString(cudaStat));
        free(a);
        cudaFree(devPtrA);
        cublasDestroy(handle);
        return EXIT_FAILURE;
    }

    // Cleanup
    cudaFree(devPtrA);
    cublasDestroy(handle);

    // Print result
    for (i = 0; i < M; i++) {
        for (j = 0; j < N; j++) {
            printf("%7.0f", a[IDX2C(i, j, M)]);
        }
        printf("\n");
    }

    free(a);
    return EXIT_SUCCESS;
}

相较官方example，修改了a的遍历打印方向。以及将cublasSetMatrix、cublasGetMatrix改为了cudaMemcpy2D。modify的作用：

  • 缩放 第 p 行，从第 q 列到最后一列（横向）alpha
  • 缩放 第 q 列，从第 p 行到最后一行（纵向）beta

output:

      1      2      3      4      5
      6      7   1536    144    160
     11     12    156     14     15
     16     17    216     19     20
     21     22    276     24     25
     26     27    336     29     30
CATALOG
    FEATURED TAGS
    Travelling 人工智能 Japan cuda
    FRIENDS