CN110598174A - 一种基于gpu架构的稀疏矩阵的回代求解方法 - Google Patents

一种基于gpu架构的稀疏矩阵的回代求解方法 Download PDF

Info

Publication number
CN110598174A
CN110598174A CN201910859642.XA CN201910859642A CN110598174A CN 110598174 A CN110598174 A CN 110598174A CN 201910859642 A CN201910859642 A CN 201910859642A CN 110598174 A CN110598174 A CN 110598174A
Authority
CN
China
Prior art keywords
matrix
nodes
sparse
substitution
back substitution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201910859642.XA
Other languages
English (en)
Other versions
CN110598174B (zh
Inventor
王晓光
程明厚
邵雪
周振亚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing CEC Huada Electronic Design Co Ltd
Original Assignee
Beijing CEC Huada Electronic Design Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing CEC Huada Electronic Design Co Ltd filed Critical Beijing CEC Huada Electronic Design Co Ltd
Priority to CN201910859642.XA priority Critical patent/CN110598174B/zh
Publication of CN110598174A publication Critical patent/CN110598174A/zh
Application granted granted Critical
Publication of CN110598174B publication Critical patent/CN110598174B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Complex Calculations (AREA)

Abstract

本发明提供一种基于GPU架构的稀疏矩阵的回代求解方法，包括以下步骤：1）计算稀疏矩阵的L矩阵和U矩阵；2）分别建立L矩阵和U矩阵的合并压缩矩阵；3）建立节点依赖关系图；4）建立节点依赖关系树；5）分层回代求解。本发明根据不同GPU的型号，确定合并矩阵的维数，可以降低矩阵维数及依赖关系树的深度。根据矩阵的规模，当前层要处理的节点数量以及稠密度等信息，选择不同的kernel函数，充分利用GPU显卡的资源，更快的完成每层节点的回代求解。

Description

一种基于GPU架构的稀疏矩阵的回代求解方法
技术领域
本发明涉及集成电路计算机辅助设计(Integrated Circuit/Computer AidedDesign)领域，尤其是EDA电路仿真技术领域，特别涉及一种基于GPU架构的稀疏矩阵的回代求解方法。
背景技术
随着微电子技术及先进制造工艺的发展，IC设计规模越来越庞大，设计中的晶体管数量急剧增加，这对于晶体管级仿真工具来说，器件的计算量和非线性方程组的维数的增加会导致完成一次迭代的时间呈超线性增长。对于在传统CPU上开发的仿真工具来说，完成一个电路设计的仿真可能需要几个月的时间。近年来NVIDIA的GPU在通用计算方面有了快速的发展，在稠密矩阵运算方面得到广泛应用，如机器学习、深度学习等方面。而工业应用中，受其GPU硬件限制，访问非连续的数据的速度要比CPU慢很多，因此GPU在稀疏矩阵的运算方面，应用的还比较少。
仿真工具主要解决的是以稀疏矩阵格式存储的非线性方程组求解的问题，目前市场上的商用仿真工具中，还没有基于GPU架构实现仿真加速的工具。本发明解决的问题是基于GPU架构的仿真加速流程中，稀疏矩阵的回代求解的加速方法。
发明内容
为了解决现有技术存在的不足，本发明的目的在于提供一种基于GPU架构的稀疏矩阵的回代求解方法，以解决稀疏矩阵的回代求解在GPU上加速的实现方法，避免频繁在CPU与GPU之间进行数据传输，将仿真流程中的多个数据处理步骤都集成在GPU架构下，保证各个流程中数据的连续性。
为实现上述目的，本发明提供的基于GPU架构的稀疏矩阵的回代求解方法，包括以下步骤：
1)对稀疏矩阵进行LU分解，得到L矩阵和U矩阵；
2)分别建立L矩阵和U矩阵的合并压缩矩阵；
3)根据所述合并压缩矩阵，建立节点依赖关系图；
4)从所述依赖关系图的叶子节点出发，保留节点最深度数，建立节点依赖关系树，并标记每个节点的深度；
5)分层回代求解。
进一步地，所述步骤2)进一步包括以下步骤：
21)确定合并压缩的维数，按合并压缩的维数进行矩阵压缩；
22)对压缩后的矩阵按照线程访问方式进行重新排列；
23)将排列后的矩阵存储在GPU的一个连续内存空间上。
进一步地，所述步骤3)将所述合并压缩矩阵的每一行或每一列作为节点，建立节点依赖关系图。
进一步地，包括，分别将压缩后的L矩阵和U矩阵的节点的非零元分为两个集合，第一集合包括对角位置的非零元，第二集合包括非对角位置的非零元。
进一步地，所述步骤5)进一步包括以下步骤：
51)计算对角位置的非零元所对应的回代结果；
52)利用所述回代结果，计算所有非对角位置的非零元；
53)将计算的数据累加至还未计算的更深层的回代结果向量。
进一步地，包括，在对角位置的回代结果确定后，将所有的非对角位置的数据更新操作全部并发执行，使用CUDA的原子操作，将数据累加到还未计算的回代结果向量中。
进一步地，包括，对于压缩后的L矩阵和压缩后的U矩阵的深度相同的节点，根据原始矩阵的规模、当前层内节点的数量和稀疏矩阵的稠密度，调用不同的kernel函数，完成该层所有节点的回代求解。
更进一步地，按照深度值从小到大，先后对L矩阵和U矩阵的每一层节点进行回代求解，最后得到稀疏矩阵的回代求解的结果向量。
为实现上述目的，本发明还提供一种计算机可读存储介质，其上存储有计算机指令，所述计算机指令运行时执行上述的基于GPU架构的稀疏矩阵的回代求解方法的步骤。
有益效果：本发明根据不同GPU的型号，确定合并矩阵的维数，可以降低矩阵维数及依赖关系树的深度。根据矩阵的规模，当前层要处理的节点数量以及稠密度等信息，选择不同的kernel函数，充分利用GPU显卡的资源，更快的完成每层节点的回代求解。
本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。
附图说明
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，并与本发明的实施例一起，用于解释本发明，并不构成对本发明的限制。在附图中：
图1为根据本发明的基于GPU架构的稀疏矩阵的回代求解方法的流程图；
图2为根据本发明的实施方式的原始16维稀疏矩阵对应的LU分解后的下三角L矩阵；
图3为根据本发明的实施方式的原始16维稀疏矩阵对应的LU分解后的上三角U矩阵；
图4为将图2中的下三角L矩阵以2为维数进行合并压缩后得到的压缩L矩阵；
图5为将图3中的上三角U矩阵以2为维数进行合并压缩后得到的压缩U矩阵；
图6为按照图4中的压缩L矩阵回代求解过程建立的依赖关系图；
图7为按照图5中的压缩U矩阵回代求解过程建立的依赖关系图；
图8为根据图6的依赖关系图构造的依赖树，并按照深度进行分层；
图9为根据图7的依赖关系图构造的依赖树，并按照深度进行分层；
图10为按照图8的分层关系，结合图4中的非零元，压缩L矩阵的回代流程图；
图11为按照图9的分层关系，结合图5中的非零元，压缩U矩阵的回代流程图；
图12为根据本发明的实施方式的在GPU上实现的稀疏矩阵的回代求解的流程图。
具体实施方式
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
图1为根据本发明的基于GPU架构的稀疏矩阵的回代求解方法的流程图，下面将参考图1，对本发明的基于GPU架构的稀疏矩阵的回代求解方法进行详细描述。
在步骤101，计算稀疏矩阵的L矩阵和U矩阵。
在该步骤中，选取的原始稀疏矩阵为一个16维的稀疏矩阵。利用LU分解，将原始16维稀疏矩阵分解为下三角L矩阵和上三角U矩阵，分别如图2和图3所示，其中，图2为根据本发明的实施方式的原始16维稀疏矩阵对应的LU分解后的下三角L矩阵，图3为根据本发明的实施方式的原始16维稀疏矩阵对应的LU分解后的上三角U矩阵。在图中，*表示非零元，0是稀疏位置，在矩阵存储时0位置的数据不需要存储。
在步骤102，建立L矩阵和U矩阵的合并压缩矩阵。
在该步骤中，对L矩阵和U矩阵进行合并压缩。具体地，根据GPU显卡硬件情况，选择合并压缩维数，以适应将压缩后的每个点对应的矩阵可以在GPU的一个block运行。本发明是在NVIDIA的GPU显卡上完成的对稀疏矩阵的回代求解方法，这里以维数2为例，对步骤101中所得的L矩阵和U矩阵分别进行压缩，压缩结果分别如图4、图5所示，其中，图4为将图2中的下三角L矩阵以2为维数进行合并压缩后得到的压缩L矩阵，图5为将图3中的上三角U矩阵以2为维数进行合并压缩后得到的压缩U矩阵。
在该步骤中进行矩阵的压缩合并，是为了提高GPU多线程的访问内存效率，因而需要重新将合并的多行与多列的稀疏矩阵的元素(数据)按照线程访问方式进行重新排列，并存储在GPU的一个连续内存空间上。
在步骤103，建立节点依赖关系图。
在该步骤中，针对压缩后的L矩阵和U矩阵，将合并压缩后的矩阵的每一行或者每一列作为节点，分别建立回代求解过程的L矩阵和U矩阵的依赖关系图，分别如图6、图7所示，其中，图6为按照图4中的压缩L矩阵回代求解过程建立的依赖关系图，图7为按照图5中的压缩U矩阵回代求解过程建立的依赖关系图。
在步骤104，建立节点依赖关系树。
在该步骤中，根据步骤103中获得的图6和图7所示的依赖关系图分别建立对应于压缩后的L矩阵和U矩阵的依赖树，并标记每个节点的深度。具体地，从叶子节点出发，只保留每个节点最深度对应的边，将得到依赖关系树，同时从0度节点出发，将具有相同深度的节点划分到同一层中，这样同时得到了L矩阵和U矩阵的回代求解过程的分层图，依赖树及分层结果分别如图8和图9所示，其中，图8为根据图6的依赖关系图构造的依赖树，并按照深度进行分层，图9为根据图7的依赖关系图构造的依赖树，并按照深度进行分层。
将合并后的矩阵的每一行或者每一列作为节点，建立依赖关系图和依赖关系树，有效降低了树中的节点数量以及树的深度。
在该步骤中，分别将每层中的L矩阵和U矩阵的节点的非零元分为两个集合，第一集合为处在对角位置的非零元的集合，第二集合为非对角位置的非零元的集合。
在步骤105，分层回代求解。
在该步骤中，对于L矩阵和U矩阵的深度相同的节点，根据原始矩阵的规模，以及当前层内节点的数量和稀疏矩阵的稠密度等信息，调用不同的kernel函数完成这一层所有节点的回代求解。
同时，为了尽量将无依赖的操作并行化，将每一个压缩矩阵的回代求解操作分为两步：第一步：计算对角位置的非零元所对应的回代结果；第二步：使用计算后的回代结果，用于计算所有非对角位置的非零元的数据，并将计算的数据累加到还未计算的更深层的回代结果向量中。其中，在对角位置的回代结果确定后，将所有的非对角位置的数据更新操作全部并发执行，使用CUDA的原子操作，将数据累加到还未计算的回代结果向量中。
现在结合图10和图11对步骤105的分层回代求解进行说明，图10为按照图8的分层关系，结合图4中的非零元，压缩L矩阵的回代流程图，图11为按照图9的分层关系，结合图5中的非零元，压缩U矩阵的回代流程图。具体地，给定输入向量b，结果输出到向量y中。第一步，执行Sub_L/Sub_U操作，计算对角矩阵对应位置的y结果。第二步，执行Sub_Update操作，将使用第一步中对应的y结果，更新部分数据到还未处理的更深层中的y的结果中，程序的具体执行顺序分别如图10和图11所示。
以L矩阵的第一层为例，在Sub_L计算确认y1，y3和y5之后，将部分结果更新到y2，y4，y6和y8。其中y2，y4，y6和y8的计算公式分别如公式(1)-(4)所示：
y2＝(b2-L21*y1)/L22 (1)
y4＝(b4-L43*y3)/L44 (2)
y6＝(b6-L65*y5)/L66 (3)
y8＝(b8-L81*y1-L82*y2-L84*y4-L86*y6-L87*y7)/L88 (4)
Sub_Update将使用CUDA的原子操作，将L21*y1，L43*y3，L65*y5和L81*y1分别更新到y2，y4，y6和y8中。
其中，Sub_L，Sub_U和Sub_Update三个kernel函数的实现有多重方法，需要根据稀疏矩阵的维数、当前层中的节点数量以及合并矩阵的稀疏度等因素，使用不同的实现方式，例如，将b使用shared memory存储，将L矩阵和U矩阵的值使用shared memory存储，使用更多线程处理一行(一列)或者非零元超长行(列)等方法，这样可以充分利用GPU显卡的资源，达到最快的加速效果。
在回代求解过程中，先执行L矩阵的所有层的回代操作，再执行U矩阵的所有层的回代操作，最后的向量y为回代求解的结果，完整的回代求解流程图如图12所示。在执行L矩阵或U矩阵的所有层的回代操作时，按照深度值，从小到大的顺序依次执行。
与CPU版本的16线程并行的回代求解方法相比，GPU的回代求解可以达到5到10倍甚至更高倍数的加速。
本发明还提供了一种计算机可读存储介质，其上存储有计算机指令，所述计算机指令运行时执行上述的基于GPU架构的稀疏矩阵的回代求解方法的步骤，所述基于GPU架构的稀疏矩阵的回代求解方法参见前述部分的介绍，不再赘述。
本领域普通技术人员可以理解：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于GPU架构的稀疏矩阵的回代求解方法，其特征在于，包括以下步骤：
1）对稀疏矩阵进行LU分解，得到L矩阵和U矩阵；
2）分别建立L矩阵和U矩阵的合并压缩矩阵；
3）根据所述合并压缩矩阵，建立节点依赖关系图；
4）从所述依赖关系图的叶子节点出发，保留节点最深度数，建立节点依赖关系树，并标记每个节点的深度；
5）分层回代求解。
2.根据权利要求1所述的基于GPU架构的稀疏矩阵的回代求解方法，其特征在于，所述步骤2）进一步包括以下步骤：
21）确定合并压缩的维数，按合并压缩的维数进行矩阵压缩；
22）对压缩后的矩阵按照线程访问方式进行重新排列；
23）将排列后的矩阵存储在GPU的一个连续内存空间上。
3.根据权利要求1所述的基于GPU架构的稀疏矩阵的回代求解方法，其特征在于，所述步骤3）将所述合并压缩矩阵的每一行或每一列作为节点，建立节点依赖关系图。
4.根据权利要求1所述的基于GPU架构的稀疏矩阵的回代求解方法，其特征在于，进一步包括，分别将压缩后的L矩阵和U矩阵的节点的非零元分为两个集合，第一集合包括对角位置的非零元，第二集合包括非对角位置的非零元。
5.根据权利要求1所述的基于GPU架构的稀疏矩阵的回代求解方法，其特征在于，所述步骤5）进一步包括以下步骤：
51）计算对角位置的非零元所对应的回代结果；
52）利用所述回代结果，计算所有非对角位置的非零元；
53）将计算的数据累加至还未计算的更深层的回代结果向量。
6.根据权利要求5所述的基于GPU架构的稀疏矩阵的回代求解方法，其特征在于，进一步包括，在对角位置的回代结果确定后，将所有的非对角位置的数据更新操作全部并发执行，使用CUDA的原子操作，将数据累加到还未计算的回代结果向量中。
7.根据权利要求5所述的基于GPU架构的稀疏矩阵的回代求解方法，其特征在于，进一步包括，对于压缩后的L矩阵和压缩后的U矩阵的深度相同的节点，根据原始矩阵的规模、当前层内节点的数量和稀疏矩阵的稠密度，调用不同的kernel函数，完成该层所有节点的回代求解。
8.根据权利要求1所述的基于GPU架构的稀疏矩阵的回代求解方法，其特征在于，所述步骤5）进一步包括，按照深度值从小到大，先后对L矩阵和U矩阵的每一层节点进行回代求解，最后得到稀疏矩阵的回代求解的结果向量。
9.一种计算机可读存储介质，其上存储有计算机指令，其特征在于，所述计算机指令运行时执行权利要求1至8任一项所述的基于GPU架构的稀疏矩阵的回代求解方法的步骤。
CN201910859642.XA 2019-09-11 2019-09-11 一种基于gpu架构的稀疏矩阵的回代求解方法 Active CN110598174B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910859642.XA CN110598174B (zh) 2019-09-11 2019-09-11 一种基于gpu架构的稀疏矩阵的回代求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910859642.XA CN110598174B (zh) 2019-09-11 2019-09-11 一种基于gpu架构的稀疏矩阵的回代求解方法

Publications (2)

Publication Number Publication Date
CN110598174A true CN110598174A (zh) 2019-12-20
CN110598174B CN110598174B (zh) 2022-06-21

Family

ID=68858918

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910859642.XA Active CN110598174B (zh) 2019-09-11 2019-09-11 一种基于gpu架构的稀疏矩阵的回代求解方法

Country Status (1)

Country Link
CN (1) CN110598174B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117311948A (zh) * 2023-11-27 2023-12-29 湖南迈曦软件有限责任公司 Cpu与gpu异构并行的自动多重子结构数据处理方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103309845A (zh) * 2013-05-16 2013-09-18 国家电网公司 一种用于电力系统动态仿真的线性方程组分块求解方法
CN103399841A (zh) * 2013-07-31 2013-11-20 清华大学 基于gpu的稀疏矩阵lu分解方法
CN108984483A (zh) * 2018-07-13 2018-12-11 清华大学 基于dag及矩阵重排的电力系统稀疏矩阵求解方法和系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103309845A (zh) * 2013-05-16 2013-09-18 国家电网公司 一种用于电力系统动态仿真的线性方程组分块求解方法
CN103399841A (zh) * 2013-07-31 2013-11-20 清华大学 基于gpu的稀疏矩阵lu分解方法
CN108984483A (zh) * 2018-07-13 2018-12-11 清华大学 基于dag及矩阵重排的电力系统稀疏矩阵求解方法和系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117311948A (zh) * 2023-11-27 2023-12-29 湖南迈曦软件有限责任公司 Cpu与gpu异构并行的自动多重子结构数据处理方法
CN117311948B (zh) * 2023-11-27 2024-03-19 湖南迈曦软件有限责任公司 Cpu与gpu异构并行的自动多重子结构数据处理方法

Also Published As

Publication number Publication date
CN110598174B (zh) 2022-06-21

Similar Documents

Publication Publication Date Title
CN110826719B (zh) 一种量子程序的处理方法、装置、存储介质和电子装置
CN102323922B (zh) 数值数据的压缩和解压缩
CN111667051A (zh) 适用边缘设备的神经网络加速器及神经网络加速计算方法
CN110516810B (zh) 一种量子程序的处理方法、装置、存储介质和电子装置
Peng et al. GLU3. 0: Fast GPU-based parallel sparse LU factorization for circuit simulation
US20200090051A1 (en) Optimization problem operation method and apparatus
WO2024027039A1 (zh) 数据处理方法、装置、设备和可读存储介质
CN116150553B (zh) 一种面向cpu+dcu异构混合架构的稀疏化amg优化方法
CN110737870B (zh) 一种用于在gpu上合并舒尔矩阵的方法及装置
CN110598174B (zh) 一种基于gpu架构的稀疏矩阵的回代求解方法
CN111640296A (zh) 交通流预测方法、系统、存储介质及终端
CN115859016B (zh) 基于处理器的运算方法、装置、计算机设备及存储介质
CN115480919A (zh) 卷积优化运算方法、装置、计算机设备及存储介质
CN115859011A (zh) 矩阵运算方法、装置及单元、电子设备
CN115293978A (zh) 卷积运算电路和方法、图像处理设备
CN114662648A (zh) 压缩系数集以供随后在神经网络中使用
Chen et al. Sparsity-oriented sparse solver design for circuit simulation
CN113312866B (zh) 一种通过合并方程变量实现约减电路的方法
Ahn et al. Common kernels and convolutions in binary-and ternary-weight neural networks
US9633147B1 (en) Power state coverage metric and method for estimating the same
CN117149446B (zh) 一种基于图形处理器的数据处理方法、装置、设备及介质
CN116805155B (zh) 一种lstm网络处理方法、装置、设备及可读存储介质
CN116301903B (zh) 一种编译器、ai网络编译方法、处理方法、执行系统
CN111274109B (zh) 一种基于请求处理模拟的系统软硬件拓扑的评估方法及系统
CN116578425B (zh) 一种基于光栅化的负载均衡方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information

Address after: 100102 floor 2, block a, No.2, lizezhong 2nd Road, Chaoyang District, Beijing

Applicant after: Beijing Huada Jiutian Technology Co.,Ltd.

Address before: 100102 floor 2, block a, No.2, lizezhong 2nd Road, Chaoyang District, Beijing

Applicant before: HUADA EMPYREAN SOFTWARE Co.,Ltd.

CB02 Change of applicant information
GR01 Patent grant
GR01 Patent grant