Mercurial > hg > octave-terminal
changeset 14773:5c269b73f467
gmres.m: Overhaul function. Fix bug #36568.
* gmres.m: Return correct exit flag (bug #36568). Return relative residual,
not norm, for Matlab compatibility. Use Octave coding conventions including
same variable names in documentation as in function. Add %!demo and %!error
tests.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Fri, 15 Jun 2012 10:15:07 -0700 |
parents | 39a2e91a246e |
children | 0d6dae0f6bc2 |
files | scripts/sparse/gmres.m |
diffstat | 1 files changed, 32 insertions(+), 15 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/sparse/gmres.m +++ b/scripts/sparse/gmres.m @@ -32,7 +32,7 @@ ## @code{min (10, numel (b) / restart)} is used. ## ## @item @var{x0} is the initial guess, -## if not given or set to [] the default value @code{zeros(size (b))} is used. +## if not given or set to [] the default value @code{zeros (size (b))} is used. ## ## @item @var{m} is the restart parameter, ## if not given or set to [] the default value @code{numel (b)} is used. @@ -57,7 +57,7 @@ ## ## @item 2 : unused, but skipped for compatibility ## -## @item 3 : algorithm reached stagnation +## @item 3 : algorithm reached stagnation (no change between iterations) ## @end table ## ## @item @var{relres} is the final value of the relative residual. @@ -72,7 +72,7 @@ ## @seealso{bicg, bicgstab, cgs, pcg} ## @end deftypefn -function [x, flag, presn, it, resids] = gmres (A, b, restart, rtol, maxit, M1, M2, x0) +function [x, flag, relres, it, resvec] = gmres (A, b, restart, rtol, maxit, M1, M2, x0) if (nargin < 2 || nargin > 8) print_usage (); @@ -142,10 +142,10 @@ ## begin loop iter = 1; restart_it = restart + 1; - resids = zeros (maxit, 1); - resids(1) = presn; + resvec = zeros (maxit, 1); + resvec(1) = presn; prec_b_norm = norm (Pm1x (b), 2); - flag = 1; + flag = 1; # Default flag is maximum # of iterations exceeded while (iter <= maxit * restart && presn > rtol * prec_b_norm) @@ -174,26 +174,36 @@ x = x_old + V(:, 1:restart_it) * Y(1:restart_it); - resids(iter) = presn; + resvec(iter) = presn; if (norm (x - x_old, inf) <= eps) - flag = 3; - break + flag = 3; # Stagnation: no change between iterations + break; endif restart_it++ ; iter++; endwhile - if (presn > rtol * prec_b_norm) - flag = 0; + if (nargout > 1) + ## Calculate extra outputs as requested + relres = presn / prec_b_norm; + if (relres <= rtol) + flag = 0; # Converged to solution within tolerance + endif + + resvec = resvec(1:iter-1); + it = [ceil(iter / restart), rem(iter, restart)]; endif - resids = resids(1:iter-1); - it = [ceil(iter / restart), rem(iter, restart)]; - endfunction +%!demo +%! dim = 20; +%! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim); +%! b = ones (dim, 1); +%! [x, flag, relres, iter, resvec] = gmres (A, b, 10, 1e-10, dim, @(x) x ./ diag (A), [], b) + %!shared A, b, dim %! dim = 100; %!test @@ -210,7 +220,7 @@ %! A = spdiags ([[1./(2:2:2*(dim-1)) 0]; 1./(1:2:2*dim-1); [0 1./(2:2:2*(dim-1))]]', -1:1, dim, dim); %! A = A'*A; %! b = rand (dim, 1); -%! [x, resids] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag (A), [], []); +%! [x, resvec] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag (A), [], []); %! assert (x, A\b, 1e-9*norm (x, Inf)); %! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag (diag (A)) \ x, [], []); %! assert (x, A\b, 1e-9*norm (x, Inf)); @@ -218,3 +228,10 @@ %! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x./diag(A), [], []); %! assert (x, A\b, 1e-7*norm (x, Inf)); + +%!error gmres (1) +%!error gmres (1,2,3,4,5,6,7,8,9) +%!error <A must be> gmres ({1},2) +%!error <A must be a function or matrix> gmres ({1},2) +%!error <M1 must be a function or matrix> gmres (1,2,3,4,5,{6}) +%!error <M2 must be a function or matrix> gmres (1,2,3,4,5,6,{7})