Mercurial > hg > octave-nkf

comparison scripts/sparse/bicgstab.m @ 9279:1673a0dc019f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix & improve preconditioning strategy in cgs & bicgstab
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 28 May 2009 07:37:13 +0200
parents 2b35cb600d50
children be55736a0783
comparison
equal deleted inserted replaced
9278:2b35cb600d50 9279:1673a0dc019f
73 maxit = min (rows (b), 20); 73 maxit = min (rows (b), 20);
74 endif 74 endif
75 75
76 ## Left preconditioner. 76 ## Left preconditioner.
77 if (nargin == 5) 77 if (nargin == 5)
78 precon = M1; 78 if (isnumeric (M1))
79 precon = @(x) M1 \ x;
80 endif
79 elseif (nargin > 5) 81 elseif (nargin > 5)
80 if (isparse(M1) && issparse(M2)) 82 if (issparse (M1) && issparse (M2))
81 precon = @(x) M1 * (M2 * x); 83 precon = @(x) M2 \ (M1 \ x);
82 else 84 else
83 precon = M1 * M2; 85 M = M1*M2;
84 endif 86 precon = @(x) M \ x;
87 endif
88 else
89 precon = @(x) x;
85 endif 90 endif
86
87 if (nargin > 4 && isnumeric(precon) )
88 precon = inv(precon);
89 endif
90
91 91
92 ## specifies initial estimate x0 92 ## specifies initial estimate x0
93 if (nargin < 7) 93 if (nargin < 7)
94 x = zeros (rows (b), 1); 94 x = zeros (rows (b), 1);
95 else 95 else
115 else 115 else
116 beta = (rho_1 / rho_2) * (alpha / omega); 116 beta = (rho_1 / rho_2) * (alpha / omega);
117 p = res + beta * (p - omega * v); 117 p = res + beta * (p - omega * v);
118 endif 118 endif
119 119
120 if (nargin > 4 && isnumeric (precon)) 120 phat = precon (p);
121 phat = precon * p;
122 elseif (nargin > 4)
123 ## Our preconditioner is a function.
124 phat = feval (precon, p);
125 else
126 phat = p;
127 endif
128 121
129 v = A * phat; 122 v = A * phat;
130 alpha = rho_1 / (rr' * v); 123 alpha = rho_1 / (rr' * v);
131 s = res - alpha * v; 124 s = res - alpha * v;
132 125
133 if (nargin > 4 && isnumeric (precon)) 126 shat = precon (s);
134 shat = precon * s;
135 elseif (nargin > 4)
136 ## Our preconditioner is a function.
137 shat = feval (precon, s);
138 else
139 shat = s;
140 endif
141 127
142 t = A * shat; 128 t = A * shat;
143 omega = (t' * s) / (t' * t); 129 omega = (t' * s) / (t' * t);
144 x = x + alpha * phat + omega * shat; 130 x = x + alpha * phat + omega * shat;
145 res = s - omega * t; 131 res = s - omega * t;