Mercurial > hg > octave-nkf
changeset 9278:2b35cb600d50
improve bicgstab and cgs
author | Radek Salac <salac.r@gmail.com> |
---|---|
date | Thu, 28 May 2009 07:03:11 +0200 |
parents | be84e9654feb |
children | 1673a0dc019f |
files | scripts/ChangeLog scripts/sparse/bicgstab.m scripts/sparse/cgs.m |
diffstat | 3 files changed, 62 insertions(+), 43 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,12 @@ +2009-05-28 Radek Salac <salac.r@gmail.com> + + * sparse/bicgstab.m: New output when calling without arguments. + Time optimization - remove certain checks linked to preconditioner which + are unacceptably slow. + * sparse/cgs.m: New output when calling without arguments. + Time optimization - remove certain checks linked to preconditioner which + are unacceptably slow. Rename internal variable to match bicgstab. + 2009-05-27 Rik Wehbring <rdrider0-list@yahoo.com> * plot/axis.m: Update documentation to reflect addition of "tight" option.
--- a/scripts/sparse/bicgstab.m +++ b/scripts/sparse/bicgstab.m @@ -1,4 +1,4 @@ -## Copyright (C) 2008, 2009 Radek Salac +## Copyright (C) 2008 Radek Salac ## ## This file is part of Octave. ## @@ -84,20 +84,10 @@ endif endif - if (nargin > 4 && isnumeric (precon)) - ## Precon can by also function. - if (det (precon) != 0) - ## We can compute inverse preconditioner and use quicker algorithm. - precon = inv (precon); - else - error ("bicgstab: preconditioner is ill conditioned"); - endif + if (nargin > 4 && isnumeric(precon) ) + precon = inv(precon); + endif - if (isinf (cond (precon))); - ## We must make test if preconditioner isn't ill conditioned. - error ("bicgstab: preconditioner is ill conditioned"); - endif - endif ## specifies initial estimate x0 if (nargin < 7) @@ -112,7 +102,7 @@ rr = res; ## Vector of the residual norms for each iteration. - resvec = [norm(res)]; + resvec = [norm(res)/norm_b]; ## Default behaviour we don't reach tolerance tol within maxit iterations. flag = 1; @@ -169,6 +159,23 @@ endif endfor + if (nargout < 2) + if (flag == 0) + printf (["bicgstab converged at iteration %i ", + "to a solution with relative residual %e\n"],iter,relres); + elseif (flag == 3) + printf (["bicgstab stopped at iteration %i ", + "without converging to the desired tolerance %e\n", + "because the method stagnated.\n", + "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); + else + printf (["bicgstab stopped at iteration %i ", + "without converging to the desired toleranc %e\n", + "because the maximum number of iterations was reached.\n", + "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); + endif + endif + endfunction %!demo
--- a/scripts/sparse/cgs.m +++ b/scripts/sparse/cgs.m @@ -1,4 +1,4 @@ -## Copyright (C) 2008, 2009 Radek Salac +## Copyright (C) 2008 Radek Salac ## ## This file is part of Octave. ## @@ -76,20 +76,8 @@ endif endif - ## Precon can also be a function. - if (nargin > 4 && isnumeric (precon)) - - ## We can compute inverse preconditioner and use quicker algorithm. - if (det (precon) != 0) - precon = inv (precon); - else - error ("cgs: preconditioner is ill conditioned"); - endif - - ## We must make test if preconditioner isn't ill conditioned. - if (isinf (cond (precon))); - error ("cgs: preconditioner is ill conditioned"); - endif + if (nargin > 4 && isnumeric(precon) ) + precon = inv(precon); endif ## Specifies initial estimate x0. @@ -99,29 +87,28 @@ x = x0; endif - relres = b - A * x; + res = b - A * x; + norm_b = norm (b); ## Vector of the residual norms for each iteration. - resvec = [norm(relres)]; + resvec = [ norm(res)/norm_b ]; ro = 0; - norm_b = norm (b); ## Default behavior we don't reach tolerance tol within maxit iterations. flag = 1; for iter = 1 : maxit if (nargin > 4 && isnumeric (precon)) - ## We have computed inverse matrix so we can use quick algorithm. - z = precon * relres; + z = precon * res; elseif (nargin > 4) ## Our preconditioner is a function. - z = feval (precon, relres); + z = feval (precon, res); else ## We don't use preconditioning. - z = relres; + z = res; endif ## Cache. ro_old = ro; - ro = relres' * z; + ro = res' * z; if (iter == 1) p = z; else @@ -133,11 +120,11 @@ alpha = ro / (p' * q); x = x + alpha * p; - relres = relres - alpha * q; - resvec = [resvec; norm(relres)]; + res = res - alpha * q; + relres = norm(res) / norm_b; + resvec = [resvec;relres]; - relres_distance = resvec (end) / norm_b; - if (relres_distance <= tol) + if (relres <= tol) ## We reach tolerance tol within maxit iterations. flag = 0; break; @@ -148,7 +135,23 @@ endif endfor; - relres = relres_distance; + if (nargout < 1) + if ( flag == 0 ) + printf (["cgs converged at iteration %i ", + "to a solution with relative residual %e\n"],iter,relres); + elseif (flag == 3) + printf (["cgs stopped at iteration %i ", + "without converging to the desired tolerance %e\n", + "because the method stagnated.\n", + "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); + else + printf (["cgs stopped at iteration %i ", + "without converging to the desired tolerance %e\n", + "because the maximum number of iterations was reached.\n", + "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); + endif + endif + endfunction