一种大型稀疏矩阵加速求解方法、系统及存储介质与流程

1.本发明属于计算电磁学技术领域，更具体地说，涉及一种大型稀疏矩阵加速求解方法、系统及存储介质。


背景技术：

2.现代eda仿真工具，无论是电路分析，还是3d电磁场分析，都面临一个随着仿真规模越来越大，计算越来越慢的问题。对于电路问题，计算速度变慢来源于复杂电路中数量巨大的电压节点和分支电流节点；对于3d电磁场问题，如果是用有限元分析，计算速度变慢则来源于稠密的网格数量。数学上来说，庞大数量的电压和电流节点，或者是稠密的网格数量，都对应着一个超大维度的稀疏矩阵。如何在保证计算精度的情况下，快速有效的求解大规模稀疏矩阵，是现代eda行业所面临的共同难题。面对这个问题，传统的做法是通过编写更高效的稀疏矩阵求解器来完成，比如一些著名的矩阵求解软件umfpack，pardiso，klu等等。使用这些矩阵求解程序来求解大规模矩阵，其计算效率远远大于自己编写的矩阵求解程序。然而，这些通用的矩阵求解软件，并不擅长处理eda行业中面临的具体问题，譬如分析带有ibis模型的非线性电路，或者是分析电路板的dc电压等等。前者由于非线性的影响，使得计算过程中每一个时间步都要面临大矩阵求逆的问题，而后者可以通过把电源平面和地平面分开考虑，使得大矩阵可以严格退耦。上述问题都是传统矩阵求解器本身所不能解决的问题，需要新的思路和方法。
3.在eda行业中普遍使用的计算方法是利用3d电磁场理论以及有限元算法，来分析chip/package/board的信号/电源完整性以及电磁兼容问题是非常有效和成熟的手段，譬如ansys公司的hfss，cadance的clerity3d等，这些工具就是使用有限元对3d电磁场进行计算和分析的。但是随着信号工作频率越来越高，芯片尺度越来越小，芯片内部结构越来越复杂。商用eda仿真工具在处理这类问题时，由于3d建模对精细网格的要求，使得计算复杂度越来越大，网格所对应的未知量的个数从几万，几十万上升到几百万不等，这种情况下，矩阵求解变得更加缓慢。一个设计，工程师通常需要等待几天甚至几周才能得到一个计算结果，耗时严重，对于eda设计人员以及芯片设计和生产企业更是非常严重的问题。
4.针对上述问题也进行了相应的改进，如中国专利申请号cn202110024925.x，公开日为2021年4月30日，该专利公开了一种稀疏矩阵加速计算方法，先读取第一读取待相乘的第一稀疏矩阵，对第一稀疏矩阵进行非零检测，并根据检测结果生成第一稀疏矩阵每行数据的第一状态信息并存储至寄存器；以及将检测到的第一稀疏矩阵的非零数据存储至ram；然后再读取待相乘的第二稀疏矩阵，对第二稀疏矩阵进行非零检测，并根据检测结果生成第二稀疏矩阵每列数据的第二状态信息并存储至寄存器；最后对第一状态信息和第二状态信息进行逻辑运算，根据逻辑运算结果读取ram中的数据并与第二稀疏矩阵的数据进行乘积运算以得到乘积矩阵的数据。该专利的不足之处在于：虽减少计算过程中数据的读取量从而加快处理速度，但是处理精度无法较好的满足。


技术实现要素：

5.1、要解决的问题
6.针对现有大型稀疏矩阵求解效率慢且精度差的问题，本发明提供一种大型稀疏矩阵加速求解方法、系统及存储介质。本发明的方法通过还原端口节点连接关系保证后续求解精度，通过将矩阵转化为无定向图之后进行分解得到最优分解，按照最优分解进行后续操作极大的提高了矩阵求解速度。本发明的系统结构简单，在满足精度的同时满足处理效率，使得精度与效率得到了较好的平衡。
7.2、技术方案
8.为解决上述问题，本发明采用如下的技术方案。
9.一种大型稀疏矩阵加速求解方法，包括以下步骤：
10.s1：对初始有限元矩阵进行还原端口节点的连接关系，得到还原后的二次有限元矩阵；
11.s2：将二次有限元矩阵转化为无定向图；
12.s3：对无定向图进行分解，并且通过评价函数选择最优分解；
13.s4：根据最优分解对二次有限元矩阵的节点重新编号；
14.s5：根据节点的重新编号对二次有限元矩阵重新排序，生成新的最终有限元矩阵；
15.s6：对最终有限元矩阵进行求解。
16.更进一步的，所述步骤s2中具体的步骤为：
17.s21：定义两个新的数据结构：点和边，点：非零元素所在列即为该行所对应的节点，该行即为非零元素所在的行；边：非零元素所在列所对应的节点形成的边；
18.s22：将二次有限元矩阵的每一行非零的元素映射成一个点和点的连接关系；
19.s23：每两个点对应一条无定向的边；
20.s24：依次重复步骤s22-s23，直至二次有限元矩阵的每一行的非零元素都存在连接关系，最终生成一个无定向图。
21.更进一步的，所述步骤s3中评价函数进行如下定义：对于一个矩阵，将其分解成n个子矩阵，n为2的倍数，子矩阵的耦合系数s为耦合矩阵的维度除以子矩阵的平均维度，最优分解即为耦合系数最小值对应的分块个数。
22.更进一步的，步骤s3中采用metis程序进行无定向图的分解。
23.更进一步的，所述步骤s1中初始有限原矩阵的来源为：通过读取网格文件，利用网格进行有限元数值建模并创建矩阵。
24.更进一步的，步骤s1中在进行初始有限元矩阵生成的同时，提前将端口节点与其它网格节点的连接关系进行备份。
25.一种应用如上述任一项大型稀疏矩阵加速求解方法的系统，包括：
26.还原模块：用于对初始有限元矩阵进行端口节点连接关系的还原，得到二次有限元矩阵；
27.转换模块：用于将二次有限元矩阵转换成无定向图；
28.分解模块：用于对无定向图进行分解，并得到最优分解；
29.重新编号模块：用于依据分解模块中的最优分解对二次有限元矩阵的节点进行重新编号；
30.重新排序模块：用于依据重新编号模块中的节点重新编号对二次有限元矩阵进行重新排序，生成最终有限元矩阵；
31.求解模块：用于对重新排序模块中的最终有限元矩阵进行求解。
32.控制模块：用于控制各个模块的工作。
33.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时执行上述任一项的方法。
34.3、有益效果
35.相比于现有技术，本发明的有益效果为：
36.(1)本发明在对大型稀疏矩阵进行求解时，通过先对初始有限元矩阵进行端口节点连接关系的还原，使得矩阵拓扑的完备性得到保障，为后续步骤提供精准度和完整度较高的支撑，从而保证整个过程的精度；其次，通过将二次有限元矩阵转化为无定向图随后进行分解，改变了以往直接对网格文件进行分解所带来的成本与效率问题；并且依据最优分解对二次有限元矩阵的节点进行重新编号随后生成新的最终有限元矩阵，避免新的最终有限元矩阵出现过大的耦合矩阵，进一步加快矩阵求解速度；最后对最终有限元矩阵利用并行求解，更深一层的提高了求解速度；整个方法在能够保证求解精度的同时大幅度提高求解效率，且操作简单，易于实施，能够在eda行业得到广泛的使用；
37.(2)本发明通过定义两个新的数据结构点和边，将二次有限元矩阵中的每一行进行扫描，得到所有点和边的关系，生成一个无定向的图，随后再对无定向的图进行分解并得到最优分解，提高效率，一改往日直接对有限元网格文件进行分解所造成的需要程序与磁盘建立io通信，并且在读取有限元网格文件之外产生额外开销；并且在评价函数中选择最小的耦合系数对应的分块个数作为最优分解，此种最优分解方案对于非线性矩阵具有更好的适用性，并且对于后续计算效率的提升帮助很大；
38.(3)本发明的系统通过将初始有限元矩阵进行还原节点连接关系后转化为无定向图，再对无定向图进行分解，得到最优分解，依据最优分解对二次有限元矩阵中的节点进行重新编号，依据重新的编号进行矩阵的重新排序，最终实现求解；在保证精度的同时提高效率，各模块之间独立工作的同时又相互协同，结构简单，易于控制。
附图说明
39.图1为本发明的流程示意图；
40.图2为原矩阵的拓扑关系图；
41.图3为分块后矩阵的拓扑关系图；
42.图4为大型系统矩阵求解示意图；
43.图5为分块后矩阵求解示意图；
44.图6为还原端口节点连接关系的示意图；
45.图7为转换成无定向图的示意图；
46.图8为最终有限元矩阵的转变示意图。
具体实施方式
47.下面结合具体实施例和附图对本发明进一步进行描述。
48.在eda行业，随着信号工作频率越来越高，芯片尺度越来越小，芯片内部结构越来越复杂，网格所对应的未知量的个数从几万，几十万上升到几百万不等，这种情况下，矩阵的求解变得更加缓慢。而本技术正是基于此环境下结合eda电路/电磁场分析的特点，利用稀疏矩阵的全部信息，还原电路/网格的拓扑以及连接关系，随后利用图分解技术对矩阵所对应的拓扑进行分块，对分块后的图进行矩阵重排，生成具有特殊模式的矩阵，一般称作带状加边的稀疏矩阵，然后利用数学知识对该系数矩阵进行并行求解，从而解决大规模矩阵计算速度慢的行业难题。因对于一个维度为n(n通常为几十万或几百万)的稀疏矩阵，通过将矩阵分成k个子矩阵，每个子矩阵的维度是p，这k个子矩阵通过一个维度是m的矩阵耦合。对于k个子矩阵，矩阵求逆运算的计算量(乘法和除法)约为k*p^3；对于耦合矩阵，其求逆的计算量约为m^3。分块后的矩阵，总的计算量约为k*p^3+m^3，其中k*p+m＝n。当矩阵维度n很大时，如果能使分块后的矩阵满足p》》m，那么分块后的计算量k*p^3+m^3约等于k*p^3，未采用分块算法的矩阵的计算量(k*p+m)^3约等于(k*p)^3。可以看出，(k*p)^3远大于k*p^3！再考虑上各个子矩阵可以采用并行计算来完成求逆运算，使得无论是计算量还是计算时间，分块后的矩阵在求解效率上都远大于未分块的矩阵。本技术具体的解决方案请参考下文。
49.实施例1
50.如图1所示，一种大型稀疏矩阵加速求解方法，包括以下步骤：
51.s1：对初始有限元矩阵进行还原端口节点的连接关系，得到还原后的二次有限元矩阵；具体的：如图6所示，图6中左边的a表示初始有限元矩阵，初始有限元矩阵中的端口节点所在行只有一个非零元1(即a中矩阵第五行为端口节点所在行)，隐藏了端口节点与其它节点之间的连接关系，使得无法通过端口节点所在行的非零元来获取的端口节点和其它网格节点的连接关系。由于本技术中后续步骤中存在分解和分块，如果不进行还原矩阵中端口节点与其它节点的连接关系，那么矩阵的连接关系是不完备的，在进行后续分解步骤时会出现伪分解或者是分解错误，伪分解会导致分块后的子矩阵条件数差，影响计算结果的精度。因此，本步骤便进行对初始有限元矩阵进行还原端口节点与其它节点的连接关系，即得到了图6中右边的b图，b即表示还原后的二次有限元矩阵，由b图可知此时的第五行端口节点左右两侧各有一个非零元1，即端口节点5与网格节点3以及网格节点7存在连接关系；并且该步骤中的还原是在进行初始有限元矩阵生成的同时，提前将端口节点与其它网格节点的连接关系进行备份，不对初始有限元矩阵进行改变，操作便捷；而初始有限元矩阵可以为通过读入并解析芯片结构的网格文件，利用网格来进行有限元数值建模并创建矩阵，得到初始有限元矩阵，初始有限元矩阵的拓扑关系如图2所示，由于利用网格构建有限元数值建模并创建矩阵是在现有技术中较为常见，具体可以参考“电磁场有限元方法”，且不属于本技术的核心改进点，因此本技术不再详细赘述；
52.s2：将二次有限元矩阵转化为无定向图如图7所示；具体的，该步骤包括如下步骤：
53.s21：定义两个新的数据结构：点和边，点：非零元素所在列即为该行所对应的节点，该行即为非零元素所在的行；边：非零元素所在列所对应的节点形成的边；
54.s22：将二次有限元矩阵(图7中左边a)的每一行非零的元素映射成一个点和点的连接关系；
55.s23：每两个点对应一条无定向的边；
56.s24：依次重复步骤s22-s23，直至二次有限元矩阵的每一行的非零元素都存在连
接关系，将这些点和边的关系，最终生成一个无定向图(如图7中右边b所示，右边b表示无定向图的数据结构，红色的数字代表被还原的与端口节点有连接关系的节点)；具体的，按照图7所示，每一行对应于一个点i，每一行的非零元素(这里假设都是1)所在位置对应于和点i有连接关系的点j，用两个矢量类型的结合，v和e就可以完整表示出矩阵或者图中有多少点，点和点之间这样连接的，因为i和j的连接以及j和i的连接记录了两次，所以是无方向的，那么v和e构成的集合就是无定向图；
57.s3：对无定向图进行分解，并且通过评价函数选择最优分解；这里需进行说明的是传统的方法是将直接对网格文件进行分解，通常来说这些文件是存储在硬盘上的，需要图分解程序读入内存来处理，然后把分解的结果保存成新的文件，这种方法是需要图分解程序和磁盘之间建立io通信，在读取网格文件之外产生额外开销；而本技术一改传统做法，该步骤根据网格文件已经建立好的矩阵出发，由于对初始有限元矩阵进行还原端口节点的连接关系得到二次有限元矩阵，二次有限元矩阵拓扑的完备性得到了保证；本技术不直接对网格文件进行分解，而是将二次有限元矩阵转换成无定向图之后再对其进行分解，避免硬盘读写网格文件所带来的额外开销；并且该步骤采用metis程序进行无定向图的分解，分解精度高且技术较为成熟；
58.该步骤中，选用了评价函数进行最优分解，评价函数的设置能够让分解程序进行时能够自动判断最优的分解结果，将其分解为最合适的数量避免了后续过程新矩阵中耦合矩阵过大的问题，从而加快矩阵平行计算效率。具体的，评价函数进行如下定义：对于一个矩阵(在本示例中即为二次有限元矩阵)，将其分解成n个子矩阵，n为2的倍数，子矩阵的耦合系数s为耦合矩阵的维度除以子矩阵的平均维度，最优分解即为耦合系数最小值对应的分块个数。采用这种方案对于非线性矩阵而言具有更好的适用性，对后续计算效率的提升具有较大的帮助。
59.s4：根据最优分解对二次有限元矩阵的节点重新编号；具体的，利用最优分解的方案对二次有限元矩阵的节点按照不同分区进行重新编号，即令最终生成的新矩阵的节点坐标(i,j)对应于原始矩阵的节点坐标(i，j)，建立一组节点映射，譬如二次有限元矩阵的原始编号为1,2,3,4,5,6,7,8,9,10,11,12；利用最优分解之后新矩阵的编号为1,2,3,11,6,7,8,9,10,12,4,5，利用新的编号可以快速重新生成分块的矩阵；
60.具体的，对步骤s4进行举例描述，譬如假设二次有限元矩阵的维度是12，即从图8进行理解的话就是存在12个节点，假设初始编号是顺序的即{1,2,3,4,5,6,7,8,9,10,11,12}，如果步骤s3中的最优分解是2，那么分解的结果可以是{111222222111}，分解结果括号中1表示该位置节点应该在第一个分块，括号里的2表示该位置节点在第二个分块。可以看出，第三个和第四个节点位于不同分块，第九个和第十个节点位于不同分块，从第三个和第四个里面，以及第九个和第十个里面任意挑出(可以挑出节点号小的节点)一个节点，认为是耦合矩阵的节点，然后重新排序，就可以得到最优分解后的新编号即{1,2,3,11,6,7,8,9,10,12,4,5}
61.s5：根据节点的重新编号即新编号对二次有限元矩阵重新排序，生成新的最终有限元矩阵如图3所示；即利用节点的重新编号对原来的二次有限元矩阵重新排序生成新的分块的矩阵(即最终有限元矩阵)，随后利用新节点重新排序激励端对应的右边项，对分块的矩阵和分块的右边项采用并行求解方案；具体的：在有限元分析和数值求解中，最重要的
就是解方程ax＝b。其中a是一个n x n的方阵，对于有限元问题，a是一个稀疏矩阵，x是一个n x 1的未知列矢量，b示一个n x 1的激励矢量，也叫做右边项，或者右端项矢量。当n的值很大时，比如超过100万，调用商业矩阵求解器直接求解未知量x就会变得很慢，图4和图5分别给出分块前后的矩阵模型，利用分块后的矩阵特性，可以加快矩阵求解速度。图4是描述ax＝b的示意图，其中a表示大型稀疏矩阵；x表示未知量；b是右边项；图4是分块后的矩阵和方程ax＝b的示意图，图中a1、a2、a3、a4表示分块后的子矩阵；a5是耦合矩阵；ei(i＝1,2,3,4)和fi(i＝1,2,3,4)用于描述子矩阵a1、a2、a3、a4和矩阵的耦合；xi(i＝1,2,3,4)对应分块后的未知量；bi(i＝1,2,3,4)对应对块后的右边项，分块后的矩阵，对于每一个子矩阵a1、a2、a3、a4，可以采用并行求解，极大的提高计算效率。
62.由图2和图3可知，分块前的矩阵，矩阵元素杂乱无章排列，分块后的矩阵，矩阵元素被分块4块，每一块的耦合被排列到图3中最右上角的子矩阵中，分块后的矩阵每一个子矩阵的矩阵元素较少，求解所需要的乘除法运算量较少，加上多线程并行技术，可以大大加快求解效率。并且依据节点的新编号进行重新生成矩阵不需要重新填写矩阵并且重排速度快，避免以往利用新的网格数据来生成矩阵；如图8所示，图8中左边a部分表示二次有限元矩阵，右边b部分表示新的最终有限元矩阵，最终有限元矩阵具有dbbd(doubleyborded block diagonal)的性质；
63.s6：对最终有限元矩阵进行求解，因最终有限元矩阵具有dbbd的性质，可以利用并行求解，具体求解方法可以参考相关文献，如：[1]:stabilizedborded block diagonal forms for parallel sparse solver.parallel computing.31(2005)275-289；[2]:parallel direct methods for block-diagonal-boardered sparse matrices.researchgate/2296092。
[0064]
本发明的方法通过对初始矩阵进行还原端口节点连接关系，为后续步骤提供精准度和完整度较高的支撑，避免后续出现条件数差的子矩阵，保证矩阵求解精度；其次，通过将二次有限元矩阵转化为无定向图随后进行分解，改变了以往直接对网格文件进行分解所带来的成本与效率问题；并且依据最优分解对二次有限元矩阵的节点进行重新编号随后生成新的最终有限元矩阵，避免新的最终有限元矩阵出现过大的耦合矩阵，进一步加快矩阵求解速度；最后对最终有限元矩阵利用并行求解，更深一层的提高了求解速度；整个方法在能够保证求解精度的同时大幅度提高求解效率，且操作简单，易于实施，能够在eda行业得到广泛的使用。
[0065]
实施例2
[0066]
一种应用如上述大型稀疏矩阵加速求解方法的系统，包括：
[0067]
还原模块：用于对初始有限元矩阵进行端口节点连接关系的还原，得到二次有限元矩阵；
[0068]
转换模块：用于将二次有限元矩阵转换成无定向图；
[0069]
分解模块：用于对无定向图进行分解，并得到最优分解；
[0070]
重新编号模块：用于依据分解模块中的最优分解对二次有限元矩阵的节点进行重新编号；
[0071]
重新排序模块：用于依据重新编号模块中的节点重新编号对二次有限元矩阵进行重新排序，生成最终有限元矩阵；
[0072]
求解模块：用于对重新排序模块中的最终有限元矩阵进行求解。
[0073]
控制模块：用于控制各个模块的工作。
[0074]
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，在计算机在执行时，可包括如上述各方法的实施例的流程。
[0075]
本发明的系统通过将初始有限元矩阵进行还原节点连接关系后转化为无定向图，再对无定向图进行分解，得到最优分解，依据最优分解对二次有限元矩阵中的节点进行重新编号，依据重新的编号进行矩阵的重新排序，最终实现求解；在保证精度的同时提高效率，使得精度与效率得到了较好的平衡，各模块之间独立工作的同时又相互协同，结构简单，易于控制。
[0076]
实施例3
[0077]
一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时执行上述的方法。所述的可读存储介质例如可以为但不限于电、磁、光、电磁、红外线或半导体的系统、装置或器件，或者任意以上的组合。
[0078]
本发明所述实例仅仅是对本发明的优选实施方式进行描述，并非对本发明构思和范围进行限定，在不脱离本发明设计思想的前提下，本领域工程技术人员对本发明的技术方案作出的各种变形和改进，均应落入本发明的保护范围。