[CUDA 性能调优实战] 稀疏矩阵向量乘法(SpMV) 的存储格式选择与性能对比

1. 稀疏矩阵与SpMV基础概念稀疏矩阵向量乘法SpMV是科学计算中最基础也最重要的操作之一几乎出现在所有涉及线性代数的应用场景中。简单来说SpMV就是将一个稀疏矩阵A与一个稠密向量x相乘得到另一个稠密向量y。这个看似简单的操作背后却隐藏着巨大的性能优化空间。我第一次接触SpMV优化是在处理一个流体力学仿真项目时。当时使用的矩阵有超过100万行但每行平均只有不到10个非零元素。如果按照传统稠密矩阵的方式存储和计算不仅会浪费大量内存计算效率也极其低下。这就是稀疏矩阵的特殊之处——它们通常具有高度的结构性非零元素分布呈现特定模式。稀疏矩阵存储格式的核心思想很简单只存储非零元素及其位置信息。常见的格式包括CSR按行压缩存储适合行访问COO坐标格式简单直观但效率不高ELL固定每行元素数适合规整矩阵HYB混合格式结合ELL和COO优势在CUDA平台上实现SpMV时选择哪种存储格式会直接影响性能。我曾经测试过一个200万阶的矩阵使用不同格式时性能差异可以达到5倍以上。这主要是因为GPU的并行架构对内存访问模式极其敏感。2. CSR格式的CUDA实现与优化CSRCompressed Sparse Row是最常用的稀疏矩阵存储格式之一。它使用三个数组data存储非零元素值col_ind存储列索引row_ptr存储每行起始位置2.1 基础实现方案最简单的CSR SpMV实现是让每个线程处理一行__global__ void csr_spmv(int n_rows, const int* row_ptr, const int* col_ind, const float* data, const float* x, float* y) { int row blockIdx.x * blockDim.x threadIdx.x; if (row n_rows) { float sum 0; for (int j row_ptr[row]; j row_ptr[row1]; j) { sum data[j] * x[col_ind[j]]; } y[row] sum; } }这种实现虽然简单但存在严重的负载不均衡问题。我在测试一个社交网络关系矩阵时发现某些行有上千个非零元而大多数行只有个位数导致GPU利用率不到30%。2.2 向量化优化方案更高级的实现是让一个warp32线程处理一行__global__ void csr_vector_spmv(int n_rows, const int* row_ptr, const int* col_ind, const float* data, const float* x, float* y) { int tid blockIdx.x * blockDim.x threadIdx.x; int warp_id tid / 32; int lane tid % 32; if (warp_id n_rows) { int row warp_id; int row_start row_ptr[row]; int row_end row_ptr[row1]; float sum 0; for (int j row_start lane; j row_end; j 32) { sum data[j] * x[col_ind[j]]; } // Warp内归约 for (int offset 16; offset 0; offset / 2) { sum __shfl_down_sync(0xffffffff, sum, offset); } if (lane 0) y[row] sum; } }这种方案在我测试的生物信息学矩阵上获得了3倍的性能提升。关键在于它更好地利用了GPU的SIMT架构减少了线程闲置。3. ELL与HYB格式的实践对比当矩阵行长度相对均匀时ELL格式往往能提供最佳性能。ELL将矩阵存储为两个二维数组data非零元素值indices列索引3.1 ELL格式实现__global__ void ell_spmv(int n_rows, int max_nnz_per_row, const int* indices, const float* data, const float* x, float* y) { int row blockIdx.x * blockDim.x threadIdx.x; if (row n_rows) { float sum 0; for (int j 0; j max_nnz_per_row; j) { int idx row j * n_rows; int col indices[idx]; sum data[idx] * x[col]; } y[row] sum; } }我在一个有限元分析矩阵上测试时ELL比CSR快了近40%。但遇到某些行特别长的情况时ELL会因为填充过多零元素而导致性能下降。3.2 HYB混合格式方案HYB结合了ELL和COO的优点。基本思路是将大部分元素用ELL格式存储超长行的额外元素用COO格式存储实现时需要特别注意原子操作的性能影响__global__ void hybrid_spmv(int n_rows, int ell_elems, const int* ell_col, const float* ell_data, int coo_nnz, const int* coo_row, const int* coo_col, const float* coo_data, const float* x, float* y) { // ELL部分 int row blockIdx.x * blockDim.x threadIdx.x; if (row n_rows) { float sum 0; for (int j 0; j ell_elems; j) { int idx row j * n_rows; sum ell_data[idx] * x[ell_col[idx]]; } atomicAdd(y[row], sum); } // COO部分 int elem blockIdx.x * blockDim.x threadIdx.x; if (elem coo_nnz) { float val coo_data[elem] * x[coo_col[elem]]; atomicAdd(y[coo_row[elem]], val); } }在一个推荐系统用的稀疏矩阵上HYB比纯ELL快了约25%比CSR快了近60%。但原子操作确实带来了额外的开销需要根据具体矩阵特性权衡。4. 性能优化实战经验经过多个项目的实践我总结了以下几点关键经验4.1 矩阵特征分析在决定使用哪种格式前必须分析矩阵的稀疏模式计算非零元分布直方图检查是否存在明显的块结构评估行长度方差我曾经遇到一个矩阵99%的行长度在10-20之间但有几十行超过1000。这种情况下使用HYB并设置ELL部分长度为20获得了最佳性能。4.2 内存访问优化GPU性能瓶颈主要在内存访问上。几个实用技巧合并访问确保相邻线程访问相邻内存地址数据预取提前加载下一批需要的数据共享内存缓存频繁访问的数据例如在CSR实现中可以将行指针缓存到共享内存__shared__ int shared_row_ptr[BLOCK_SIZE1]; if (threadIdx.x BLOCK_SIZE) { shared_row_ptr[threadIdx.x] row_ptr[blockIdx.x * BLOCK_SIZE threadIdx.x]; } __syncthreads();4.3 自动格式选择策略对于通用库开发可以实现自动选择器def select_format(matrix): nnz_per_row matrix.nnz / matrix.shape[0] max_nnz max(len(row) for row in matrix.rows) if max_nnz 2 * nnz_per_row: return ELL elif nnz_per_row 10: return COO else: return HYB if max_nnz 5 * nnz_per_row else CSR在实际项目中我开发了一个动态分析工具可以自动测试不同格式的性能并选择最优方案。对于批处理的SpMV操作这个工具平均能提升35%的性能。5. 不同应用场景下的选择建议根据我的项目经验不同领域的最佳选择往往不同5.1 科学计算领域典型特征矩阵来自物理方程离散化通常具有规则的稀疏模式非零元分布相对均匀推荐方案优先尝试ELL格式如果存在少量超长行再考虑HYB。例如在计算流体力学中ELL通常表现最佳。5.2 社交网络分析典型特征矩阵表示用户关系图存在幂律分布少数节点有大量连接稀疏模式高度不规则推荐方案CSR-Vector或Merge-based方案。在处理Twitter社交图谱时Merge-based方案比传统CSR快了近3倍。5.3 推荐系统典型特征矩阵表示用户-物品交互通常非常稀疏密度0.1%行长度差异较大推荐方案HYB格式。在电商推荐系统项目中HYB比纯CSR节省了40%的计算时间。6. 现代GPU上的新特性利用新一代GPU架构如Ampere、Hopper提供了许多有用的特性6.1 张量核心加速虽然SpMV主要是内存受限型操作但对于某些块稀疏矩阵可以使用张量核心wmma::load_matrix_sync(a_frag, a_ptr, stride); wmma::load_matrix_sync(b_frag, b_ptr, stride); wmma::mma_sync(c_frag, a_frag, b_frag, c_frag);6.2 异步数据拷贝利用CUDA 11的异步拷贝减少内存延迟__shared__ extern float smem[]; cuda::memcpy_async(smem, global_ptr, size, thread_block()); cuda::wait(thread_block());6.3 持久线程块对于超大稀疏矩阵可以配置持久线程块提高GPU利用率cudaOccupancyMaxPersistentBlocks(num_blocks, kernel, block_size, smem_size);在A100上测试表明这些新技术可以带来额外的15-20%性能提升但需要特别注意不同架构的兼容性问题。我在移植一个HPC应用到新GPU平台时就遇到了张量核心使用不当导致的精度问题最终通过混合精度计算解决了这个问题。