优化稀疏矩阵向量乘提升不可压缩管流模拟效率的方法

技术领域

本发明属于流体力学数值模拟，尤其涉及一种优化稀疏矩阵向量乘提升不可压缩管流模拟效率的方法。

背景技术

不可压缩管流问题作为流体力学的重要研究对象，对此进行的研究成果被广泛应用于等离子体物理学、磁流体力学和天体物理学等相关领域的科学研究和工业技术中。随着计算机科学技术的发展，需经过现场试验才能解决的不可压缩管流问题已可以由计算机数值模拟方法解决。不可压缩管流模拟需求解大型线性方程组，即Ax=b形式(A为矩阵，x、b为向量)，而求解方程组的矩阵A往往是稀疏矩阵，其中非零元素数目占矩阵总容量不足5%。因此可以通过优化稀疏向量矩阵乘过程来提升不可压缩管流模拟的效率。

稀疏矩阵向量乘运算因矩阵稀疏模式与存储格式，矩阵元素缺乏时间局部性，非零元素的间接内存引用和向量的不规则访问等因素导致计算性能并不高。当前研究主要从降低时间复杂度和计算机体系结构两方面改进稀疏矩阵向量乘运算的效果。降低稀疏矩阵向量乘的时间复杂度的优化方法在大规模矩阵运算过程中因浮点运算操作会导致误差累积，主要适用于对精度要求低的领域。针对计算机处理器体系结构特点进行稀疏矩阵向量乘优化的方法主要从数据结构的改进、底层指令优化、并行优化技术和自动调优等四方面提高了运算效率。

这些方法受矩阵非零元素分布特性的影响大，并不具有一般性。对于总体稀疏但局部存在较多稠密子矩阵特点的矩阵，基于传统稀疏矩阵存储格式的稀疏矩阵向量乘运算的空间和时间局部性较差，为此进行的优化也会引入额外的开销，不可压缩管流的模拟效率仍有待提升。

发明内容

针对上述方法存在的问题，本发明提出一种优化稀疏矩阵向量乘提升不可压缩管流模拟效率的方法，通过结合四叉树结构和CSR存储结构的组合优势对稀疏矩阵进行递归式分解和重排列实现稀疏矩阵的存储，并在CPU/GPU异构并行系统上实现基于QCSR存储结构的高效稀疏矩阵向量乘方法。本发明所公开的方法可在稀疏矩阵向量乘数值计算过程中更有效提高数据局部性和缓存命中率，并且更具有普适性，取得了更好的计算加速和整体加速效果，进而提升了不可压缩管流的效率。

本发明方法的具体步骤如下：

步骤1、读取方程组，利用QCSR存储结构存储稀疏矩阵。首先按照行主序依次读入稀疏矩阵的非零元素，然后按照分块方式对矩阵进行行平均划分，并根据列主序对每个划分块内的元素重新排序存储，在当前划分块区间内，以第一个未在块区间内的非零元素所在列为起始列，按分块方式对矩阵进行列平均划分，如此递归进行划分。最后，当子矩阵块大小小于设定的分块大小时，递归划分矩阵结束，得到叶子结点分块矩阵，依次保存各叶子结点边界内的元素信息即可完成存储。此时叶子结点分块矩阵的大小与Cache容量相适应。

步骤2、将转换后的稀疏矩阵数据传输到GPU。采用分页锁定的传输方式申请一块写结合类型的内存，再将稀疏矩阵的非零元素拷贝到全局存储器，降低了数据传输的时间。

步骤3、数据存取优化。全局存储器访问方式采用硬件优化的合并访问模式，采用线程块中连续的线程访问连续元素，计算出非零元素的部分乘积。为了满足合并访问模式的要求，对稀疏矩阵进行数组尾部补零对齐，使每行元素个数为16的倍数，起始全局存储器地址对齐到访问数据长度的16倍。

步骤4、线程映射优化。一个线程块对应一个行分块，线程块中的线程对应一个叶子节点分块的一行，为每个矩阵行分配一个线程串完成矩阵中相对应的一行非零元素的计算。为了降低稀疏矩阵非零元素分布的影响，在已经进行补零对齐的基础上，各线程块先计算含有叶子节点分块数目较多的行分块。最后将各个线程计算的结果并行求和规约得到最终的计算结果。

步骤5、数据复用优化。利用共享存储器存储中间结果元素和行索引，利用纹理存储器存储输入向量。

步骤6、执行稀疏矩阵向量乘数值计算。利用GPU的编程模型CUDA实现稀疏矩阵向量乘的并行运算。

步骤7、将计算结果传输到CPU。通过并行运算得到稀疏矩阵向量乘的结果，并将结果向量传输到CPU。

本发明提供的优化稀疏矩阵向量乘提升不可压缩管流模拟效率的方法，具有以下优势：

1、采用高度压缩及灵活方便的QCSR存储结构，可在稀疏矩阵向量乘数值计算过程中更有效提高数据局部性和缓存命中率，同时通过对矩阵存储的重排序，降低了运算过程中矩阵内非零元素分布的影响因素，更具有普适性。

2、采用CPU/GPU协同运算方法，并结合线程映射优化、数据存取优化、数据传输优化和数据复用优化四个策略，提高了稀疏矩阵向量乘的加速效率，进而提升了不可压缩管流的模拟效率。

附图说明

图1为本发明流程图。

具体实施方式

下面结合附图和实施方法对本发明作进一步的详细说明。

如图1所示，本发明方法包括以下步骤：

步骤1、读取方程组，利用QCSR存储结构存储稀疏矩阵。首先按照行主序依次读入稀疏矩阵的非零元素，然后按照分块方式对矩阵进行行平均划分，并根据列主序对每个划分块内的元素重新排序存储，在当前划分块区间内，以第一个未在块区间内的非零元素所在列为起始列，按分块方式对矩阵进行列平均划分，如此递归进行划分。最后，当子矩阵块大小小于设定的分块大小时，递归划分矩阵结束，得到叶子结点分块矩阵，依次保存各叶子结点边界内的元素信息即可完成存储。此时叶子结点分块矩阵的大小与Cache容量相适应。

步骤2、将转换后的稀疏矩阵数据传输到GPU。采用分页锁定的传输方式申请一块写结合类型的内存，再将稀疏矩阵的非零元素拷贝到全局存储器，降低了数据传输的时间。

步骤3、数据存取优化。全局存储器访问方式采用硬件优化的合并访问模式，采用线程块中连续的线程访问连续元素，计算出非零元素的部分乘积。为了满足合并访问模式的要求，对稀疏矩阵进行数组尾部补零对齐，使每行元素个数为16的倍数，起始全局存储器地址对齐到访问数据长度的16倍。

步骤4、线程映射优化。一个线程块对应一个行分块，线程块中的线程对应一个叶子节点分块的一行，为每个矩阵行分配一个线程串完成矩阵中相对应的一行非零元素的计算。为了降低稀疏矩阵非零元素分布的影响，在已经进行补零对齐的基础上，各线程块先计算含有叶子节点分块数目较多的行分块。最后将各个线程计算的结果并行求和规约得到最终的计算结果。

步骤5、数据复用优化。利用共享存储器存储中间结果元素和行索引，利用纹理存储器存储输入向量。

步骤6、执行稀疏矩阵向量乘数值计算。利用GPU的编程模型CUDA实现稀疏矩阵向量乘的并行运算。

步骤7、将计算结果传输到CPU。通过并行运算得到稀疏矩阵向量乘的结果，并将结果向量传输到CPU。