cuda编程笔记（16）--GPU上的稀疏矩阵

稀疏矩阵（Sparse Matrix）是大多数元素为 0 的矩阵。直接存储所有元素会浪费内存，所以常用压缩存储格式：

CSR（Compressed Sparse Row）

结构：

• values：存储非零元素（按行展开）

• col_indices：对应列索引

• row_ptr：每行的起始位置在 values 中的索引，长度 = 行数 + 1，最后一个位置存储values 的总长度

A = [10  0  0  0
      0 20  0 30
      0  0 40  0]
      
values   = [10, 20, 30, 40]
col_idx  = [0, 1, 3, 2]
row_ptr  = [0, 1, 3, 4]   // row_ptr[i+1]-row_ptr[i] = 第 i 行非零元素个数

• row_ptr：每行在 values 中的起始位置索引

  • 第一行：非零元素 10 在 values[0] → row_ptr[0] = 0

  • 第二行：非零元素 20, 30 在 values[1] 和 values[2] → row_ptr[1] = 1

  • 第三行：非零元素 40 在 values[3] → row_ptr[2] = 3

row_ptr 的长度 = 行数 + 1，最后一个元素是 values 的总长度

为何说row_ptr[i+1]-row_ptr[i] = 第 i 行非零元素个数？

用具体例子：

• 第 0 行：row_ptr[0]=0, row_ptr[1]=1 → 只有 values[0] = 10

• 第 1 行：row_ptr[1]=1, row_ptr[2]=3 → values[1]=20, values[2]=30

• 第 2 行：row_ptr[2]=3, row_ptr[3]=4 → values[3]=40

这样就明确了每行非零元素在 values 中的区间。

优点：

• 行遍历高效

• 内存紧凑

缺点：

• 列访问不方便

ELLPACK / ELL

结构：

• 每行存储固定数量的非零元素（假设 max_nnz_per_row = 最大行非零数）

• values[row][k] 和 col_indices[row][k]

A = [10 0 0 0
     0 20 0 30
     0 0 40 0]

max_nnz_per_row = 2

values    = [[10, 0],
             [20, 30],
             [40, 0]]

col_idx   = [[0, -1],
             [1, 3],
             [2, -1]]    // -1 表示无效

优点：

• GPU 并行处理好（每行一个线程）

• 内存访问对齐

缺点：

• 如果行非零元素差异大，内存浪费

SpMV 操作

目标：
给定稀疏矩阵 A 和向量 x，计算 y = A * x。

CSR SpMV：

用伪代码表示

for i in 0..num_rows-1:
    y[i] = 0
    //一行乘向量的结果存到y[i]
    for j in row_ptr[i] .. row_ptr[i+1]-1:
        y[i] += values[j] * x[col_idx[j]]

ELL SpMV：

for i in 0..num_rows-1:
    y[i] = 0
    for k in 0..max_nnz_per_row-1:
        if col_idx[i][k] != -1:
            y[i] += values[i][k] * x[col_idx[i][k]]

CSR 更适合 CPU，ELL 对 GPU 友好（内存访问规整）。

GPU 实现（cuSPARSE）

cuSPARSE 提供高性能稀疏矩阵运算接口。

//初始化句柄
cusparseHandle_t handle;
cusparseCreate(&handle);

创建 CSR 矩阵描述符

cusparseStatus_t cusparseCreateCsr(
    cusparseSpMatDescr_t *spMat,
    int64_t rows, int64_t cols, int64_t nnz,
    void *csrRowPtr, void *csrColInd, void *csrValues,
    cusparseIndexType_t rowType,
    cusparseIndexType_t colType,
    cusparseIndexBase_t idxBase,
    cudaDataType valueType
);

cusparseIndexType_t — 索引数据类型

• 用于描述 CSR/COO 等稀疏矩阵的 行列索引数组类型

cusparseCreateCsr 中 rowType 和 colType 都要用这个枚举。

cusparseIndexBase_t — 索引基

• 用于描述 CSR/COO 索引是从 0 还是从 1 开始

• 典型值：

示例：row_ptr 第一个元素是 0 或 1。

cudaDataType — 数据类型

• 用于描述稀疏矩阵或向量元素的数据类型

• 常见值：

创建稠密向量描述符

cusparseStatus_t cusparseCreateDnVec(
    cusparseDnVecDescr_t *vecDescr,
    int64_t size,
    void *values,
    cudaDataType valueType
);

查询 buffer 大小

cusparseStatus_t cusparseSpMV_bufferSize(
    cusparseHandle_t handle,
    cusparseOperation_t op,
    const void *alpha,
    cusparseSpMatDescr_t matA,
    cusparseDnVecDescr_t vecX,
    const void *beta,
    cusparseDnVecDescr_t vecY,
    cudaDataType computeType,
    cusparseSpMVAlg_t alg,
    size_t *bufferSize
);

作用：先查询 buffer 大小（查询的时候alpha和beta指针不可以是空的），方便在 GPU 上分配临时存储。

cusparseOperation_t — 矩阵操作

• 用于指定稀疏矩阵是否转置

• 常用值：

cusparseSpMVAlg_t — SpMV 算法类型

• 用于选择 稀疏矩阵-向量乘法的执行算法

• 常用值：

对小矩阵一般用 DEFAULT 就够；大矩阵可根据性能需求选特定算法。

执行 SpMV

cusparseStatus_t cusparseSpMV(
    cusparseHandle_t handle,
    cusparseOperation_t op,
    const void *alpha,
    cusparseSpMatDescr_t matA,
    cusparseDnVecDescr_t vecX,
    const void *beta,
    cusparseDnVecDescr_t vecY,
    cudaDataType computeType,
    cusparseSpMVAlg_t alg,
    void *externalBuffer
);

相当于执行：

销毁描述符

cusparseDestroySpMat(cusparseSpMatDescr_t matA);
cusparseDestroyDnVec(cusparseDnVecDescr_t vec);

释放 GPU 内存描述符资源

使用流程总结

• cusparseCreate → 初始化 handle

• cusparseCreateCsr → 创建 CSR 矩阵描述符

• cusparseCreateDnVec → 创建输入/输出向量描述符

• cusparseSpMV_bufferSize → 查询 buffer 大小

• cudaMalloc → 分配 buffer

• cusparseSpMV → 执行 SpMV

• cusparseDestroySpMat / cusparseDestroyDnVec → 清理描述符

• cusparseDestroy → 清理 handle

完整示例

#ifndef __CUDACC__
#define __CUDACC__
#endif
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cublas_v2.h>
#include<cub/cub.cuh>
#include<cusparse_v2.h>
#include <iostream>
#include<cstdio>
#include <vector>


void error_handling(cudaError_t res) {
    if (res !=cudaSuccess) {
        std::cout << "error!" << std::endl;
    }
}

int main() {
    // -----------------------
    // 1. 初始化数据
    // -----------------------
    // 矩阵 A (3x4)
    // [10 0 0 0
    //  0 20 0 30
    //  0 0 40 0]
    int rows = 3, cols = 4, nnz = 4;
    int h_csrRowPtr[] = { 0, 1, 3, 4 };
    int h_csrColInd[] = { 0, 1, 3, 2 };
    float h_csrValues[] = { 10.f, 20.f, 30.f, 40.f };
    float h_x[] = { 1.f, 2.f, 3.f, 4.f };
    float h_y[3] = { 0.f, 0.f, 0.f };
    // -----------------------
    // 2. 申请 GPU 内存
    // -----------------------
    int* d_csrRowPtr, * d_csrColInd;
    float* d_csrValues, * d_x, * d_y;
    cudaMalloc(&d_csrRowPtr, (rows + 1) * sizeof(int));
    cudaMalloc(&d_csrColInd, nnz * sizeof(int));
    cudaMalloc(&d_csrValues, nnz * sizeof(float));
    cudaMalloc(&d_x, cols * sizeof(float));
    cudaMalloc(&d_y, rows * sizeof(float));

    cudaMemcpy(d_csrRowPtr, h_csrRowPtr, (rows + 1) * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_csrColInd, h_csrColInd, nnz * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_csrValues, h_csrValues, nnz * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_x, h_x, cols * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_y, h_y, rows * sizeof(float), cudaMemcpyHostToDevice);

    // -----------------------
    // 3. 初始化 cuSPARSE handle
    // -----------------------
    cusparseHandle_t handle;
    cusparseCreate(&handle);

    // -----------------------
    // 4. 创建 CSR 矩阵描述符
    // -----------------------
    cusparseSpMatDescr_t matA;
    cusparseCreateCsr(&matA,
                        rows,cols,nnz,
                        d_csrRowPtr,d_csrColInd,d_csrValues,
                        CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I,
                        CUSPARSE_INDEX_BASE_ZERO,CUDA_R_32F);
    // -----------------------
    // 5. 创建向量描述符
    // -----------------------
    cusparseDnVecDescr_t vecX, vecY;
    cusparseCreateDnVec(&vecX, cols, d_x, CUDA_R_32F);
    cusparseCreateDnVec(&vecY, rows, d_y, CUDA_R_32F);

    // -----------------------
    // 6. 查询 buffer 大小并分配
    // -----------------------
    float alpha = 1.0f, beta = 0.0f;
    size_t bufferSize = 0;
    cusparseSpMV_bufferSize(handle,
        CUSPARSE_OPERATION_NON_TRANSPOSE,
        &alpha, matA, vecX, &beta,  vecY,
        CUDA_R_32F, CUSPARSE_SPMV_ALG_DEFAULT, &bufferSize);

    void* dBuffer;
    cudaMalloc(&dBuffer, bufferSize);
    // -----------------------
   // 7. 执行 SpMV
   // -----------------------

    cusparseSpMV(handle, 
        CUSPARSE_OPERATION_NON_TRANSPOSE,
        &alpha,matA,vecX,&beta,vecY,
        CUDA_R_32F,CUSPARSE_SPMV_ALG_DEFAULT,dBuffer);
    // -----------------------
   // 8. 拷贝结果回 CPU 并打印
   // -----------------------
    cudaMemcpy(h_y, d_y, rows * sizeof(float), cudaMemcpyDeviceToHost);
    std::cout << "y = [ ";
    for (int i = 0; i < rows; i++) std::cout << h_y[i] << " ";
    std::cout << "]" << std::endl;

    // -----------------------
    // 9. 清理
    // -----------------------
    cusparseDestroySpMat(matA);
    cusparseDestroyDnVec(vecX);
    cusparseDestroyDnVec(vecY);
    cusparseDestroy(handle);
    cudaFree(d_csrRowPtr);
    cudaFree(d_csrColInd);
    cudaFree(d_csrValues);
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(dBuffer);

    return 0;
}

输出：

y = [ 10 160 120 ]