(19)国家知识产权局 (12)发明 专利申请 (10)申请公布号 (43)申请公布日 (21)申请 号 202210522339.2 (22)申请日 2022.05.13 (71)申请人 吉林大学 地址 130012 吉林省长 春市前进大街269 9 号 (72)发明人 张雪松　闫昭　赫枫龄　雷新丽　 马琳　 (74)专利代理 机构 北京京万通知识产权代理有 限公司 1 1440 专利代理师 许天易 (51)Int.Cl. G06F 9/50(2006.01) G06F 12/0877(2016.01) G06T 1/20(2006.01) (54)发明名称 一种用于仿真系统的基于GPU的中小型稠 密 矩阵的批量并行LU分解方法 (57)摘要 本发明针对仿真系统的小中型稠密矩阵实 现一种基于GPU的高性能批量并行LU分解方法， 该方法包括： S1、 计算仿真系统的模型单元在每 个时间步的雅可比矩阵， 并将这些矩阵作为子矩 阵， 并拼接为大矩阵H； S2、 对每个子矩阵内部分 块， 并按照分块的大小为每个子矩阵分配GPU线 程组； S3、 所述 GPU线程组按照对角顺序从左上到 右下， 依次完成每个对角块内所有列的主元选 取、 对角块所在行右侧块的更新、 对角块内部列 与主元列交换， 以及对角块所在列下方和右下方 的所有分块数据更新， 从而完成对大矩阵H内的 所有子矩阵的LU分解。 本发明能随着批量数量的 增加， 性能以近似线性的趋势增加， 并且在隐式 求解包含大批量小型常微分方程组的模型时， 能 够加快求 解速度， 有效缩短仿真时间。 权利要求书2页 说明书7页 附图4页 CN 114911619 A 2022.08.16 CN 114911619 A 1.一种用于仿真系统的基于GPU的中小型稠密矩阵的批量LU分解方法， 其特征在于， 包 括： S1、 计算仿真系统 的模型单元在每个时间步的雅可比矩阵， 并将这些矩阵作为子矩阵， 并拼接为大矩阵H； S2、 对每个子矩阵内部分块， 并按照分块的大小为每 个子矩阵分配GPU 线程组； S3、 所述GPU线程组按照对角顺序从左上到右下， 依次完成每个对角块 内所有列的主元 选取、 对角块所在行右侧块的更新、 对角块内部列与主 元列交换， 以及 对角块所在列下方和 右下方的所有分块数据更新， 从而完成对大矩阵H内的所有子矩阵的LU分解。 2.根据权利要求1所述的批量LU分解方法， 其特征在于， 在步骤S1中， 所述仿真系统具 有n个模型单元， 所述雅可比矩阵为m阶中小型稠密矩阵Hi， i的范围是[0， n ‑1]， 其中， m的获 取方法为： 将当前待求解的所有模型单元的维度向上扩充到2的整数次幂m(即雅可比矩阵 扩展到m阶)， 对于增 加的行和列， 除了对角线元 素的值为1外， 其 他全部用0填充。 3.根据权利要求1所述的批量 LU分解方法， 其特 征在于， 在步骤S2中， 包括： S21、 将子矩阵Hi分为若干子块， 以保证后续的处理 中， 线程组能够一次性从全局内存读 入整个分块内的数据， 减少全局内存的访问次数； S22、 按照分块的大小为子矩阵Hi分配对应数量的GPU 线程， 构成线程组Groupi。 4.根据权利要求3所述的批量LU分解方法， 其特征在于， 将子矩阵Hi分为若干子块的方 法包括： 如果m大于等于16时， 设置块大小bl ocksize为8或16； 如果m小于16， 则将整个子矩阵作为 一块， 此时块大小bl ocksize＝m。 5.根据权利要求2所述的批量 LU分解方法， 其特 征在于， 在步骤S3中， 具体步骤 包括： S31、 每个线程组Groupi选定子矩阵Hi内索引为j的对角块Bj,j， 确定Bj,j内所有列对应的 主元列， 并对位于Bj,j右侧所在行 上的所有分块进行 更新； 其中， j初始化 为0； S32、 将得到的所有 主元列与对角块Bj,j内部对应的列互换位置； S33、 对位于对角块Bj,j所在列下 方的所有分块进行 更新； S34、 对位于对角块Bj,j右下方的所有分块进行 更新； S35、 如果对角块Bj,j是子矩阵Hi中的右下角分块， 批量并行LU分解结束， 输出大矩阵H， 其中， 每个子矩阵Hi都由Li和Ui构成， 分别是单位下三角矩阵和上三角矩阵， 用于后续的线 性方程组快速求 解； 否则， 令j＝j+1， 跳转到步骤S31， 继续处 理下一个分块。 6.根据权利要求5所述的批量 LU分解方法， 其特 征在于， 步骤S31包括： S311： 初始化维度为m的一维矩阵Toi， 令Toi[l]＝m,0 ≤l<m， 其中， 数组索引l对应子矩阵Hi的列序号， Toi[l]用于记录序号为l的列， 其最终的主元列序 号； S312： 线程组从当前块Bj,j开始， 并行读取块 内第p行的所有元素， 存入线程 内部的缓存 中； 其中， p按照子矩阵Hi的行索引计算， 初始化为j ×blocksize； 在此步骤中， 由于块Bj,j的 行元素之间在全局内存中紧密排列， GPU仅需发射一次指令即可实现线程组的全部数据读 取， 达到了内存合并访问的效果； 其中， j ×blocksize≤p<(j+1) ×blocksize)； S313： 线程组继续并行读取Bj,j右侧相邻的下一分块中p行所有的数据， 并与缓存中的 数据比较， 将具有较大绝对值的元素及其所在列序号保存在内部缓存中； 此过程持续迭代，权　利　要　求　书 1/2 页 2 CN 114911619 A 2直到子矩阵最后一个块数据被读取； S314： 比较线程组内缓存的所有元素， 取绝对值最大的元素所对应的列q作为列p对应 的主元列， 更新元 素Toi[p]＝q； S315： 对子矩阵Hi中的当前行p， 从列索引j ×blocksiz e开始， 遍历行中后续的所有元素 如果Toi[l]>q， 则执 行如下运 算： 其中， 为元素 对应的真正主元； S316： 对于P1区域的每 个元素 的列索引y， 如果Toi[y]>q， 则执 行如下消元运 算： 其中， P1区域为位于子矩阵Hi的[p+1,(j+1) ×blocksize ‑1]行和[j×blocksize,m ‑1] 列之间的矩阵； p+1≤x<(j+1) ×blocksize,j ×blocksize≤y<m； S317： 如果p是块Bj,j内的最后一行， 步骤S31结束， 此时Bj,j中包含了当前分块内的LU分 解结果Lj,j和Uj,j； 否则， 令p＝p+1， 跳转到步骤S312， 继续处 理下一行。 7.根据权利要求5所述的批量LU分解方法， 其特征在于， 在步骤S32中， 对于Bj,j内的任 意一列p， 线程组内的第p号线程执 行列colp与列 之间的元 素交换。 8.根据权利 要求5所述的批量LU分解方法， 其特征在于， 在步骤S33中， 设P2区域为位于 子矩阵Hi的[(j+1) ×blocksize,m ‑1]行和[j ×blocksize,(j+1) ×blocksize ‑1]列之间的 矩阵， 对于其中的每 个元素 执行如下运 算： 其中， (j+1) ×blocksize≤x<m,j ×blocksize≤y<(j+1) ×blocksize。 9.根据权利 要求5所述的批量LU分解方法， 其特征在于， 在步骤S34中， 设P3区域为位于 子矩阵Hi的[(j+1) ×blocksize,m ‑1]行和[(j+1) ×blocksize,m ‑1]列之间的矩阵； P1和P3 区域依据blocksize进一步划分为多个子区域， 线程组Groupi依次对其中的任意子区域P1t、 P3t,0≤t<k‑j， 进行如下处 理： 线程组Groupi中的每个GPU线程Ty,0≤y<blocksize， 执 行如下运 算： 读入子区域P1t中序号为的y列向量A到线程的缓存中； 计算子区域P3t中序号为的y列向量D中的每 个元素： D[x]＝D[x]‑CxA   (4) 其中， (j+1) ×blocksize≤x<m。 10.根据权利要求5所述的批量LU分解方法， 其特征在于， 在步骤S35中， 判断当前处理 的分块Bj,j是否子矩阵Hi中的右下角分块， 如果是， 批量并行LU分解结束； 否则， 令j＝j+1， 跳转到步骤S31， 继续处 理下一个分块。权　利　要　求　书 2/2 页 3 CN 114911619 A 3