--- a/scripts/sparse/pcr.m
+++ b/scripts/sparse/pcr.m
@@ -100,9 +100,9 @@
 ## 
 ## @example
 ## @group
-## 	n = 10; 
-## 	a = sparse (diag (1:n));
-## 	b = rand (N, 1);
+##      n = 10; 
+##      a = sparse (diag (1:n));
+##      b = rand (N, 1);
 ## @end group
 ## @end example
 ## 
@@ -131,10 +131,10 @@
 ## 
 ## @example
 ## @group
-##   function y = apply_m (x)		
+##   function y = apply_m (x)           
 ##     k = floor (length(x)-2); 
 ##     y = x; 
-##     y(1:k) = x(1:k)./[1:k]';	
+##     y(1:k) = x(1:k)./[1:k]'; 
 ##   endfunction
 ## 
 ##   [x, flag, relres, iter, resvec] = ...
@@ -150,7 +150,7 @@
 ## @group
 ##   function y = apply_m (x, varargin)
 ##     k = varargin@{1@}; 
-##     y = x; y(1:k) = x(1:k)./[1:k]';	 
+##     y = x; y(1:k) = x(1:k)./[1:k]';   
 ##   endfunction
 ## 
 ##   [x, flag, relres, iter, resvec] = ...
@@ -160,8 +160,8 @@
 ## 
 ## @sc{References}
 ## 
-## 	[1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of
-##  	Equations", section 9.5.4; Springer, 1994
+##      [1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of
+##      Equations", section 9.5.4; Springer, 1994
 ##
 ## @seealso{sparse, pcg}
 ## @end deftypefn
@@ -197,19 +197,19 @@
   endif
 
   ##  init
-  if (isnumeric (a))		# is A a matrix?
+  if (isnumeric (a))            # is A a matrix?
     r = b - a*x;
-  else				# then A should be a function!
+  else                          # then A should be a function!
     r = b - feval (a, x, varargin{:});
   endif
 
-  if (isnumeric (m))		# is M a matrix?
-    if (isempty (m))		# if M is empty, use no precond
+  if (isnumeric (m))            # is M a matrix?
+    if (isempty (m))            # if M is empty, use no precond
       p = r;
-    else			# otherwise, apply the precond
+    else                        # otherwise, apply the precond
       p = m \ r;
     endif
-  else				# then M should be a function!
+  else                          # then M should be a function!
     p = feval (m, r, varargin{:});
   endif
 
@@ -218,58 +218,58 @@
   b_bot_old = 1;
   q_old = p_old = s_old = zeros (size (x));
 
-  if (isnumeric (a))		# is A a matrix?
+  if (isnumeric (a))            # is A a matrix?
     q = a * p;
-  else				# then A should be a function!
+  else                          # then A should be a function!
     q = feval (a, p, varargin{:});
   endif
-	
+        
   resvec(1) = abs (norm (r)); 
 
   ## iteration
   while (resvec(iter-1) > tol*resvec(1) && iter < maxit)
 
-    if (isnumeric (m))		# is M a matrix?
-      if (isempty (m))		# if M is empty, use no precond
-	s = q;
-      else			# otherwise, apply the precond
-	s = m \ q;
+    if (isnumeric (m))          # is M a matrix?
+      if (isempty (m))          # if M is empty, use no precond
+        s = q;
+      else                      # otherwise, apply the precond
+        s = m \ q;
       endif
-    else			# then M should be a function!
+    else                        # then M should be a function!
       s = feval (m, q, varargin{:});
     endif
     b_top = r' * s;
     b_bot = q' * s;
-	
+        
     if (b_bot == 0.0)
       breakdown = true;
       break;
     endif
     lambda = b_top / b_bot;
-	
+        
     x += lambda*p;
     r -= lambda*q;
-	
-    if (isnumeric(a))		# is A a matrix?
+        
+    if (isnumeric(a))           # is A a matrix?
       t = a*s;
-    else			# then A should be a function!
+    else                        # then A should be a function!
       t = feval (a, s, varargin{:});
     endif
-	
+        
     alpha0 = (t'*s) / b_bot;
     alpha1 = (t'*s_old) / b_bot_old;
-	
+        
     p_temp = p;
     q_temp = q;
 
     p = s - alpha0*p - alpha1*p_old;
     q = t - alpha0*q - alpha1*q_old;
-	
+        
     s_old = s;
     p_old = p_temp;
     q_old = q_temp;
     b_bot_old = b_bot;
-	
+        
     resvec(iter) = abs (norm (r));
     iter++;
   endwhile
@@ -286,7 +286,7 @@
   elseif (nargout < 2 && ! breakdown)
     fprintf (stderr, "pcr: converged in %d iterations. \n", iter);
     fprintf (stderr, "the initial residual norm was reduced %g times.\n",
-	     1.0 / relres);
+             1.0 / relres);
   endif
 
   if (breakdown)
@@ -301,130 +301,130 @@
 
 %!demo
 %!
-%!	# Simplest usage of PCR (see also 'help pcr')
+%!      # Simplest usage of PCR (see also 'help pcr')
 %!
-%!	N = 20; 
-%!	A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution
-%!  	x = pcr(A,b);
-%!	printf('The solution relative error is %g\n', norm(x-y)/norm(y));
+%!      N = 20; 
+%!      A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution
+%!      x = pcr(A,b);
+%!      printf('The solution relative error is %g\n', norm(x-y)/norm(y));
 %!
-%!	# You shouldn't be afraid if PCR issues some warning messages in this
-%!	# example: watch out in the second example, why it takes N iterations 
-%!	# of PCR to converge to (a very accurate, by the way) solution
+%!      # You shouldn't be afraid if PCR issues some warning messages in this
+%!      # example: watch out in the second example, why it takes N iterations 
+%!      # of PCR to converge to (a very accurate, by the way) solution
 %!demo
 %!
-%!	# Full output from PCR
-%!	# We use this output to plot the convergence history  
+%!      # Full output from PCR
+%!      # We use this output to plot the convergence history  
 %!
-%!	N = 20; 
-%!	A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution
-%!  	[x, flag, relres, iter, resvec] = pcr(A,b);
-%!	printf('The solution relative error is %g\n', norm(x-X)/norm(X));
-%!	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
-%!	semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;');
+%!      N = 20; 
+%!      A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution
+%!      [x, flag, relres, iter, resvec] = pcr(A,b);
+%!      printf('The solution relative error is %g\n', norm(x-X)/norm(X));
+%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
+%!      semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;');
 %!demo
 %!
-%!	# Full output from PCR
-%!	# We use indefinite matrix based on the Hilbert matrix, with one
-%!	# strongly negative eigenvalue
-%!	# Hilbert matrix is extremely ill conditioned, so is ours, 
-%!	# and that's why PCR WILL have problems
+%!      # Full output from PCR
+%!      # We use indefinite matrix based on the Hilbert matrix, with one
+%!      # strongly negative eigenvalue
+%!      # Hilbert matrix is extremely ill conditioned, so is ours, 
+%!      # and that's why PCR WILL have problems
 %!
-%!	N = 10; 
-%!	A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution
-%!	printf('Condition number of A is   %g\n', cond(A));
-%!  	[x, flag, relres, iter, resvec] = pcr(A,b,[],200);
-%!	if (flag == 3)
-%!	  printf('PCR breakdown. System matrix is [close to] singular\n');
-%!	end
-%!	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
-%!	semilogy([0:iter],resvec,'o-g;absolute residual;');
+%!      N = 10; 
+%!      A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution
+%!      printf('Condition number of A is   %g\n', cond(A));
+%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],200);
+%!      if (flag == 3)
+%!        printf('PCR breakdown. System matrix is [close to] singular\n');
+%!      end
+%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
+%!      semilogy([0:iter],resvec,'o-g;absolute residual;');
 %!demo
 %!
-%!	# Full output from PCR
-%!	# We use an indefinite matrix based on the 1-D Laplacian matrix for A, 
-%!	# and here we have cond(A) = O(N^2)
-%!	# That's the reason we need some preconditioner; here we take
-%!	# a very simple and not powerful Jacobi preconditioner, 
-%!	# which is the diagonal of A
+%!      # Full output from PCR
+%!      # We use an indefinite matrix based on the 1-D Laplacian matrix for A, 
+%!      # and here we have cond(A) = O(N^2)
+%!      # That's the reason we need some preconditioner; here we take
+%!      # a very simple and not powerful Jacobi preconditioner, 
+%!      # which is the diagonal of A
 %!
-%!	# Note that we use here indefinite preconditioners!
+%!      # Note that we use here indefinite preconditioners!
 %!
-%!	N = 100; 
-%!	A = zeros(N,N);
-%!	for i=1:N-1 # form 1-D Laplacian matrix
-%!		A(i:i+1,i:i+1) = [2 -1; -1 2];
-%!	endfor
-%!	A = [A, zeros(size(A)); zeros(size(A)), -A];
-%!	b = rand(2*N,1); X = A\b; #X is the true solution
-%!	maxit = 80;
-%!	printf('System condition number is %g\n',cond(A));
-%!	# No preconditioner: the convergence is very slow!
+%!      N = 100; 
+%!      A = zeros(N,N);
+%!      for i=1:N-1 # form 1-D Laplacian matrix
+%!              A(i:i+1,i:i+1) = [2 -1; -1 2];
+%!      endfor
+%!      A = [A, zeros(size(A)); zeros(size(A)), -A];
+%!      b = rand(2*N,1); X = A\b; #X is the true solution
+%!      maxit = 80;
+%!      printf('System condition number is %g\n',cond(A));
+%!      # No preconditioner: the convergence is very slow!
 %!
-%!  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit);
-%!	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
-%!	semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;');
+%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit);
+%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
+%!      semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;');
 %!
-%!	pause(1);
-%!	# Test Jacobi preconditioner: it will not help much!!!
+%!      pause(1);
+%!      # Test Jacobi preconditioner: it will not help much!!!
 %!
-%!	M = diag(diag(A)); # Jacobi preconditioner
-%!  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
-%!	hold on;
-%!	semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;');
+%!      M = diag(diag(A)); # Jacobi preconditioner
+%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
+%!      hold on;
+%!      semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;');
 %!
-%!	pause(1);
-%!	# Test nonoverlapping block Jacobi preconditioner: this one should give
-%!	# some convergence speedup!
+%!      pause(1);
+%!      # Test nonoverlapping block Jacobi preconditioner: this one should give
+%!      # some convergence speedup!
 %!
-%!	M = zeros(N,N);k=4;
-%!	for i=1:k:N # get k x k diagonal blocks of A
-%!		M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
-%!	endfor
-%!	M = [M, zeros(size(M)); zeros(size(M)), -M];
-%!  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
-%!	semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;');
-%!	hold off;
+%!      M = zeros(N,N);k=4;
+%!      for i=1:k:N # get k x k diagonal blocks of A
+%!              M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
+%!      endfor
+%!      M = [M, zeros(size(M)); zeros(size(M)), -M];
+%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
+%!      semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;');
+%!      hold off;
 %!test
 %!
-%!	#solve small indefinite diagonal system
+%!      #solve small indefinite diagonal system
 %!
-%!	N = 10; 
-%!	A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution
-%!  	[x, flag] = pcr(A,b,[],N+1);
-%!	assert(norm(x-X)/norm(X)<1e-10);
-%!	assert(flag,0);
+%!      N = 10; 
+%!      A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution
+%!      [x, flag] = pcr(A,b,[],N+1);
+%!      assert(norm(x-X)/norm(X)<1e-10);
+%!      assert(flag,0);
 %!
 %!test
 %!
-%!	#solve tridiagonal system, do not converge in default 20 iterations
-%!	#should perform max allowable default number of iterations
+%!      #solve tridiagonal system, do not converge in default 20 iterations
+%!      #should perform max allowable default number of iterations
 %!
-%!	N = 100; 
-%!	A = zeros(N,N);
-%!	for i=1:N-1 # form 1-D Laplacian matrix
-%!		A(i:i+1,i:i+1) = [2 -1; -1 2];
-%!	endfor
-%!	b = ones(N,1); X = A\b; #X is the true solution
-%!  	[x, flag, relres, iter, resvec] = pcr(A,b,1e-12);
-%!	assert(flag,1);
-%!	assert(relres>0.6);
-%!	assert(iter,20);
+%!      N = 100; 
+%!      A = zeros(N,N);
+%!      for i=1:N-1 # form 1-D Laplacian matrix
+%!              A(i:i+1,i:i+1) = [2 -1; -1 2];
+%!      endfor
+%!      b = ones(N,1); X = A\b; #X is the true solution
+%!      [x, flag, relres, iter, resvec] = pcr(A,b,1e-12);
+%!      assert(flag,1);
+%!      assert(relres>0.6);
+%!      assert(iter,20);
 %!
 %!test
 %!
-%!	#solve tridiagonal system with 'prefect' preconditioner
-%!	#converges in one iteration
+%!      #solve tridiagonal system with 'prefect' preconditioner
+%!      #converges in one iteration
 %!
-%!	N = 100; 
-%!	A = zeros(N,N);
-%!	for i=1:N-1 # form 1-D Laplacian matrix
-%!		A(i:i+1,i:i+1) = [2 -1; -1 2];
-%!	endfor
-%!	b = ones(N,1); X = A\b; #X is the true solution
-%!  	[x, flag, relres, iter] = pcr(A,b,[],[],A,b);
-%!	assert(norm(x-X)/norm(X)<1e-6);
-%!	assert(relres<1e-6);
-%!	assert(flag,0);
-%!	assert(iter,1); #should converge in one iteration
+%!      N = 100; 
+%!      A = zeros(N,N);
+%!      for i=1:N-1 # form 1-D Laplacian matrix
+%!              A(i:i+1,i:i+1) = [2 -1; -1 2];
+%!      endfor
+%!      b = ones(N,1); X = A\b; #X is the true solution
+%!      [x, flag, relres, iter] = pcr(A,b,[],[],A,b);
+%!      assert(norm(x-X)/norm(X)<1e-6);
+%!      assert(relres<1e-6);
+%!      assert(flag,0);
+%!      assert(iter,1); #should converge in one iteration
 %!