--- a/scripts/sparse/bicg.m
+++ b/scripts/sparse/bicg.m
@@ -15,30 +15,30 @@
 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-##  
+##
 ## @deftypefn {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
 ## @deftypefnx {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicg (@var{A}, @var{b}, ...)
 ## Solve @code{A x = b} using the Bi-conjugate gradient iterative method.
 ##
 ## @itemize @minus
-## @item @var{rtol} is the relative tolerance, if not given 
+## @item @var{rtol} is the relative tolerance, if not given
 ## or set to [] the default value 1e-6 is used.
-## @item @var{maxit} the maximum number of outer iterations, 
-## if not given or set to [] the default value 
+## @item @var{maxit} the maximum number of outer iterations,
+## if not given or set to [] the default value
 ## @code{min (20, numel (b))} is used.
-## @item @var{x0} the initial guess, if not given or set to [] 
-## the default value @code{zeros (size (b))} is used. 
+## @item @var{x0} the initial guess, if not given or set to []
+## the default value @code{zeros (size (b))} is used.
 ## @end itemize
 ##
-## @var{A} can be passed as a matrix or as a function handle or 
+## @var{A} can be passed as a matrix or as a function handle or
 ## inline function @code{f} such that @code{f(x, "notransp") = A*x}
 ## and @code{f(x, "transp") = A'*x}.
 ##
-## The preconditioner @var{P} is given as @code{P = M1 * M2}. 
-## Both @var{M1} and @var{M2} can be passed as a matrix or as 
-## a function handle or inline function @code{g} such that 
-## @code{g(x, 'notransp') = M1 \ x} or @code{g(x, 'notransp') = M2 \ x} and 
+## The preconditioner @var{P} is given as @code{P = M1 * M2}.
+## Both @var{M1} and @var{M2} can be passed as a matrix or as
+## a function handle or inline function @code{g} such that
+## @code{g(x, 'notransp') = M1 \ x} or @code{g(x, 'notransp') = M2 \ x} and
 ## @code{g(x, 'transp') = M1' \ x} or @code{g(x, 'transp') = M2' \ x}.
 ##
 ## If colled with more than one output parameter
@@ -52,7 +52,7 @@
 ## @end itemize
 ## (the value 2 is unused but skipped for compatibility).
 ## @item @var{relres} is the final value of the relative residual.
-## @item @var{iter} is the number of iterations performed. 
+## @item @var{iter} is the number of iterations performed.
 ## @item @var{resvec} is a vector containing the relative residual at each iteration.
 ## @end itemize
 ##
@@ -64,7 +64,7 @@
 function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0)
 
   if (nargin >= 2 && isvector (full (b)))
-    
+
     if (ischar (A))
       fun = str2func (A);
       Ax  = @(x) feval (fun, x, "notransp");
@@ -105,7 +105,7 @@
       error (["bicg: preconditioner is expected to " ...
               "be a function or matrix"]);
     endif
-    
+
     if (nargin < 6 || isempty (M2))
       M2m1x = @(x, ignore) x;
       M2tm1x = M2m1x;
@@ -144,44 +144,44 @@
     res0  = Inf;
 
     if (any (r0 != 0))
-      
+
       for k = 1:maxit
 
-	a  = (s0' * Pm1x (r0)) ./ (f' * Ax (d));
+        a  = (s0' * Pm1x (r0)) ./ (f' * Ax (d));
 
-	x += a * d;
-	y += conj (a) * f;
+        x += a * d;
+        y += conj (a) * f;
 
-	r1 = r0 - a * Ax (d);
-	s1 = s0 - conj (a) * Atx (f);
+        r1 = r0 - a * Ax (d);
+        s1 = s0 - conj (a) * Atx (f);
 
-	beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0));
+        beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0));
 
-	d = Pm1x (r1) + beta * d;
-	f = Ptm1x (s1) + conj (beta) * f;
+        d = Pm1x (r1) + beta * d;
+        f = Ptm1x (s1) + conj (beta) * f;
 
-	r0 = r1;
-	s0 = s1;
+        r0 = r1;
+        s0 = s1;
 
-	res1 = norm (b - Ax (x)) / bnorm;
-	if (res1 < tol)
+        res1 = norm (b - Ax (x)) / bnorm;
+        if (res1 < tol)
           flag = 0;
           if (nargout < 2)
-	    printf ("bicg converged at iteration %i ", k);
+            printf ("bicg converged at iteration %i ", k);
             printf ("to a solution with relative residual %e\n", res1);
           endif
-	  break;
-	endif
-	
-	if (res0 <= res1)
+          break;
+        endif
+
+        if (res0 <= res1)
           flag = 3;
-	  printf ("bicg stopped at iteration %i ", k);
+          printf ("bicg stopped at iteration %i ", k);
           printf ("without converging to the desired tolerance %e\n", tol);
           printf ("because the method stagnated.\n");
           printf ("The iterate returned (number %i) ", k-1);
           printf ("has relative residual %e\n", res0);
-	  break
-	endif
+          break
+        endif
         res0 = res1;
         if (nargout > 4)
           resvec(k) = res0;
@@ -190,7 +190,7 @@
 
       if (k == maxit)
         flag = 1;
-	printf ("bicg stopped at iteration %i ", maxit);
+        printf ("bicg stopped at iteration %i ", maxit);
         printf ("without converging to the desired tolerance %e\n", tol);
         printf ("because the maximum number of iterations was reached. ");
         printf ("The iterate returned (number %i) has ", maxit);
@@ -203,33 +203,33 @@
         printf ("bicg converged after 0 interations\n");
       endif
     endif
-  
+
   else
     print_usage ();
   endif
-  
+
 endfunction;
 
 
 %!test
-%! n = 100; 
+%! n = 100;
 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
-%! b = sum (A, 2); 
-%! tol = 1e-8; 
+%! b = sum (A, 2);
+%! tol = 1e-8;
 %! maxit = 15;
-%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); 
-%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); 
+%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
+%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, maxit, M1, M2);
 %! assert (x, ones (size (b)), 1e-7);
 %!
 %!test
-%! n = 100; 
+%! n = 100;
 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
-%! b = sum (A, 2); 
-%! tol = 1e-8; 
+%! b = sum (A, 2);
+%! tol = 1e-8;
 %! maxit = 15;
-%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); 
-%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); 
+%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
+%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
 %!
 %! function y = afun (x, t, a)
 %!  switch t
@@ -240,15 +240,15 @@
 %!  endswitch
 %! endfunction
 %!
-%! [x, flag, relres, iter, resvec] = bicg (@(x, t) afun (x, t, A), 
+%! [x, flag, relres, iter, resvec] = bicg (@(x, t) afun (x, t, A),
 %!                                         b, tol, maxit, M1, M2);
 %! assert (x, ones (size (b)), 1e-7);
 
 %!test
-%! n = 100; 
-%! tol = 1e-8; 
+%! n = 100;
+%! tol = 1e-8;
 %! a = sprand (n, n, .1);
 %! A = a' * a + 100 * eye (n);
-%! b = sum (A, 2); 
+%! b = sum (A, 2);
 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, [], diag (diag (A)));
 %! assert (x, ones (size (b)), 1e-7);