Mercurial > hg > octave-lojdl
changeset 5838:376e02b2ce70
[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))