CN113239591B - 面向dcu集群的大规模有限元网格并行分区的方法及装置 - Google Patents
面向dcu集群的大规模有限元网格并行分区的方法及装置 Download PDFInfo
- Publication number
- CN113239591B CN113239591B CN202110541398.XA CN202110541398A CN113239591B CN 113239591 B CN113239591 B CN 113239591B CN 202110541398 A CN202110541398 A CN 202110541398A CN 113239591 B CN113239591 B CN 113239591B
- Authority
- CN
- China
- Prior art keywords
- stack
- graph
- tid
- matrix
- ptr
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 214
- 238000005192 partition Methods 0.000 title claims description 35
- 230000009977 dual effect Effects 0.000 claims abstract description 55
- 238000001228 spectrum Methods 0.000 claims abstract description 26
- 238000000638 solvent extraction Methods 0.000 claims abstract description 15
- 238000009826 distribution Methods 0.000 claims abstract description 10
- 230000008569 process Effects 0.000 claims description 166
- 239000011159 matrix material Substances 0.000 claims description 114
- 238000004891 communication Methods 0.000 claims description 64
- 238000004364 calculation method Methods 0.000 claims description 24
- 230000000903 blocking effect Effects 0.000 claims description 21
- 238000012545 processing Methods 0.000 claims description 20
- 230000015654 memory Effects 0.000 claims description 17
- 230000006870 function Effects 0.000 claims description 11
- 238000005516 engineering process Methods 0.000 claims description 8
- 238000013507 mapping Methods 0.000 claims description 7
- 238000003860 storage Methods 0.000 claims description 6
- 238000005457 optimization Methods 0.000 claims description 5
- 230000003595 spectral effect Effects 0.000 claims description 3
- 238000003491 array Methods 0.000 claims 1
- 238000012360 testing method Methods 0.000 description 11
- 238000010586 diagram Methods 0.000 description 10
- 238000004422 calculation algorithm Methods 0.000 description 6
- 230000006872 improvement Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000003287 optical effect Effects 0.000 description 2
- 241000668842 Lepidosaphes gloverii Species 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000003825 pressing Methods 0.000 description 1
- UPPMZCXMQRVMME-UHFFFAOYSA-N valethamate Chemical compound CC[N+](C)(CC)CCOC(=O)C(C(C)CC)C1=CC=CC=C1 UPPMZCXMQRVMME-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- 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)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Analysis (AREA)
- Computing Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Algebra (AREA)
- Geometry (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Complex Calculations (AREA)
Abstract
本发明涉及高性能计算技术领域，特别是指一种面向DCU集群的大规模有限元网格并行分区的方法及装置，该方法包括：将求解区域离散为网格，通过带权对偶图描述离散所得的网格的分布情况，其中，所述带权对偶图中的顶点表示网格，所述带权对偶图中的边表示网格间的连接关系，通过子图划分任务的并行，可以极大提高递归谱二分法划分的速度，而且，在子图划分并行之上，还可以进行谱二分法计算热点的并行，进一步缩短了有限元网格划分的处理时间，解决了串行递归谱二分法存在着划分大规模网格时间需求过长的问题，使得解决问题的效率大大提高。
Description
技术领域
本发明涉及技术领域，特别是指一种面向DCU集群的大规模有限元网格并行分区的方法及装置。
背景技术
有限元方法的原理和基本概念是用较简单的问题代替复杂问题后进行求解，它将求解域分解为许多称为有限元网格的小的互联子域，对每一单元假定一个合适的近似解，然后根据求解域整体需要满足的条件得到问题的最终解。对求解域的离散化使得有限元方法不仅具有较好的计算精度，而且能适应各种复杂形状。依托超级计算机的海量计算吞吐能力和大规模并行计算技术，有限元方法成为解决复杂工程分析计算问题的有效途径，应用范围包含汽车、航空、机械制造、材料加工电子电器等各种科学技术领域。
使用有限元法进行大规模并行模拟的两个重要步骤是将离散得到的网格进行网格分区和处理器映射，使得多个网格聚合形成分区并与实际处理器数量相匹配，并做到分区间通信量少且负载均衡。常见的网格划分方法有递归二分法、启发式算法等。在各种区域划分算法中，递归谱二分法相对于传统算法与其他递归二分法而言，具有负载均衡、拓扑结构好的特点，非常适合用于处理高精度模拟的情况，但缺点就在于计算耗时非常长，并随着网格数量的提升大幅增长。
以基于递归谱二分法建立的串行网格分区模块为例，测试结果可以如图1所示，在处理1千万规模的网格分区时就需要消耗近11小时，而亿级网格甚至可能需要数天时间，由此可见，串行递归谱二分法存在着划分大规模网格时间需求过长的问题，导致解决问题的效率大大降低。
发明内容
本发明实施例提供了一种面向DCU集群的大规模有限元网格并行分区的方法及装置。所述技术方案如下：
一方面，提供了一种面向DCU集群的大规模有限元网格并行分区的方法，该方法应用于电子设备，该方法包括：
S1、将求解区域离散为网格，通过带权对偶图描述离散所得的网格的分布情况，其中，所述带权对偶图中的顶点表示网格，所述带权对偶图中的边表示网格间的连接关系；
S2、若当前进程编号为0，则输入带权对偶图G0，将所述带权对偶图G0压入堆栈GraphStack；若当前进程编号非0，则进入阻塞通信状态，直到完成与其他进程的通信获得输入所述带权对偶图G0及其进程编号PRE_PROC，将所述带权对偶图G0压入堆栈GraphStack；
S3、从所述堆栈GraphStack的栈顶取出图Gtop，如果图Gtop的划分深度达到预期值，则将所述图Gtop压入堆栈SubGraphStack；如果图Gtop的划分深度未达到预期值，则使用谱二分法将图Gtop划分为子图G1、G2，子图的划分深度为Gtop→depth+1；
设定同时运行的进程数量为P，根据当前进程编号得到通信目的进程的编号，如果所述通信目的进程的编号NEXT_PROC未达到最大进程数量P，则与进程NEXT_PROC进行通信，发送所述子图G2，将所述进程NEXT_PROC压入堆栈ProcStack中，将所述子图G1压入堆栈GraphStack中；如果所述通信目的进程的编号达到最大进程数量，则将子图G1、G2都压入堆栈GraphStack中；
重复执行S3，直到所述堆栈GraphStack为空时，执行S4；
S4、当所述堆栈ProcStack非空时，从所述堆栈ProcStack栈顶取出编号Ptop，进入阻塞通信状态直到完成与进程Ptop的通信，接收得到进程Ptop返回的子图集合GraphSet，将集合中的子图压入所述堆栈SubGraphStack中；
重复执行S4，直到ProcStack为空时，执行S5；
S5、使用所述堆栈SubGraphStack中的所有子图创建GraphSet，若当前进程为非0号进程，则进入阻塞通信状态，当完成与进程PRE_PROC的通信时，发送子图集合GraghSet；若当前进程为0号进程，则为所述GraghSet中的子图进行分区编号，输出网格与分区编号的对应结果。
可选地，所述S3中的所述使用谱二分法将图Gtop划分为子图G1、G2，包括下述S31-S33：
S31、计算得到所述图Gtop的拉普拉斯矩阵A；
S32、使用Lanczos迭代方法，将所述拉普拉斯矩阵A转换为三对角矩阵T；
S33、计算所述三对角矩阵T的费德勒向量，根据各分量的正负将所述图Gtop中对应的顶点划分为两个集合，得到子图G1、G2。
可选地，所述S32的使用Lanczos迭代方法，将所述拉普拉斯矩阵A转换为三对角矩阵T，包括：
Lanczos迭代的数学定义如下：
若A为对称矩阵，则存在正交矩阵Q＝[q1,q2,…,qn]使得T＝QTAQ，其中qi为N维列向量，T为三对角阵；
给定任意列向量q1(||q1||2＝1)，通过迭代得到正交矩阵Q以及三对角阵T；
其中，迭代格式如下：
④、μ1＝Aq1；i＝1
⑤、ri＝μi-αiqi；βi＝||ri||2
⑥、若βi＝0或达到迭代次数限制则退出，否则：
qi+1＝ri/βi；μi+1＝Aqi+1-βiqi；i＝i+1
转到执行②；
其中，所述Aq1为稀疏矩阵向量乘法。
可选地，所述稀疏矩阵向量乘法基于CSR格式存储；
所述稀疏矩阵向量乘法采用HIP编程的方式进行并行优化。
可选地，设定稀疏矩阵的尺寸为NxN,向量的尺寸为Nx1，所述稀疏矩阵中所有非零元素的个数为NotZeroNums，设定数组IA[N+1]、JA[NotZeroNums]、VA[NotZeroNums]存储稀疏矩阵，设定数组x[N]存储向量；
所述稀疏矩阵向量乘法采用HIP编程的方式进行并行优化，包括下述步骤S321-S323：
S321、对于向量规模小于4096的稀疏矩阵向量乘法，不采取辅助优化技术；
S322、引入中间数组v[NotZeroNums],将矩阵向量乘运算分解为乘运算与加运算两个核心执行；
在计算核心一中完成乘运算部分，设置每个Block中的线程数量为threadsPerBlock，Block数量为ceil(NotZeroNums/threadsPerBlock),ceil为向上取整函数，计算核心一的步骤包括下述S3221-S3224：
S3221、计算线程全局编号：
tid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x；
S3222、计算线程总数icr＝hipBlockDim_x*hipGridDim_x；
S3223、根据tid取出VA[tid]、JA[tid]、x[JA[tid]]，计算VA[tid]*x[JA[tid]]存入v[tid]；
S3224、计算tid＝tid+icr，重复执行S73-S74，直到tid超过NotZeroNums，则执行S8；
S323、在计算核心二中完成加运算部分，设定Block数量为ceil(N/threadsPerBlock),ceil为向上取整函数，计算核心二的步骤包括下述S81-S89：
S3231、计算线程全局编号：
gid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x
S3232、取得线程局部编号tid＝hipThreadIdx_x；
S3233、设定共享内存空间ptr_s[threadsPerBlock+1]、v_s[sizeSharedMemory]；
S3234、根据编号从数组IA中取得元素下标范围，存入ptr_s中：
ptr_s[tid]＝IA[gid]
S3235、若tid为0，即该线程是Block内的0号线程，则计算：
MaxIndex＝gid+threadsPerBlock
其中，MaxIndex为所述Block负责的最大行编号，若MaxIndex大于N，则取MaxIndex为N；取出最大元素下标IA[MaxIndex]存入ptr_s[threadsPerBlock]中；
S3236、计算变量
temp＝(ptr_s[threadsPerBlock]-ptr_s[0])/threadsPerBlock+1；nlen＝min(temp*threadsPerBlock,sizeSharedMemory)。
得到一次循环过程中加载的元素数量nlen，设定循环变量i的初始值为ptr_s[0]，中间结果sum的初始值为0；
S3237、计算index＝i+tid，以v_s[tid]、v[index]为起点，threadsPerBlock为步长，从v数组中加载数据到v_s数组；
加载完成后，若ptr_s[tid]<＝i+nlen-1且ptr_s[tid+1]>i,即加载到共享内存空间中的元素属于第gid行，则计算元素下标的起始范围：
row_st＝max(ptr_s[tid]–i,0)
row_ed＝min(ptr_s[tid+1]–i,nlen)
累加v_s[row_st]至v_s[row_ed-1]的所有元素到变量sum上；
S3238、i＝i+nlen；重复执行S87直到i达到ptr_s[threadsPerBlock]时，执行S89；
S3239、将sum存入y[gid]。
可选地，所述S3中的根据当前进程编号得到通信目的进程的编号，包括：
设当前进程编号为i，则通过下述公式(1)计算得到通信目的进程的编号：
另一方面，提供了一种面向DCU集群的大规模有限元网格并行分区的装置，该装置应用于电子设备，该装置包括：
离散单元，被配置为执行S1：将求解区域离散为网格，通过带权对偶图描述离散所得的网格的分布情况，其中，所述带权对偶图中的顶点表示网格，所述带权对偶图中的边表示网格间的连接关系；
压栈单元，被配置为执行S2：若当前进程编号为0，则输入带权对偶图G0，将所述带权对偶图G0压入堆栈GraphStack；若当前进程编号非0，则进入阻塞通信状态，直到完成与其他进程的通信获得输入所述带权对偶图G0及其进程编号PRE_PROC，将所述带权对偶图G0压入堆栈GraphStack；
处理单元，被配置为执行S3：从所述堆栈GraphStack的栈顶取出图Gtop，如果图Gtop的划分深度达到预期值，则将所述图Gtop压入堆栈SubGraphStack；如果图Gtop的划分深度未达到预期值，则使用谱二分法将图Gtop划分为子图G1、G2，子图的划分深度为Gtop→depth+1；
设定同时运行的进程数量为P，根据当前进程编号得到通信目的进程的编号，如果所述通信目的进程的编号NEXT_PROC未达到最大进程数量P，则与进程NEXT_PROC进行通信，发送所述子图G2，将所述进程NEXT_PROC压入堆栈ProcStack中，将所述子图G1压入堆栈GraphStack中；如果所述通信目的进程的编号达到最大进程数量，则将子图G1、G2都压入堆栈GraphStack中；
重复执行S3，直到所述堆栈GraphStack为空时，执行S4；
接收单元，被配置为执行S4：当所述堆栈ProcStack非空时，从所述堆栈ProcStack栈顶取出编号Ptop，进入阻塞通信状态直到完成与进程Ptop的通信，接收得到进程Ptop返回的子图集合GraphSet，将集合中的子图压入所述堆栈SubGraphStack中；
重复执行S4，直到ProcStack为空时，执行S5；
输出单元，被配置为执行S5：使用所述堆栈SubGraphStack中的所有子图创建GraphSet，若当前进程为非0号进程，则进入阻塞通信状态，当完成与进程PRE_PROC的通信时，发送子图集合GraghSet；若当前进程为0号进程，则为所述GraghSet中的子图进行分区编号，输出网格与分区编号的对应结果。
可选地，所述被配置为执行S3的处理单元，进一步被配置为执行下述S31-S33：
S31、计算得到所述图Gtop的拉普拉斯矩阵A；
S32、使用Lanczos迭代方法，将所述拉普拉斯矩阵A转换为三对角矩阵T；
S33、计算所述三对角矩阵T的费德勒向量，根据各分量的正负将所述图Gtop中对应的顶点划分为两个集合，得到子图G1、G2。
可选地，所述被配置为执行S32的处理单元，进一步被配置为执行：
Lanczos迭代的数学定义如下：
若A为对称矩阵，则存在正交矩阵Q＝[q1,q2,…,qn]使得T＝QTAQ，其中qi为N维列向量，T为三对角阵；
给定任意列向量q1(||q1||2＝1)，通过迭代得到正交矩阵Q以及三对角阵T；
其中，迭代格式如下：
⑦、μ1＝Aq1；i＝1
⑧、ri＝μi-αiqi；βi＝||ri||2
⑨、若βi＝0或达到迭代次数限制则退出，否则：
qi+1＝ri/βi；μi+1＝Aqi+1-βiqi；i＝i+1
转到执行②；
其中，所述Aq1为稀疏矩阵向量乘法。
可选地，所述稀疏矩阵向量乘法基于CSR格式存储；
所述稀疏矩阵向量乘法采用HIP编程的方式进行并行优化。
可选地，设定稀疏矩阵的尺寸为NxN,向量的尺寸为Nx1，所述稀疏矩阵中所有非零元素的个数为NotZeroNums，设定数组IA[N+1]、JA[NotZeroNums]、VA[NotZeroNums]存储稀疏矩阵，设定数组x[N]存储向量；
所述被配置为执行S32的处理单元，进一步被配置为执行下述S321-S323：
S321、对于向量规模小于4096的稀疏矩阵向量乘法，不采取辅助优化技术；
S322、引入中间数组v[NotZeroNums],将矩阵向量乘运算分解为乘运算与加运算两个核心执行；
在计算核心一中完成乘运算部分，设置每个Block中的线程数量为threadsPerBlock，Block数量为ceil(NotZeroNums/threadsPerBlock),ceil为向上取整函数，计算核心一的步骤包括下述S3221-S3224：
S3221、计算线程全局编号：
tid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x；
S3222、计算线程总数icr＝hipBlockDim_x*hipGridDim_x；
S3223、根据tid取出VA[tid]、JA[tid]、x[JA[tid]]，计算VA[tid]*x[JA[tid]]存入v[tid]；
S3224、计算tid＝tid+icr，重复执行S73-S74，直到tid超过NotZeroNums，则执行S8；
S323、在计算核心二中完成加运算部分，设定Block数量为ceil(N/threadsPerBlock),ceil为向上取整函数，计算核心二的步骤包括下述S81-S89：
S3231、计算线程全局编号：
gid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x
S3232、取得线程局部编号tid＝hipThreadIdx_x；
S3233、设定共享内存空间ptr_s[threadsPerBlock+1]、v_s[sizeSharedMemory]；
S3234、根据编号从数组IA中取得元素下标范围，存入ptr_s中：
ptr_s[tid]＝IA[gid]
S3235、若tid为0，即该线程是Block内的0号线程，则计算：
MaxIndex＝gid+threadsPerBlock
其中，MaxIndex为所述Block负责的最大行编号，若MaxIndex大于N，则取MaxIndex为N；取出最大元素下标IA[MaxIndex]存入ptr_s[threadsPerBlock]中；
S3236、计算变量
temp＝(ptr_s[threadsPerBlock]-ptr_s[0])/threadsPerBlock+1；nlen＝min(temp*threadsPerBlock,sizeSharedMemory)。
得到一次循环过程中加载的元素数量nlen，设定循环变量i的初始值为ptr_s[0]，中间结果sum的初始值为0；
S3237、计算index＝i+tid，以v_s[tid]、v[index]为起点，threadsPerBlock为步长，从v数组中加载数据到v_s数组；
加载完成后，若ptr_s[tid]<＝i+nlen-1且ptr_s[tid+1]>i,即加载到共享内存空间中的元素属于第gid行，则计算元素下标的起始范围：
row_st＝max(ptr_s[tid]–i,0)
row_ed＝min(ptr_s[tid+1]–i,nlen)
累加v_s[row_st]至v_s[row_ed-1]的所有元素到变量sum上；
S3238、i＝i+nlen；重复执行S87直到i达到ptr_s[threadsPerBlock]时，执行S89；
S3239、将sum存入y[gid]。
可选地，所述被配置为执行S3的处理单元，进一步被配置为执行：
设当前进程编号为i，则通过下述公式(1)计算得到通信目的进程的编号：
一方面，提供了一种电子设备，所述电子设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由所述处理器加载并执行以实现上述面向DCU集群的大规模有限元网格并行分区的方法。
一方面，提供了一种计算机可读存储介质，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现上述面向DCU集群的大规模有限元网格并行分区的方法。
本发明实施例提供的技术方案带来的有益效果至少包括：
上述方案中，通过子图划分任务的并行，可以极大提高递归谱二分法划分的速度，而且，在子图划分并行之上，还可以进行谱二分法计算热点的并行，进一步缩短了有限元网格划分的处理时间，解决了串行递归谱二分法存在着划分大规模网格时间需求过长的问题，使得解决问题的效率大大提高。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的现有技术中递归谱二分法的串行网格划分模块运行时间测试结果的示意图；
图2是本发明实施例提供的一种实施环境图；
图3是本发明实施例提供的一种面向DCU集群的大规模有限元网格并行分区的方法流程图；
图4是本发明实施例提供的一种面向DCU集群的大规模有限元网格并行分区的网格的对偶图表示的示意图；
图5(a)是本发明实施例提供的一种面向DCU集群的大规模有限元网格并行分区的测试结果与现有技术测试结果的对比示意图；
图5(b)是本发明实施例提供的一种子图划分并行之上进行谱二分法计算热点的并行的测试结果与单独进行子图划分并行的测试结果的对比示意图；
图6是本发明实施例提供的一种面向DCU集群的大规模有限元网格并行分区的装置框图；
图7是本发明实施例提供的一种电子设备的结构示意图。
具体实施方式
为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。
本发明实施例提供了一种面向DCU集群的大规模有限元网格并行分区的方法，如图2所示，该实施环境可以包括至少一个终端101、以及用于为该终端101提供服务的服务器102。至少一个终端101通过无线或者有线网络和服务器102连接，该终端101可以为能够访问服务器102的计算机设备或智能终端等。对于有限元网格并行分区的过程，可以是终端101独自进行处理，也可以是服务器102独自进行处理，还可以是终端101与服务器102共同进行处理，下面例举其中一种可行的方式：终端101中可以设置有限元网格并行分区处理程序，服务器102中存储有稀疏矩阵向量乘法并行算法，服务器102通过稀疏矩阵向量乘法并行算法进行计算后，将计算结果发送至终端101，终端101通过计算结果对有限元网格进行并行分区。
本发明实施例提供了一种面向DCU集群的大规模有限元网格并行分区的方法，该方法可以由电子设备实现，该电子设备可以是终端或服务器。参照图3所示的面向DCU集群的大规模有限元网格并行分区的方法流程图，该方法的处理流程可以包括如下的步骤：
S1、将求解区域离散为网格，通过带权对偶图描述离散所得的网格的分布情况，其中，带权对偶图中的顶点表示网格，带权对偶图中的边表示网格间的连接关系。
一种可行的实施方式中，在三维数值模拟问题中，一般将求解区域离散为一系列四面体或六面体网格；而二维问题中一般离散为三角形或四边形。可以使用一个带权对偶图G来描述离散所得网格的分布情况，图中的顶点代表网格，图中的边代表网格间的连接关系，如图4所示，该图中给出了二维四边形网格的带权对偶图转换示例，其中点连接的边权值为1，线连接的边权值为2，无连接的边权值为0。三维情况下还可以存在面连接，对应边权重为4。
使用如下的数据结构描述由网格转换而来的图G，其中nel为网格数量。ELIST数组存储子图包含的所有网格的编号；PMAP数组存储网格编号对应的分区编号；IA、JA、VA数组以CSR格式存储图的拓扑结构。
对图G递归使用谱二分法就可以完成网格分区任务，若设定划分的最大深度为N，则最终可以得到2N个分区。子图划分任务的并行基于MPI进行实现，设定同时运行的进程数量为P。在读取网格文件并转换得到描述网格分布信息的带权对偶图G0后启动并行程序，设定图G0的划分深度G0→depth为0。
S2、若当前进程编号为0，则输入带权对偶图G0，将带权对偶图G0压入堆栈GraphStack；若当前进程编号非0，则进入阻塞通信状态，直到完成与其他进程的通信获得输入带权对偶图G0及其进程编号PRE_PROC，将带权对偶图G0压入堆栈GraphStack。
S3、从堆栈GraphStack的栈顶取出图Gtop，如果图Gtop的划分深度达到预期值，则将图Gtop压入堆栈SubGraphStack；如果图Gtop的划分深度未达到预期值，则使用谱二分法将图Gtop划分为子图G1、G2，子图的划分深度为Gtop→depth+1。
其中，谱二分法是计算图的拉普拉斯矩阵的费德勒向量的方法，费德勒向量是图的拉普拉斯矩阵的次小特征值对应的特征向量，可以根据各分量的正负符号对图进行划分。根据费德勒向量各分量的正负符号将分量对应的图中顶点划分到两个不同的子图中。当矩阵为大型稀疏矩阵时一般采用Lanczos迭代法将其转换为较低阶数的三对角矩阵进行费德勒向量的近似求解。
一种可行的实施方式中，根据当前进程编号得到通信目的进程的编号，设当前进程编号为i，则通过下述公式(1)计算得到通信目的进程的编号：
如果通信目的进程的编号NEXT_PROC未达到最大进程数量P，则与进程NEXT_PROC进行通信，发送子图G2，将进程NEXT_PROC压入堆栈ProcStack中，将子图G1压入堆栈GraphStack中；如果通信目的进程的编号达到最大进程数量，则将子图G1、G2都压入堆栈GraphStack中。
可选地，S3中使用谱二分法将图Gtop划分为子图G1、G2，具体可以包括下述S31-S33：
S31、计算得到图Gtop的拉普拉斯矩阵A。
其中，拉普拉斯矩阵为无向图的度矩阵与邻接矩阵的差，它是图论函数的拉普拉斯算子，所有的特征值均为正数，且可以描述图的相关信息，如连通度等
S32、使用Lanczos迭代方法，将拉普拉斯矩阵A转换为三对角矩阵T。
其中，Lanczos迭代为利用正交矩阵Q将对称矩阵A转换为三对角矩阵T的一种算法。方法为构建等式AQ＝QT，通过展开矩阵乘得到递推关系。通过给定初始向量q递推构造矩阵Q和T。
一种可行的实施方式中，Lanczos迭代的数学定义如下：
若A为对称矩阵，则存在正交矩阵Q＝[q1，q2，…，qn]使得T＝QTAQ，其中qi为N维列向量，T为三对角阵；
给定任意列向量q1(||q1||2＝1)，通过迭代得到正交矩阵Q以及三对角阵T；
其中，迭代格式如下：
①、μ1＝Aq1；i＝1
②、ri＝μi-αiqi；βi＝||ri||2
③、若βi＝0或达到迭代次数限制则退出，否则：
qi+1＝ri/βi；μi+1＝Aqi+1-βiqi；i＝i+1
转到执行②；
其中，Aq1为稀疏矩阵向量乘法。
可选地，在图划分过程中计算热点主要集中在Lanczos迭代过程中的稀疏矩阵向量乘运算，因此给出了一种DCU上的稀疏矩阵向量乘法并行方案，以CSR格式为例进行实现，谱二分法的计算热点为稀疏矩阵向量乘法，稀疏矩阵向量乘法可以基于CSR格式存储，为了进一步提升运算速度，稀疏矩阵向量乘法可以采用HIP编程的方式进行并行优化。
设定稀疏矩阵的尺寸为NxN,向量的尺寸为Nx1，稀疏矩阵中所有非零元素的个数为NotZeroNums，设定数组IA[N+1]、JA[NotZeroNums]、VA[NotZeroNums]存储稀疏矩阵，设定数组x[N]存储向量。
稀疏矩阵向量乘法采用HIP编程的方式进行并行优化，包括下述步骤S321-323：
S321、对于向量规模小于4096的稀疏矩阵向量乘法，不采取辅助优化技术。这样，可以减少效果不好的并行计算，真正达到计算速度的提升。
S322、引入中间数组v[NotZeroNums],将矩阵向量乘运算分解为乘运算与加运算两个核心执行。
在计算核心一中完成乘运算部分，设置每个Block中的线程数量为threadsPerBlock，Block数量为ceil(NotZeroNums/threadsPerBlock),ceil为向上取整函数，计算核心一的步骤包括下述S3221-S3224：
S3221、计算线程全局编号：
tid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x。
S3222、计算线程总数icr＝hipBlockDim_x*hipGridDim_x。
S3223、根据tid取出VA[tid]、JA[tid]、x[JA[tid]]，计算VA[tid]*x[JA[tid]]存入v[tid]。
S3224、计算tid＝tid+icr，重复执行S73-S74，直到tid超过NotZeroNums，则执行S323。
S323、在计算核心二中完成加运算部分，设定Block数量为ceil(N/threadsPerBlock),ceil为向上取整函数，计算核心二的步骤包括下述S3231-S3239：
S3231、计算线程全局编号：
gid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x
S3232、取得线程局部编号tid＝hipThreadIdx_x。
S3233、设定共享内存空间ptr_s[threadsPerBlock+1]、v_s[sizeSharedMemory]。
S3234、根据编号从数组IA中取得元素下标范围，存入ptr_s中：
ptr_s[tid]＝IA[gid]
S3235、若tid为0，即该线程是Block内的0号线程，则计算：
MaxIndex＝gid+threadsPerBlock
其中，MaxIndex为Block负责的最大行编号，若MaxIndex大于N，则取MaxIndex为N。取出最大元素下标IA[MaxIndex]存入ptr_s[threadsPerBlock]中。
S3236、计算变量
temp＝(ptr_s[threadsPerBlock]-ptr_s[0])/threadsPerBlock+1。nlen＝min(temp*threadsPerBlock,sizeSharedMemory)。
得到一次循环过程中加载的元素数量nlen，设定循环变量i的初始值为ptr_s[0]，中间结果sum的初始值为0。
S3237、计算index＝i+tid，以v_s[tid]、v[index]为起点，threadsPerBlock为步长，从v数组中加载数据到v_s数组。
加载完成后，若ptr_s[tid]<＝i+nlen-1且ptr_s[tid+1]>i,即加载到共享内存空间中的元素属于第gid行，则计算元素下标的起始范围：
row_st＝max(ptr_s[tid]–i,0)
row_ed＝min(ptr_s[tid+1]–i,nlen)
累加v_s[row_st]至v_s[row_ed-1]的所有元素到变量sum上。
S3238、i＝i+nlen。重复执行S87直到i达到ptr_s[threadsPerBlock]时，执行S89。
S3239、将sum存入y[gid]。
下面对上述步骤S321-323通过代码进行说明：
(1)设定数组IA，IA[0]＝0，对矩阵中某行i，扫描稀疏矩阵A中该行所有的非零元，将非零元的数量记作N，则IA[i+1]＝IA[i]+N。可以通过IA[i+1]-IA[i]计算得到第i行的非零元素数量N。IA[i]记录了第i行元素在JA、VA中存储的起始地址。
(2)对矩阵中某行i，扫描稀疏矩阵A中该行所有的非零元，将非零元素的列号按顺序存入数组JA中，将其数据顺序存入数组VA中。
/>
S33、计算三对角矩阵T的费德勒向量，根据各分量的正负将图Gtop中对应的顶点划分为两个集合，得到子图G1、G2。
重复执行S3，直到堆栈GraphStack为空时，执行S4。
S4、当堆栈ProcStack非空时，从堆栈ProcStack栈顶取出编号Ptop，进入阻塞通信状态直到完成与进程Ptop的通信，接收得到进程Ptop返回的子图集合GraphSet，将集合中的子图压入堆栈SubGraphStack中。
重复执行S4，直到ProcStack为空时，执行S5。
S5、使用堆栈SubGraphStack中的所有子图创建GraphSet，若当前进程为非0号进程，则进入阻塞通信状态，当完成与进程PRE_PROC的通信时，发送子图集合GraghSet；若当前进程为0号进程，则为GraghSet中的子图进行分区编号，输出网格与分区编号的对应结果。
下面通过代码对上述步骤S2-S5进行说明：
/>
/>
通过上述面向DCU集群的大规模有限元网格并行分区的方法进行试验，其测试结果可以如图5(a)所示，由测试结果可以看出，通过应用基于子图任务的并行技术，在8进程的并行规模下，与同样基于谱二分法的串行分区软件相比，本申请的方法可以获得近5.6倍的速度提升。
本申请实施例中，通过子图划分任务的并行，可以极大提高递归谱二分法划分的速度，而且，在子图划分并行之上，还可以进行谱二分法计算热点的并行，其测试结果如图5(b)所示，由测试结果可以看出，通过应用基于计算热点的并行技术，可以获得近30％的速度提升，进一步缩短了有限元网格划分的处理时间，解决了串行递归谱二分法存在着划分大规模网格时间需求过长的问题，使得解决问题的效率大大提高。
图6是根据一示例性实施例示出的一种面向DCU集群的大规模有限元网格并行分区的装置框图。参照图6，该装置包括离散单元610、压栈单元620、处理单元630、接收单元640以及输出单元650。
离散单元610，被配置为执行S1：将求解区域离散为网格，通过带权对偶图描述离散所得的网格的分布情况，其中，所述带权对偶图中的顶点表示网格，所述带权对偶图中的边表示网格间的连接关系；
压栈单元620，被配置为执行S2：若当前进程编号为0，则输入带权对偶图G0，将所述带权对偶图G0压入堆栈GraphStack；若当前进程编号非0，则进入阻塞通信状态，直到完成与其他进程的通信获得输入所述带权对偶图G0及其进程编号PRE_PROC，将所述带权对偶图G0压入堆栈GraphStack；
处理单元630，被配置为执行S3：从所述堆栈GraphStack的栈顶取出图Gtop，如果图Gtop的划分深度达到预期值，则将所述图Gtop压入堆栈SubGraphStack；如果图Gtop的划分深度未达到预期值，则使用谱二分法将图Gtop划分为子图G1、G2，子图的划分深度为Gtop→depth+1；
设定同时运行的进程数量为P，根据当前进程编号得到通信目的进程的编号，如果所述通信目的进程的编号NEXT_PROC未达到最大进程数量P，则与进程NEXT_PROC进行通信，发送所述子图G2，将所述进程NEXT_PROC压入堆栈ProcStack中，将所述子图G1压入堆栈GraphStack中；如果所述通信目的进程的编号达到最大进程数量，则将子图G1、G2都压入堆栈GraphStack中；
重复执行S3，直到所述堆栈GraphStack为空时，执行S4；
接收单元640，被配置为执行S4：当所述堆栈ProcStack非空时，从所述堆栈ProcStack栈顶取出编号Ptop，进入阻塞通信状态直到完成与进程Ptop的通信，接收得到进程Ptop返回的子图集合GraphSet，将集合中的子图压入所述堆栈SubGraphStack中；
重复执行S4，直到ProcStack为空时，执行S5；
输出单元650，被配置为执行S5：使用所述堆栈SubGraphStack中的所有子图创建GraphSet，若当前进程为非0号进程，则进入阻塞通信状态，当完成与进程PRE_PROC的通信时，发送子图集合GraghSet；若当前进程为0号进程，则为所述GraghSet中的子图进行分区编号，输出网格与分区编号的对应结果。
可选地，所述被配置为执行S3的处理单元630，进一步被配置为执行下述S31-S33：
S31、计算得到所述图Gtop的拉普拉斯矩阵A；
S32、使用Lanczos迭代方法，将所述拉普拉斯矩阵A转换为三对角矩阵T；
S33、计算所述三对角矩阵T的费德勒向量，根据各分量的正负将所述图Gtop中对应的顶点划分为两个集合，得到子图G1、G2。
可选地，所述被配置为执行S32的处理单元630，进一步被配置为执行：
Lanczos迭代的数学定义如下：
若A为对称矩阵，则存在正交矩阵Q＝[q1，q2，…，qn]使得T＝QTAQ，其中qi为N维列向量，T为三对角阵；
给定任意列向量q1(||q1||2＝1)，通过迭代得到正交矩阵Q以及三对角阵T；
其中，迭代格式如下：
⑩、μ1＝Aq1；i＝1
ri＝μi-αiqi；βi＝||ri||2
若βi＝0或达到迭代次数限制则退出，否则：
qi+1＝ri/βi；μi+1＝Aqi+1-βiqi；i＝i+1
转到执行②；
其中，所述Aq1为稀疏矩阵向量乘法。
可选地，所述稀疏矩阵向量乘法基于CSR格式存储；
所述稀疏矩阵向量乘法采用HIP编程的方式进行并行优化。
可选地，设定稀疏矩阵的尺寸为NxN,向量的尺寸为Nx1，所述稀疏矩阵中所有非零元素的个数为NotZeroNums，设定数组IA[N+1]、JA[NotZeroNums]、VA[NotZeroNums]存储稀疏矩阵，设定数组x[N]存储向量；
所述被配置为执行S32的处理单元630，进一步被配置为执行下述S321-S323：
S321、对于向量规模小于4096的稀疏矩阵向量乘法，不采取辅助优化技术；
S322、引入中间数组v[NotZeroNums],将矩阵向量乘运算分解为乘运算与加运算两个核心执行；
在计算核心一中完成乘运算部分，设置每个Block中的线程数量为threadsPerBlock，Block数量为ceil(NotZeroNums/threadsPerBlock),ceil为向上取整函数，计算核心一的步骤包括下述S3221-S3224：
S3221、计算线程全局编号：
tid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x；
S3222、计算线程总数icr＝hipBlockDim_x*hipGridDim_x；
S3223、根据tid取出VA[tid]、JA[tid]、x[JA[tid]]，计算VA[tid]*x[JA[tid]]存入v[tid]；
S3224、计算tid＝tid+icr，重复执行S73-S74，直到tid超过NotZeroNums，则执行S8；
S323、在计算核心二中完成加运算部分，设定Block数量为ceil(N/threadsPerBlock),ceil为向上取整函数，计算核心二的步骤包括下述S81-S89：
S3231、计算线程全局编号：
gid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x
S3232、取得线程局部编号tid＝hipThreadIdx_x；
S3233、设定共享内存空间ptr_s[threadsPerBlock+1]、v_s[sizeSharedMemory]；
S3234、根据编号从数组IA中取得元素下标范围，存入ptr_s中：
ptr_s[tid]＝IA[gid]
S3235、若tid为0，即该线程是Block内的0号线程，则计算：
MaxIndex＝gid+threadsPerBlock
其中，MaxIndex为所述Block负责的最大行编号，若MaxIndex大于N，则取MaxIndex为N；取出最大元素下标IA[MaxIndex]存入ptr_s[threadsPerBlock]中；
S3236、计算变量
temp＝(ptr_s[threadsPerBlock]-ptr_s[0])/threadsPerBlock+1；nlen＝min(temp*threadsPerBlock,sizeSharedMemory)。
得到一次循环过程中加载的元素数量nlen，设定循环变量i的初始值为ptr_s[0]，中间结果sum的初始值为0；
S3237、计算index＝i+tid，以v_s[tid]、v[index]为起点，threadsPerBlock为步长，从v数组中加载数据到v_s数组；
加载完成后，若ptr_s[tid]<＝i+nlen-1且ptr_s[tid+1]>i,即加载到共享内存空间中的元素属于第gid行，则计算元素下标的起始范围：
row_st＝max(ptr_s[tid]–i,0)
row_ed＝min(ptr_s[tid+1]–i,nlen)
累加v_s[row_st]至v_s[row_ed-1]的所有元素到变量sum上；
S3238、i＝i+nlen；重复执行S87直到i达到ptr_s[threadsPerBlock]时，执行S89；
S3239、将sum存入y[gid]。
可选地，所述被配置为执行S3的处理单元630，进一步被配置为执行：
设当前进程编号为i，则通过下述公式(1)计算得到通信目的进程的编号：
图7是本发明实施例提供的一种电子设备700的结构示意图，该电子设备700可因配置或性能不同而产生比较大的差异，可以包括一个或一个以上处理器(centralprocessing units，CPU)701和一个或一个以上的存储器702，其中，所述存储器702中存储有至少一条指令，所述至少一条指令由所述处理器701加载并执行以实现下述面向DCU集群的大规模有限元网格并行分区的方法的步骤：
S1、将求解区域离散为网格，通过带权对偶图描述离散所得的网格的分布情况，其中，所述带权对偶图中的顶点表示网格，所述带权对偶图中的边表示网格间的连接关系；
S2、若当前进程编号为0，则输入带权对偶图G0，将所述带权对偶图G0压入堆栈GraphStack；若当前进程编号非0，则进入阻塞通信状态，直到完成与其他进程的通信获得输入所述带权对偶图G0及其进程编号PRE_PROC，将所述带权对偶图G0压入堆栈GraphStack；
S3、从所述堆栈GraphStack的栈顶取出图Gtop，如果图Gtop的划分深度达到预期值，则将所述图Gtop压入堆栈SubGraphStack；如果图Gtop的划分深度未达到预期值，则使用谱二分法将图Gtop划分为子图G1、G2，子图的划分深度为Gtop→depth+1；
设定同时运行的进程数量为P，根据当前进程编号得到通信目的进程的编号，如果所述通信目的进程的编号NEXT_PROC未达到最大进程数量P，则与进程NEXT_PROC进行通信，发送所述子图G2，将所述进程NEXT_PROC压入堆栈ProcStack中，将所述子图G1压入堆栈GraphStack中；如果所述通信目的进程的编号达到最大进程数量，则将子图G1、G2都压入堆栈GraphStack中；
重复执行S3，直到所述堆栈GraphStack为空时，执行S4；
S4、当所述堆栈ProcStack非空时，从所述堆栈ProcStack栈顶取出编号Ptop，进入阻塞通信状态直到完成与进程Ptop的通信，接收得到进程Ptop返回的子图集合GraphSet，将集合中的子图压入所述堆栈SubGraphStack中；
重复执行S4，直到ProcStack为空时，执行S5；
S5、使用所述堆栈SubGraphStack中的所有子图创建GraphSet，若当前进程为非0号进程，则进入阻塞通信状态，当完成与进程PRE_PROC的通信时，发送子图集合GraghSet；若当前进程为0号进程，则为所述GraghSet中的子图进行分区编号，输出网格与分区编号的对应结果。
本申请实施例中，通过子图划分任务的并行，可以极大提高递归谱二分法划分的速度，而且，在子图划分并行之上，还可以进行谱二分法计算热点的并行，进一步缩短了有限元网格划分的处理时间，解决了串行递归谱二分法存在着划分大规模网格时间需求过长的问题，使得解决问题的效率大大提高。
在示例性实施例中，还提供了一种计算机可读存储介质，例如包括指令的存储器，上述指令可由终端中的处理器执行以完成上述面向DCU集群的大规模有限元网格并行分区的方法。例如，所述计算机可读存储介质可以是ROM、随机存取存储器(RAM)、CD-ROM、磁带、软盘和光数据存储设备等。
本领域普通技术人员可以理解实现上述实施例的全部或部分步骤可以通过硬件来完成，也可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，上述提到的存储介质可以是只读存储器，磁盘或光盘等。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种面向DCU集群的大规模有限元网格并行分区的方法，其特征在于，所述方法包括：
S1、将求解区域离散为网格，通过带权对偶图描述离散所得的网格的分布情况，其中，所述带权对偶图中的顶点表示网格，所述带权对偶图中的边表示网格间的连接关系；
S2、若当前进程编号为0，则输入带权对偶图G0，将所述带权对偶图G0压入堆栈GraphStack；若当前进程编号非0，则进入阻塞通信状态，直到完成与其他进程的通信获得输入所述带权对偶图G0及其进程编号PRE_PROC，将所述带权对偶图G0压入堆栈GraphStack；
S3、从所述堆栈GraphStack的栈顶取出图Gtop，如果图Gtop的划分深度达到预期值，则将所述图Gtop压入堆栈SubGraphStack；如果图Gtop的划分深度未达到预期值，则使用谱二分法将图Gtop划分为子图G1、G2，子图的划分深度为Gtop→depth+1；
设定同时运行的进程数量为P，根据当前进程编号得到通信目的进程的编号，如果所述通信目的进程的编号NEXT_PROC未达到最大进程数量P，则与进程NEXT_PROC进行通信，发送所述子图G2，将所述进程NEXT_PROC压入堆栈ProcStack中，将所述子图G1压入堆栈GraphStack中；如果所述通信目的进程的编号达到最大进程数量，则将子图G1、G2都压入堆栈GraphStack中；
重复执行S3，直到所述堆栈GraphStack为空时，执行S4；
S4、当所述堆栈ProcStack非空时，从所述堆栈ProcStack栈顶取出编号Ptop，进入阻塞通信状态直到完成与进程Ptop的通信，接收得到进程Ptop返回的子图集合GraphSet，将集合中的子图压入所述堆栈SubGraphStack中；
重复执行S4，直到ProcStack为空时，执行S5；
S5、使用所述堆栈SubGraphStack中的所有子图创建GraphSet，若当前进程为非0号进程，则进入阻塞通信状态，当完成与进程PRE_PROC的通信时，发送子图集合GraghSet；若当前进程为0号进程，则为所述GraghSet中的子图进行分区编号，输出网格与分区编号的对应结果。
2.根据权利要求1所述的方法，其特征在于，所述S3中的所述使用谱二分法将图Gtop划分为子图G1、G2，包括下述S31-S33：
S31、计算得到所述图Gtop的拉普拉斯矩阵A；
S32、使用Lanczos迭代方法，将所述拉普拉斯矩阵A转换为三对角矩阵T；
S33、计算所述三对角矩阵T的费德勒向量，根据各分量的正负将所述图Gtop中对应的顶点划分为两个集合，得到子图G1、G2。
3.根据权利要求2所述的方法，其特征在于，所述S32的使用Lanczos迭代方法，将所述拉普拉斯矩阵A转换为三对角矩阵T，包括：
Lanczos迭代的数学定义如下：
若A为对称矩阵，则存在正交矩阵Q＝[q1,q2,…,qn]使得T＝QTAQ，其中qi为N维列向量，T为三对角阵；
给定任意列向量q1(||q1||2＝1)，通过迭代得到正交矩阵Q以及三对角阵T；
其中，迭代格式如下：
①、μ1＝Aq1；i＝1
②、ri＝μi-αiqi；βi＝||ri||2
③、若βi＝0或达到迭代次数限制则退出，否则：
qi+1＝ri/βi；μi+1＝Aqi+1-βiqi；i＝i+1
转到执行②；
其中，所述Aq1为稀疏矩阵向量乘法。
4.根据权利要求3所述的方法，其特征在于，所述稀疏矩阵向量乘法基于CSR格式存储；
所述稀疏矩阵向量乘法采用HIP编程的方式进行并行优化。
5.根据权利要求3所述的方法，其特征在于，设定稀疏矩阵的尺寸为NxN,向量的尺寸为Nx1，所述稀疏矩阵中所有非零元素的个数为NotZeroNums，设定数组IA[N+1]、JA[NotZeroNums]、VA[NotZeroNums]存储稀疏矩阵，设定数组x[N]存储向量；
所述稀疏矩阵向量乘法采用HIP编程的方式进行并行优化，包括下述步骤S321-S323：
S321、对于向量规模小于4096的稀疏矩阵向量乘法，不采取辅助优化技术；
S322、引入中间数组v[NotZeroNums],将矩阵向量乘运算分解为乘运算与加运算两个核心执行；
在计算核心一中完成乘运算部分，设置每个Block中的线程数量为threadsPerBlock，Block数量为ceil(NotZeroNums/threadsPerBlock),ceil为向上取整函数，计算核心一的步骤包括下述S3221-S3224：
S3221、计算线程全局编号：
tid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x；
S3222、计算线程总数icr＝hipBlockDim_x*hipGridDim_x；
S3223、根据tid取出VA[tid]、JA[tid]、x[JA[tid]]，计算VA[tid]*x[JA[tid]]存入v[tid]；
S3224、计算tid＝tid+icr，重复执行S73-S74，直到tid超过NotZeroNums，则执行S8；
S323、在计算核心二中完成加运算部分，设定Block数量为ceil(N/threadsPerBlock),ceil为向上取整函数，计算核心二的步骤包括下述S81-S89：
S3231、计算线程全局编号：
gid＝hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x
S3232、取得线程局部编号tid＝hipThreadIdx_x；
S3233、设定共享内存空间ptr_s[threadsPerBlock+1]、v_s[sizeSharedMemory]；
S3234、根据编号从数组IA中取得元素下标范围，存入ptr_s中：
ptr_s[tid]＝IA[gid]
S3235、若tid为0，即该线程是Block内的0号线程，则计算：
MaxIndex＝gid+threadsPerBlock
其中，MaxIndex为所述Block负责的最大行编号，若MaxIndex大于N，则取MaxIndex为N；取出最大元素下标IA[MaxIndex]存入ptr_s[threadsPerBlock]中；
S3236、计算变量
temp＝(ptr_s[threadsPerBlock]-ptr_s[0])/threadsPerBlock+1；nlen＝min(temp*threadsPerBlock,sizeSharedMemory)；
得到一次循环过程中加载的元素数量nlen，设定循环变量i的初始值为ptr_s[0]，中间结果sum的初始值为0；
S3237、计算index＝i+tid，以v_s[tid]、v[index]为起点，threadsPerBlock为步长，从v数组中加载数据到v_s数组；
加载完成后，若ptr_s[tid]<＝i+nlen-1且ptr_s[tid+1]>i,即加载到共享内存空间中的元素属于第gid行，则计算元素下标的起始范围：
row_st＝max(ptr_s[tid]–i,0)
row_ed＝min(ptr_s[tid+1]–i,nlen)
累加v_s[row_st]至v_s[row_ed-1]的所有元素到变量sum上；
S3238、i＝i+nlen；重复执行S87直到i达到ptr_s[threadsPerBlock]时，执行S89；
S3239、将sum存入y[gid]。
6.根据权利要求1所述的方法，其特征在于，所述S3中的根据当前进程编号得到通信目的进程的编号，包括：
设当前进程编号为i，则通过下述公式(1)计算得到通信目的进程的编号：
7.一种面向DCU集群的大规模有限元网格并行分区的装置，其特征在于，所述装置包括：
离散单元，被配置为执行S1：将求解区域离散为网格，通过带权对偶图描述离散所得的网格的分布情况，其中，所述带权对偶图中的顶点表示网格，所述带权对偶图中的边表示网格间的连接关系；
压栈单元，被配置为执行S2：若当前进程编号为0，则输入带权对偶图G0，将所述带权对偶图G0压入堆栈GraphStack；若当前进程编号非0，则进入阻塞通信状态，直到完成与其他进程的通信获得输入所述带权对偶图G0及其进程编号PRE_PROC，将所述带权对偶图G0压入堆栈GraphStack；
处理单元，被配置为执行S3：从所述堆栈GraphStack的栈顶取出图Gtop，如果图Gtop的划分深度达到预期值，则将所述图Gtop压入堆栈SubGraphStack；如果图Gtop的划分深度未达到预期值，则使用谱二分法将图Gtop划分为子图G1、G2，子图的划分深度为Gtop→depth+1；
设定同时运行的进程数量为P，根据当前进程编号得到通信目的进程的编号，如果所述通信目的进程的编号NEXT_PROC未达到最大进程数量P，则与进程NEXT_PROC进行通信，发送所述子图G2，将所述进程NEXT_PROC压入堆栈ProcStack中，将所述子图G1压入堆栈GraphStack中；如果所述通信目的进程的编号达到最大进程数量，则将子图G1、G2都压入堆栈GraphStack中；
重复执行S3，直到所述堆栈GraphStack为空时，执行S4；
接收单元，被配置为执行S4：当所述堆栈ProcStack非空时，从所述堆栈ProcStack栈顶取出编号Ptop，进入阻塞通信状态直到完成与进程Ptop的通信，接收得到进程Ptop返回的子图集合GraphSet，将集合中的子图压入所述堆栈SubGraphStack中；
重复执行S4，直到ProcStack为空时，执行S5；
输出单元，被配置为执行S5：使用所述堆栈SubGraphStack中的所有子图创建GraphSet，若当前进程为非0号进程，则进入阻塞通信状态，当完成与进程PRE_PROC的通信时，发送子图集合GraghSet；若当前进程为0号进程，则为所述GraghSet中的子图进行分区编号，输出网格与分区编号的对应结果。
8.根据权利要求7所述的装置，其特征在于，所述被配置为执行S3的处理单元进一步被配置为执行下述S31-S33：
S31、计算得到所述图Gtop的拉普拉斯矩阵A；
S32、使用Lanczos迭代方法，将所述拉普拉斯矩阵A转换为三对角矩阵T；
S33、计算所述三对角矩阵T的费德勒向量，根据各分量的正负将所述图Gtop中对应的顶点划分为两个集合，得到子图g1、G2。
9.一种电子设备，其特征在于，所述电子设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由所述处理器加载并执行以实现如权利要求1至6任一所述的一种面向DCU集群的大规模有限元网格并行分区的方法。
10.一种计算机可读存储介质，其特征在于，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现如权利要求1至6任一所述的一种面向DCU集群的大规模有限元网格并行分区的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110541398.XA CN113239591B (zh) | 2021-05-18 | 2021-05-18 | 面向dcu集群的大规模有限元网格并行分区的方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110541398.XA CN113239591B (zh) | 2021-05-18 | 2021-05-18 | 面向dcu集群的大规模有限元网格并行分区的方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113239591A CN113239591A (zh) | 2021-08-10 |
| CN113239591B true CN113239591B (zh) | 2023-10-27 |
Family
ID=77135160
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110541398.XA Active CN113239591B (zh) | 2021-05-18 | 2021-05-18 | 面向dcu集群的大规模有限元网格并行分区的方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113239591B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115938494B (zh) * | 2022-11-24 | 2024-01-09 | 中国科学院大气物理研究所 | 气相化学模块的dcu加速计算方法、设备及存储介质 |
| CN115618827B (zh) * | 2022-12-20 | 2023-03-10 | 西安葡萄城软件有限公司 | 防止电子表格系统堆栈溢出的计算方法及装置 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6718290B1 (en) * | 1998-12-10 | 2004-04-06 | Georgia Tech Research Corporation | Systems and methods for encoding tetrahedral meshes |
| CN104363650A (zh) * | 2014-09-19 | 2015-02-18 | 西北大学 | 一种野外条件下无线传感器网络定位优化方法 |
| CN104461467A (zh) * | 2013-09-25 | 2015-03-25 | 广州中国科学院软件应用技术研究所 | 针对SMP集群系统采用MPI和OpenMP混合并行提高计算速度的方法 |
| CN106250102A (zh) * | 2015-06-12 | 2016-12-21 | 中国石油化工股份有限公司 | 交错网格有限差分正演模拟优化的方法 |
| CN109960739A (zh) * | 2019-03-27 | 2019-07-02 | 中南大学 | 一种基于子图生长的图分割方法 |
| CN111125949A (zh) * | 2019-12-06 | 2020-05-08 | 北京科技大学 | 一种有限元分析的大规模并行网格划分系统及方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20180342030A1 (en) * | 2017-05-24 | 2018-11-29 | The Research Foundation For The State University Of New York | Neutral radistricting using a multi-level weighted graph partitioning algorithm |
-
2021
- 2021-05-18 CN CN202110541398.XA patent/CN113239591B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6718290B1 (en) * | 1998-12-10 | 2004-04-06 | Georgia Tech Research Corporation | Systems and methods for encoding tetrahedral meshes |
| CN104461467A (zh) * | 2013-09-25 | 2015-03-25 | 广州中国科学院软件应用技术研究所 | 针对SMP集群系统采用MPI和OpenMP混合并行提高计算速度的方法 |
| CN104363650A (zh) * | 2014-09-19 | 2015-02-18 | 西北大学 | 一种野外条件下无线传感器网络定位优化方法 |
| CN106250102A (zh) * | 2015-06-12 | 2016-12-21 | 中国石油化工股份有限公司 | 交错网格有限差分正演模拟优化的方法 |
| CN109960739A (zh) * | 2019-03-27 | 2019-07-02 | 中南大学 | 一种基于子图生长的图分割方法 |
| CN111125949A (zh) * | 2019-12-06 | 2020-05-08 | 北京科技大学 | 一种有限元分析的大规模并行网格划分系统及方法 |
Non-Patent Citations (1)
| Title |
|---|
| 流动数值模拟中一种并行自适应有限元算法;周春华;;计算物理(第04期);412-418 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113239591A (zh) | 2021-08-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Ballard et al. | Communication lower bounds and optimal algorithms for numerical linear algebra | |
| CN113239591B (zh) | 面向dcu集群的大规模有限元网格并行分区的方法及装置 | |
| Koanantakool et al. | Communication-avoiding parallel sparse-dense matrix-matrix multiplication | |
| Grigori et al. | Parallel symbolic factorization for sparse LU with static pivoting | |
| JP6784780B2 (ja) | 大規模再生可能エネルギーのデータについて確率モデルを構築する方法 | |
| CN109726441B (zh) | 体和面混合gpu并行的计算电磁学dgtd方法 | |
| US20200159810A1 (en) | Partitioning sparse matrices based on sparse matrix representations for crossbar-based architectures | |
| WO2018027706A1 (zh) | Fft处理器及运算方法 | |
| Margaris et al. | Parallel implementations of the jacobi linear algebraic systems solve | |
| CN116258042B (zh) | 一种基于ddm的大规模传热异构并行仿真方法 | |
| Danalis et al. | Efficient quality threshold clustering for parallel architectures | |
| CN107247686B (zh) | 一种基于并行算法的fetd仿真模拟方法 | |
| Hussain et al. | Communication-avoiding and memory-constrained sparse matrix-matrix multiplication at extreme scale | |
| Zayer et al. | Sparse matrix assembly on the GPU through multiplication patterns | |
| Blelloch et al. | Improved Parallel Cache-Oblivious Algorithms for Dynamic Programming [Extend Abstract] | |
| Kopysov et al. | Hybrid Multi-GPU solver based on Schur complement method | |
| Jamal et al. | A hybrid CPU/GPU approach for the parallel algebraic recursive multilevel solver pARMS | |
| Wang et al. | Accelerating on-line training of LS-SVM with run-time reconfiguration | |
| CN109948253B (zh) | 薄板无网格Galerkin结构模态分析的GPU加速方法 | |
| Kumar et al. | Multi-threaded nested filtering factorization preconditioner | |
| Abu-Sufah et al. | On implementing sparse matrix multi-vector multiplication on GPUs | |
| Alberdi-Rodriguez et al. | Recent memory and performance improvements in Octopus code | |
| Boman et al. | A nested dissection partitioning method for parallel sparse matrix-vector multiplication | |
| Reisner et al. | Scalable line and plane relaxation in a parallel structured multigrid solver | |
| Bienz | Reducing communication in sparse solvers |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |