【讨论】预调件共轭梯度法（PCG） - 数学 - 计算数学 - 第4页小木虫论坛-学术科研互动平台

31楼2011-11-24 13:43:59

已阅申请数学EPI 回复此楼编辑查看我的主页

saladin983

铁杆木虫 (正式写手)

数学EPI: 4
应助: 6 (幼儿园)
贵宾: 0.2
金币: 9197.9
红花: 3
帖子: 687
在线: 115.7小时
虫号: 448168
注册: 2007-11-01
专业: 计算数学与科学工程计算

★
小木虫(金币+0.5):给个红包，谢谢回帖

引用回帖:

31楼: Originally posted by mastergxm at 2011-11-24 07:43:59:

上面是Iterative Methods for Sparse Linear Systems 书中关于BICGSTAB的算法步骤，我后来按照上面的步骤编了程序，但是总是有问题，后来仔细看了步骤，发现每一 ...

:= 在伪程序中表示赋值，就是通常意义下的等号。至于这个记号怎么来的，我就不清楚了。程序的问题，应该跟这个无关。如果用matlab的话，自带BiCGStab的函数。其它语言的网上应该也都能找到吧，你可以找个参照一下。这个算法实现起来简单，debug应该也不会太难吧。

32楼2011-11-25 02:59:18

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

匿名

用户注销 (知名作家)

应助: 1 (幼儿园)
金币: 3228.2
散金: 24
红花: 1
帖子: 6711
在线: 384.6小时
虫号: 0
注册: 2009-11-17
专业: 大气环境与全球气候变化

★
小木虫(金币+0.5):给个红包，谢谢回帖

本帖仅楼主可见

33楼2012-01-02 21:41:35

已阅申请数学EPI 回复此楼编辑查看我的主页

saladin983

铁杆木虫 (正式写手)

数学EPI: 4
应助: 6 (幼儿园)
贵宾: 0.2
金币: 9197.9
红花: 3
帖子: 687
在线: 115.7小时
虫号: 448168
注册: 2007-11-01
专业: 计算数学与科学工程计算

★
小木虫(金币+0.5):给个红包，谢谢回帖

引用回帖:

33楼: Originally posted by mastergxm at 2012-01-02 15:41:35:
我发现这个BiCGSTAB在有些时候似乎不管用。
我就碰到这样一种情况，以上面图中求解过程中的各变量为例。
初始值我是这样设置的：X0＝0；r*o因为可以任意设置，所以我在程序中将r*o＝X0；并且的我的方程组 ...

对于Krylov子空间类的迭代算法的效率，最重要的自然是系数矩阵的特征值分布，但是初始向量的确也会有很大的影响。比如说，在系数矩阵的（绝对值）最小的特征值接近于0的情况下，若初始迭代向量是相应的特征值的话，则最开始的迭代过程非常耗时。事实上在很多情况下，CG的迭代过程表现出开始时残量下降较慢，而最后一段则下降较大，这是一个在实际运算中很麻烦的状况。通常需要通过条件预处理来克服。这一点，我相信GMRES不会比BiCGStab更好。

回到你的问题，我比较困惑的是，如果右端向量是0的话，你的解不应该就是0么？难道是不满秩的方程组或者是过定的？但这些情况下都是不能用这些迭代法直接求解的。

34楼2012-01-03 22:55:26

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

匿名

用户注销 (知名作家)

应助: 1 (幼儿园)
金币: 3228.2
散金: 24
红花: 1
帖子: 6711
在线: 384.6小时
虫号: 0
注册: 2009-11-17
专业: 大气环境与全球气候变化

★
小木虫(金币+0.5):给个红包，谢谢回帖

本帖仅楼主可见

35楼2012-01-04 09:08:41

已阅申请数学EPI 回复此楼编辑查看我的主页

saladin983

铁杆木虫 (正式写手)

数学EPI: 4
应助: 6 (幼儿园)
贵宾: 0.2
金币: 9197.9
红花: 3
帖子: 687
在线: 115.7小时
虫号: 448168
注册: 2007-11-01
专业: 计算数学与科学工程计算

★
小木虫(金币+0.5):给个红包，谢谢回帖

引用回帖:

35楼: Originally posted by mastergxm at 2012-01-04 03:08:41:
解的确应为0，但是第7步求出的X_j+1 中的绝大部分x都为0，但是个别x就非常大，接近我所声明变量类型的最大值，看来还是我的程序中没有考虑到这种情况。没有考虑到在计算出如出现NaN，负无穷等情况下的处理办法。 ...

其实实际算法实现中，经常会对每一步的残量r进行判断，甚至直接用\|r\|
ILU是不完全的LU，主旨就是节省计算量。建议还是直接读算法来对比一下吧。

36楼2012-01-06 00:57:33

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

匿名

用户注销 (知名作家)

应助: 1 (幼儿园)
金币: 3228.2
散金: 24
红花: 1
帖子: 6711
在线: 384.6小时
虫号: 0
注册: 2009-11-17
专业: 大气环境与全球气候变化

★
小木虫: 金币+0.5, 给个红包，谢谢回帖

本帖仅楼主可见

37楼2012-08-03 14:23:46

已阅申请数学EPI 回复此楼编辑查看我的主页

saladin983

铁杆木虫 (正式写手)

数学EPI: 4
应助: 6 (幼儿园)
贵宾: 0.2
金币: 9197.9
红花: 3
帖子: 687
在线: 115.7小时
虫号: 448168
注册: 2007-11-01
专业: 计算数学与科学工程计算

★ ★ ★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小雨萌萌: 金币+2, 谢谢关注~ 2012-09-03 12:58:34

引用回帖:

2405247楼: Originally posted by mastergxm at 2012-08-03 08:23:46
对于PCG方法有个疑问，不知道是我没有弄情楚，还是实际上就是如此。
81/2a/905677_1343974829_332.jpg
这是PCG的计算步骤。其中用到的P就是preconditioner。书中介绍说把A作ILU分解后，用LU的乘积矩阵作为这 ...

ILU分解得到的是上三角矩阵和下三角矩阵，应用P的时候，实际上是求解相关的两个方程组，由于是三角矩阵，计算量并不大。事实上，在应用中，几乎不会出现使用P的逆矩阵的情况，而是以求解相关线性方程组，亦即Py=r的方式实现。

38楼2012-08-05 04:03:27

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

匿名

用户注销 (知名作家)

应助: 1 (幼儿园)
金币: 3228.2
散金: 24
红花: 1
帖子: 6711
在线: 384.6小时
虫号: 0
注册: 2009-11-17
专业: 大气环境与全球气候变化

★ ★ ★ ★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小雨萌萌: 金币+3, 3Q 2012-09-03 13:01:05

本帖仅楼主可见

39楼2012-08-05 10:51:58

已阅申请数学EPI 回复此楼编辑查看我的主页

saladin983

铁杆木虫 (正式写手)

数学EPI: 4
应助: 6 (幼儿园)
贵宾: 0.2
金币: 9197.9
红花: 3
帖子: 687
在线: 115.7小时
虫号: 448168
注册: 2007-11-01
专业: 计算数学与科学工程计算

★ ★ ★ ★ ★ ★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小雨萌萌: 金币+5, 3Q 2012-09-03 13:01:48

引用回帖:

2405363楼: Originally posted by mastergxm at 2012-08-05 04:51:58
f1/77/905677_1344167710_647.jpg
先求解(5)式，得到y之后，再求解(4)式，就可得到x，但是(4)、(5)式的系数矩阵不对称，CG法显然又不能求解，所以感觉这也就不是PCG了。第二种理解是直接用ILU分解后得到的L和U ...

就是你所说的先求解(5)式，得到y之后，再求解(4)式，这里因为L和U都是三角矩阵，直接求解的计算量很小，求解y和x的时候，按顺序计算各个分量即可，这里是不需要使用CG法。之前说的两个相关的方程组，就是指可以轻易地直接求解的（5）和（4）式。

事实上，P的应用，就是求解一个形如Py=r的方程组，这里也可以用其它的办法，ILU是一种比较廉价的做法，也经常见用（代数）多重网格法。还有些时候，构造出来的P可能会有特殊结构，使得Py=r非常易解，比如P是三对角矩阵。

你的理解有些偏差。CG法里应用preconditioner所需要的不一定是显式的P，算法里事实上只用到Py=r里的y，这个y总是用尽可能低廉的算法达到合适的精度，这里如何求解y会比较需要经验。CG求解的是Ax=b这个方程组，而如果利用precondioner的话，CG算法的每一步都要求解一个Py=r的线性方程组，这是两个不同的层级。在Py=r这个级别上，对于算法的高效性要求更高，而对精度的要求相对要低点，所以才会容许使用ILU分解这种近似算法。当然，你也可以直接用CG方法求解Py=r到适当的精度，但是这个时候，这里的CG法和求解Ax=b的CG法已经是一个嵌套的关系了。