Mercurial > hg > octave-max

changeset 5838:376e02b2ce70

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-06-01 20:23:53 by jwe]
author jwe
date Thu, 01 Jun 2006 20:23:54 +0000
parents 55404f3b0da1
children 0bc78d9ddcd4
files scripts/general/bicubic.m scripts/general/gradient.m scripts/general/interp1.m scripts/general/interp2.m scripts/general/interpft.m scripts/general/polyarea.m scripts/general/quadl.m scripts/miscellaneous/inputname.m scripts/plot/__plt3__.m scripts/plot/ndgrid.m scripts/plot/plot3.m scripts/polynomial/pchip.m scripts/polynomial/spline.m scripts/sparse/pcg.m scripts/sparse/pcr.m scripts/strings/mat2str.m src/ChangeLog src/DLD-FUNCTIONS/__pchip_deriv__.cc src/Makefile.in
diffstat 19 files changed, 531 insertions(+), 535 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/general/bicubic.m
+++ b/scripts/general/bicubic.m
@@ -20,49 +20,50 @@
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
 ##
-## Bicubic interpolation method. Returns a matrix @var{zi}, corresponding the
-## the bicubic interpolations at @var{xi} and @var{yi} of the data supplied
+## Return a matrix @var{zi} corresponding to the the bicubic
+## interpolations at @var{xi} and @var{yi} of the data supplied
 ## as @var{x}, @var{y} and @var{z}. 
 ##
 ## For further information please see bicubic.pdf available at
 ## @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic}
+## @seealso{interp2}
 ## @end deftypefn
-## @seealso{interp2}
 
 ## Bicubic interpolation method.
 ## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn>
 
-function F = bicubic(X, Y, Z, XI, YI, spline_alpha)
+function F = bicubic (X, Y, Z, XI, YI, spline_alpha)
 
   if (nargin < 1 || nargin > 6)
     print_usage ();
   endif
 
-  if (nargin == 6 && prod(size(spline_alpha))==1)
+  if (nargin == 6 && prod (size (spline_alpha)) == 1)
     a = spline_alpha
   else
     a = 0.5;
   endif
 
-  if (nargin <=2) # bicubic(Z) or bicubic(Z, 2)
+  if (nargin <= 2)
+    ## bicubic (Z) or bicubic (Z, 2)
     if (nargin == 1) 
       n = 1;
     else
       n = Y;
     endif
-    Z= X;
+    Z = X;
     X = [];
-    [rz,cz] = size(Z);
-    s = linspace(1, cz, (cz-1)*pow2(n)+1);
-    t = linspace(1, rz, (rz-1)*pow2(n)+1);
+    [rz, cz] = size (Z);
+    s = linspace (1, cz, (cz-1)*pow2(n)+1);
+    t = linspace (1, rz, (rz-1)*pow2(n)+1);
   elseif (nargin == 3)
-    if (!isvector (X) || !isvector (Y))
+    if (! isvector (X) || ! isvector (Y))
       error ("XI and YI must be vector");
     endif
     s = Y;
     t = Z;
     Z = X;
-    [rz, cz] = size(Z);
+    [rz, cz] = size (Z);
   elseif (nargin == 5 || nargin == 6)
     [rz, cz] = size (Z) ; 
     if (isvector (X) && isvector (Y))
@@ -88,24 +89,24 @@
     YI(ylast_ind) = Y(rz);
 
 
-    X = reshape(X, 1, cz);
-    X(cz) *= (1+(sign(X(cz)))*eps);
+    X = reshape (X, 1, cz);
+    X(cz) *= 1 + sign (X(cz))*eps;
     if (X(cz) == 0) 
       X(cz) = eps;
     endif; 
-    XI = reshape(XI,1,length(XI));
-    [m, i] = sort([X, XI]);
-    o = cumsum ( i<= cz);
-    xidx = o([find( i> cz)]);
+    XI = reshape (XI, 1, length (XI));
+    [m, i] = sort ([X, XI]);
+    o = cumsum (i <= cz);
+    xidx = o(find (i > cz));
     
-    Y = reshape(Y, rz, 1);
-    Y(rz) *= (1+(sign(Y(rz)))*eps);
+    Y = reshape (Y, rz, 1);
+    Y(rz) *= 1 + sign (Y(rz))*eps;
     if (Y(rz) == 0) 
       Y(rz) = eps;
     endif; 
-    YI = reshape(YI,length(YI),1);
-    [m, i] = sort([Y; YI]);
-    o = cumsum ( i<= rz);
+    YI = reshape (YI, length (YI), 1);
+    [m, i] = sort ([Y; YI]);
+    o = cumsum (i <= rz);
     yidx = o([find( i> rz)]);
     
     ## set s and t used follow codes
@@ -119,21 +120,21 @@
     error ("Z at least a 3 by 3 matrices");
   endif
 
-  inds = floor(s);
-  d = find(s==cz);
-  s = s - floor(s);
+  inds = floor (s);
+  d = find (s == cz);
+  s = s - floor (s);
   inds(d) = cz-1;
   s(d) = 1.0;
   
   d = [];
-  indt = floor(t);
-  d = find(t==rz);
-  t = t - floor(t);
+  indt = floor (t);
+  d = find (t == rz);
+  t = t - floor (t);
   indt(d) = rz-1;
   t(d) = 1.0;
   d = [];
 
-  p = zeros(size(Z)+2);
+  p = zeros (size (Z) + 2);
   p(2:rz+1,2:cz+1) = Z;
   p(1,:)      =    (6*(1-a))*p(2,:)    -3*p(3,:)   + (6*a-2)*p(4,:);
   p(rz+2,:)   =    (6*(1-a))*p(rz+1,:) -3*p(rz,:)  + (6*a-2)*p(rz-1,:);
@@ -164,16 +165,16 @@
   cs2 = cs2([1,1,1,1],:);
   cs3 = cs3([1,1,1,1],:);
 
-  lent = length(ct0);
-  lens = length(cs0);
-  F = zeros(lent,lens);
+  lent = length (ct0);
+  lens = length (cs0);
+  F = zeros (lent, lens);
   
   for i = 1:lent
     it = indt(i);
     int = [it, it+1, it+2, it+3];
     F(i,:) = [ct0(i),ct1(i),ct2(i),ct3(i)] * ...
-	(p(int,inds) .* cs0 + p(int, inds+1) .* cs1 + ...
-	 p(int, inds+2) .* cs2 + p(int, inds+3) .* cs3);
+	(p(int,inds) .* cs0 + p(int,inds+1) .* cs1 + ...
+	 p(int,inds+2) .* cs2 + p(int,inds+3) .* cs3);
   endfor
 
   ## set points outside the table to NaN
--- a/scripts/general/gradient.m
+++ b/scripts/general/gradient.m
@@ -54,14 +54,14 @@
 
 function [varargout] = gradient (M, varargin)
   
-  if ((nargin < 1))
+  if (nargin < 1)
     print_usage ()
   endif
 
   transposed = false;
   if (isvector (M))
     ## make a column vector
-    transposed = (size(M,2) == 1);
+    transposed = (size (M, 2) == 1);
     M = M(:)';
   endif
 
@@ -71,25 +71,25 @@
     print_usage ()
   endif
     
-  d = cell(1,nd);
+  d = cell (1, nd);
   if (nargin == 1)
     for i=1:nd
-      d{i} = ones(sz(i), 1);
+      d{i} = ones (sz(i), 1);
     endfor
   elseif (nargin == 2)
     if (isscalar (varargin{1}))
-      for i=1:nd
-	d{i} = varargin{1} * ones(sz(i), 1);
+      for i = 1:nd
+	d{i} = varargin{1} * ones (sz(i), 1);
       endfor
     else
-      for i=1:nd
+      for i = 1:nd
 	d{i} = varargin{1};
       endfor
     endif
   else
     for i=1:nd
       if (isscalar (varargin{1}))
-	d{i} = varargin{i} * ones(sz(i), 1);
+	d{i} = varargin{i} * ones (sz(i), 1);
       else
 	d{i} = varargin{i};
       endif
@@ -101,21 +101,21 @@
     d{2} = tmp;
   endif
 
-  for i = 1:max(2,min(nd,nargout))
+  for i = 1:max (2, min (nd, nargout))
     mr = sz(i);
-    mc = prod([sz(1:i-1),sz(i+1:nd)]);
-    Y = zeros (size(M), class(M));
+    mc = prod ([sz(1:i-1), sz(i+1:nd)]);
+    Y = zeros (size (M), class (M));
 
     if (mr > 1)
       ## top and bottom boundary
-      Y(1, :) = diff (M(1:2, :)) / d{i}(1);
-      Y(mr, :) = diff (M(mr - 1:mr, :)) / d{i}(mr - 1);
+      Y(1,:) = diff (M(1:2,:)) / d{i}(1);
+      Y(mr,:) = diff (M(mr-1:mr,:)) / d{i}(mr-1);
     endif
 
     if (mr > 2)
       ## interior points
-      Y(2:mr-1, :) = (M(3:mr, :) .- M(1:mr - 2, :)) ./ ...
-          kron (d{i}(1:mr - 2) .+ d{i}(2:mr - 1), ones(1, mc));
+      Y(2:mr-1,:) = (M(3:mr,:) .- M(1:mr-2,:)) ...
+	  ./ kron (d{i}(1:mr-2) .+ d{i}(2:mr-1), ones (1, mc));
     endif
     varargout{i} = ipermute (Y, [i:nd,1:i-1]);
     M = permute (M, [2:nd,1]);
--- a/scripts/general/interp1.m
+++ b/scripts/general/interp1.m
@@ -78,7 +78,8 @@
 ## @seealso{interpft}
 ## @end deftypefn
 
-## 2000-03-25 Paul Kienzle
+## Author: Paul Kienzle
+## Date: 2000-03-25
 ##    added 'nearest' as suggested by Kai Habel
 ## 2000-07-17 Paul Kienzle
 ##    added '*' methods and matrix y
@@ -86,9 +87,9 @@
 ## 2002-01-23 Paul Kienzle
 ##    fixed extrapolation
 
-function yi = interp1(x, y, varargin)
+function yi = interp1 (x, y, varargin)
 
-  if ( nargin < 3 || nargin > 6)
+  if (nargin < 3 || nargin > 6)
     print_usage ();
   endif
 
@@ -99,13 +100,13 @@
   firstnumeric = true;
 
   if (nargin > 2)
-    for i = 1:length(varargin)
+    for i = 1:length (varargin)
       arg = varargin{i};
-      if (ischar(arg))
+      if (ischar (arg))
 	arg = tolower (arg);
-	if (strcmp("extrap",arg))
+	if (strcmp ("extrap", arg))
 	  extrap = "extrap";
-	elseif (strcmp("pp",arg))
+	elseif (strcmp ("pp", arg))
 	  pp = true;
 	else
 	  method = arg;
@@ -123,21 +124,21 @@
 
   ## reshape matrices for convenience
   x = x(:);
-  nx = size(x,1);
+  nx = size (x, 1);
   if (isvector(y) && size (y, 1) == 1)
-    y = y (:);
+    y = y(:);
   endif
   ndy = ndims (y);
-  szy = size(y);
+  szy = size (y);
   ny = szy(1);
   nc = prod (szy(2:end));
   y = reshape (y, ny, nc);
-  szx = size(xi);
+  szx = size (xi);
   xi = xi(:);
 
   ## determine sizes
   if (nx < 2 || ny < 2)
-     error ("interp1: table too short");
+    error ("interp1: table too short");
   endif
 
   ## determine which values are out of range and set them to extrap,
@@ -150,15 +151,17 @@
   else
      maxx = x(nx);
   endif
-  if (!pp)
-    if ischar(extrap) && strcmp(extrap,"extrap")
-      range=1:size(xi,1);
-      yi = zeros(size(xi,1), size(y,2));
+
+  if (! pp)
+    if (ischar (extrap) && strcmp (extrap, "extrap"))
+      range = 1:size (xi, 1);
+      yi = zeros (size (xi, 1), size (y, 2));
     else
-      range = find(xi >= minx & xi <= maxx);
-      yi = extrap*ones(size(xi,1), size(y,2));
-      if isempty(range), 
-	if (!isvector(y) && length(szx) == 2 && (szx(1) == 1 || szx(2) == 1))
+      range = find (xi >= minx & xi <= maxx);
+      yi = extrap * ones (size (xi, 1), size (y, 2));
+      if (isempty (range))
+	if (! isvector (y) && length (szx) == 2
+	    && (szx(1) == 1 || szx(2) == 1))
 	  if (szx(1) == 1)
 	    yi = reshape (yi, [szx(2), szy(2:end)]);
 	  else
@@ -173,38 +176,38 @@
     endif
   endif
 
-  if strcmp(method, "nearest")
+  if (strcmp (method, "nearest"))
     if (pp)
-      yi = mkpp ([x(1);(x(1:end-1)+x(2:end))/2;x(end)],y,szy(2:end));
+      yi = mkpp ([x(1); (x(1:end-1)+x(2:end))/2; x(end)], y, szy(2:end));
     else
-      idx = lookup(0.5*(x(1:nx-1)+x(2:nx)), xi)+1;
+      idx = lookup (0.5*(x(1:nx-1)+x(2:nx)), xi) + 1;
       yi(range,:) = y(idx,:);
     endif
-  elseif strcmp(method, "*nearest")
+  elseif (strcmp (method, "*nearest"))
     if (pp)
-      yi = mkpp ([minx; minx + [0.5:(ny-1)]'*dx; maxx],y,szy(2:end));
+      yi = mkpp ([minx; minx+[0.5:(ny-1)]'*dx; maxx], y, szy(2:end));
     else
-      idx = max(1,min(ny,floor((xi-minx)/dx+1.5)));
+      idx = max (1, min (ny, floor((xi-minx)/dx+1.5)));
       yi(range,:) = y(idx,:);
     endif
-  elseif strcmp(method, "linear")
+  elseif (strcmp (method, "linear"))
     dy = y(2:ny,:) - y(1:ny-1,:);
     dx = x(2:nx) - x(1:nx-1);
     if (pp)
-      yi = mkpp(x, [dy./dx, y(1:end-1)],szy(2:end));
+      yi = mkpp (x, [dy./dx, y(1:end-1)], szy(2:end));
     else
       ## find the interval containing the test point
-      idx = lookup (x(2:nx-1), xi)+1; 
+      idx = lookup (x(2:nx-1), xi)+1;
 				# 2:n-1 so that anything beyond the ends
 				# gets dumped into an interval
       ## use the endpoints of the interval to define a line
       s = (xi - x(idx))./dx(idx);
       yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:);
     endif
-  elseif strcmp(method, "*linear")
+  elseif (strcmp (method, "*linear"))
     if (pp)
       dy = [y(2:ny,:) - y(1:ny-1,:)];
-      yi = mkpp (minx + [0:ny-1]*dx, [dy ./ dx, y(1:end-1)], szy(2:end));
+      yi = mkpp (minx + [0:ny-1]*dx, [dy./dx, y(1:end-1)], szy(2:end));
     else
       ## find the interval containing the test point
       t = (xi - minx)/dx + 1;
@@ -215,43 +218,43 @@
       s = t - idx;
       yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); 
     endif
-  elseif strcmp(method, "pchip") || strcmp(method, "*pchip")
+  elseif (strcmp (method, "pchip") || strcmp (method, "*pchip"))
     if (nx == 2 || method(1) == "*") 
-      x = linspace(minx, maxx, ny); 
+      x = linspace (minx, maxx, ny);
     endif
     ## Note that pchip's arguments are transposed relative to interp1
     if (pp)
-      yi = pchip(x.', y.');
+      yi = pchip (x.', y.');
       yi.d = szy(2:end);
     else
-      yi(range,:) = pchip(x.', y.', xi.').';
+      yi(range,:) = pchip (x.', y.', xi.').';
     endif
 
-  elseif strcmp(method, "cubic") || (strcmp(method, "*cubic") && pp)
+  elseif (strcmp (method, "cubic") || (strcmp (method, "*cubic") && pp))
     ## FIXME Is there a better way to treat pp return return and *cubic
-    if (method(1) == "*") 
-      x = linspace(minx, maxx, ny).'; 
+    if (method(1) == "*")
+      x = linspace (minx, maxx, ny).'; 
       nx = ny;
     endif
 
     if (nx < 4 || ny < 4)
       error ("interp1: table too short");
     endif
-    idx = lookup(x(3:nx-2), xi) + 1;
+    idx = lookup (x(3:nx-2), xi) + 1;
 
     ## Construct cubic equations for each interval using divided
     ## differences (computation of c and d don't use divided differences
     ## but instead solve 2 equations for 2 unknowns). Perhaps
     ## reformulating this as a lagrange polynomial would be more efficient.
-    i=1:nx-3;
-    J = ones(1,nc);
-    dx = diff(x);
+    i = 1:nx-3;
+    J = ones (1, nc);
+    dx = diff (x);
     dx2 = x(i+1).^2 - x(i).^2;
     dx3 = x(i+1).^3 - x(i).^3;
-    a=diff(y,3)./dx(i,J).^3/6;
-    b=(diff(y(1:nx-1,:),2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2;
-    c=(diff(y(1:nx-2,:),1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J);
-    d=y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J);
+    a = diff (y, 3)./dx(i,J).^3/6;
+    b = (diff (y(1:nx-1,:), 2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2;
+    c = (diff (y(1:nx-2,:), 1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J);
+    d = y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J);
 
     if (pp)
       xs = [x(1);x(3:nx-2)];
@@ -264,7 +267,7 @@
       yi(range,:) = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ...
 		     + c(idx,:)).*xi(:,J) + d(idx,:);
     endif
-  elseif strcmp(method, "*cubic")
+  elseif (strcmp (method, "*cubic"))
     if (nx < 4 || ny < 4)
       error ("interp1: table too short");
     endif
@@ -272,7 +275,7 @@
     ## From: Miloje Makivic 
     ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html
     t = (xi - minx)/dx + 1;
-    idx = max(min(floor(t), ny-2), 2);
+    idx = max (min (floor (t), ny-2), 2);
     t = t - idx;
     t2 = t.*t;
     tp = 1 - 0.5*t;
@@ -280,28 +283,28 @@
     b = (t2 + t).*tp;
     c = (t2 - t).*tp/3;
     d = (t2 - 1).*t/6;
-    J = ones(1,nc);
+    J = ones (1, nc);
 
     yi(range,:) = a(:,J) .* y(idx,:) + b(:,J) .* y(idx+1,:) ...
 		  + c(:,J) .* y(idx-1,:) + d(:,J) .* y(idx+2,:);
 
-  elseif strcmp(method, "spline") || strcmp(method, "*spline")
+  elseif (strcmp (method, "spline") || strcmp (method, "*spline"))
     if (nx == 2 || method(1) == "*") 
       x = linspace(minx, maxx, ny); 
     endif
     ## Note that spline's arguments are transposed relative to interp1
     if (pp)
-      yi = spline(x.', y.');
+      yi = spline (x.', y.');
       yi.d = szy(2:end);
     else
-      yi(range,:) = spline(x.', y.', xi.').';
+      yi(range,:) = spline (x.', y.', xi.').';
     endif
   else
-    error(["interp1 doesn't understand method '", method, "'"]);
+    error ("interp1: invalid method '%s'", method);
   endif
 
-  if (!pp)
-    if (!isvector(y) && length(szx) == 2 && (szx(1) == 1 || szx(2) == 1))
+  if (! pp)
+    if (! isvector (y) && length (szx) == 2 && (szx(1) == 1 || szx(2) == 1))
       if (szx(1) == 1)
 	yi = reshape (yi, [szx(2), szy(2:end)]);
       else
--- a/scripts/general/interp2.m
+++ b/scripts/general/interp2.m
@@ -82,109 +82,109 @@
 ##     * Modified demo and test for new gnuplot interface
 ## 2005-09-07 Hoxide <hoxide_dirac@yahoo.com.cn>
 ##     * Add bicubic interpolation method
-##     * Fix the eat line bug when the last element of XI or YI is negative or zero.
+##     * Fix the eat line bug when the last element of XI or YI is
+##       negative or zero.
 ## 2005-11-26 Pierre Baldensperger <balden@libertysurf.fr>
 ##     * Rather big modification (XI,YI no longer need to be
 ##       "meshgridded") to be consistent with the help message
 ##       above and for compatibility.
 
-
 function ZI = interp2 (varargin)
-  Z = X = Y = XI = YI = [];
-  n = [];
+  Z = X = Y = XI = YI = n = [];
   method = "linear";
   extrapval = NaN;
 
-  switch nargin
-  case 1
-    Z = varargin{1};
-  case 2
-    if (ischar(varargin{2}))
-      [Z,method] = deal(varargin{:});
-    else
-      [Z,n] = deal(varargin{:});
-    endif
-  case 3
-    if (ischar(varargin{3}))
-      [Z,n,method] = deal(varargin{:});
-    else
-      [Z,XI,YI] = deal(varargin{:});
-    endif
-  case 4
-    if (ischar(varargin{4}))
-      [Z,XI,YI,method] = deal(varargin{:});
-    else
-      [Z,n,method,extrapval] = deal(varargin{:});
-    endif
-  case 5
-    if (ischar(varargin{4}))
-      [Z,XI,YI,method, extrapval] = deal(varargin{:});
-    else
-      [X,Y,Z,XI,YI] = deal(varargin{:});
-    endif
-  case 6 
-      [X,Y,Z,XI,YI,method] = deal(varargin{:});
-  case 7
-      [X,Y,Z,XI,YI,method,extrapval] = deal(varargin{:});
-  otherwise
-    print_usage ();
+  switch (nargin)
+    case 1
+      Z = varargin{1};
+    case 2
+      if (ischar (varargin{2}))
+	[Z, method] = deal (varargin{:});
+      else
+	[Z, n] = deal (varargin{:});
+      endif
+    case 3
+      if (ischar (varargin{3}))
+	[Z, n, method] = deal (varargin{:});
+      else
+	[Z, XI, YI] = deal (varargin{:});
+      endif
+    case 4
+      if (ischar (varargin{4}))
+	[Z, XI, YI, method] = deal (varargin{:});
+      else
+	[Z, n, method, extrapval] = deal (varargin{:});
+      endif
+    case 5
+      if (ischar (varargin{4}))
+	[Z, XI, YI, method, extrapval] = deal (varargin{:});
+      else
+	[X, Y, Z, XI, YI] = deal (varargin{:});
+      endif
+    case 6 
+	[X, Y, Z, XI, YI, method] = deal (varargin{:});
+    case 7
+	[X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:});
+    otherwise
+      print_usage ();
   endswitch
 
   ## Type checking.
-  if (!ismatrix(Z))
-    error("interp2 expected matrix Z"); 
+  if (!ismatrix (Z))
+    error ("interp2 expected matrix Z"); 
   endif
-  if (!isempty(n) && !isscalar(n))
-    error("interp2 expected scalar n"); 
+  if (!isempty (n) && !isscalar (n))
+    error ("interp2 expected scalar n"); 
   endif
-  if (!ischar(method))
-    error("interp2 expected string 'method'"); 
+  if (!ischar (method))
+    error ("interp2 expected string 'method'"); 
   endif
-  if (!isscalar(extrapval))
-    error("interp2 expected n extrapval");
+  if (!isscalar (extrapval))
+    error ("interp2 expected n extrapval");
   endif
 
-  ## Define X,Y,XI,YI if needed
+  ## Define X, Y, XI, YI if needed
   [zr, zc] = size (Z);
-  if (isempty(X))
-    X=[1:zc]; 
-    Y=[1:zr]; 
+  if (isempty (X))
+    X = 1:zc; 
+    Y = 1:zr;
   endif
-  if (!isnumeric(X) || !isnumeric(Y))
-    error("interp2 expected numeric X,Y"); 
+  if (! isnumeric (X) || ! isnumeric (Y))
+    error ("interp2 expected numeric X, Y"); 
   endif
-  if (!isempty(n))
-    p=2^n; 
-    XI=[p:p*zc]/p; 
-    YI=[p:p*zr]'/p; 
+  if (! isempty (n))
+    p = 2^n; 
+    XI = (p:p*zc)/p; 
+    YI = (p:p*zr)'/p; 
   endif
-  if (!isnumeric(XI) || !isnumeric(YI))
-    error("interp2 expected numeric XI,YI"); 
+  if (! isnumeric (XI) || ! isnumeric (YI))
+    error ("interp2 expected numeric XI, YI"); 
   endif
 
   ## If X and Y vectors produce a grid from them
   if (isvector (X) && isvector (Y))
     [X, Y] = meshgrid (X, Y);
-  elseif (! all(size (X) == size (Y)))
+  elseif (! all (size (X) == size (Y)))
     error ("X and Y must be matrices of same size");
   endif
-  if (any(size (X) != size (Z)))
+  if (any (size (X) != size (Z)))
     error ("X and Y size must match Z dimensions");
   endif
 
   ## If Xi and Yi are vectors of different orientation build a grid
-  if ((rows(XI)==1 && columns(YI)==1) || (columns(XI)==1 && rows(YI)==1))
+  if ((rows (XI) == 1 && columns (YI) == 1)
+      || (columns (XI) == 1 && rows (YI) == 1))
     [XI, YI] = meshgrid (XI, YI);
-  elseif (any(size(XI) != size(YI)))
+  elseif (any (size (XI) != size (YI)))
     error ("XI and YI must be matrices of same size");
   endif
 
-  shape = size(XI);
-  XI = reshape(XI, 1, prod(shape));
-  YI = reshape(YI, 1, prod(shape));
+  shape = size (XI);
+  XI = reshape (XI, 1, prod (shape));
+  YI = reshape (YI, 1, prod (shape));
 
-  xidx = lookup(X(1, 2:end-1), XI) + 1;
-  yidx = lookup(Y(2:end-1, 1), YI) + 1;
+  xidx = lookup (X(1, 2:end-1), XI) + 1;
+  yidx = lookup (Y(2:end-1, 1), YI) + 1;
 
   if (strcmp (method, "linear"))
     ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
@@ -197,9 +197,9 @@
     c = Z(2:zr, 1:(zc - 1)) - a;
     d = Z(2:zr, 2:zc) - a - b - c;
 
-    idx = sub2ind(size(a),yidx,xidx);
+    idx = sub2ind (size (a), yidx, xidx);
 
-    ## scale XI,YI values to a 1-spaced grid
+    ## scale XI, YI values to a 1-spaced grid
     Xsc = (XI - X(1, xidx)) ./ (X(1, xidx + 1) - X(1, xidx));
     Ysc = (YI - Y(yidx, 1)') ./ (Y(yidx + 1, 1) - Y(yidx, 1))';
 
@@ -211,12 +211,12 @@
     ytable = Y(:, 1)';
     ii = (XI - xtable(xidx) > xtable(xidx + 1) - XI);
     jj = (YI - ytable(yidx) > ytable(yidx + 1) - YI);
-    idx = sub2ind(size(Z),yidx+jj,xidx+ii);
+    idx = sub2ind (size (Z), yidx+jj, xidx+ii);
     ZI = Z(idx);
 
   elseif (strcmp (method, "cubic"))
     ## FIXME bicubic doesn't handle arbitrary XI, YI
-    ZI = bicubic(X, Y, Z, XI(1,:), YI(:,1));
+    ZI = bicubic (X, Y, Z, XI(1,:), YI(:,1));
 
   elseif (strcmp (method, "spline"))
     ## FIXME Implement 2-D (or in fact ND) spline interpolation
@@ -227,8 +227,8 @@
   endif
 
   ## set points outside the table to NaN
-  ZI( XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1) ) = extrapval;
-  ZI = reshape(ZI,shape);
+  ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1)) = extrapval;
+  ZI = reshape (ZI, shape);
 
 endfunction
 
--- a/scripts/general/interpft.m
+++ b/scripts/general/interpft.m
@@ -48,48 +48,48 @@
   endif
 
   if (nargin == 2)
-    if (isvector(x) && size(x,1) == 1)
+    if (isvector (x) && size (x, 1) == 1)
       dim = 2;
     else
       dim = 1;
     endif
   endif
 
-  if (!isscalar (n))
+  if (! isscalar (n))
     error ("interpft: n must be an integer scalar");
   endif
 
-  nd = ndims(x);
+  nd = ndims (x);
 
   if (dim < 1 || dim > nd)
     error ("interpft: integrating over invalid dimension");
   endif
 
-  perm = [dim:nd,1:(dim-1)];
-  x = permute(x, perm);
+  perm = [dim:nd, 1:(dim-1)];
+  x = permute (x, perm);
   m = size (x, 1);
 
   inc = 1;
-  while (inc * n < m)
+  while (inc*n < m)
     inc++;
   endwhile
   y = fft (x) / m;
   k = floor (m / 2);
-  sz = size(x);
+  sz = size (x);
   sz(1) = n * inc - m;
-  idx = cell(nd,1);
+  idx = cell (nd, 1);
   for i = 2:nd
     idx{i} = 1:sz(i);
   endfor
   idx{1} = 1:k;
-  z = cat (1, y(idx{:}), zeros(sz));
+  z = cat (1, y(idx{:}), zeros (sz));
   idx{1} = k+1:m;
   z = cat (1, z, y(idx{:}));
   z = n * ifft (z);
 
   if (inc != 1)
     sz(1) = n;
-    z = inc * reshape(z(1:inc:end),sz);
+    z = inc * reshape (z(1:inc:end), sz);
   endif
 
   z = ipermute (z, perm);
--- a/scripts/general/polyarea.m
+++ b/scripts/general/polyarea.m
@@ -40,7 +40,7 @@
 ##        (Is simple closed curve the technical definition of this?). 
 
 ## Author: David M. Doolin <doolin@ce.berkeley.edu>
-## Date: 1999/04/14 15:52:20 $
+## Date: 1999-04-14
 ## Modified-by: 
 ##    2000-01-15 Paul Kienzle <pkienzle@kienzle.powernet.co.uk>
 ##    * use matlab compatible interface
@@ -51,14 +51,14 @@
 function a = polyarea (x, y, dim)
   if (nargin != 2 && nargin != 3)
     print_usage ();
-  elseif any (size(x) != size(y))
-    error ("polyarea: x and y must have the same shape");
-  else
+  elseif (ndims (x) == ndims (y) && size (x) == size (y))
     if (nargin == 2)
       a = abs (sum (x .* (shift (y, -1) - shift (y, 1)))) / 2;
     else
       a = abs (sum (x .* (shift (y, -1, dim) - shift (y, 1, dim)), dim)) / 2;
     endif
+  else
+    error ("polyarea: x and y must have the same shape");
   endif
 endfunction
 
--- a/scripts/general/quadl.m
+++ b/scripts/general/quadl.m
@@ -56,19 +56,19 @@
 ##   * replace global variable terminate2 with local function need_warning
 ##   * add paper ref to docs
 
-function Q = quadl(f,a,b,tol,trace,varargin)
-  need_warning(1);
+function Q = quadl (f, a, b, tol, trace, varargin)
+  need_warning (1);
   if (nargin < 4)
-    tol=[]; 
+    tol = []; 
   endif
   if (nargin < 5)
     trace = []; 
   endif
-  if (isempty(tol))
+  if (isempty (tol))
     tol = eps; 
   endif
-  if (isempty(trace))
-    trace=0; 
+  if (isempty (trace))
+    trace = 0; 
   endif
   if (tol < eps)
     tol = eps;
@@ -78,25 +78,38 @@
   h = (b-a)/2;
   alpha = sqrt(2/3); 
   beta = 1/sqrt(5);
+
   x1 = .942882415695480; 
   x2 = .641853342345781;
   x3 = .236383199662150;
-  x = [a,m-x1*h,m-alpha*h,m-x2*h,m-beta*h,m-x3*h,m,m+x3*h,...
-       m+beta*h,m+x2*h,m+alpha*h,m+x1*h,b];
-  y = feval(f,x,varargin{:});
+
+  x = [a, m-x1*h, m-alpha*h, m-x2*h, m-beta*h, m-x3*h, m, m+x3*h, ...
+       m+beta*h, m+x2*h, m+alpha*h, m+x1*h, b];
+
+  y = feval (f, x, varargin{:});
+
   fa = y(1); 
   fb = y(13);
-  i2 = (h/6)*(y(1)+y(13)+5*(y(5)+y(9)));
-  i1 = (h/1470)*(77*(y(1)+y(13))+432*(y(3)+y(11))+ ...
-		 625*(y(5)+y(9))+672*y(7));
-  is = h*(.0158271919734802*(y(1)+y(13))+.0942738402188500 ...
-	  *(y(2)+y(12))+.155071987336585*(y(3)+y(11))+ ...
-	  .188821573960182*(y(4)+y(10))+.199773405226859 ...
-	  *(y(5)+y(9))+.224926465333340*(y(6)+y(8))...
-	  +.242611071901408*y(7));
+
+  i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9)));
+
+  i1 = (h/1470)*(77*(y(1)+y(13))
+		 + 432*(y(3)+y(11))
+		 + 625*(y(5)+y(9))
+		 + 672*y(7));
+
+  is = h*(.0158271919734802*(y(1)+y(13))
+	  +.0942738402188500*(y(2)+y(12))
+	  + .155071987336585*(y(3)+y(11))
+	  + .188821573960182*(y(4)+y(10))
+	  + .199773405226859*(y(5)+y(9))
+	  + .224926465333340*(y(6)+y(8))
+	  + .242611071901408*y(7));
+
   s = sign(is); 
-  if (s == 0),
-    s=1;
+
+  if (s == 0)
+    s = 1;
   endif
   erri1 = abs(i1-is);
   erri2 = abs(i2-is);
@@ -105,18 +118,18 @@
     R = erri1/erri2; 
   endif
   if (R > 0 && R < 1)
-    tol=tol/R; 
+    tol = tol/R; 
   endif
   is = s*abs(is)*tol/eps;
   if (is == 0)
     is = b-a;
   endif
-  Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin{:});
+  Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:});
 endfunction
 
 ## ADAPTLOBSTP  Recursive function used by QUADL.
 ##
-##   Q = ADAPTLOBSTP('F',A,B,FA,FB,IS,TRACE) tries to
+##   Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to
 ##   approximate the integral of F(X) from A to B to
 ##   an appropriate relative error. The argument 'F' is
 ##   a string containing the name of f.  The remaining
@@ -124,7 +137,7 @@
 ##
 ##   Walter Gautschi, 08/03/98
 
-function Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin)
+function Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
   h = (b-a)/2; 
   m = (a+b)/2;
   alpha = sqrt(2/3); 
@@ -133,37 +146,36 @@
   ml = m-beta*h; 
   mr = m+beta*h; 
   mrr = m+alpha*h;
-  x = [mll,ml,m,mr,mrr];
-  y = feval(f,x,varargin{:});
+  x = [mll, ml, m, mr, mrr];
+  y = feval(f, x, varargin{:});
   fmll = y(1); 
   fml = y(2); 
   fm = y(3); 
   fmr = y(4); 
   fmrr = y(5);
-  i2 = (h/6)*(fa+fb+5*(fml+fmr));
-  i1 = (h/1470)*(77*(fa+fb)+432*(fmll+fmrr)+625*(fml+fmr) ...
-		 +672*fm);
-  if ((is+(i1-i2) == is) || (mll <= a) || (b <= mrr))
-    if (((m <= a) || (b <= m)) && need_warning())
-      warning(['Interval contains no more machine number. ',...
-               'Required tolerance may not be met.']);
-      need_warning(0);
+  i2 = (h/6)*(fa + fb + 5*(fml+fmr));
+  i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm);
+  if (is+(i1-i2) == is || mll <= a || b <= mrr)
+    if ((m <= a || b <= m) && need_warning ())
+      warning ("quadl: interval contains no more machine number");
+      warning ("quadl: required tolerance may not be met");
+      need_warning (0);
     endif
     Q = i1;
     if (trace)
-      disp([a b-a Q]);
+      disp ([a, b-a, Q]);
     endif
   else
-    Q = adaptlobstp(f,a,mll,fa,fmll,is,trace,varargin{:})+...
-	adaptlobstp(f,mll,ml,fmll,fml,is,trace,varargin{:})+...
-	adaptlobstp(f,ml,m,fml,fm,is,trace,varargin{:})+...
-	adaptlobstp(f,m,mr,fm,fmr,is,trace,varargin{:})+...
-	adaptlobstp(f,mr,mrr,fmr,fmrr,is,trace,varargin{:})+...
-	adaptlobstp(f,mrr,b,fmrr,fb,is,trace,varargin{:});
+    Q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:})
+	 + adaptlobstp (f, mll, ml, fmll, fml, is, trace, varargin{:})
+	 + adaptlobstp (f, ml, m, fml, fm, is, trace, varargin{:})
+	 + adaptlobstp (f, m, mr, fm, fmr, is, trace, varargin{:})
+	 + adaptlobstp (f, mr, mrr, fmr, fmrr, is, trace, varargin{:})
+	 + adaptlobstp (f, mrr, b, fmrr, fb, is, trace, varargin{:}));
   endif
 endfunction
 
-function r = need_warning(v)
+function r = need_warning (v)
   persistent w = [];
   if (nargin == 0)
     r = w;
--- a/scripts/miscellaneous/inputname.m
+++ b/scripts/miscellaneous/inputname.m
@@ -6,13 +6,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} inputname (@var{n})
-## 
-##    Return the text defining n-th input to the function.
-##
+## Return the text defining @var{n}-th input to the function.
 ## @end deftypefn
 
 function s = inputname (n)
-  s = evalin ('caller', sprintf('deblank(argn(%d,:));', n));
+  s = evalin ("caller", sprintf ("deblank (argn(%d,:));", n));
 endfunction
 
 ## Warning: heap big magic in the following tests!!!
--- a/scripts/plot/__plt3__.m
+++ b/scripts/plot/__plt3__.m
@@ -35,30 +35,31 @@
 
 function __plt3__ (x, y, z, fmt)
 
-  if (isvector(x) && isvector(y))
-    if (isvector(z))
-      x = x(:); y=y(:); z=z(:);
-    elseif (length(x) == rows(z) && length(y) == columns(z))
-      error("plot3: [length(x), length(y)] must match size(z)");
+  if (isvector (x) && isvector (y))
+    if (isvector (z))
+      x = x(:);
+      y = y(:);
+      z = z(:);
+    elseif (length (x) == rows (z) && length (y) == columns (z))
+      error ("plot3: [length(x), length(y)] must match size(z)");
     else
-      [x,y] = meshgrid(x,y);
+      [x, y] = meshgrid (x, y);
     endif
   endif
 
-  if (any(size(x) != size(y)) || any(size(x) != size(z)))
-    error("plot3: x, y, and z must have the same shape");
+  if (any (size (x) != size (y)) || any (size (x) != size (z)))
+    error ("plot3: x, y, and z must have the same shape");
   endif
 
   unwind_protect
-    __gnuplot_set__  parametric;
+    __gnuplot_set__ parametric;
     __gnuplot_raw__ ("set nohidden3d;\n");
 
     tmp = [([x; NaN*ones(1,size(x,2))])(:), ...
 	   ([y; NaN*ones(1,size(y,2))])(:), ...
 	   ([z; NaN*ones(1,size(z,2))])(:)];
 
-    cmd =  ["__gnuplot_splot__ tmp ", fmt, ";\n"];
-    eval (cmd);
+    eval (sprintf ("__gnuplot_splot__ tmp %s\n", fmt));
   unwind_protect_cleanup
     __gnuplot_set__ noparametric; 
   end_unwind_protect
--- a/scripts/plot/ndgrid.m
+++ b/scripts/plot/ndgrid.m
@@ -32,43 +32,40 @@
 
 ## Author: Alexander Barth <abarth@marine.usf.edu>
 
-function varargout = ndgrid(varargin)
+function varargout = ndgrid (varargin)
 
   if (nargin == 1)
-    n = max([nargout 2]);  
+    n = max ([nargout, 2]);  
     ## If only one input argument is given, repeat it n-times
     varargin{1:n} = varargin{1};
   elseif (nargin >= nargout)
-    n = max([nargin 2]);  
+    n = max ([nargin, 2]);  
   else
     error ("ndgrid: wrong number of input arguments");
   endif
 
-
   ## Determine the size of the output arguments
   
-  shape = zeros(1,n);
+  shape = zeros (1, n);
 
-  for i=1:n
-    if (~isvector(varargin{i}))
+  for i = 1:n
+    if (! isvector (varargin{i}))
       error ("ndgrid: arguments must be vectors");
     endif
 
-    shape(i) = length(varargin{i});
+    shape(i) = length (varargin{i});
   endfor
 
-  for i=1:n
+  for i = 1:n
     ## size for reshape
-    r = ones(1,n);
+    r = ones (1, n);
     r(i) = shape(i);
 
     ## size for repmat
     s = shape;
     s(i) = 1;
 
-    varargout{i} = repmat(reshape(varargin{i},r) ,s);
+    varargout{i} = repmat (reshape (varargin{i}, r), s);
   endfor
 
-
-
 endfunction
--- a/scripts/plot/plot3.m
+++ b/scripts/plot/plot3.m
@@ -137,7 +137,7 @@
 ## Author: Paul Kienzle
 ##         (modified from __plt__.m)
 
-function plot3(varargin)
+function plot3 (varargin)
 
   hold_state = ishold ();
   
@@ -148,7 +148,7 @@
     z_set = 0;
     
     ## Gather arguments, decode format, and plot lines.
-    for arg = 1:length(varargin)
+    for arg = 1:nargin
       new = varargin{arg};
       
       if (ischar (new))
@@ -161,13 +161,13 @@
 	x_set = 0;
 	y_set = 0;
 	z_set = 0;
-      elseif (!x_set)
+      elseif (! x_set)
 	x = new;
 	x_set = 1;
-      elseif (!y_set)
+      elseif (! y_set)
 	y = new;
 	y_set = 1;
-      elseif (!z_set)
+      elseif (! z_set)
 	z = new;
 	z_set = 1;
       else
--- a/scripts/polynomial/pchip.m
+++ b/scripts/polynomial/pchip.m
@@ -91,16 +91,16 @@
     y = reshape (y, [prod(szy(1:end-1)), szy(end)])';
   endif
 
-  h = diff(x);
-  if all(h<0)
-    x = flipud(x);
-    h = diff(x);
-    y = flipud(y);
-  elseif (any(h <= 0))
+  h = diff (x);
+  if (all (h < 0))
+    x = flipud (x);
+    h = diff (x);
+    y = flipud (y);
+  elseif (any (h <= 0))
     error("pchip: x must be strictly monotonic")
   endif
 
-  if (rows(y) != n)
+  if (rows (y) != n)
     error("pchip: size of x and y must match");
   endif
 
@@ -112,7 +112,7 @@
   dy = diff (y) ./ h;
 
   a = y;
-  b = __pchip_deriv__(x,y);
+  b = __pchip_deriv__ (x, y);
   c = - (b(2:n, :) + 2 * b(1:n - 1, :)) ./ h + 3 * diff (a) ./ h .^ 2;
   d = (b(1:n - 1, :) + b(2:n, :)) ./ h.^2 - 2 * diff (a) ./ h.^3;
 
--- a/scripts/polynomial/spline.m
+++ b/scripts/polynomial/spline.m
@@ -112,88 +112,93 @@
 
   b = c = zeros (size (a));
   h = diff (x);
-  idx = ones (columns (a),1);
+  idx = ones (columns (a), 1);
 
   if (complete)
 
     if (n == 3)
       dg = 1.5 * h(1) - 0.5 * h(2);
-      c(2:n - 1,:) = 1/dg(1);
+      c(2:n-1,:) = 1/dg(1);
     else
-      dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
+      dg = 2 * (h(1:n-2) .+ h(2:n-1));
       dg(1) = dg(1) - 0.5 * h(1);
-      dg(n - 2) = dg(n-2) - 0.5 * h(n - 1);
+      dg(n-2) = dg(n-2) - 0.5 * h(n-1);
 
-      e = h(2:n - 2);
+      e = h(2:n-2);
 
       size(a)
       size(h)
       n
 
-      g = 3 * diff (a(2:n,:)) ./ h(2:n - 1,idx)\
-        - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2,idx);
-      g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) \
+      g = 3 * diff (a(2:n,:)) ./ h(2:n-1,idx) ...
+          - 3 * diff (a(1:n-1,:)) ./ h(1:n-2,idx);
+      g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) ...
           - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs);
-      g(n - 2,:) = 3 / 2 * (3 * (a(n,:) - a(n - 1,:)) / h(n - 1) - dfe)\
-          - 3 * (a(n - 1,:) - a(n - 2,:)) / h(n - 2);
+      g(n-2,:) = 3 / 2 * (3 * (a(n,:) - a(n-1,:)) / h(n-1) - dfe) ...
+          - 3 * (a(n-1,:) - a(n-2,:)) / h(n-2);
 
-      c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g;
+      c(2:n-1,:) = spdiags ([[e(:); 0], dg, [0; e(:)]],
+			      [-1, 0, 1], n-2, n-2) \ g;
     endif
 
     c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs
-	      - c(2,:) * h(1)) / (2 * h(1)); 
-    c(n,:) = - (3 / h(n - 1) * (a(n,:) - a(n - 1,:)) - 3 * dfe
-		+ c(n - 1,:) * h(n - 1)) / (2 * h(n - 1));
-    b(1:n - 1,:) = diff (a) ./ h(1:n - 1, idx)\
-      - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
-    d = diff (c) ./ (3 * h(1:n - 1, idx));
+	      - c(2,:) * h(1)) / (2 * h(1));
+    c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe
+		+ c(n-1,:) * h(n-1)) / (2 * h(n-1));
+    b(1:n-1,:) = diff (a) ./ h(1:n-1, idx) ...
+	- h(1:n-1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
+    d = diff (c) ./ (3 * h(1:n-1, idx));
 
   else
 
-    g = zeros(n - 2,columns(a));
-    g(1,:) = 3 / (h(1) + h(2)) * (a(3,:) - a(2,:)\
-				  - h(2) / h(1) * (a(2,:) - a(1,:)));
-    g(n - 2,:) = 3 / (h(n - 1) + h(n - 2)) *\
-	(h(n - 2) / h(n - 1) * (a(n,:) - a(n - 1,:)) -\
-	 (a(n - 1,:) - a(n - 2,:)));
+    g = zeros (n-2, columns (a));
+    g(1,:) = 3 / (h(1) + h(2)) ...
+	* (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:)));
+    g(n-2,:) = 3 / (h(n-1) + h(n-2)) ...
+	* (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:)));
 
     if (n > 4)
 
-      g(2:n - 3,:) = 3 * diff (a(3:n - 1,:)) ./ h(3:n - 2,idx)\
-          - 3 * diff (a(2:n - 2,:)) ./ h(2:n - 3,idx);
+      g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ...
+          - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx);
 
-      dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
+      dg = 2 * (h(1:n-2) .+ h(2:n-1));
       dg(1) = dg(1) - h(1);
-      dg(n - 2) = dg(n-2) - h(n - 1);
+      dg(n-2) = dg(n-2) - h(n-1);
 
-      ldg = udg = h(2:n - 2);
+      ldg = udg = h(2:n-2);
       udg(1) = udg(1) - h(1);
-      ldg(n - 3) = ldg(n-3) - h(n - 1);
-      c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
+      ldg(n - 3) = ldg(n-3) - h(n-1);
+      c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]],
+			      [-1, 0, 1], n-2, n-2) \ g;
 
     elseif (n == 4)
 
       dg = [h(1) + 2 * h(2), 2 * h(2) + h(3)];
       ldg = h(2) - h(3);
       udg = h(2) - h(1);
-      c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
+      c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]],
+			      [-1, 0, 1], n-2, n-2) \ g;
       
     else # n == 3
 	    
-      dg= [h(1) + 2 * h(2)];
-      c(2:n - 1,:) = g/dg(1);
+      dg = h(1) + 2 * h(2);
+      c(2:n-1,:) = g/dg(1);
 
     endif
 
     c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
-    c(n,:) = c(n - 1,:) + h(n - 1) / h(n - 2) * (c(n - 1,:) - c(n - 2,:));
-    b = diff (a) ./ h(1:n - 1, idx)\
-	- h(1:n - 1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
-    d = diff (c) ./ (3 * h(1:n - 1, idx));
+    c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:));
+    b = diff (a) ./ h(1:n-1, idx) ...
+	- h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
+    d = diff (c) ./ (3 * h(1:n-1, idx));
 
   endif
 
-  d = d(1:n-1,:); c=c(1:n-1,:); b=b(1:n-1,:); a=a(1:n-1,:);
+  d = d(1:n-1,:);
+  c = c(1:n-1,:);
+  b = b(1:n-1,:);
+  a = a(1:n-1,:);
   coeffs = [d(:), c(:), b(:), a(:)];
   ret = mkpp (x, coeffs, szy(1:end-1));
 
--- 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
--- a/scripts/sparse/pcr.m
+++ b/scripts/sparse/pcr.m
@@ -18,20 +18,20 @@
 ## 02110-1301, USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x} =} pcr (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M}, @var{x0}, @dots{})
+## @deftypefn {Function File} {@var{x} =} pcr (@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}] =} pcr (@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 Residuals 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 non-singular; if @code{pcr}
-## finds @var{A} to be numerically singular, you will get a warning
+## of a function which computes @code{@var{a} * @var{x}}. In principle
+## @var{a} should be symmetric and non-singular; if @code{pcr}
+## finds @var{a} to be numerically singular, 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{pcr} @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{pcr}. 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
@@ -93,7 +93,7 @@
 ## @var{resvec} describes the convergence history of the method,
 ## so that @code{@var{resvec} (i)} contains the Euclidean norms of the 
 ## residualafter the (@var{i}-1)-th iteration, @code{@var{i} =
-## 1,2,...@var{iter}+1}.
+## 1,2, @dots{}, @var{iter}+1}.
 ## @end itemize
 ## 
 ## Let us consider a trivial problem with a diagonal matrix (we exploit the
@@ -114,7 +114,7 @@
 ## @end example
 ## 
 ## @sc{Example 2:} @code{pcr} with a function which computes
-## @code{@var{A} * @var{x}}.
+## @code{@var{a} * @var{x}}.
 ##
 ## @example
 ## @group
@@ -128,7 +128,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
@@ -165,22 +165,14 @@
 ## @seealso{sparse, pcg}
 ## @end deftypefn
 
-## REVISION HISTORY
-##
-## 2004-08-14, Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
-## 
-## 	Added 4 demos and 4 tests
-##  
-## 2004-08-01, Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
-## 
-## 	First version of pcr code
+## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
 
-function [x, flag, relres, iter, resvec] = pcr(A,b,tol,maxit,M,x0,varargin)
+function [x, flag, relres, iter, resvec] = pcr (A, b, tol, maxit, M, x0, varargin)
 
   breakdown = false;
 
-  if ((nargin < 6) || isempty(x0))
-    x = zeros(size(b));
+  if (nargin < 6 || isempty (x0))
+    x = zeros (size (b));
   else
     x = x0;
   endif
@@ -189,84 +181,86 @@
     M = [];
   endif
 
-  if ((nargin < 4) || isempty(maxit))
+  if (nargin < 4 || isempty (maxit))
     maxit = 20;
   endif
 
-  maxit = maxit + 2;
+  maxit += 2;
 
-  if ((nargin < 3) || isempty(tol))
+  if (nargin < 3 || isempty (tol))
     tol = 1e-6;
   endif
 
   if (nargin < 2)
-    print_usage();
+    print_usage ();
   endif
 
   ##  init
-  if (isnumeric(A))		# is A a matrix?
-    r = b-A*x;
+  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
 
-  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
-      p = M\r;
+      p = M \ r;
     endif
   else				# then M should be a function!
-    p = feval(M,r,varargin{:});
+    p = feval (M, r, varargin{:});
   endif
 
   iter = 2;
 
   b_bot_old = 1;
-  q_old = p_old = s_old = zeros(size(x));
+  q_old = p_old = s_old = zeros (size (x));
 
-  if (isnumeric(A))		# is A a matrix?
-    q = A*p;
+  if (isnumeric (A))		# is A a matrix?
+    q = A * p;
   else				# then A should be a function!
-    q = feval(A,p,varargin{:});
+    q = feval (A, p, varargin{:});
   endif
 	
-  resvec(1) = abs(norm(r)); 
+  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
+  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;
+	s = M \ q;
       endif
     else			# then M should be a function!
-      s = feval(M,q,varargin{:});
+      s = feval (M, q, varargin{:});
     endif
-    b_top = r'*s;
-    b_bot = q'*s;
+    b_top = r' * s;
+    b_bot = q' * s;
 	
     if (b_bot == 0.0)
       breakdown = true;
       break;
     endif
-    lambda = b_top/b_bot;
+    lambda = b_top / b_bot;
 	
-    x = x + lambda*p;
-    r = r - lambda*q;
+    x += lambda*p;
+    r -= lambda*q;
 	
     if (isnumeric(A))		# is A a matrix?
       t = A*s;
     else			# then A should be a function!
-      t = feval(A,s,varargin{:});
+      t = feval (A, s, varargin{:});
     endif
 	
-    alpha0 = (t'*s)/b_bot;
-    alpha1 = (t'*s_old)/b_bot_old;
+    alpha0 = (t'*s) / b_bot;
+    alpha1 = (t'*s_old) / b_bot_old;
 	
-    p_temp = p; q_temp = q;
+    p_temp = p;
+    q_temp = q;
+
     p = s - alpha0*p - alpha1*p_old;
     q = t - alpha0*q - alpha1*q_old;
 	
@@ -275,32 +269,30 @@
     q_old = q_temp;
     b_bot_old = b_bot;
 	
-	
-    resvec(iter) = abs(norm(r));
-    iter = iter + 1;
+    resvec(iter) = abs (norm (r));
+    iter++;
   endwhile
 
   flag = 0;
-  relres = resvec(iter-1)./resvec(1);
-  iter = iter - 2;
-  if (iter >= (maxit-2))
+  relres = resvec(iter-1) ./ resvec(1);
+  iter -= 2;
+  if (iter >= maxit-2)
     flag = 1;
     if (nargout < 2)
-      warning("PCR: maximum number of iterations (%d) reached\n", iter);
-      warning("The initial residual norm was reduced %g times.\n", 1.0/relres);
+      warning ("PCR: 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) && (~breakdown))
-      fprintf(stderr, "PCR: converged in %d iterations. \n", iter);
-      fprintf(stderr, "The initial residual norm was reduced %g times.\n",...
-	      1.0/relres);
-    endif
+  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);
   endif
+
   if (breakdown)
     flag = 3;
     if (nargout < 2)
-      warning("PCR: breakdown occured.\n");
-      warning("System matrix singular or preconditioner indefinite?\n");
+      warning ("PCR: breakdown occured.\n");
+      warning ("System matrix singular or preconditioner indefinite?\n");
     endif
   endif
 
--- a/scripts/strings/mat2str.m
+++ b/scripts/strings/mat2str.m
@@ -49,62 +49,67 @@
 ## @seealso{sprintf, int2str}
 ## @end deftypefn
 
-function s = mat2str(x,n,cls)
+function s = mat2str (x, n, cls)
 
-  if (nargin < 2 || isempty(n))
-    n = 17;		   # default precision
+  if (nargin < 2 || isempty (n))
+    ## Default precision
+    n = 17;
   endif
 
   if (nargin < 3)
-    if (ischar(n))
+    if (ischar (n))
       cls = n;
       n = 17;
     else
-      cls = '';
+      cls = "";
     endif
   endif
 
-  if (nargin < 1 || nargin > 3 || ischar(x) || isstruct(x) ||
-     ischar(n) || isstruct(n) || isstruct(cls))
-    usage ("mat2str");
+  if (nargin < 1 || nargin > 3 || ! isnumeric (x))
+    print_usage ();
   endif
   
   if (ndims (x) > 2)
-    error ("mat2str: x must be two dimensional");
+    error ("mat2str: X must be two dimensional");
   endif
 
-  if (!(COMPLEX = is_complex(x)))
-    FMT = sprintf("%%.%dg", n(1));
+  x_is_complex = is_complex (x);
+
+  if (! x_is_complex)
+    fmt = sprintf ("%%.%dg", n(1));
   else
-    if (length(n) == 1 )
+    if (length (n) == 1 )
       n = [n, n];
     endif
-    FMT = sprintf("%%.%dg%%+.%dgi", n(1), n(2));
+    fmt = sprintf ("%%.%dg%%+.%dgi", n(1), n(2));
   endif
 
-  [nr, nc] = size(x);
+  nel = numel (x);
 
-  if (nr*nc == 0)         # empty .. only print brackets
+  if (nel == 0)
+    ## Empty, only print brackets
     s = "[]";
-  elseif (nr*nc == 1)	# scalar x .. don't print brackets
-    if (!COMPLEX)
-      s = sprintf( FMT, x );
+  elseif (nel == 1)
+    ## Scalar X, don't print brackets
+    if (! x_is_complex)
+      s = sprintf (fmt, x);
     else
-      s = sprintf( FMT, real(x), imag(x) );
+      s = sprintf (fmt, real (x), imag (x));
     endif
-  else			# non-scalar x .. print brackets
-    FMT = [FMT, ','];
-    if (!COMPLEX)
-      s = sprintf( FMT, x.' );
+  else
+    ## Non-scalar X, print brackets
+    fmt = [fmt, ","];
+    if (! x_is_complex)
+      s = sprintf (fmt, x.');
     else
       x = x.';
-      s = sprintf( FMT, [ real(x(:))'; imag(x(:))' ] );
+      s = sprintf (fmt, [real(x(:))'; imag(x(:))']);
     endif
 
     s = ["[", s];
-    s (length(s)) = "]";
-    IND = find(s == ",");
-    s (IND(nc:nc:length(IND)) ) = ";";
+    s(end) = "]";
+    ind = find (s == ",");
+    s(ind(nc:nc:end)) = ";";
   endif
 
   if (strcmp ("class", cls))
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,7 @@
+2006-06-01  John W. Eaton  <jwe@octave.org>
+
+	* Makefile.in (DLD_XSRC): Add __pchip_deriv__.cc to the list.
+
 2006-06-01  David Bateman  <dbateman@free.fr>
 
 	* DLD-FUNCTIONS/__pchip_deriv__.cc: New file.
--- a/src/DLD-FUNCTIONS/__pchip_deriv__.cc
+++ b/src/DLD-FUNCTIONS/__pchip_deriv__.cc
@@ -32,61 +32,60 @@
 #include "utils.h"
 #include "f77-fcn.h"
 
-extern "C" {
+extern "C"
+{
   extern int F77_FUNC (dpchim, DPCHIM)
-    (const int &n, double *x, double *f, 
-	  double *d, const int &incfd, int *ierr);
+    (const int &n, double *x, double *f, double *d, const int &incfd,
+     int *ierr);
 }
 
 DEFUN_DLD(__pchip_deriv__, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {} __pchip_deriv__ (@var{x}, @var{y})\n\
 Wrapper for SLATEC/PCHIP function DPCHIM to calculate the derivates for\n\
-piecewise polynomials. You should be using @code{pchip} function instead.\n\
+piecewise polynomials.  You should be using @code{pchip} function instead.\n\
 @end deftypefn")
 {
   octave_value retval;
-  const int nargin = args.length();
+  const int nargin = args.length ();
 
   if (nargin == 2)
     {
-      ColumnVector xvec(args(0).vector_value());
-      Matrix ymat(args(1).matrix_value());
-      int nx = xvec.length();
-      int nyr = ymat.rows();
-      int nyc = ymat.columns();
+      ColumnVector xvec (args(0).vector_value ());
+      Matrix ymat (args(1).matrix_value ());
+
+      int nx = xvec.length ();
+      int nyr = ymat.rows ();
+      int nyc = ymat.columns ();
 
       if (nx != nyr)
         {
-          error("number of rows for x and y must match");
+          error ("number of rows for x and y must match");
           return retval;
         }
 
-      ColumnVector dvec(nx),yvec(nx);
-      Matrix dmat(nyr,nyc);
+      ColumnVector dvec (nx), yvec (nx);
+      Matrix dmat (nyr, nyc);
 
       int ierr;
       const int incfd = 1;
       for (int c = 0; c < nyc; c++)
         {
-          for (int r=0; r<nx; r++)
-            {
-              yvec(r) = ymat(r,c);
-            }
+          for (int r = 0; r < nx; r++)
+	    yvec(r) = ymat(r,c);
 
-          F77_FUNC (dpchim, DPCHIM) (nx, xvec.fortran_vec(), 
-				     yvec.fortran_vec(), 
-				     dvec.fortran_vec(), incfd, &ierr);
+          F77_FUNC (dpchim, DPCHIM) (nx, xvec.fortran_vec (), 
+				     yvec.fortran_vec (), 
+				     dvec.fortran_vec (), incfd, &ierr);
 
 	  if (ierr < 0)
             {
-	      error ("DPCHIM error: %i\n",ierr);
+	      error ("DPCHIM error: %i\n", ierr);
               return retval;
             }
+
           for (int r=0; r<nx; r++)
-            {
-              dmat(r,c) = dvec(r);
-            }
+	    dmat(r,c) = dvec(r);
         }
 
       retval = dmat;
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -53,7 +53,7 @@
 	quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \
 	spchol.cc spdet.cc spkron.cc splu.cc spparms.cc spqr.cc \
 	sqrtm.cc svd.cc syl.cc time.cc \
-	__gnuplot_raw__.l __glpk__.cc __qp__.cc
+	__gnuplot_raw__.l __glpk__.cc __pchip_deriv__.cc __qp__.cc
 
 DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))