Mercurial > hg > octave-lojdl
changeset 9279:1673a0dc019f
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 | 40fb718a2e67 |
files | scripts/ChangeLog scripts/sparse/bicgstab.m scripts/sparse/cgs.m |
diffstat | 3 files changed, 27 insertions(+), 44 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2009-05-28 Jaroslav Hajek <highegg@gmail.com> + + * sparse/bicgstab.m: Improve preconditioning; avoid explicit inverse. + * sparse/cgs.m: Improve preconditioning; avoid explicit inverse. + 2009-05-28 Radek Salac <salac.r@gmail.com> * sparse/bicgstab.m: New output when calling without arguments.
--- a/scripts/sparse/bicgstab.m +++ b/scripts/sparse/bicgstab.m @@ -75,20 +75,20 @@ ## Left preconditioner. if (nargin == 5) - precon = M1; + if (isnumeric (M1)) + precon = @(x) M1 \ x; + endif elseif (nargin > 5) - if (isparse(M1) && issparse(M2)) - precon = @(x) M1 * (M2 * x); + if (issparse (M1) && issparse (M2)) + precon = @(x) M2 \ (M1 \ x); else - precon = M1 * M2; - endif + M = M1*M2; + precon = @(x) M \ x; + endif + else + precon = @(x) x; endif - if (nargin > 4 && isnumeric(precon) ) - precon = inv(precon); - endif - - ## specifies initial estimate x0 if (nargin < 7) x = zeros (rows (b), 1); @@ -117,27 +117,13 @@ p = res + beta * (p - omega * v); endif - if (nargin > 4 && isnumeric (precon)) - phat = precon * p; - elseif (nargin > 4) - ## Our preconditioner is a function. - phat = feval (precon, p); - else - phat = p; - endif + phat = precon (p); v = A * phat; alpha = rho_1 / (rr' * v); s = res - alpha * v; - if (nargin > 4 && isnumeric (precon)) - shat = precon * s; - elseif (nargin > 4) - ## Our preconditioner is a function. - shat = feval (precon, s); - else - shat = s; - endif + shat = precon (s); t = A * shat; omega = (t' * s) / (t' * t);
--- a/scripts/sparse/cgs.m +++ b/scripts/sparse/cgs.m @@ -65,19 +65,19 @@ endif ## Left preconditioner. - precon = []; if (nargin == 5) - precon = M1; + if (isnumeric (M1)) + precon = @(x) M1 \ x; + endif elseif (nargin > 5) - if (isparse(M1) && issparse(M2)) - precon = @(x) M1 * (M2 * x); + if (issparse (M1) && issparse (M2)) + precon = @(x) M2 \ (M1 \ x); else - precon = M1 * M2; + M = M1*M2; + precon = @(x) M \ x; endif - endif - - if (nargin > 4 && isnumeric(precon) ) - precon = inv(precon); + else + precon = @(x) x; endif ## Specifies initial estimate x0. @@ -96,15 +96,7 @@ flag = 1; for iter = 1 : maxit - if (nargin > 4 && isnumeric (precon)) - z = precon * res; - elseif (nargin > 4) - ## Our preconditioner is a function. - z = feval (precon, res); - else - ## We don't use preconditioning. - z = res; - endif + z = precon (res); ## Cache. ro_old = ro;