第３５卷第３期 物探化探计算技术 ２０１３年５月 文章编号：１ＯＯ１—１７４９（２０１３）０３—０２９３—０４ 稀疏矩阵快速回代的Ｃｈｏｌｅｓｋｙ分解法 宋 滔，王绪本 （成都理工大学摘地球物理学院，成都６１００５９） 要：采用一维压缩存储正演计算中的对称稀疏矩阵，进行Ｃｈｏｌｅｓｋｙ分解，利用分解后二个 矩阵的对称性和稀疏性，对占用时间较多的回代过程采用先消去列的方法，实现快速回代。算例 表明，采用该方法，对于点源场的求解与传统顺代回代求解法对比可以提高五倍的速度，对于大 地电磁的正演问题，提高了二倍的速度。 关键词：对称稀疏矩阵；Ｃｈｏｌｅｓｋｙ分解；快速回代 中图分类号：Ｏ １５１．２１ 文献标识码：Ａ ＤＯＩ： ０前言 在使用有限单元法进行地球物理的数值模拟 时，最终需要求解一个大型稀疏线型方程组Ａｘ— Ｂ，而这个线性方程组往往是对称正定的。在解点 源问题时，如果采用齐次边界条件，那么矩阵Ａ只 和地下电性参数有关，不包含源信息，这个时候只 有列向量Ｂ包含有源信息。如果使用高斯消元法 求解，那么对于每一个电源点都要进行消元和迭 ５１ 代，这大大增加了计算量。阮百尧Ⅲ使用Ｃｈｏｌ— ｅｓｋｙ分解法，首先将矩阵Ａ分解为两个对称矩阵， 对于不同的列向量Ｂ只需要进行一次顺代和一次 回代便可以求解方程组。但是通过对计算过程的 ＿１ ） （２） ２ Ｃｈｏｌｅｓｋｙ分解法 矩阵Ａ是正定对称的，那么可以使用Ｃｈｏｌ— ｅｓｋｙ分解法进行计算，即可以得到 Ａ＝ＬＬ 分析发现，其算法的回代过程几乎占用了整个求解 过程时间的一半，没有充分利用分解后矩阵的对称 性和稀疏性。在此基础上，作者对该算法进行了一 定的改进。 其中 １ 大型对称稀疏矩阵的变带宽存储 一ｒ Ｚ１ Ｏ Ｏ ＩＬＩ ｚ２ Ｚ２２０ 对于变带宽的对称稀疏矩阵，采用二个一维数 组来存储Ⅲ，用一个一维数组Ｇ存储矩阵的下三 ，一ｌ Ｊ／ ２ … 基金项目：国家高技术研究发展计划（８６３计划）资助（２Ｏ０９ＡＡ０６Ｚ１０８） 收稿日期：２Ｏ１２一Ｏ９—１９ 改回日期：２Ｏ１２—１２—１６ ２９４ ｉ １ 物探化探计算技术 ３５卷 ｉ＝ｊ （４） 高的 如下矩阵 所示。 ，广ｌ （ａｉｊ－－ ｌｉｋｌｊｋ） ∥＞ Ｉ ２ １ Ｌ—Ｊ ０ ０ １ １ ２ ｌ ０ Ｉｏ， ＜ 从式（３）与式（４）中可以看出，下三角矩阵Ｌ 中的元素对应矩阵Ａ中下三角的元素，所以分解 Ｏ （８） １ １ １ ０ １ ２ Ｌ １ １ １ ， ２ １ ０ ０ Ｏ ］●● ，● ０ ０ 的结果矩阵Ｌ可以使用存储Ａ的二个一维数组来 存储索引。由于 是Ｌ的转置，所以 中的元 素也可以使用这二个一维数组来访问。 由于Ａ是稀疏的，含有大量的零元素，所以Ｌ 和Ｌ 中也含有大量零元素，在顺代和回代的过程 中应该充分利用这一特性。 上式（２）中矩阵Ａ分解得到： Ｌ—Ｉ 『 Ｌ —．．．．．．．．．．．．．．．．．．．．。．。．．．一．． ．．．．． Ｇ＝＝＝（１ ２ ２ １ Ｃ ＩＤ一（１，３，４，７，９） 对于线型方程组Ａｘ—Ｂ，即ＬＬＴ３２＝＝＝Ｂ，进行顺 代和回代便可求解。］ ３顺代和回代 顺代的求解过程很简单，对于Ｌｙ—Ｂ，Ｌ是下 三角矩阵，只需要依次求解Ｙ 、Ｙ 、…、Ｙ 即可。下 一步是求解Ｌ 一 ，由于Ｌ 是Ｌ的转置，同时共 用一个存储空间： ０ ０ ０］ Ｉｚ Ｌ＝ ０ ０ １ （６） ‘’‘ ‘‘‘ ’Ｉ Ｚ ２ … Ｚ ｌ 那么 厂Ｚ１１ Ｚ２１ ０ Ｚ １＿ Ｌｒ￣－－ｉＬＬ ０：０ ０１２ｚ … Ｊ ０ ｌ （７） Ｌ 是上三角矩阵，通过回代便可完成方程求 解。方法（１）：按行依次求解ｚ 、ｚ 一 、…、 ，矩阵 元素均可通过矩阵Ｌ的索引方式得到；方法（２）： 回带的同时，消去方程组中包含当前求解的ＪＣ 。 这两个方法在Ｌ非稀疏时，计算时间是一样的，但 是当Ｌ是稀疏矩阵时，方法（２）的计算效率是非常 Ｌ ２ ） ０ １ １ ２ （９） １ 使用方法（２）回带时，首先解出ｚ ，然后消去 方程组中的ｚｓ，即将第５列中的元素消去，因为 ／ ５ 中的列是Ｌ中的行，在数组Ｇ中是按行存储的，所 以这一步骤是很容易实现的。同时，Ｇ是压缩存储 的，所以在此例中，消去最后一列，仅需要计算一步 （消去 中的ｌ ）。在大型稀疏矩阵中，这种消去 列的实现方法，对计算速度的提升是非常可观的。 ４ Ｃｈｏｌｅｓｋｙ分解法的求解程序 ！Ｎ整变量，输入参数，方程组的阶数。 ！Ｎ Ｐ整变量，输入参数，为系数矩阵Ａ压缩存 储的元素个数。 ！Ａ输入、输出参数，Ｎ Ｐ个元素的一维实数组， 输入时存放系数矩阵的压缩存储元素；输出时存 放。 ！Ｃｈｏ ｌｅｓｋｙ分解得到的下三角阵中变带宽内的 元素。 ！Ｂ输入、输出参数，Ｎ个元素的一维实数组。 输入时存放方程组右端的ｎ维常向量；输出时，存 放解。 ！向量。 ！ＩＤ输入参数，Ｎ个元素的一维整数组。存放 系数矩阵Ａ的各个对角线元素在压缩的一维数组 中的。 ！位置坐标。 ＳＵＢＲＯＵＴ１ＮＥ ＣＨＯＬＥＳＫＹ（Ｎ，ＮＰ，Ａ，Ｂ，ＩＤ） ＩＭＰＬＩＣＩＴ ＮＯＮＥ ＩＮＴＥＧＥＲ Ｎ，ＮＰ，ＩＤ（Ｎ） ＲＥＡＬ Ａ（ＮＰ），Ｂ（Ｎ） ＩＮＴＥＧＥＲ １０，Ｉ，ＩＭ，ＪＭ，Ｊ０，ＪＴ，Ｊ，ＫＴ，Ｋ， ＪＫ，ＩＮ，Ｋ０，ＫＭ，ＫＩ，ＮＮ，ＩＪ ！分解 １０＝１ ３期 宋滔等：稀疏矩阵快速回代的Ｃｈｏｌｅｓｋｙ分解法 ２９５ ＩＭ—ＩＤ（１） Ａ（ＩＭ）一ＳＱＲＴ（Ａ（ＩＭ）） １０＝：：ＩＭ＋１ Ｄｏ Ｉ一２．Ｎ ＩＭ—ＩＤ（Ｉ） ＪＭ—Ｉ—ＩＭ＋ＩＯ一１ ＩＦ（ＪＭ．ＬＴ．１）ＴＨＥＮ Ｊｏ＝＝１ ＥＬＳＥ Ｊ０＝ＩＤ（ＪＭ）＋１ ＥＮＤ ＩＦ ＤＯ ＪＴ—１０，ＩＭ一１ Ｊ＝＝：Ｉ—ＩＭ＋ＪＴ ＪＭ—ＩＤ（Ｊ） Ｄｏ ＫＴ—１０．ＪＴ一１ Ｋ—Ｉ—ｌＭ＋ＫＴ ＪＫ—ＪＭ一（Ｊ—Ｋ） ＩＦ（ＪＫ．ＬＴ．Ｊ０）ＣＹＣＬＥ Ａ（ＪＴ）一Ａ（ＪＴ）一Ａ（Ｋ］１） Ａ（ＪＫ） ＥＮＤ ＤＯ Ｊｏ＝ＪＭ＋１ Ａ（ＪＴ）一Ａ（ＪＴ）／Ａ（ＪＭ） ＥＮＤ ＤＯ Ｄｏ Ｋ—１０。ＩＭ一１ Ａ（ＩＭ）一Ａ（ＩＭ）一Ａ（Ｋ） Ａ（Ｋ） ＥＮＤ Ｄｏ Ａ（ＩＭ）一ＳＱＲＴ（Ａ（ＩＭ）） １０—１Ｍ＋１ ＥＮＤ Ｄｏ ！顺代 １０＝１ Ｄｏ Ｉ一１．Ｎ ＩＭ—ＩＤ（Ｉ） Ｊ一１ Ｄ０ ＫＴ—ＩＭ一１，ＩＯ，一１ Ｂ（Ｉ）一Ｂ（Ｉ）一Ａ（ＫＴ） Ｂ（Ｉ—Ｊ） Ｊ＝Ｊ＋１ ＥＮＤ Ｄｏ Ｂ（Ｉ）一Ｂ（Ｉ）／Ａ（ＩＭ） １０＝ＩＭ＋１ ＥＮＤ ＤＯ ！回代 ＩＮ—ＩＤ（Ｎ） Ｂ（Ｎ）＝＝：Ｂ（Ｎ）／Ａ（ＩＮ） ＤＯ ｉ＿－＿Ｎ，２，一１ ＮＮ：＝：ＩＤ（Ｉ）一ＩＤ（１—１）一１ Ｄ０ Ｊ：：＝１．ＮＮ ＩＪ—Ｉ—Ｊ Ｂ（Ｕ）一Ｂ（ＩＪ）一Ａ（Ⅱ×１）一Ｊ） Ｂ（Ｉ） ＥＮＤ ＤＯ Ｂ（Ｉ—１）一Ｂ（Ｉ一１）／Ａ（ＩＤ（１—１）） ＥＮＤ ＤＯ ＲＥＴＵＲＮ ＥＮＤ ５算例分析 作者将本算法分别用于点源二维电场和二维 大地电磁的计算，与传统的ｃｈｏｌｅｓｋｙ分解算法还 有不带平方根的ｃｈｏｌｅｓｋｙ算法进行对比。在以下 计算分析中，所使用的计算机为ＣＰＵ：Ｔ６６００， ２．２ ＧＨｚ，内存２ Ｇ。 ５．１点源二维电场的计算 作者采用矩形网格剖分，电导率分块均匀变 化，双线性插值有限元进行模拟。ｚ方向（水平沿测 线方向）包括１００个测区网格，两边分别有１２个稀 疏网格；Ｙ方向（垂直方向）一共设置２０个网格，点 源点移动１００次，即分解后，进行１００次回代和顺 带。地下空间设置为均匀半空间，分别使用本方法 和ｃｈｏｌｅｓｋｙ分解法求解，如下页表１所示。 从计算结果中可以看出，快速回代的ｃｈｏｌｅｓｋｙ 方法对于点电源场求解速度的提升是非常大的，与 传统的顺带回代方法相对比，提升了五倍左右，而 且单元数越多，其提升越明显。 ５．２二维大地电磁正演模拟 地面测线长４ ｋｍ，测点间距１００ １ＴＩ，频率范围 为１ ０００ Ｈｚ～Ｏ．０１ Ｈｚ，对数等间隔采样，一共４１ 个频点。 作者采用矩形网格剖分，电导率分块均匀变 化，双二次插值有限元进行模拟。ｚ方向（水平沿测 线方向）包括４１个测区网格，两边分别有１８个稀 疏网格；Ｙ方向（垂直方向）一共设置５６个网格， ＴＥ模式下空气网格设置１４个。地下空间设置为 均匀半空间，分别使用本方法和不带平方根的 ｃｈｏｌｅｓｋｙ分解法求解，如下页表２、表３所示。 从表２、表３可以看出，使用快速回代的ｃｈｏｌｅｓｋｙ 算法的速度，要比传统的ｃｈｏｌｅｓｋｙ分解法 借。 ６结论 通过算例分析，快速回代的方法充分利用了 ２９６ 物探化探计算技术 表１点源场求解时间对比 Ｔａｂ．１ Ｔｈｅ ｓｏｌｕｔｉｏｎ ｔｉｍｅ ｏｆ ｐｏｉｎｔ ｓｏｕｒｃｅ ｆｉｅｌｄ ３５卷 求解方法 快速回代的ｃｈｏｌｅｓｋｙ Ｃｈｏｌｅｓｋｙ分解 单元数 ２４８０ ２４８０ 节点数 ２６２５ ２６２５ 回代次数 １００ １Ｏ０ 平均计算时间／ｓ ０．０１２ ０．０６６ 总时间／ｓ １．２１ ６．６４ 表２ ＴＥ模式求解时间对比 Ｔａｂ．２ Ｔｈｅ ｓｏｌｕｔｉｏｎ ｔｉｍｅ ｏｆ ＴＥ ｍｏｄｅ ＴＥ模式 快速回代的ｃｈｏｌｅｓｋｙ 单元数 ５３９０ 节点数 １６４６５ 频点数 ４１ 平均计算时间／ｓ ２．３ 总时间／ｓ ９３．８ Ｃｈｏｌｅｓｋｙ分解 ５３９０ １６４６５ ４１ ４．７ １９３．５ 表３ ＴＭ模式求解时间对比 Ｔａｂ．３ Ｔｈｅ ｓｏｌｕｔｉｏｎ ｔｉｍｅ ｏｆ ＴＭ ｍｏｄｅ ＴＭ模式 快速回代的ｃｈｏｌｅｓｋｙ Ｃｈｏｌｅｓｋｙ分解 单元数 ４３１２ ４３１２ 节点数 １３２０３ １３２０３ 频点数 ４１ ４１ 平均计算时间／ｓ １．２ ２．５ 总时间／ｓ ４９．７ １０１．９ Ｃｈｏｌｅｓｋｙ分解后矩阵的对称性和稀疏性，大大加 元法正演模拟ＥＪ］．地球物理学报，１９９７，４Ｏ（３）：４２１． Ｅ８３刘云，王绪本．大地电磁二维自适应地形有限元正演 模拟ＥＪ３．地震地质，２０１０，３２（３）：３８２． 快了整体求解的速度。该方法对于使用有限元进 行点源场和大地电磁模拟，所形成的大型对称稀疏 线性方程组的求解是有效的。 参考文献： ［１］阮百尧，熊彬．大型对称变带宽方程组的Ｃｈｏｌｅｓｋｙ分 ［９］徐世浙．地球物理中的有限单元法［Ｍ］．北京：科学 出版社，１９９４． ［１Ｏ］罗延钟，张桂青．电子计算机在电法勘探中的应用 ＥＭ３．武汉：武汉地质学院出版社，１９８７． ［１１］朴化荣．电磁测深法原理［Ｍ］．北京：地质出版社， １９９０． 解法［Ｊ］．物探化探计算技术，２０００，２２（４）：３６１． Ｅ２］刘德贵，费景高，于江永，等．Ｆ（）ＲＴＲＡＮ算法汇编第 一分册ＦＭ］．北京：国防工业出版社，１９８０． 国电力出版社，２００２． ［１２］曾国．大地电磁二维有限元正演数值模拟［Ｄ］．湖 南：中南大学，２００８． Ｅ１３３王绪本，李永年，高永才．大地电磁测深二维地形影 ［３］彭国伦．Ｆｏｒｔｒａｎ ９５程序设计（第一版）ＥＭ３．北京：中 Ｅ４３徐士良．计算机常用算法（第二版）［Ｍ］．北京：清华大 学出版社，１９９５． 响及其校正方法研究［Ｊ］．物探化探计算技术，１９９９， ２１（４）：３２７． ［５］阮百尧，徐世浙．电导率分块线性变化二维地电断面 ［１４］徐世浙．点电源二维电场问题的付氏变换的波数ｋ 的选择口］．物探化探计算技术，１９８８，１０（３）：２３５． Ｅ１５３强建科，罗延钟．三维地形直流电阻率有限元法模拟 ［Ｊ］．地球物理学报，２００７，５０（５）：１６０６． 电阻率测深有限元数值模拟ＥＪ３．地球科学一中国地 质大学学报，１９９８，２３（３）：３０３． ［６］熊彬，阮百尧．电位双二次变化二维地电断面电阻率 测深有限元数值模拟［Ｊ］．地球物理学报，２００２，４５ （２）：２８５． 作者简介：宋滔（１９８８一），男（汉），四川什邡人，硕 士，研究方向为地球物理电磁法的数值模拟。 Ｅ７３史明娟，徐世浙，刘斌．大地电磁二次函数插值的有限
