代码：https://github.com/mrzhuzhe/riven/solver

1. 概述

动机

迭代法解线性方程组，在稀疏矩阵时比直接法快几十倍，所以自己实现一下迭代法，加深理解

• 1. 测试一下不收敛的原因
• 2. 测试一下不同方法收敛的速度
• 3. 实现时留意是否易于向量化并行
• 4. 测试一些预条件加速法，比如Multigrid或者preconditioner，或者矩阵重排
• 5. 测试求解器嵌套的情况

2.相关工作

相比较于直接法

并行化：迭代法相对直接法，并行化更直观

速度：迭代法中途有很多乘法计算在矩阵稀疏时乘以0，所以矩阵稀疏时计算起来会非常快（几十倍）
相对直接法中间过程都会变成dense矩阵，所以稀疏和不稀疏计算时间一样

误差：直接法只有浮点误差，迭代法是近似收敛法，会有收敛误差

不收敛：迭代法对矩阵本身的性质：对称正定（SPD），对称，矩阵稳定性（最大特征值/最小特征值）有要求，矩阵性质不好时会很难收敛

• 1. Painless conjugate gradient https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
• 2. A Multigrid Tutorial, 2nd Edition https://www.math.hkust.edu.hk/~mamu/courses/531/tutorial_with_corrections.pdf
• 3. Iterative Methods for Sparse Linear Systems, 2nd ed. https://www-users.cse.umn.edu/~saad/IterMethBook_2ndEd.pdf
• 4. Cosmol的各种迭代求解器用法和对比 For Usage and compare https://doc.comsol.com/5.5/doc/com.comsol.help.comsol/comsol_ref_solver.27.123.html

3.实验

代码还是都在 https://github.com/mrzhuzhe/riven/solver 这个目录下

https://github.com/mrzhuzhe/riven 现在在这个repos项目中增加了Eigen作为子项目
Clone 时需要 --recursive

3.1 Jacobian iteration

其中第一个参数："32"是矩阵大小，可以从外部传入默认是64

其中第二个参数：1 是柏松矩阵 2 是稠密对称矩阵 0 是随机稠密矩阵

后面都保持这个约定

jacobian 有个小trick就是计算时，因为都是根据对角线来，如果矩阵是基于对角线的带状矩阵，可以把一行的计算简化成几个固定乘法的计算

3.2 Guassion-seidel iteration

高斯赛德二和jacobian都需要矩阵近似对称正定，且收敛需要的步数较多，但是好处是非常容易并行化

3.3 过于简单的multigrid

你只需要写一个f2c和c2f函数

再 16 - 8 - 4 - 8 - 16 每个尺度上跑几个iteration 给下一个尺度当初始条件就可以了

3.4 Conjugate gradient

根据painless conjugate gradient文档共轭梯度就是向梯度方向前进，但是不是垂直于梯度方向（sgd），而是共轭于自己上一次的迭代方向，每次迭代的方向都是已经确定的，每次只需计算步长即可

共轭梯度为krylov子空间法，虽然共轭梯度也需要正定，但是似乎比jacobian等迭代法有更好的泛化性能

3.5 Biconjugate gradient

双共轭梯度泛化性比共轭梯度更好一些，但是需要额外再计算一个residuel和方向

还有一个计算量更大但数值更稳定的bicg-stablized 在 solver/inc/bicg_stab.h

3.6 Preconditioner conjugate gradient

预条件法改变了优化的目标residuel，且在迭代过程中要解一个子线性方程，但是会让迭代收敛大幅加快

3.7 GMRES

GMRES在计算过程中需要求特征值

4. 结论

1. 共轭梯度为什么时候会fail , fail和不收敛有关系吗?

Fail是因为共轭梯度计算alpha和beta时有除法，判断到分母为0时算法就会当作异常

• p.T*A*p 除以 r.T * r 带来的奇点

2. 为什么迭代法都是针对一列b和x计算的，多列x和b如何处理？

• 每列单独求解
• Blcok conjugate gradient https://www.sciencedirect.com/science/article/pii/0024379580902475 每次求 rows*(n)列，求cols/n组, 似乎对收敛性有影响

5. 后续工作

后面会根据实际使用时遇到的问题来持续优化

还是研究具体场景

• 1. openfoam 流体，电磁
• 2. openroad vlsi电路仿真