求助:A为大型非对称矩阵时,AX=b的迭代算法

hao1982

我需要AX=B迭代求解过程只涉及A和向量乘积的计算，不破坏A的完整性,列如共厄梯度法CG和广义极小残差法GMRES.

现在我试过GMRES法中的一种，效果不是很好 ,系数矩阵A维数增大或者A中有一列的数比起A中其他数小几个数量级时根本算不对
不知道有没有做过这方面的人,给点建议或者推荐一种方法,谢谢

如有人清楚请加我QQ68836285或者发EMAIL   hao.19820125@163.com

[ 本帖最后由 咕噜噜 于 2007-6-15 08:52 编辑 ]

hao1982

要是哪能找到GMRES的原程序及相关资料也好呀!

[ 本帖最后由 hao1982 于 2006-7-6 16:52 编辑 ]

xj2070

(1) 将方程变换 A'*Ax=A'b;

    A‘*A 是对称的，可以用处理对称方程的方法求解
   但是这种处理可能增大方程的条件数

（2）用matlb中的GMRES。。。。。。help GMRES

hao1982

hailiangwww

看看这个，在网上搜到的

1. //*****************************************************************
   
2. // Iterative template routine -- GMRES
   
3. //
   
4. // GMRES solves the unsymmetric linear system Ax = b using the
   
5. // Generalized Minimum Residual method
   
6. //
   
7. // GMRES follows the algorithm described on p. 20 of the
   
8. // SIAM Templates book.
   
9. //
   
10. // The return value indicates convergence within max_iter (input)
    
11. // iterations (0), or no convergence within max_iter iterations (1).
    
12. //
    
13. // Upon successful return, output arguments have the following values:
    
14. //  
    
15. //        x  --  approximate solution to Ax = b
    
16. // max_iter  --  the number of iterations performed before the
    
17. //               tolerance was reached
    
18. //      tol  --  the residual after the final iteration
    
19. //  
    
20. //*****************************************************************
    
21. 
    
22. 
    
23. template < class Matrix, class Vector >
    
24. void
    
25. Update(Vector &x, int k, Matrix &h, Vector &s, Vector v[])
    
26. {
    
27.   Vector y(s);
    
28. 
    
29.   // Backsolve:  
    
30.   for (int i = k; i >= 0; i--) {
    
31.     y(i) /= h(i,i);
    
32.     for (int j = i - 1; j >= 0; j--)
    
33.       y(j) -= h(j,i) * y(i);
    
34.   }
    
35. 
    
36.   for (int j = 0; j <= k; j++)
    
37.     x += v[j] * y(j);
    
38. }
    
39. 
    
40. 
    
41. template < class Real >
    
42. Real
    
43. abs(Real x)
    
44. {
    
45.   return (x > 0 ? x : -x);
    
46. }
    
47. 
    
48. 
    
49. template < class Operator, class Vector, class Preconditioner,
    
50.            class Matrix, class Real >
    
51. int
    
52. GMRES(const Operator &A, Vector &x, const Vector &b,
    
53.       const Preconditioner &M, Matrix &H, int &m, int &max_iter,
    
54.       Real &tol)
    
55. {
    
56.   Real resid;
    
57.   int i, j = 1, k;
    
58.   Vector s(m+1), cs(m+1), sn(m+1), w;
    
59.   
    
60.   Real normb = norm(M.solve(b));
    
61.   Vector r = M.solve(b - A * x);
    
62.   Real beta = norm(r);
    
63.   
    
64.   if (normb == 0.0)
    
65.     normb = 1;
    
66.   
    
67.   if ((resid = norm(r) / normb) <= tol) {
    
68.     tol = resid;
    
69.     max_iter = 0;
    
70.     return 0;
    
71.   }
    
72. 
    
73.   Vector *v = new Vector[m+1];
    
74. 
    
75.   while (j <= max_iter) {
    
76.     v[0] = r * (1.0 / beta);    // ??? r / beta
    
77.     s = 0.0;
    
78.     s(0) = beta;
    
79.    
    
80.     for (i = 0; i < m && j <= max_iter; i++, j++) {
    
81.       w = M.solve(A * v[i]);
    
82.       for (k = 0; k <= i; k++) {
    
83.         H(k, i) = dot(w, v[k]);
    
84.         w -= H(k, i) * v[k];
    
85.       }
    
86.       H(i+1, i) = norm(w);
    
87.       v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
    
88. 
    
89.       for (k = 0; k < i; k++)
    
90.         ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
    
91.       
    
92.       GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
    
93.       ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
    
94.       ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
    
95.       
    
96.       if ((resid = abs(s(i+1)) / normb) < tol) {
    
97.         Update(x, i, H, s, v);
    
98.         tol = resid;
    
99.         max_iter = j;
    
100.         delete [] v;
     
101.         return 0;
     
102.       }
     
103.     }
     
104.     Update(x, m - 1, H, s, v);
     
105.     r = M.solve(b - A * x);
     
106.     beta = norm(r);
     
107.     if ((resid = beta / normb) < tol) {
     
108.       tol = resid;
     
109.       max_iter = j;
     
110.       delete [] v;
     
111.       return 0;
     
112.     }
     
113.   }
     
114.   
     
115.   tol = resid;
     
116.   delete [] v;
     
117.   return 1;
     
118. }
     
119. 
     
120. 
     
121. #include <math.h>
     
122. 
     
123. 
     
124. template<class Real>
     
125. void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
     
126. {
     
127.   if (dy == 0.0) {
     
128.     cs = 1.0;
     
129.     sn = 0.0;
     
130.   } else if (abs(dy) > abs(dx)) {
     
131.     Real temp = dx / dy;
     
132.     sn = 1.0 / sqrt( 1.0 + temp*temp );
     
133.     cs = temp * sn;
     
134.   } else {
     
135.     Real temp = dy / dx;
     
136.     cs = 1.0 / sqrt( 1.0 + temp*temp );
     
137.     sn = temp * cs;
     
138.   }
     
139. }
     
140. 
     
141. 
     
142. template<class Real>
     
143. void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
     
144. {
     
145.   Real temp  =  cs * dx + sn * dy;
     
146.   dy = -sn * dx + cs * dy;
     
147.   dx = temp;
     
148. }

复制代码

closest

CG一般只用于对称正定的矩阵。GMRES理论上不依赖于矩阵的特性，但常常收敛很慢。这时需要高质量的预处理器，这方面的算法很多，可根据具体问题采用。
如果A是稀疏矩阵，A‘*A还会破坏稀疏性。
A的阶数有多大？

gghhjj

原帖由 closest 于 2006-11-15 13:39 发表
CG一般只用于对称正定的矩阵。GMRES理论上不依赖于矩阵的特性，但常常收敛很慢。这时需要高质量的预处理器，这方面的算法很多，可根据具体问题采用。
如果A是稀疏矩阵，A‘*A还会破坏稀疏性。
A的阶数有多大？

能否介绍一下你说的高质量的预处理器？

closest

预处理算法至少有十几种，以Saad(迭代法权威)的一系列论文为代表，许多算法的复杂程度远超GMRES本身。这些算法更多的是基于一些假设和测试，理论上很难分析，也还没有通用的预处理器。对于具体的问题来说，有兴趣和时间的话可以逐个尝试这些算法，作一些改进也不是很难。
主要的问题在于对于很多矩阵，目前还没有可靠的预处理器。预处理器的表现往往两极化，对于某些矩阵收敛特别快，对于其它的矩阵则没什么作用.

gghhjj

原帖由 closest 于 2006-11-16 09:47 发表
预处理算法至少有十几种，以Saad(迭代法权威)的一系列论文为代表，许多算法的复杂程度远超GMRES本身。这些算法更多的是基于一些假设和测试，理论上很难分析，也还没有通用的预处理器。对于具体的问题来说，有兴趣 ...

谢谢你的介绍，有时间我看看