GeoFEA的直接和迭代线性求解器
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
装，消去过程便可执行，而这时方程右端的 aij和 bi
并不需要完全组装。那么，问题是：该如何判断方 程（2）右端的三项乘积中所有三个矩阵元素都被完 全组装。在波前求解进行之前，程序必须确定两个
中国土木工程学会第十届土力学及岩土工程学术会议论文集
信息： （1）整个求解域中单元的编号； （2）按当前单元编号对最后出现的结点号进行
摘 要：有限元法离散所产生的大规模线性方程组的求解是有限元模拟过程中最耗计算时间和内存的部分，因此，需要使用 高效的线性求解器来快速精确地进行有限元模拟。GeoFEA 的线性求解器分为直接和 Krylov 子空间迭代求解器两大类。直接 求解器的特点是可以用固定计算量精确地得到线性方程系统的解，但对于涉及十万甚至上百万个自由度的大规模岩土工程问 题，现有的直接求解器所需要的内存和计算时间是难以接受的；迭代求解器的特点是内存需求量小，当结合高效的预处理子 时，对三维问题所产生大规模线性方程系统的计算时间要显著少于直接求解器。 关 键 词：GeoFEA；线性方程系统；直接求解器；迭代求解器
在 Krylov 子空间迭代法中，矩阵与矢量的乘积 构成了所需计算量的主要部分。因此，对于有限元 法或有限差分法得到的稀疏矩阵，迭代法中的矩阵 与矢量乘积要优越于直接求解法的 LU 分解。一方 面，LU 分解通常使被分解稀疏矩阵的零值元素变 为非零而产生稠密的 L，U 分解因子，从而导致难 以承受的计算量和存储量。然而，现有的岩土工程 有限元软件中，迭代求解法并没有完全取代直接求 解法主要基于下面两个原因：
{ } Kk ( A, r0 ) = span r0 , Ar0 ,K Ak−1r0 , k = 1, 2,K (3)
这样，近似解可以通过这个子空间和初始假设解 x0 根据下式产生：
xk ∈ x0 + Kkk ( A, r0 ), k = 1, 2,K
(4)
因此，Krylov 子空间迭代法的迭代框架可以概括如 下：
一般来说，一个高效的预处理子具有下面性质： (1) 对于某一 Krylov 子空间迭代法，预处理后 的线性方程系统具有更快的收敛性质，达到收敛时 所需要的迭代步数要明显低于原始线性方程系统达 到收敛时的迭代步数； (2) 每一迭代步中，预处理所带来的计算量有
限，可以用被减少的迭代步数来弥补使得最终收敛 时所需要的求解时间明显少于原始线形方程系统的 求解时间；
⎧⎨⎩ abii**j
= =
aij − aik ak−k1akj bi − aik ak−k1bk
(2)
在传统的高斯消去求解法中，方程（2）右端的 aij 和 bi 都是完全组装的，并且高斯消去或者是 LU
分解都是基于这些已经完全组装的矩阵元素上的。 而 Irons 观察到对于高斯消去过程，只要方程（2） 右端的三项乘积中所有三个矩阵元素都被完全组
3 GeoFEA 的对称 Krylov 子空间求解 器
求解线性方程系统（1）有多种方法，主要分为 定常迭代和非定常迭代两大类。而属于非定常迭代
的 Krylov 子空间迭代法最为流行。Krylov 子空间迭 代法的两个关键成分：一个是如何构筑子空间，另
一个是如何构筑近似解。因此，根据构筑方法不同， 我们可以得到不同的 Krylov 子空间迭代法（参见 Chen[7]; Chen et al[8]）。Krylov 子空间主要通过矩阵 与矢量的乘积来构筑，可以具体表示为：
对于对称正定线性方程系统，共轭剃度法是目 前最为有效的方法。但对于对称非正定的线性系统， 比较常用的是最小残量 (MINRES)方法和对称 LQ (SYMMLQ)方法，但最小残量法和对称 LQ 这两种 方法的限制是所使用的预处理子只能是对称正定 的。 Freund 发展了简化的或对称的拟最小残量方 法 (SQMR) 来求解对称非正定线性系统并可以结 合任意对称预处理子。近年来，Smith 和 Griffiths[2] 也成功地将预处理共轭剃度法用于比奥固结产生的 对称非正定线性系统，并且很少出现问题。 下式给 出来比奥固结方程离散得到的 2×2 的耦合对称非 正定线性方程
(3) 预处理子所需要的内存存储量有限。 方程（1）的预处理形式可以表达为：
PL−1 APR−1PR x = PL−1b
(6)
这里，一个广义的预处理子可以表示为 P = PL PR。 当 PL = I, PR ≠ I 表示预处理方式为右端处理；当 PL≠ I, PR =I 表示预处理方式为左端处理；当 PL ≠ I, PR ≠ I 表示为左右端处理方式。GeoFEA 提供了几种对 称预处理子，包括修正雅克比预处理子（MJ）、广 义雅克比预处理子（GJ）和修正对称连续超松弛预 处理子（MSSOR）。观察到标准雅克比预处理子在 求解比奥固结线性方程时性能变差，Phoon et al[3] 提出了广义雅克比预处理子来改善了对角预处理子
Ax
=
b，这里
A
=
⎡ ⎢ ⎢⎣
K BT
B⎤ ⎥
−C ⎥⎦ N×N
(5)
这里，N = m + n; A ∈ ℜ(m+n)×(m+n) 是 2×2 耦合矩阵; b ∈ ℜm+n；K ∈ ℜm×m，是对应于土体骨架刚度的对称 正定矩阵；C ∈ ℜn×n 是对应于流体刚度的对称半正 定矩阵； B ∈ ℜm×n 是土体骨架刚度与流体刚度的 耦合矩阵；m，n 分别对应结点位移自由度总数和超 静孔隙水压力自由度总数。
常见，其基本特点可以概括如下：
(1) 单元接单元的矩阵存储指的是有限元离散的单 元刚度矩阵不需要组装。 这样，Krylov 子空间迭 代法中的整体刚度矩阵与矢量的乘积 w = Av 可以
∑ 通过单元接单元地进行，即
Av
=
⎛ ⎜⎝
e
LTe
Ae
Le
⎞ ⎟⎠
v
，
其伪代码如下：
for i = 1 to nel do
End For Krylov 子空间迭代法之所以流行有几个很重要的 原因：一方面，预处理技术的使用可以有效地降低 病态线性方程系统的条件数，或者说，预处理技术 可以改善系数矩阵的特征值分布而使其更为聚集， 因为对称子空间迭代法的收敛性与线性方程的系数 矩阵的特征值分布紧密相关。 通过应用预处理技
术，Krylov 子空间迭代法收敛时所需要的迭代步数 要远小于线性方程系统的总自由度数。另一方面， 用有限元、有限差分和有限体积等方法离散所产生 的线性系统的刚度矩阵通常都是稀疏的，而稀疏矩 阵又尤其适合矩阵与矢量乘积的计算。而且矩阵矢 量的乘积可以在不组装整体系数或刚度矩阵的条件 下得到。 3.1 对称 Krylov 子空间迭代求解器
w = 0; for j = 1 to n do
if v(j) ≠ 0 then for k = jebe(j) to jebe(j+1) – 1 do
r = iebe(k)
w(r) = w(r) + v(j)*ebea(k) end for end if end for 这里，iebe(:)，jebe(:)和 ebea(:)是整体矩阵 A 的压缩 稀疏列（CSC）的存储格式。对于对称矩阵的上三 角或下三角的存储格式，上面的伪代码需要进行相 应调整。尽管从某种程度来讲，两种存储格式都是 稀疏存储。但对于单元接单元的存储格式，仍有大 量的重复的整体刚度矩阵位置 (i, j)的元素， 而对 于已经完全组装的完全稀疏矩阵格式，(i, j)位置的 元素是唯一的。这样对于稀疏矩阵存储的格式，矩 阵与矢量乘积的所需要的浮点运算次数要明显少于 单元接单元的存储格式。
一个线性方程系统可用下面的通用形式表示： Ax = b, 这里 A= [aij] ∈ ℜN×N , 并且 x, b ∈ ℜN (1) 求解线性方程系统（1）的各种直接求解法都是传统 高斯消去算法的变体。其中，较著名的有 Irons[6]提 出的用于求解对称有限元矩阵方程的波前直接求解 法。波前直接求解法采用单元刚度矩阵组装与高斯 消去相互交织的求解方式，基本求解原理可以概括 为
（1）对所有的耦合、非耦合问题，迭代求解法
并不能确保都收敛； （2）直接求解法能够用固定的浮点运算数得到
精确的、可靠的解。 基 于 上 面 考 虑 ， GeoFEA 既 提 供 了 波 前
（Frontal）直接求解器，也提供了多种可供选择的 对称和非对称 Krylov 子空间迭代求解器。
2 GeoFEA 的波前直接求解器
选用迭代初始解 x0，通常可设 x0 = 0； For k = 1 to max_it, Do
（a）如果必要，应用预处理子（Preconditioner）； （b）计算当前步的近似解，xk； （c）如果必要，计算当前步残量，rk = b–Axk ； （d）使用适当的收敛准则来检查收敛情况，如果 收敛，则跳出循环；
(2) 单元接单元的存储格式的一个显著的优势就是 内存需求量少。因为尽管稀疏求解法只需存储整体 刚度的非零元素，但需要组装工作空间，所以最后 所需的峰值存储量还是高于单元接单元所需要的总 存储量。如果考虑对称性质，所需存储量可缩减约 一半。如何形成高效的刚度矩阵的稀疏存储也是一 个非常重要的问题，在 GeoFEA 中，我们采用 Chen[7] 所描述的快速排序形成压缩的稀疏行（或者是压缩 的稀疏列）的矩阵形成、存储方式。 3.3 预处理子
标记。 这里，所谓最后出现的结点号指的是该结点在随后 的所有单元组装中不再出现，这样便可对该结点对 应的未知量按（2）式进行消去。值得提及的是，对 于波前直接求解法，不同的有限单元编号顺序会对 波前求解法的求解性能具有显著的影响。因此， GeoFEA 自动采用单元编号优化算法对单元编号进 行优化，使得最大波前宽尽可能小，从而节省计算 时间和存储量。从波前宽角度考虑，狭长三维问题 很适合采用波前直接求解法。
然而，问题求解方法的稳定性和可靠性仍是一个需
要慎重的方面，所以 Phoon et al[3]主张使用拟最小
残量算法来求解比奥固结方程离散所产生的大规模
对称非正定线性方程系统。
3.2 刚度矩阵的存储 通常有限元计算中所指的单元接单元迭代求解
和稀疏迭代求解实质上是指线性系统刚度矩阵的存
储。这两种矩阵存储格式在目前有限元软件中最为
中国土木工程学会第十届土力学及岩土工程学术会议论文集
GeoFEA 的直接和迭代线性求解器
1
2
2
2
陈曦 ，方思瀚 ，王效宁 ，罗章权
（1. 北京伟思德克科技有限责任公司，上地信息路甲 28 号科实大厦 B 座 10 层 D， 北京市海淀区，100085； 2. 吉奥软件公司，2 Kaki Bukit Place, #07-00 Tritech Building, 新加坡, 416180）
1 引言
岩土工程问题的大规模有限元分析中，空间离 散和时间插分所产生的线性方程系统经常含有上 万、数十万、甚至成百万个未知数。求解这样大规 模的线性方程系统，常规的高斯消去基础上的直接 解法需要消耗大量计算时间和内存，以至于基于个 人计算机的模拟是十分困难的。自从 Hestenes 和 Stiefel[1]发现共轭梯度法（CG–Conjugate Gradient） 以来，求解线性方程的 Krylov 子空间迭代法有了长 Байду номын сангаас的发展，尤其是近二十年的发展尤为显著。基于 个人计算机的 Krylov 子空间迭代计算法已经成功 地应用于大规模二维，尤其是三维岩土工程有限元 分析中 （见 Smith and Griffiths[2]，Phoon et al[3]， Chen[4]）。大型三维有限元分析软件 GeoFEA[5]是由 新加坡吉奥软件公司开发的新一代岩土工程软件， 本文对 GeoFEA 的波前直接求解器、对称和非对称 Krylov 子空间迭代求解器的内在算法及适用范围做 了较为详细的阐述。
中国土木工程学会第十届土力学及岩土工程学术会议论文集
Ae ; ve v ; we = Aeve ; w we end for 这里，nel 为模型中的单元总数，Ae 为单元刚度矩 阵， Le 是一个 nel×N 的布尔矩阵，nel 而为模型中单 元刚度矩阵的维度。ve 是通过布尔矩阵 Le 从整体 矢量 v 出来对应的单元矢量，而 we 为组装到整体矢 量 w 的对应的单元矢量。稀疏矩阵存储需要组装并 只存储整体刚度矩阵的非零元素，所以迭代法中的 整体刚度矩阵与矢量的乘积 w = Av 是基于已经组装 的稀疏矩阵 A，其伪代码如下：
并不需要完全组装。那么，问题是：该如何判断方 程（2）右端的三项乘积中所有三个矩阵元素都被完 全组装。在波前求解进行之前，程序必须确定两个
中国土木工程学会第十届土力学及岩土工程学术会议论文集
信息： （1）整个求解域中单元的编号； （2）按当前单元编号对最后出现的结点号进行
摘 要：有限元法离散所产生的大规模线性方程组的求解是有限元模拟过程中最耗计算时间和内存的部分，因此，需要使用 高效的线性求解器来快速精确地进行有限元模拟。GeoFEA 的线性求解器分为直接和 Krylov 子空间迭代求解器两大类。直接 求解器的特点是可以用固定计算量精确地得到线性方程系统的解，但对于涉及十万甚至上百万个自由度的大规模岩土工程问 题，现有的直接求解器所需要的内存和计算时间是难以接受的；迭代求解器的特点是内存需求量小，当结合高效的预处理子 时，对三维问题所产生大规模线性方程系统的计算时间要显著少于直接求解器。 关 键 词：GeoFEA；线性方程系统；直接求解器；迭代求解器
在 Krylov 子空间迭代法中，矩阵与矢量的乘积 构成了所需计算量的主要部分。因此，对于有限元 法或有限差分法得到的稀疏矩阵，迭代法中的矩阵 与矢量乘积要优越于直接求解法的 LU 分解。一方 面，LU 分解通常使被分解稀疏矩阵的零值元素变 为非零而产生稠密的 L，U 分解因子，从而导致难 以承受的计算量和存储量。然而，现有的岩土工程 有限元软件中，迭代求解法并没有完全取代直接求 解法主要基于下面两个原因：
{ } Kk ( A, r0 ) = span r0 , Ar0 ,K Ak−1r0 , k = 1, 2,K (3)
这样，近似解可以通过这个子空间和初始假设解 x0 根据下式产生：
xk ∈ x0 + Kkk ( A, r0 ), k = 1, 2,K
(4)
因此，Krylov 子空间迭代法的迭代框架可以概括如 下：
一般来说，一个高效的预处理子具有下面性质： (1) 对于某一 Krylov 子空间迭代法，预处理后 的线性方程系统具有更快的收敛性质，达到收敛时 所需要的迭代步数要明显低于原始线性方程系统达 到收敛时的迭代步数； (2) 每一迭代步中，预处理所带来的计算量有
限，可以用被减少的迭代步数来弥补使得最终收敛 时所需要的求解时间明显少于原始线形方程系统的 求解时间；
⎧⎨⎩ abii**j
= =
aij − aik ak−k1akj bi − aik ak−k1bk
(2)
在传统的高斯消去求解法中，方程（2）右端的 aij 和 bi 都是完全组装的，并且高斯消去或者是 LU
分解都是基于这些已经完全组装的矩阵元素上的。 而 Irons 观察到对于高斯消去过程，只要方程（2） 右端的三项乘积中所有三个矩阵元素都被完全组
3 GeoFEA 的对称 Krylov 子空间求解 器
求解线性方程系统（1）有多种方法，主要分为 定常迭代和非定常迭代两大类。而属于非定常迭代
的 Krylov 子空间迭代法最为流行。Krylov 子空间迭 代法的两个关键成分：一个是如何构筑子空间，另
一个是如何构筑近似解。因此，根据构筑方法不同， 我们可以得到不同的 Krylov 子空间迭代法（参见 Chen[7]; Chen et al[8]）。Krylov 子空间主要通过矩阵 与矢量的乘积来构筑，可以具体表示为：
对于对称正定线性方程系统，共轭剃度法是目 前最为有效的方法。但对于对称非正定的线性系统， 比较常用的是最小残量 (MINRES)方法和对称 LQ (SYMMLQ)方法，但最小残量法和对称 LQ 这两种 方法的限制是所使用的预处理子只能是对称正定 的。 Freund 发展了简化的或对称的拟最小残量方 法 (SQMR) 来求解对称非正定线性系统并可以结 合任意对称预处理子。近年来，Smith 和 Griffiths[2] 也成功地将预处理共轭剃度法用于比奥固结产生的 对称非正定线性系统，并且很少出现问题。 下式给 出来比奥固结方程离散得到的 2×2 的耦合对称非 正定线性方程
(3) 预处理子所需要的内存存储量有限。 方程（1）的预处理形式可以表达为：
PL−1 APR−1PR x = PL−1b
(6)
这里，一个广义的预处理子可以表示为 P = PL PR。 当 PL = I, PR ≠ I 表示预处理方式为右端处理；当 PL≠ I, PR =I 表示预处理方式为左端处理；当 PL ≠ I, PR ≠ I 表示为左右端处理方式。GeoFEA 提供了几种对 称预处理子，包括修正雅克比预处理子（MJ）、广 义雅克比预处理子（GJ）和修正对称连续超松弛预 处理子（MSSOR）。观察到标准雅克比预处理子在 求解比奥固结线性方程时性能变差，Phoon et al[3] 提出了广义雅克比预处理子来改善了对角预处理子
Ax
=
b，这里
A
=
⎡ ⎢ ⎢⎣
K BT
B⎤ ⎥
−C ⎥⎦ N×N
(5)
这里，N = m + n; A ∈ ℜ(m+n)×(m+n) 是 2×2 耦合矩阵; b ∈ ℜm+n；K ∈ ℜm×m，是对应于土体骨架刚度的对称 正定矩阵；C ∈ ℜn×n 是对应于流体刚度的对称半正 定矩阵； B ∈ ℜm×n 是土体骨架刚度与流体刚度的 耦合矩阵；m，n 分别对应结点位移自由度总数和超 静孔隙水压力自由度总数。
常见，其基本特点可以概括如下：
(1) 单元接单元的矩阵存储指的是有限元离散的单 元刚度矩阵不需要组装。 这样，Krylov 子空间迭 代法中的整体刚度矩阵与矢量的乘积 w = Av 可以
∑ 通过单元接单元地进行，即
Av
=
⎛ ⎜⎝
e
LTe
Ae
Le
⎞ ⎟⎠
v
，
其伪代码如下：
for i = 1 to nel do
End For Krylov 子空间迭代法之所以流行有几个很重要的 原因：一方面，预处理技术的使用可以有效地降低 病态线性方程系统的条件数，或者说，预处理技术 可以改善系数矩阵的特征值分布而使其更为聚集， 因为对称子空间迭代法的收敛性与线性方程的系数 矩阵的特征值分布紧密相关。 通过应用预处理技
术，Krylov 子空间迭代法收敛时所需要的迭代步数 要远小于线性方程系统的总自由度数。另一方面， 用有限元、有限差分和有限体积等方法离散所产生 的线性系统的刚度矩阵通常都是稀疏的，而稀疏矩 阵又尤其适合矩阵与矢量乘积的计算。而且矩阵矢 量的乘积可以在不组装整体系数或刚度矩阵的条件 下得到。 3.1 对称 Krylov 子空间迭代求解器
w = 0; for j = 1 to n do
if v(j) ≠ 0 then for k = jebe(j) to jebe(j+1) – 1 do
r = iebe(k)
w(r) = w(r) + v(j)*ebea(k) end for end if end for 这里，iebe(:)，jebe(:)和 ebea(:)是整体矩阵 A 的压缩 稀疏列（CSC）的存储格式。对于对称矩阵的上三 角或下三角的存储格式，上面的伪代码需要进行相 应调整。尽管从某种程度来讲，两种存储格式都是 稀疏存储。但对于单元接单元的存储格式，仍有大 量的重复的整体刚度矩阵位置 (i, j)的元素， 而对 于已经完全组装的完全稀疏矩阵格式，(i, j)位置的 元素是唯一的。这样对于稀疏矩阵存储的格式，矩 阵与矢量乘积的所需要的浮点运算次数要明显少于 单元接单元的存储格式。
一个线性方程系统可用下面的通用形式表示： Ax = b, 这里 A= [aij] ∈ ℜN×N , 并且 x, b ∈ ℜN (1) 求解线性方程系统（1）的各种直接求解法都是传统 高斯消去算法的变体。其中，较著名的有 Irons[6]提 出的用于求解对称有限元矩阵方程的波前直接求解 法。波前直接求解法采用单元刚度矩阵组装与高斯 消去相互交织的求解方式，基本求解原理可以概括 为
（1）对所有的耦合、非耦合问题，迭代求解法
并不能确保都收敛； （2）直接求解法能够用固定的浮点运算数得到
精确的、可靠的解。 基 于 上 面 考 虑 ， GeoFEA 既 提 供 了 波 前
（Frontal）直接求解器，也提供了多种可供选择的 对称和非对称 Krylov 子空间迭代求解器。
2 GeoFEA 的波前直接求解器
选用迭代初始解 x0，通常可设 x0 = 0； For k = 1 to max_it, Do
（a）如果必要，应用预处理子（Preconditioner）； （b）计算当前步的近似解，xk； （c）如果必要，计算当前步残量，rk = b–Axk ； （d）使用适当的收敛准则来检查收敛情况，如果 收敛，则跳出循环；
(2) 单元接单元的存储格式的一个显著的优势就是 内存需求量少。因为尽管稀疏求解法只需存储整体 刚度的非零元素，但需要组装工作空间，所以最后 所需的峰值存储量还是高于单元接单元所需要的总 存储量。如果考虑对称性质，所需存储量可缩减约 一半。如何形成高效的刚度矩阵的稀疏存储也是一 个非常重要的问题，在 GeoFEA 中，我们采用 Chen[7] 所描述的快速排序形成压缩的稀疏行（或者是压缩 的稀疏列）的矩阵形成、存储方式。 3.3 预处理子
标记。 这里，所谓最后出现的结点号指的是该结点在随后 的所有单元组装中不再出现，这样便可对该结点对 应的未知量按（2）式进行消去。值得提及的是，对 于波前直接求解法，不同的有限单元编号顺序会对 波前求解法的求解性能具有显著的影响。因此， GeoFEA 自动采用单元编号优化算法对单元编号进 行优化，使得最大波前宽尽可能小，从而节省计算 时间和存储量。从波前宽角度考虑，狭长三维问题 很适合采用波前直接求解法。
然而，问题求解方法的稳定性和可靠性仍是一个需
要慎重的方面，所以 Phoon et al[3]主张使用拟最小
残量算法来求解比奥固结方程离散所产生的大规模
对称非正定线性方程系统。
3.2 刚度矩阵的存储 通常有限元计算中所指的单元接单元迭代求解
和稀疏迭代求解实质上是指线性系统刚度矩阵的存
储。这两种矩阵存储格式在目前有限元软件中最为
中国土木工程学会第十届土力学及岩土工程学术会议论文集
GeoFEA 的直接和迭代线性求解器
1
2
2
2
陈曦 ，方思瀚 ，王效宁 ，罗章权
（1. 北京伟思德克科技有限责任公司，上地信息路甲 28 号科实大厦 B 座 10 层 D， 北京市海淀区，100085； 2. 吉奥软件公司，2 Kaki Bukit Place, #07-00 Tritech Building, 新加坡, 416180）
1 引言
岩土工程问题的大规模有限元分析中，空间离 散和时间插分所产生的线性方程系统经常含有上 万、数十万、甚至成百万个未知数。求解这样大规 模的线性方程系统，常规的高斯消去基础上的直接 解法需要消耗大量计算时间和内存，以至于基于个 人计算机的模拟是十分困难的。自从 Hestenes 和 Stiefel[1]发现共轭梯度法（CG–Conjugate Gradient） 以来，求解线性方程的 Krylov 子空间迭代法有了长 Байду номын сангаас的发展，尤其是近二十年的发展尤为显著。基于 个人计算机的 Krylov 子空间迭代计算法已经成功 地应用于大规模二维，尤其是三维岩土工程有限元 分析中 （见 Smith and Griffiths[2]，Phoon et al[3]， Chen[4]）。大型三维有限元分析软件 GeoFEA[5]是由 新加坡吉奥软件公司开发的新一代岩土工程软件， 本文对 GeoFEA 的波前直接求解器、对称和非对称 Krylov 子空间迭代求解器的内在算法及适用范围做 了较为详细的阐述。
中国土木工程学会第十届土力学及岩土工程学术会议论文集
Ae ; ve v ; we = Aeve ; w we end for 这里，nel 为模型中的单元总数，Ae 为单元刚度矩 阵， Le 是一个 nel×N 的布尔矩阵，nel 而为模型中单 元刚度矩阵的维度。ve 是通过布尔矩阵 Le 从整体 矢量 v 出来对应的单元矢量，而 we 为组装到整体矢 量 w 的对应的单元矢量。稀疏矩阵存储需要组装并 只存储整体刚度矩阵的非零元素，所以迭代法中的 整体刚度矩阵与矢量的乘积 w = Av 是基于已经组装 的稀疏矩阵 A，其伪代码如下：