求助:A为大型非对称矩阵时,AX=b的迭代算法
我需要AX=B迭代求解过程只涉及A和向量乘积的计算，不破坏A的完整性,列如共厄梯度法CG和广义极小残差法GMRES.现在我试过GMRES法中的一种，效果不是很好 ,系数矩阵A维数增大或者A中有一列的数比起A中其他数小几个数量级时根本算不对
不知道有没有做过这方面的人,给点建议或者推荐一种方法,谢谢
如有人清楚请加我QQ68836285或者发EMAIL hao.19820125@163.com
[ 本帖最后由 咕噜噜 于 2007-6-15 08:52 编辑 ]
帮忙呀
要是哪能找到GMRES的原程序及相关资料也好呀![ 本帖最后由 hao1982 于 2006-7-6 16:52 编辑 ]
re:
(1) 将方程变换 A'*Ax=A'b;A‘*A 是对称的，可以用处理对称方程的方法求解
但是这种处理可能增大方程的条件数
（2）用matlb中的GMRES。。。。。。help GMRES
THANKS
看看这个，在网上搜到的//*****************************************************************
// Iterative template routine -- GMRES
//
// GMRES solves the unsymmetric linear system Ax = b using the
// Generalized Minimum Residual method
//
// GMRES follows the algorithm described on p. 20 of the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x--approximate solution to Ax = b
// max_iter--the number of iterations performed before the
// tolerance was reached
// tol--the residual after the final iteration
//
//*****************************************************************
template < class Matrix, class Vector >
void
Update(Vector &x, int k, Matrix &h, Vector &s, Vector v[])
{
Vector y(s);
// Backsolve:
for (int i = k; i >= 0; i--) {
y(i) /= h(i,i);
for (int j = i - 1; j >= 0; j--)
y(j) -= h(j,i) * y(i);
}
for (int j = 0; j <= k; j++)
x += v * y(j);
}
template < class Real >
Real
abs(Real x)
{
return (x > 0 ? x : -x);
}
template < class Operator, class Vector, class Preconditioner,
class Matrix, class Real >
int
GMRES(const Operator &A, Vector &x, const Vector &b,
const Preconditioner &M, Matrix &H, int &m, int &max_iter,
Real &tol)
{
Real resid;
int i, j = 1, k;
Vector s(m+1), cs(m+1), sn(m+1), w;
Real normb = norm(M.solve(b));
Vector r = M.solve(b - A * x);
Real beta = norm(r);
if (normb == 0.0)
normb = 1;
if ((resid = norm(r) / normb) <= tol) {
tol = resid;
max_iter = 0;
return 0;
}
Vector *v = new Vector;
while (j <= max_iter) {
v = r * (1.0 / beta); // ??? r / beta
s = 0.0;
s(0) = beta;
for (i = 0; i < m && j <= max_iter; i++, j++) {
w = M.solve(A * v);
for (k = 0; k <= i; k++) {
H(k, i) = dot(w, v);
w -= H(k, i) * v;
}
H(i+1, i) = norm(w);
v = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
for (k = 0; k < i; k++)
ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
if ((resid = abs(s(i+1)) / normb) < tol) {
Update(x, i, H, s, v);
tol = resid;
max_iter = j;
delete [] v;
return 0;
}
}
Update(x, m - 1, H, s, v);
r = M.solve(b - A * x);
beta = norm(r);
if ((resid = beta / normb) < tol) {
tol = resid;
max_iter = j;
delete [] v;
return 0;
}
}
tol = resid;
delete [] v;
return 1;
}
#include <math.h>
template<class Real>
void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
if (dy == 0.0) {
cs = 1.0;
sn = 0.0;
} else if (abs(dy) > abs(dx)) {
Real temp = dx / dy;
sn = 1.0 / sqrt( 1.0 + temp*temp );
cs = temp * sn;
} else {
Real temp = dy / dx;
cs = 1.0 / sqrt( 1.0 + temp*temp );
sn = temp * cs;
}
}
template<class Real>
void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
Real temp=cs * dx + sn * dy;
dy = -sn * dx + cs * dy;
dx = temp;
} CG一般只用于对称正定的矩阵。GMRES理论上不依赖于矩阵的特性，但常常收敛很慢。这时需要高质量的预处理器，这方面的算法很多，可根据具体问题采用。
如果A是稀疏矩阵，A‘*A还会破坏稀疏性。
A的阶数有多大？ 原帖由 closest 于 2006-11-15 13:39 发表
CG一般只用于对称正定的矩阵。GMRES理论上不依赖于矩阵的特性，但常常收敛很慢。这时需要高质量的预处理器，这方面的算法很多，可根据具体问题采用。
如果A是稀疏矩阵，A‘*A还会破坏稀疏性。
A的阶数有多大？
能否介绍一下你说的高质量的预处理器？ 预处理算法至少有十几种，以Saad(迭代法权威)的一系列论文为代表，许多算法的复杂程度远超GMRES本身。这些算法更多的是基于一些假设和测试，理论上很难分析，也还没有通用的预处理器。对于具体的问题来说，有兴趣和时间的话可以逐个尝试这些算法，作一些改进也不是很难。
主要的问题在于对于很多矩阵，目前还没有可靠的预处理器。预处理器的表现往往两极化，对于某些矩阵收敛特别快，对于其它的矩阵则没什么作用. 原帖由 closest 于 2006-11-16 09:47 发表
预处理算法至少有十几种，以Saad(迭代法权威)的一系列论文为代表，许多算法的复杂程度远超GMRES本身。这些算法更多的是基于一些假设和测试，理论上很难分析，也还没有通用的预处理器。对于具体的问题来说，有兴趣 ...
谢谢你的介绍，有时间我看看
页:
[1]