--- a/scripts/sparse/gmres.m	Thu Jun 14 22:11:04 2012 -0700
+++ b/scripts/sparse/gmres.m	Fri Jun 15 10:15:07 2012 -0700
@@ -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})