--- a/scripts/sparse/pcg.m
+++ b/scripts/sparse/pcg.m
@@ -18,20 +18,20 @@
 ## 02110-1301, USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M}, @var{x0}, @dots{})
+## @deftypefn {Function File} {@var{x} =} pcg (@var{a}, @var{b}, @var{tol}, @var{maxit}, @var{m}, @var{x0}, @dots{})
 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@dots{})
 ##
-## Solves the linear system of equations @code{@var{A} * @var{x} =
+## Solves the linear system of equations @code{@var{a} * @var{x} =
 ## @var{b}} by means of the  Preconditioned Conjugate Gradient iterative
 ## method. The input arguments are
 ##
 ## @itemize
 ## @item
-## @var{A} can be either a square (preferably sparse) matrix or a
+## @var{a} can be either a square (preferably sparse) matrix or a
 ## function handle, inline function or string containing the name
-## of a function which computes @code{@var{A} * @var{x}}. In principle
-## @var{A} should be symmetric and positive definite; if @code{pcg}
-## finds @var{A} to not be positive definite, you will get a warning
+## of a function which computes @code{@var{a} * @var{x}}. In principle
+## @var{a} should be symmetric and positive definite; if @code{pcg}
+## finds @var{a} to not be positive definite, you will get a warning
 ## message and the @var{flag} output parameter will be set.
 ## 
 ## @item
@@ -39,8 +39,8 @@
 ## 
 ## @item
 ## @var{tol} is the required relative tolerance for the residual error,
-## @code{@var{b} - @var{A} * @var{x}}. The iteration stops if @code{norm
-## (@var{b} - @var{A} * @var{x}) <= @var{tol} * norm (@var{b} - @var{A} *
+## @code{@var{b} - @var{a} * @var{x}}. The iteration stops if @code{norm
+## (@var{b} - @var{a} * @var{x}) <= @var{tol} * norm (@var{b} - @var{a} *
 ## @var{x0})}. If @var{tol} is empty or is omitted, the function sets
 ## @code{@var{tol} = 1e-6} by default.
 ## 
@@ -50,15 +50,15 @@
 ## arguments, a default value equal to 20 is used.
 ## 
 ## @item
-## @var{M} is the (left) preconditioning matrix, so that the iteration is
+## @var{m} is the (left) preconditioning matrix, so that the iteration is
 ## (theoretically) equivalent to solving by @code{pcg} @code{@var{P} *
-## @var{x} = @var{M} \ @var{b}}, with @code{@var{P} = @var{M} \ @var{A}}.
+## @var{x} = @var{m} \ @var{b}}, with @code{@var{P} = @var{m} \ @var{a}}.
 ## Note that a proper choice of the preconditioner may dramatically
 ## improve the overall performance of the method. Instead of matrix
-## @var{M}, the user may pass a function which returns the results of 
-## applying the inverse of @var{M} to a vector (usually this is the
+## @var{m}, the user may pass a function which returns the results of 
+## applying the inverse of @var{m} to a vector (usually this is the
 ## preferred way of using the preconditioner). If @code{[]} is supplied
-## for @var{M}, or @var{M} is omitted, no preconditioning is applied.
+## for @var{m}, or @var{m} is omitted, no preconditioning is applied.
 ## 
 ## @item
 ## @var{x0} is the initial guess. If @var{x0} is empty or omitted, the 
@@ -66,14 +66,14 @@
 ## @end itemize
 ## 
 ## The arguments which follow @var{x0} are treated as parameters, and
-## passed in a proper way to any of the functions (@var{A} or @var{M})
+## passed in a proper way to any of the functions (@var{a} or @var{m})
 ## which are passed to @code{pcg}. See the examples below for further
 ## details. The output arguments are
 ##
 ## @itemize
 ## @item
 ## @var{x} is the computed approximation to the solution of
-## @code{@var{A} * @var{x} = @var{b}}.
+## @code{@var{a} * @var{x} = @var{b}}.
 ## 
 ## @item
 ## @var{flag} reports on the convergence. @code{@var{flag} = 0} means
@@ -95,23 +95,23 @@
 ## @code{@var{resvec} (i,2)} is the preconditioned residual norm,
 ## after the (@var{i}-1)-th iteration, @code{@var{i} =
 ## 1,2,...@var{iter}+1}. The preconditioned residual norm is defined as
-## @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{M} \ @var{r})} where
-## @code{@var{r} = @var{b} - @var{A} * @var{x}}, see also the
-## description of @var{M}. If @var{eigest} is not required, only
+## @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{m} \ @var{r})} where
+## @code{@var{r} = @var{b} - @var{a} * @var{x}}, see also the
+## description of @var{m}. If @var{eigest} is not required, only
 ## @code{@var{resvec} (:,1)} is returned.
 ## 
 ## @item
 ## @var{eigest} returns the estimate for the smallest @code{@var{eigest}
 ## (1)} and largest @code{@var{eigest} (2)} eigenvalues of the
-## preconditioned matrix @code{@var{P} = @var{M} \ @var{A}}. In 
+## preconditioned matrix @code{@var{P} = @var{m} \ @var{a}}. In 
 ## particular, if no preconditioning is used, the extimates for the
-## extreme eigenvalues of @var{A} are returned. @code{@var{eigest} (1)}
+## extreme eigenvalues of @var{a} are returned. @code{@var{eigest} (1)}
 ## is an overestimate and @code{@var{eigest} (2)} is an underestimate, 
 ## so that @code{@var{eigest} (2) / @var{eigest} (1)} is a lower bound
 ## for @code{cond (@var{P}, 2)}, which nevertheless in the limit should
 ## theoretically be equal to the actual value of the condition number. 
 ## The method which computes @var{eigest} works only for symmetric positive
-## definite @var{A} and @var{M}, and the user is responsible for
+## definite @var{a} and @var{m}, and the user is responsible for
 ## verifying this assumption. 
 ## @end itemize
 ## 
@@ -133,7 +133,7 @@
 ## @end example
 ## 
 ## @sc{Example 2:} @code{pcg} with a function which computes
-## @code{@var{A} * @var{x}}
+## @code{@var{a} * @var{x}}
 ## 
 ## @example
 ## @group
@@ -147,7 +147,7 @@
 ## 
 ## @sc{Example 3:} Preconditioned iteration, with full diagnostics. The
 ## preconditioner (quite strange, because even the original matrix
-## @var{A} is trivial) is defined as a function
+## @var{a} is trivial) is defined as a function
 ## 
 ## @example
 ## @group
@@ -190,31 +190,12 @@
 ## @seealso{sparse, pcr}
 ## @end deftypefn
 
-## REVISION HISTORY
-##
-## 2004-05-21, Piotr Krzyzanowski:
-##	Added 4 demos and 4 tests
-##
-## 2004-05-18, Piotr Krzyzanowski:
-##	Warnings use warning() function now
-##
-## 2004-04-29, Piotr Krzyzanowski:
-##	Added more warning messages when FLAG is not required
-##
-## 2004-04-28, Piotr Krzyzanowski:
-## 	When eigest is required, resvec returns both the Euclidean and the
-##	preconditioned residual norm convergence history
-##
-## 2004-04-20, Piotr Krzyzanowski: 
-## 	Corrected eigenvalue estimator. Changed the tridiagonal matrix
-##	eigenvalue solver to regular eig
-## 
+## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
 
-function [x, flag, relres, iter, resvec, eigest] = ...
-		pcg( A, b, tol, maxit, M, x0, varargin )
+function [x, flag, relres, iter, resvec, eigest] = pcg (A, b, tol, maxit, M, x0, varargin)
 
-  if ((nargin < 6) || isempty(x0))
-    x = zeros(size(b));
+  if (nargin < 6 || isempty (x0))
+    x = zeros (size (b));
   else
     x = x0;
   endif
@@ -223,125 +204,123 @@
     M = [];
   endif
 
-  if ((nargin < 4) || isempty(maxit))
-    maxit = min(size(b,1),20);
+  if (nargin < 4 || isempty (maxit))
+    maxit = min (size (b, 1), 20);
   endif
 
-  maxit = maxit + 2;
+  maxit += 2;
 
-  if ((nargin < 3) || isempty(tol))
+  if (nargin < 3 || isempty (tol))
     tol = 1e-6;
   endif
 
   preconditioned_residual_out = false;
   if (nargout > 5)
-    T = zeros(maxit,maxit);
+    T = zeros (maxit, maxit);
     preconditioned_residual_out = true;
   endif
 
   matrix_positive_definite = true;	# assume A is positive definite
 
-  p = zeros(size(b)); 
+  p = zeros (size (b));
   oldtau = 1; 
-  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!
-    r = b - feval(A,x,varargin{:});
+    r = b - feval (A, x, varargin{:});
   endif
 
-  resvec(1,1) = norm(r); 
+  resvec(1,1) = norm (r);
   alpha = 1;
   iter = 2;
 
-  while ((resvec(iter-1,1) > tol*resvec(1,1)) && (iter < maxit))
-    if (isnumeric(M))		# is M a matrix?
-      if isempty(M)		# if M is empty, use no precond
+  while (resvec(iter-1,1) > tol*resvec(1,1) && iter < maxit)
+    if (isnumeric (M))		# is M a matrix?
+      if (isempty (M))		# if M is empty, use no precond
 	z = r;
       else			# otherwise, apply the precond
 	z = M \ r;
       endif
     else			# then M should be a function!
-      z = feval(M,r,varargin{:});
+      z = feval (M, r, varargin{:});
     endif
-    tau = z'*r; 
-    resvec(iter-1,2) = sqrt(tau);
-    beta = tau/oldtau;
+    tau = z' * r; 
+    resvec(iter-1,2) = sqrt (tau);
+    beta = tau / oldtau;
     oldtau = tau;
     p = z + beta*p;
-    if (isnumeric(A))		# is A a matrix?
-      w = A*p;
+    if (isnumeric (A))		# is A a matrix?
+      w = A * p;
     else			# then A should be a function!
-      w = feval(A,p,varargin{:});
+      w = feval (A, p, varargin{:});
     endif
     oldalpha = alpha; 		# needed only for eigest
-    alpha = tau/(p'*w);
+    alpha = tau / (p'*w);
     if (alpha <= 0.0) # negative matrix?
       matrix_positive_definite = false;
     endif
-    x = x + alpha*p;
-    r = r - alpha*w;
-    if ((nargout > 5) && (iter > 2))
+    x += alpha*p;
+    r -= alpha*w;
+    if (nargout > 5 && iter > 2)
       T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
 	  [1 sqrt(beta); sqrt(beta) beta]./oldalpha;
       ## EVS = eig(T(2:iter-1,2:iter-1));
       ## fprintf(stderr,"PCG condest: %g (iteration: %d)\n", max(EVS)/min(EVS),iter);
     endif
-    resvec(iter,1) = norm(r);
-    iter = iter + 1;
+    resvec(iter,1) = norm (r);
+    iter++;
   endwhile
 
   if (nargout > 5)
-    if (matrix_positive_definite )
+    if (matrix_positive_definite)
       if (iter > 3)
 	T = T(2:iter-2,2:iter-2);
 	l = eig(T);
-	eigest = [min(l) max(l)];
-	## fprintf(stderr, "PCG condest: %g\n",eigest(2)/eigest(1));
+	eigest = [min(l), max(l)];
+	## fprintf (stderr, "PCG condest: %g\n", eigest(2)/eigest(1));
       else
-	eigest = [NaN NaN];
-	warning("PCG: eigenvalue estimate failed: iteration converged too fast.");
+	eigest = [NaN, NaN];
+	warning ("PCG: eigenvalue estimate failed: iteration converged too fast.");
       endif
     else
-      eigest = [NaN NaN];
+      eigest = [NaN, NaN];
     endif
 
     ## apply the preconditioner once more and finish with the precond
     ## residual
-    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
 	z = r;
       else			# otherwise, apply the precond
-	z = M\r;
+	z = M \ r;
       endif
     else			# then M should be a function!
-      z = feval(M,r,varargin{:});
+      z = feval (M, r, varargin{:});
     endif
-    resvec(iter-1,2) = sqrt(r'*z);
+    resvec(iter-1,2) = sqrt (r'*z);
   else
-    resvec = resvec(:,1);	
+    resvec = resvec(:,1);
   endif
 
   flag = 0;
   relres = resvec(iter-1,1)./resvec(1,1);
-  iter = iter - 2;
-  if (iter >= (maxit-2))
+  iter -= 2;
+  if (iter >= maxit-2)
     flag = 1;
     if (nargout < 2)
-      warning("PCG: maximum number of iterations (%d) reached\n", iter);
-      warning("The initial residual norm was reduced %g times.\n", 1.0/relres);
+      warning ("PCG: maximum number of iterations (%d) reached\n", iter);
+      warning ("The initial residual norm was reduced %g times.\n", 1.0/relres);
     endif
-  else
-    if (nargout < 2)
-      fprintf(stderr, "PCG: converged in %d iterations. ", iter);
-      fprintf(stderr, "The initial residual norm was reduced %g times.\n",...
-	      1.0/relres);
-    endif
+  elseif (nargout < 2)
+    fprintf (stderr, "PCG: converged in %d iterations. ", iter);
+    fprintf (stderr, "The initial residual norm was reduced %g times.\n",...
+	     1.0/relres);
   endif
 
-  if (!matrix_positive_definite)
+  if (! matrix_positive_definite)
     flag = 3;
     if (nargout < 2)
-      warning("PCG: matrix not positive definite?\n");
+      warning ("PCG: matrix not positive definite?\n");
     endif
   endif
 endfunction