Mercurial > forge
changeset 2312:21b6809d768c octave-forge
Remove files just added to octave 2.9.5+
| author | adb014 |
|---|---|
| date | Fri, 02 Jun 2006 06:36:59 +0000 |
| parents | 3833cc92119e |
| children | b1d9540811e7 |
| files | main/general/bicubic.m main/general/gradient.m main/general/interp1.m main/general/interp2.m main/general/interpft.m main/general/polyarea.m main/general/quadl.m main/miscellaneous/inputname.m main/plot/__plt3__.m main/plot/plot3.m main/sparse/pcg.m main/sparse/pcr.m main/splines/Makefile main/splines/dpchim.f main/splines/dpchst.f main/splines/pchip.m main/splines/pchip_deriv.cc main/strings/mat2str.m |
| diffstat | 18 files changed, 0 insertions(+), 2808 deletions(-) [+] |
line wrap: on
line diff
--- a/main/general/bicubic.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,194 +0,0 @@ -## Copyright (C) 2005 Hoxide Ma -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; if not, write to the Free Software -## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -## -*- 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 -## 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} -## @end deftypefn -## @seealso{interp2} - -## Bicubic interpolation method. -## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn> - -function F = bicubic(X, Y, Z, XI, YI) - global spline_alpha; - if nargin<1 || nargin>6 - usage ("bicubic (X, Y, Z, XI, YI)") - endif - - if exist("spline_alpha","var") && prod(size(spline_alpha))==1 - a = spline_alpha - else - a = 0.5; - endif - - if nargin <=2 # bicubic(Z) or bicubic(Z, 2) - if nargin==1 - n = 1; - else - n = Y; - endif - 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); - elseif nargin == 3 - if !isvector (X) || !isvector (Y) - error ("XI and YI must be vector"); - endif - s = Y; - t = Z; - Z = X; - [rz, cz] = size(Z); - elseif nargin == 5 - [rz, cz] = size (Z) ; - if isvector (X) && isvector (Y) - if(rz != length (Y) || cz != length (X)) - error ("length of X and Y must match the size of Z"); - endif - elseif size(X) == size(Y) && size(X) == size(Z) - X = X(1,:); - Y = Y(:,1); - else - error ("X, Y and Z must be martrices of same size"); - endif - - ## mark values outside the lookup table - xfirst_ind = find (XI < X(1)); - xlast_ind = find (XI > X(cz)); - yfirst_ind = find (YI < Y(1)); - ylast_ind = find (YI > Y(rz)); - ## set value outside the table preliminary to min max index - XI(xfirst_ind) = X(1); - XI(xlast_ind) = X(cz); - YI(yfirst_ind) = Y(1); - YI(ylast_ind) = Y(rz); - - - 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)]); - - 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); - yidx = o([find( i> rz)]); - - ## set s and t used follow codes - s = xidx + ((XI .- X(xidx))./(X(xidx+1) .- X(xidx))); - t = yidx + ((YI - Y(yidx))./(Y(yidx+1) - Y(yidx))); - else - usage ("bicubic (X, Y, Z, XI, YI)"); - endif - - if rz < 3 || cz < 3 - error ("Z at least a 3 by 3 matrices"); - endif - - 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(d) = rz-1; - t(d) = 1.0; - d = []; - - 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,:); - p(:,1) = (6*(1-a))*p(:,2) -3*p(:,3) + (6*a-2)*p(:,4); - p(:,cz+2) = (6*(1-a))*p(:,cz+1) -3*p(:,cz) + (6*a-2)*p(:,cz-1); - - ## calculte the C1(t) C2(t) C3(t) C4(t) and C1(s) C2(s) C3(s) C4(s) - t2= t.*t; - t3= t2.*t; - - ct0= -a .* t3 + (2 * a) .* t2 - a .* t ; # -a G0 - ct1 = (2-a) .* t3 + (-3+a) .* t2 + 1 ; # F0 - a G1 - ct2 = (a-2) .* t3 + (-2 *a + 3) .* t2 + a .* t ; # F1 + a G0 - ct3 = a .* t3 - a .* t2; # a G1 - t = [];t2=[]; t3=[]; - - s2= s.*s; - s3= s2.*s; - - cs0= -a .* s3 + (2 * a) .* s2 - a .*s ; # -a G0 - cs1 = (2-a) .* s3 + (-3 + a) .* s2 + 1 ; # F0 - a G1 - cs2 = (a-2) .* s3 + (-2 *a + 3) .* s2 + a .*s ; # F1 + a G0 - cs3 = a .* s3 - a .* s2; # a G1 - s=[] ; s2 = []; s3 = []; - - cs0 = cs0([1,1,1,1],:); - cs1 = cs1([1,1,1,1],:); - cs2 = cs2([1,1,1,1],:); - cs3 = cs3([1,1,1,1],:); - - 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); - endfor - - ## set points outside the table to NaN - if (! (isempty (xfirst_ind) && isempty (xlast_ind))) - F(:, [xfirst_ind, xlast_ind]) = NaN; - endif - if (! (isempty (yfirst_ind) && isempty (ylast_ind))) - F([yfirst_ind; ylast_ind], :) = NaN; - endif - -endfunction - -%!demo -%! A=[13,-1,12;5,4,3;1,6,2]; -%! x=[0,1,4]+10; y=[-10,-9,-8]; -%! xi=linspace(min(x),max(x),17); -%! yi=linspace(min(y),max(y),26); -%! mesh(xi,yi,bicubic(x,y,A,xi,yi)); -%! [x,y] = meshgrid(x,y); -%! __gnuplot_raw__ ("set nohidden3d;\n") -%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
--- a/main/general/gradient.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,115 +0,0 @@ -## Copyright (C) 2000 Kai Habel -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; if not, write to the Free Software -## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{X} = } gradient (@var{M}) -## @deftypefnx {Function File} {[@var{X},@var{Y}] = }gradient (@var{M}) -## @deftypefnx {Function File} {[...] = }gradient (...,@var{dx},@var{dy}) -## -## @var{X} = gradient (@var{M}) calculates the one dimensional -## gradient if M is a vector. Is M a Matrix the gradient is calculated -## for each row. -## [@var{X},@var{Y}] = gradient (@var{M}) calculates the one dimensinal -## gradient for each direction if @var{M} is a matrix. -## Spacing values between two points can be provided by the -## @var{dx}, @var{dy} parameters, otherwise the spacing is set to 1. -## A scalar value specifies an equidistant spacing. -## A vector value can be used to specify a variable spacing. The length -## must match their respective dimension of @var{M}. -## -## At boundary points a linear extrapolation is applied. Interior points -## are calculated with the first approximation of the numerical gradient -## @example -## y'(i) = 1/(x(i+1)-x(i-1)) *(y(i-1)-y(i+1)). -## @end example -## -## @end deftypefn - -## Author: Kai Habel <kai.habel@gmx.de> - -function [varargout] = gradient (M, dx, dy) - - if ((nargin < 1) || (nargin > 3)) - usage ("gradient(M,dx,dy)"); - elseif (nargin == 1) - dx = dy = 1; - elseif (nargin == 2) - dy = 1; - endif - - if (isvector (M)) - ## make a row vector - M = M(:)'; - endif - - if !(is_matrix (M)) - error ("first argument must be a vector or matrix"); - else - if (is_scalar (dx)) - dx = dx * ones (1, columns (M) - 1); - else - if (length (dx) != columns (M)) - error ("columns of M must match length of dx"); - else - dx = diff (dx); - dx = dx(:)'; - endif - endif - - if (is_scalar (dy)) - dy = dy * ones (rows (M) - 1, 1); - else - if (length (dy) != rows (M)) - error ("rows of M must match length of dy") - else - dy = diff (dy); - endif - endif - endif - [mr, mc] = size (M); - X = zeros (size (M)); - - if (mc > 1) - ## left and right boundary - X(:, 1) = diff (M(:, 1:2)')' / dx(1); - X(:,mc) = diff (M(:, mc - 1:mc)')' / dx(mc - 1); - endif - - if (mc > 2) - ## interior points - X(:, 2:mc - 1) = (M(:, 3:mc) .- M(:,1:mc - 2))\ - ./ kron (dx(1:mc - 2) .+ dx(2:mc - 1), ones (mr, 1)); - # ./ (ones (mr, 1) * (dx(1:mc - 2) .+ dx(2:mc - 1))); - endif - vr_val_cnt = 1; varargout{vr_val_cnt++} = X; - - if (nargout == 2) - Y = zeros (size (M)); - if (mr > 1) - ## top and bottom boundary - Y(1, :) = diff (M(1:2, :)) / dy(1); - Y(mr, :) = diff (M(mr - 1:mr, :)) / dy(mr - 1); - endif - - if (mr > 2) - ## interior points - Y(2:mr-1, :) = (M(3:mr, :) .- M(1:mr - 2, :))\ - ./kron (dy(1:mr - 2) .+ dy(2:mr - 1), ones(1, mc)); - endif - varargout{vr_val_cnt++} = Y; - - endif -endfunction
--- a/main/general/interp1.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,301 +0,0 @@ -## Copyright (C) 2000 Paul Kienzle -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; if not, write to the Free Software -## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -## usage: yi = interp1(x, y, xi [, 'method' [, 'extrap']]) -## -## Interpolate the function y=f(x) at the points xi. The sample -## points x must be strictly monotonic. If y is a matrix with -## length(x) rows, yi will be a matrix of size rows(xi) by columns(y), -## or its transpose if xi is a row vector. -## -## Method is one of: -## 'nearest': return nearest neighbour. -## 'linear': linear interpolation from nearest neighbours -## 'pchip': piece-wise cubic hermite interpolating polynomial -## 'cubic': cubic interpolation from four nearest neighbours -## 'spline': cubic spline interpolation--smooth first and second -## derivatives throughout the curve -## ['*' method]: same as method, but assumes x is uniformly spaced -## only uses x(1) and x(2); usually faster, never slower -## -## Method defaults to 'linear'. -## -## If extrap is the string 'extrap', then extrapolate values beyond -## the endpoints. If extrap is a number, replace values beyond the -## endpoints with that number. If extrap is missing, assume NaN. -## -## Example: -## xf=[0:0.05:10]; yf = sin(2*pi*xf/5); -## xp=[0:10]; yp = sin(2*pi*xp/5); -## lin=interp1(xp,yp,xf); -## spl=interp1(xp,yp,xf,'spline'); -## cub=interp1(xp,yp,xf,'cubic'); -## near=interp1(xp,yp,xf,'nearest'); -## plot(xf,yf,';original;',xf,lin,';linear;',xf,spl,';spline;',... -## xf,cub,';cubic;',xf,near,';nearest;',xp,yp,'*;;'); -## -## See also: interp - -## 2000-03-25 Paul Kienzle -## added 'nearest' as suggested by Kai Habel -## 2000-07-17 Paul Kienzle -## added '*' methods and matrix y -## check for proper table lengths -## 2002-01-23 Paul Kienzle -## fixed extrapolation - -function yi = interp1(x, y, xi, method, extrap) - - if nargin<3 || nargin>5 - usage("yi = interp1(x, y, xi [, 'method' [, 'extrap']])"); - endif - - if nargin < 4, - method = 'linear'; - else - method = tolower(method); - endif - - if nargin < 5 - extrap = NaN; - endif - - ## reshape matrices for convenience - x = x(:); - if size(y,1)==1, y=y(:); endif - transposed = (size(xi,1)==1); - xi = xi(:); - - ## determine sizes - nx = size(x,1); - [ny, nc] = size(y); - if (nx < 2 || ny < 2) - error ("interp1: table too short"); - endif - - ## determine which values are out of range and set them to extrap, - ## unless extrap=='extrap' in which case, extrapolate them like we - ## should be doing in the first place. - minx = x(1); - if (method(1) == '*') - dx = x(2) - x(1); - maxx = minx + (ny-1)*dx; - else - maxx = x(nx); - endif - 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 transposed, yi = yi.'; endif - return; - endif - xi = xi(range); - endif - - if strcmp(method, 'nearest') - idx = lookup(0.5*(x(1:nx-1)+x(2:nx)), xi)+1; - yi(range,:) = y(idx,:); - - elseif strcmp(method, '*nearest') - idx = floor((xi-minx)/dx+1.5); - yi(range,:) = y(idx,:); - - elseif strcmp(method, 'linear') - ## find the interval containing the test point - 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 - dy = y(2:ny,:) - y(1:ny-1,:); - dx = x(2:nx) - x(1:nx-1); - s = (xi - x(idx))./dx(idx); - yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); - - elseif strcmp(method, '*linear') - ## find the interval containing the test point - t = (xi - minx)/dx + 1; - idx = floor(t); - - ## use the endpoints of the interval to define a line - dy = [y(2:ny,:) - y(1:ny-1,:); zeros(1,nc)]; - s = t - idx; - yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); - - elseif strcmp(method, 'pchip') || strcmp(method, '*pchip') - if (nx == 2) x = linspace(minx, maxx, ny); endif - yi(range,:) = pchip(x, y, xi); - - elseif strcmp(method, 'cubic') - if (nx < 4 || ny < 4) - error ("interp1: table too short"); - endif - 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); - 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); - yi(range,:) = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ... - + c(idx,:)).*xi(:,J) + d(idx,:); - - elseif strcmp(method, '*cubic') - if (nx < 4 || ny < 4) - error ("interp1: table too short"); - endif - - ## 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); - t = t - idx; - t2 = t.*t; - tp = 1 - 0.5*t; - a = (1 - t2).*tp; - b = (t2 + t).*tp; - c = (t2 - t).*tp/3; - d = (t2 - 1).*t/6; - 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') - if (nx == 2) x = linspace(minx, maxx, ny); endif - yi(range,:) = spline(x, y, xi); - - else - error(["interp1 doesn't understand method '", method, "'"]); - endif - if transposed, yi=yi.'; endif - -endfunction - -%!demo -%! xf=0:0.05:10; yf = sin(2*pi*xf/5); -%! xp=0:10; yp = sin(2*pi*xp/5); -%! lin=interp1(xp,yp,xf,'linear'); -%! spl=interp1(xp,yp,xf,'spline'); -%! cub=interp1(xp,yp,xf,'pchip'); -%! near=interp1(xp,yp,xf,'nearest'); -%! plot(xf,yf,';original;',xf,near,';nearest;',xf,lin,';linear;',... -%! xf,cub,';pchip;',xf,spl,';spline;',xp,yp,'*;;'); -%! %-------------------------------------------------------- -%! % confirm that interpolated function matches the original - -%!demo -%! xf=0:0.05:10; yf = sin(2*pi*xf/5); -%! xp=0:10; yp = sin(2*pi*xp/5); -%! lin=interp1(xp,yp,xf,'*linear'); -%! spl=interp1(xp,yp,xf,'*spline'); -%! cub=interp1(xp,yp,xf,'*cubic'); -%! near=interp1(xp,yp,xf,'*nearest'); -%! plot(xf,yf,';*original;',xf,near,';*nearest;',xf,lin,';*linear;',... -%! xf,cub,';*cubic;',xf,spl,';*spline;',xp,yp,'*;;'); -%! %-------------------------------------------------------- -%! % confirm that interpolated function matches the original - -%!shared xp, yp, xi, style -%! xp=0:5; yp = sin(2*pi*xp/5); -%! xi = sort([-1, max(xp)*rand(1,6), max(xp)+1]); - -%!test style = 'nearest'; -%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1]), [NaN, NaN]); -%!assert (interp1(xp,yp,xp,style), yp, 100*eps); -%!assert (interp1(xp,yp,xp',style), yp', 100*eps); -%!assert (interp1(xp',yp',xp',style), yp', 100*eps); -%!assert (interp1(xp',yp',xp,style), yp, 100*eps); -%!assert (isempty(interp1(xp',yp',[],style))); -%!assert (isempty(interp1(xp,yp,[],style))); -%!assert (interp1(xp,[yp',yp'],xi(:),style),... -%! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]); -%!assert (interp1(xp,[yp',yp'],xi,style), -%! interp1(xp,[yp',yp'],xi,['*',style])); - -%!test style = 'linear'; -%!assert (interp1(xp, yp, [-1, max(xp)+1]), [NaN, NaN]); -%!assert (interp1(xp,yp,xp,style), yp, 100*eps); -%!assert (interp1(xp,yp,xp',style), yp', 100*eps); -%!assert (interp1(xp',yp',xp',style), yp', 100*eps); -%!assert (interp1(xp',yp',xp,style), yp, 100*eps); -%!assert (isempty(interp1(xp',yp',[],style))); -%!assert (isempty(interp1(xp,yp,[],style))); -%!assert (interp1(xp,[yp',yp'],xi(:),style),... -%! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]); -%!assert (interp1(xp,[yp',yp'],xi,style), -%! interp1(xp,[yp',yp'],xi,['*',style]),100*eps); - -%!test style = 'cubic'; -%!assert (interp1(xp, yp, [-1, max(xp)+1]), [NaN, NaN]); -%!assert (interp1(xp,yp,xp,style), yp, 100*eps); -%!assert (interp1(xp,yp,xp',style), yp', 100*eps); -%!assert (interp1(xp',yp',xp',style), yp', 100*eps); -%!assert (interp1(xp',yp',xp,style), yp, 100*eps); -%!assert (isempty(interp1(xp',yp',[],style))); -%!assert (isempty(interp1(xp,yp,[],style))); -%!assert (interp1(xp,[yp',yp'],xi(:),style),... -%! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]); -%!assert (interp1(xp,[yp',yp'],xi,style), -%! interp1(xp,[yp',yp'],xi,['*',style]),1000*eps); - -%!test style = 'spline'; -%!assert (interp1(xp, yp, [-1, max(xp) + 1]), [NaN, NaN]); -%!assert (interp1(xp,yp,xp,style), yp, 100*eps); -%!assert (interp1(xp,yp,xp',style), yp', 100*eps); -%!assert (interp1(xp',yp',xp',style), yp', 100*eps); -%!assert (interp1(xp',yp',xp,style), yp, 100*eps); -%!assert (isempty(interp1(xp',yp',[],style))); -%!assert (isempty(interp1(xp,yp,[],style))); -%!assert (interp1(xp,[yp',yp'],xi(:),style),... -%! [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]); -%!assert (interp1(xp,[yp',yp'],xi,style), -%! interp1(xp,[yp',yp'],xi,['*',style]),10*eps); - -%!# test linear extrapolation -%!assert (interp1([1:5],[3:2:11],[0,6],'linear','extrap'), [1, 13], eps); -%!assert (interp1(xp, yp, [-1, max(xp)+1],'linear',5), [5, 5]); - -%!error interp1 -%!error interp1(1:2,1:2,1,'bogus') - -%!error interp1(1,1,1, 'nearest'); -%!assert (interp1(1:2,1:2,1.4,'nearest'),1); -%!error interp1(1,1,1, 'linear'); -%!assert (interp1(1:2,1:2,1.4,'linear'),1.4); -%!error interp1(1:3,1:3,1, 'cubic'); -%!assert (interp1(1:4,1:4,1.4,'cubic'),1.4); -%!error interp1(1:2,1:2,1, 'spline'); -%!assert (interp1(1:3,1:3,1.4,'spline'),1.4); - -%!error interp1(1,1,1, '*nearest'); -%!assert (interp1(1:2:4,1:2:4,1.4,'*nearest'),1); -%!error interp1(1,1,1, '*linear'); -%!assert (interp1(1:2:4,1:2:4,[0,1,1.4,3,4],'*linear'),[NaN,1,1.4,3,NaN]); -%!error interp1(1:3,1:3,1, '*cubic'); -%!assert (interp1(1:2:8,1:2:8,1.4,'*cubic'),1.4); -%!error interp1(1:2,1:2,1, '*spline'); -%!assert (interp1(1:2:6,1:2:6,1.4,'*spline'),1.4);
--- a/main/general/interp2.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,249 +0,0 @@ -## Copyright (C) 2000 Kai Habel -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; if not, write to the Free Software -## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{zi}=} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}) -## @deftypefnx {Function File} {@var{zi}=} interp2 (@var{Z}, @var{xi}, @var{yi}) -## @deftypefnx {Function File} {@var{zi}=} interp2 (@var{Z}, @var{n}) -## @deftypefnx {Function File} {@var{zi}=} interp2 (... , '@var{method}') -## Two-dimensional interpolation. @var{X},@var{Y} and @var{Z} describe a -## surface function. If @var{X} and @var{Y} are vectors their length -## must correspondent to the size of @var{Z}. @var{X} and @var{Y} must be -## monotonic. If they are matrices they must have the 'meshgrid' format. -## -## ZI = interp2 (X, Y, Z, XI, YI, ...) returns a matrix corresponding -## to the points described by the matrices @var{XI}, @var{YI}. -## -## If the last argument is a string, the interpolation method can -## be specified. The method can be 'linear', 'nearest' or 'cubic'. -## If it is omitted 'linear' interpolation is assumed. -## -## ZI = interp2 (Z, XI, YI) assumes X = 1:rows(Z) and Y = 1:columns(Z) -## -## ZI = interp2 (Z, n) interleaves the Matrix Z n-times. If n is ommited -## n = 1 is assumed -## -## @seealso{interp1} -## @end deftypefn - -## Author: Kai Habel <kai.habel@gmx.de> -## 2005-03-02 Thomas Weber <weber@num.uni-sb.de> -## * Add test cases -## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net> -## * Simplify -## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com> -## * 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. -## 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 = []; - method = "linear"; - - 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 - [Z,XI,YI,method] = deal(varargin{:}); - case 5 - [X,Y,Z,XI,YI] = deal(varargin{:}); - case 6 - [X,Y,Z,XI,YI,method] = deal(varargin{:}); - otherwise - usage ("interp2 (X, Y, Z, XI, YI, method)"); - endswitch - - % Type checking. - if !ismatrix(Z), error("interp2 expected matrix Z"); endif - if !isempty(n) && !isscalar(n), error("interp2 expected scalar n"); endif - if !ischar(method), error("interp2 expected string 'method'"); endif - - % Define X,Y,XI,YI if needed - [zr, zc] = size (Z); - if isempty(X), X=[1:zc]; Y=[1:zr]; endif - 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; endif - 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)) - error ("X and Y must be matrices of same size"); - endif - 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) - [XI, YI] = meshgrid (XI, 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)); - - 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 - ## - ## a-b - ## | | - ## c-d - a = Z(1:(zr - 1), 1:(zc - 1)); - b = Z(1:(zr - 1), 2:zc) - a; - c = Z(2:zr, 1:(zc - 1)) - a; - d = Z(2:zr, 2:zc) - a - b - c; - - idx = sub2ind(size(a),yidx,xidx); - - ## 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))'; - - ## apply plane equation - ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc; - - elseif strcmp (method, "nearest") - xtable = X(1, :); - 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); - ZI = Z(idx); - - elseif strcmp (method, "cubic") - ## FIXME bicubic doesn't handle arbitrary XI, YI - ZI = bicubic(X, Y, Z, XI(1,:), YI(:,1)); - - else - error ("interpolation method not (yet) supported"); - 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) ) = NaN; - - ZI = reshape(ZI,shape); - - -endfunction - -%!demo -%! A=[13,-1,12;5,4,3;1,6,2]; -%! x=[0,1,4]; y=[10,11,12]; -%! xi=linspace(min(x),max(x),17); -%! yi=linspace(min(y),max(y),26)'; -%! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear')); -%! [x,y] = meshgrid(x,y); -%! __gnuplot_raw__ ("set nohidden3d;\n") -%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; - -%!demo -%! A=[13,-1,12;5,4,3;1,6,2]; -%! x=[0,1,4]; y=[10,11,12]; -%! xi=linspace(min(x),max(x),17); -%! yi=linspace(min(y),max(y),26)'; -%! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest')); -%! [x,y] = meshgrid(x,y); -%! __gnuplot_raw__ ("set nohidden3d;\n") -%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; - -%!#demo -%! A=[13,-1,12;5,4,3;1,6,2]; -%! x=[0,1,2]; y=[10,11,12]; -%! xi=linspace(min(x),max(x),17); -%! yi=linspace(min(y),max(y),26); -%! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic')); -%! [x,y] = meshgrid(x,y); -%! __gnuplot_raw__ ("set nohidden3d;\n") -%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; - -%!test % simple test -%! x = [1,2,3]; -%! y = [4,5,6,7]; -%! [X, Y] = meshgrid(x,y); -%! Orig = X.^2 + Y.^3; -%! xi = [1.2,2, 1.5]; -%! yi = [6.2, 4.0, 5.0]'; -%! -%! Expected = ... -%! [243, 245.4, 243.9; -%! 65.6, 68, 66.5; -%! 126.6, 129, 127.5]; -%! Result = interp2(x,y,Orig, xi, yi); -%! -%! assert(Result, Expected, 1000*eps); - -%!test % 2^n form -%! x = [1,2,3]; -%! y = [4,5,6,7]; -%! [X, Y] = meshgrid(x,y); -%! Orig = X.^2 + Y.^3; -%! xi = [1:0.25:3]; yi = [4:0.25:7]'; -%! Expected = interp2(x,y,Orig, xi, yi); -%! Result = interp2(Orig,2); -%! -%! assert(Result, Expected, 10*eps); - -%!test % matrix slice -%! A = eye(4); -%! assert(interp2(A,[1:4],[1:4]),[1,1,1,1]); - -%!test % non-gridded XI,YI -%! A = eye(4); -%! assert(interp2(A,[1,2;3,4],[1,3;2,4]),[1,0;0,1]); - -%!test % for values outside of boundaries -%! x = [1,2,3]; -%! y = [4,5,6,7]; -%! [X, Y] = meshgrid(x,y); -%! Orig = X.^2 + Y.^3; -%! xi = [0,4]; -%! yi = [3,8]'; -%! assert(interp2(x,y,Orig, xi, yi),[nan,nan;nan,nan]); - -%!test % for values at boundaries -%! A=[1,2;3,4]; -%! x=[0,1]; -%! y=[2,3]'; -%! assert(interp2(x,y,A,x,y,'linear'), A); -%! assert(interp2(x,y,A,x,y,'nearest'), A); -
--- a/main/general/interpft.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,128 +0,0 @@ -## Copyright (C) 2001 Paul Kienzle -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; if not, write to the Free Software -## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -## y = interpft (x, n) -## Resample x with n points using Fourier interpolation. -## The data x must be equally spaced and n >= length(x). -## y = interpft (x, n, m) -## Resample matrix x to NxM useing Fourier interpolation - -## Author: Paul Kienzle -## 2001-02-11 -## * initial version -## 2002-03-17 aadler -## * added code to work on matrices as well -## $Id$ - -function z = interpft (x, n, m) - if (nargin < 2 || nargin > 3) - usage ("y=interpft (x, n)"); - endif - - transpose = ( rows (x) == 1 ); - if (transpose) x = x (:); endif - -if nargin==2 - [nr, nc] = size (x); - if n > nr - y = fft (x) / nr; - k = floor (nr / 2); - z = n * ifft ( [ y(1:k,:); zeros(n-nr,nc); y(k+1:nr) ] ); - elseif n < nr - warning ("interpft: Poor results possible: n should be bigger than x"); - ## XXX FIXME XXX the following doesn't work so well - y = fft (x) / nr; - k = ceil (n / 2); - z = n * ifft ( [ y(1:k,:); y(nr-k+2:nr) ] ); - else - z = x; - endif - -else - # this code is for matrices as well as vectors, - # it could replace the first code except - # 1) it depends on sparse being available - # 2) it's slightly slower - try save_empty_list_elements_ok= empty_list_elements_ok; - catch save_empty_list_elements_ok= 0; - end - try save_warn_empty_list_elements= warn_empty_list_elements; - catch save_warn_empty_list_elements= 0; - end - unwind_protect - empty_list_elements_ok= 1; - warn_empty_list_elements= 0; - - ny= n; - if nargin==2; my=1; else; my= m; endif - - [nx,mx]= size(x); - - if nx < ny - nhx= ceil( nx/2 ); - yconv= sparse( [ 1:nhx, ny-(nx-1-nhx:-1:0) ], 1:nx, 1, ny, nx ); - elseif nx > ny - nhy= ceil( ny/2 ); - yconv= sparse( 1:ny, [ 1:nhy, nx-(ny-1-nhy:-1:0) ], 1, ny, nx ); - elseif nx == ny - yconv= sparse( 1:ny, 1:nx, 1, ny, nx); - end - - if mx < my - mhx= ceil( mx/2 ); - xconv= sparse( 1:mx, [ 1:mhx, my-(mx-1-mhx:-1:0) ], 1, mx, my ); - elseif mx > my - mhy= ceil( my/2 ); - xconv= sparse( [ 1:mhy, mx-(my-1-mhy:-1:0) ], 1:my, 1, mx, my ); - elseif mx == my - xconv= sparse( 1:mx, 1:my, 1, mx, my); - end - - z= ifft2( yconv*fft2(x)*xconv) * (ny*my)/(nx*mx); - - unwind_protect_cleanup - empty_list_elements_ok= save_empty_list_elements_ok; - warn_empty_list_elements= save_warn_empty_list_elements; - end_unwind_protect - -endif - - if isreal (x), z = real (z); endif - if transpose, z = z.'; endif - -endfunction - -%!demo -%! t = 0 : 0.3 : pi; dt = t(2)-t(1); -%! n = length (t); k = 100; -%! ti = t(1) + [0 : k-1]*dt*n/k; -%! y = sin (4*t + 0.3) .* cos (3*t - 0.1); -%! yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1); -%! plot (ti, yp, 'g;sin(4t+0.3)cos(3t-0.1);', ... -%! ti, interp1 (t, y, ti, 'spline'), 'b;spline;', ... -%! ti, interpft (y ,k), 'c;interpft;', ... -%! t, y, 'r+;data;'); - -%!#demo -%! t = 0:0.3:pi; dt = t(2)-t(1); -%! n = length(t); k = 100; -%! ti = t(1)+[0:k-1]*dt*n/k; -%! y = sin (4*t+0.3).*cos(3*t-0.1); -%! yp = sin (4*ti+0.3).*cos(3*ti-0.1); -%! plot (ti, yp, 'g;sin(4t+0.3)cos(3t-0.1);', ... -%! t, interp1(ti,yp,t,'spline'), 'bx;spline;', ... -%! t, interpft(yp, n), 'co;interpft;', ... -%! t, y, 'r+;data;')
--- a/main/general/polyarea.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -## Copyright (C) 1999 David M. Doolin <doolin@ce.berkeley.edu> -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; if not, write to the Free Software -## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -## A = polyarea(X,Y) -## -## Determines area of a polygon by triangle method. -## -## input: -## polygon: x, y form vertex pairs, one polygon per column. -## -## output: -## area: Area of polygons. -## -## todo: Add moments for centroid, etc. -## -## bugs and limitations: -## Probably ought to be an optional check to make sure that -## traversing the vertices doesn't make any sides cross -## (Is simple closed curve the technical definition of this?). - -## Author: David M. Doolin <doolin@ce.berkeley.edu> -## Date: 1999/04/14 15:52:20 $ -## Modified-by: -## 2000-01-15 Paul Kienzle <pkienzle@kienzle.powernet.co.uk> -## * use matlab compatible interface -## * return absolute value of area so traversal order doesn't matter -## 2005-10-13 Torsten Finke -## * optimization saving half the sums and multiplies - -function a = polyarea (x, y) - if nargin != 2 - usage ("polyarea (x, y)"); - elseif any (size(x) != size(y)) - error ("polyarea: x and y must have the same shape"); - else - a = abs (sum (x .* (shift (y, -1) - shift (y, 1)))) / 2; - endif -endfunction - -%!assert (polyarea([1;1;3;3;1],[1;3;3;1;1]), 4, eps)
--- a/main/general/quadl.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,131 +0,0 @@ -% Copyright (C) 1998 Walter Gautschi -% -% This program is free software; you can redistribute it and/or modify it -% under the terms of the GNU General Public License as published by the -% Free Software Foundation; either version 2, or (at your option) any -% later version. -% -% This program is distributed in the hope that it will be useful, but WITHOUT -% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -% for more details. -% -% You should have received a copy of the GNU General Public License -% along with Octave; see the file COPYING. If not, write to the Free -% Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -%QUADL Numerically evaluate integral using adaptive -% Lobatto rule. -% -% Q=QUADL('F',A,B) approximates the integral of -% F(X) from A to B to machine precision. 'F' is a -% string containing the name of the function. The -% function F must return a vector of output values if -% given a vector of input values. -% -% Q=QUADL('F',A,B,TOL) integrates to a relative -% error of TOL. -% -% Q=QUADL('F',A,B,TOL,TRACE) displays the left -% end point of the current interval, the interval -% length, and the partial integral. -% -% Q=QUADL('F',A,B,TOL,TRACE,P1,P2,...) allows -% coefficients P1, ... to be passed directly to the -% function F: G=F(X,P1,P2,...). To use default values -% for TOL or TRACE, one may pass the empty matrix ([]). -% -% W. Gander and W. Gautschi, Adaptive Quadrature - Revisited, -% BIT Vol. 40, No. 1, March 2000, pp. 84--101. -% http://www.inf.ethz.ch/personal/gander/ - -% Author: Walter Gautschi -% Date: 08/03/98 -% Reference: Gander, Computermathematik, Birkhaeuser, 1992. - -% 2003-08-05 Shai Ayal -% * permission from author to release as GPL -% 2004-02-10 Paul Kienzle -% * renamed to quadl for compatibility -% * 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); - if(nargin<4), tol=[]; end; - if(nargin<5), trace=[]; end; - if(isempty(tol)), tol=eps; end; - if(isempty(trace)), trace=0; end; - if tol < eps - tol = eps; - end - m=(a+b)/2; 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{:}); - 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)); - s=sign(is); if(s==0), s=1; end; - erri1=abs(i1-is); - erri2=abs(i2-is); - R=1; if(erri2~=0), R=erri1/erri2; end; - if(R>0 & R<1), tol=tol/R; end; - is=s*abs(is)*tol/eps; - if(is==0), is=b-a, end; - Q=adaptlobstp(f,a,b,fa,fb,is,trace,varargin{:}); - - - -function Q=adaptlobstp(f,a,b,fa,fb,is,trace,varargin) -%ADAPTLOBSTP Recursive function used by ADAPTLOB. -% -% 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 -% arguments are generated by ADAPTLOB or by recursion. -% -% See also ADAPTLOB. - -% Walter Gautschi, 08/03/98 - - h=(b-a)/2; m=(a+b)/2; - alpha=sqrt(2/3); beta=1/sqrt(5); - mll=m-alpha*h; ml=m-beta*h; mr=m+beta*h; mrr=m+alpha*h; - 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); - end; - Q=i1; - if(trace), disp([a b-a Q]), end; - 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{:}); - end - -function r = need_warning(v) - persistent w = []; - if nargin == 0, r = w; - else w=v; end
--- a/main/miscellaneous/inputname.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -## s=inputname(n) -## Return the text defining nth input to the function. -function s=inputname(n) - s=evalin('caller',sprintf('deblank(argn(%d,:));',n)); - -## Warning: heap big magic in the following tests!!! -## The test function builds a private context for each -## test, with only the specified values shared between -## them. It does this using the following template: -## -## function [<shared>] = testfn(<shared>) -## <test> -## -## To test inputname, I need a function context invoked -## with known parameter names. So define a couple of -## shared parameters, et voila!, the test is trivial. -%!shared hello,worldly -%!assert(inputname(1),'hello'); -%!assert(inputname(2),'worldly');
--- a/main/plot/__plt3__.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -## Copyright (C) 1996 John W. Eaton -## -## This program is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2, or (at your option) -## any later version. -## -## This program is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; see the file COPYING. If not, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## __plt3__ is not a user callable function - -## Author: Paul Kienzle <kienzle.powernet.co.uk> -## 2001-04-06 Paul Kienzle <kienzle.powernet.co.uk> -## * __gnuplot_set__ nohidden3d; vector X,Y, matrix Z => meshgrid(X,Y) - -## Modified to use new gnuplot interface in octave > 2.9.0 -## Dmitri A. Sergatskov <dasergatskov@gmail.com> -## April 18, 2005 - -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)"); - else - [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"); - endif - - unwind_protect - __gnuplot_set__ parametric; - __gnuplot_raw__ ("set nohidden3d;\n"); - for i=1:columns(x) - tmp = [x(:,i), y(:,i), z(:,i)]; - cmd = ["__gnuplot_splot__ tmp ", fmt, ";\n"]; - eval (cmd); - endfor - unwind_protect_cleanup - __gnuplot_set__ noparametric; - end_unwind_protect -endfunction
--- a/main/plot/plot3.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,111 +0,0 @@ -## Copyright (C) 1996 John W. Eaton -## -## This program is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2, or (at your option) -## any later version. -## -## This program is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; see the file COPYING. If not, write to the Free -## Software Foundation, 59 Temple Place - Suite 330, Boston, MA -## 02111-1307, USA. - -## usage: plot3 (x, y, z) -## plot3 (x1, y1, z1, x2, y2, z2, ...) -## plot3 (x, y, z, fmt) -## -## If all arguments are vectors of the same length, a single line is -## drawn in three space. -## -## If all arguments are matrices, each column is drawn as a separate -## line. In this case, all matrices must have the same number of rows -## and columns and no attempt is made to transpose the arguments to make -## the number of rows match. -## -## To see possible options for FMT please see __pltopt__. -## -## Example -## -## z = [0:0.05:5]; -## plot3(cos(2*pi*z), sin(2*pi*z), z, ";helix;"); - -## Author: Paul Kienzle -## (modified from __plt__.m) - -function plot3(varargin) - - ## use the following for 2.1.38 and below - # varargin = list(varargin, all_va_args); - - hold_state = ishold (); - - unwind_protect - - x_set = 0; - y_set = 0; - z_set = 0; - - ## Gather arguments, decode format, and plot lines. - for arg = 1:length(varargin) - new = nth (varargin, arg); - - if (ischar (new)) - if (! z_set) - error ("plot3: needs x, y, z"); - endif - fmt = __pltopt__ ("plot3", new); - __plt3__(x, y, z, fmt); - hold on; - x_set = 0; - y_set = 0; - z_set = 0; - elseif (!x_set) - x = new; - x_set = 1; - elseif (!y_set) - y = new; - y_set = 1; - elseif (!z_set) - z = new; - z_set = 1; - else - __plt3__ (x, y, z, ""); - hold on; - x = new; - y_set = 0; - z_set = 0; - endif - -### Code copied from __plt__; don't know if it is needed -### -### ## Something fishy is going on. I don't think this should be -### ## necessary, but without it, sometimes not all the lines from a -### ## given plot command appear on the screen. Even with it, the -### ## delay might not be long enough for some systems... -### -### usleep (1e5); - - endfor - - ## Handle last plot. - - if (z_set) - __plt3__ (x, y, z, ""); - elseif (x_set) - error ("plot3: needs x, y, z"); - endif - - unwind_protect_cleanup - - if (! hold_state) - hold off; - endif - - end_unwind_protect - -endfunction
--- a/main/sparse/pcg.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,440 +0,0 @@ -function [x, flag, relres, iter, resvec, eigest] = ... - pcg( A, b, tol, maxit, M, x0, varargin ) -%[X, FLAG, RELRES, ITER, RESVEC, EIGEST] = pcg( A, B, TOL, MAXIT, M, X0, ... ) -% -%Solves the linear system of equations A*x = B by means of the Preconditioned -%Conjugate Gradient iterative method. -% -%INPUT ARGUMENTS -%--------------- -% -%A can be either a square (preferably sparse) matrix or a name of a function -%which computes A*x. In principle A should be symmetric and positive -%definite; if PCG finds A not positive definite, you will get a warning message -%and the FLAG output parameter will be set. -% -%B is the right hand side vector. -% -%TOL is the required relative tolerance for the residual error, B-A*X. The -%iteration stops if ||B-A*X|| <= TOL*||B-A*X0||, where ||.|| denotes the -%Euclidean norm. If TOL is empty or is omitted, the function sets TOL = 1e-6 by -%default. -% -%MAXIT is the maximum allowable number of iterations; if [] is supplied for -%MAXIT, or PCG has less arguments, a default value equal to 20 is used. -% -%M is the (left) preconditioning matrix, so that the iteration is -%(theoretically) equivalent to solving by PCG P*x = M\B, with P = M\A. -%Note that a proper choice of the preconditioner may -%dramatically improve the overall performance of the method! Instead of matrix M, -%the user may pass a function which returns the results of applying the inverse of -%M to a vector (usually this is the preferred way of using the preconditioner). -%If [] is supplied for M, or M is omitted, no preconditioning is -%applied. -% -%X0 is the initial guess. If X0 is empty or omitted, the function sets X0 to a -%zero vector by default. -% -%The arguments which follow X0 are treated as parameters, and passed in a proper -%way to any of the functions (A or M) which are passed to pcg. See the EXAMPLES -%for details. -% -%OUTPUT ARGUMENTS -%---------------- -% -%X is the computed approximation to the solution x of A*x=B. -% -%FLAG reports on the convergence. FLAG = 0 means the solution converged -%and the tolerance criterion given by TOL is satisfied. FLAG = 1 means that the -%MAXIT limit for the iteration count was reached. FLAG = 3 reports that the -%(preconditioned) matrix was found not positive definite. -% -%RELRES is the ratio of the final residual to its initial value, measured in the -%Euclidean norm. -% -%ITER is the actual number of iterations performed. -% -%RESVEC describes the convergence history of the method. RESVEC(i,1) is the -%Euclidean norm of the residual, and RESVEC(i,2) is the preconditioned residual -%norm, after the (i-1)-th iteration, i = 1,2,...ITER+1. The preconditioned -%residual norm is defined as |||r|||^2 = r'*(M\r) where r = B-A*x, see also the -%description of M. If EIGEST is not required, only RESVEC(:,1) is returned. -% -%EIGEST returns the estimate for the smallest (EIGEST(1)) and largest -%(EIGEST(2)) eigenvalues of the preconditioned matrix P=M\A. -%In particular, if no preconditioning is used, the -%estimates for the extreme eigenvalues of A are returned. EIGEST(1) is an -%overestimate and EIGEST(2) is an underestimate, so that EIGEST(2)/EIGEST(1) is -%a lower bound for cond(P,2), which nevertheless in the limit should -%theoretically be equal to the actual value of the condition number. -%The method which computes EIGEST works only for symmetric positive -%definite A and M, and the user is responsible for verifying this assumption. -% -%EXAMPLES -%-------- -% -%Let us consider a trivial problem with a diagonal matrix (we exploit the -%sparsity of A) -% -% N = 10; -% A = diag([1:N]); A = sparse(A); -% b = rand(N,1); -% -%EX. 1. Simplest use of PCG -% -% x = pcg(A,b) -% -%EX. 2. PCG with a function which computes A*x -% -% function y = applyA(x) y = [1:N]'.*x; end -% -% x = pcg('applyA',b) -% -%EX. 3. Preconditioned iteration, with full diagnostics. The preconditioner (quite -%strange, because even the original matrix A is trivial) is defined as a -%function: -% -% function y = applyM(x) -% K = floor(length(x)-2); -% y = x; -% y(1:K) = x(1:K)./[1:K]'; -% end -% -% [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],[],'applyM') -% semilogy([1:iter+1], resvec); -% -%EX. 4. Finally, a preconditioner which depends on a parameter K: -% -% function y = applyM(x, varargin) -% K = varargin{1}; -% y = x; y(1:K) = x(1:K)./[1:K]'; -% end -% -% [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],[],'applyM',[],3) -% -%You can also run -% -% demo('pcg') -% -%from the command line to see more simple examples of how the pcg works. -% -%SEE ALSO: sparse, pcr, gmres -% -%REFERENCES -% -% [1] C.T.Kelley, 'Iterative methods for linear and nonlinear equations', -% SIAM, 1995 (the base PCG algorithm) -% -% [2] Y.Saad, 'Iterative methods for sparse linear systems', PWS 1996 -% (condition number estimate from PCG) Revised version of this book is -% available online at http://www-users.cs.umn.edu/~saad/books.html -% - -%% Copyright (C) 2004 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl> -%% -%% This program is free software; you can redistribute it and/or modify -%% it under the terms of the GNU General Public License as published by -%% the Free Software Foundation; either version 2 of the License, or -%% (at your option) any later version. -%% -%% This program is distributed in the hope that it will be useful, -%% but WITHOUT ANY WARRANTY; without even the implied warranty of -%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -%% GNU General Public License for more details. -%% -%% You should have received a copy of the GNU General Public License -%% along with this program; if not, write to the Free Software -%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -%% -%% 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 -%% -if ((nargin < 6) || isempty(x0)) - x = zeros(size(b)); -else - x = x0; -end - -if (nargin < 5) - M = []; -end - -if ((nargin < 4) || isempty(maxit)) - maxit = min(size(b,1),20); -end - -maxit = maxit + 2; - -if ((nargin < 3) || isempty(tol)) - tol = 1e-6; -end - -preconditioned_residual_out = false; -if (nargout > 5) - T = zeros(maxit,maxit); - preconditioned_residual_out = true; -end - -matrix_positive_definite = true; % assume A is positive definite - -if (~exist('OCTAVE_VERSION')) % a trick to maintain source - % compatibility with M*TLAB - stderr = 2; -end - -p = zeros(size(b)); -oldtau = 1; -if (isnumeric(A)) % is A a matrix? - r = b - A*x; -else % then A should be a function! - r = b - feval(A,x,varargin{:}); -end - -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 - z = r; - else % otherwise, apply the precond - z = M\r; - end - else % then M should be a function! - z = feval(M,r,varargin{:}); - end - 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; - else % then A should be a function! - w = feval(A,p,varargin{:}); - end - oldalpha = alpha; % needed only for eigest - alpha = tau/(p'*w); - if (alpha <= 0.0) % negative matrix? - matrix_positive_definite = false; - end - x = x + alpha*p; - r = 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); - end; - resvec(iter,1) = norm(r); - iter = iter +1; -end - - -if (nargout > 5) - 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)); - else - eigest = [NaN NaN]; - warning('PCG: eigenvalue estimate failed: iteration converged too fast.'); - end - else - eigest = [NaN NaN]; - end - - % 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 - z = r; - else % otherwise, apply the precond - z = M\r; - end - else % then M should be a function! - z = feval(M,r,varargin{:}); - end - resvec(iter-1,2) = sqrt(r'*z); -else - resvec = resvec(:,1); -end - -flag = 0; -relres = resvec(iter-1,1)./resvec(1,1); -iter = 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); - end -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); - end -end -if (~matrix_positive_definite) - flag = 3; - if (nargout < 2) - warning('PCG: matrix not positive definite?\n'); - end -end - -end -%!demo -%! -%! # Simplest usage of pcg (see also 'help pcg') -%! -%! N = 10; -%! A = diag([1:N]); b = rand(N,1); y = A\b; #y is the true solution -%! x = pcg(A,b); -%! printf('The solution relative error is %g\n', norm(x-y)/norm(y)); -%! -%! # You shouldn't be afraid if pcg issues some warning messages in this -%! # example: watch out in the second example, why it takes N iterations -%! # of pcg to converge to (a very accurate, by the way) solution -%!demo -%! -%! # Full output from pcg, except for the eigenvalue estimates -%! # We use this output to plot the convergence history -%! -%! N = 10; -%! A = diag([1:N]); b = rand(N,1); X = A\b; #X is the true solution -%! [x, flag, relres, iter, resvec] = pcg(A,b); -%! printf('The solution relative error is %g\n', norm(x-X)/norm(X)); -%! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)'); -%! semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;'); -%!demo -%! -%! # Full output from pcg, including the eigenvalue estimates -%! # Hilbert matrix is extremely ill conditioned, so pcg WILL have problems -%! -%! N = 10; -%! A = hilb(N); b = rand(N,1); X = A\b; #X is the true solution -%! [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],200); -%! printf('The solution relative error is %g\n', norm(x-X)/norm(X)); -%! printf('Condition number estimate is %g\n', eigest(2)/eigest(1)); -%! printf('Actual condition number is %g\n', cond(A)); -%! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); -%! semilogy([0:iter],resvec,['o-g;absolute residual;';'+-r;absolute preconditioned residual;']); -%!demo -%! -%! # Full output from pcg, including the eigenvalue estimates -%! # We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2) -%! # and that's the reasone we need some preconditioner; here we take -%! # a very simple and not powerful Jacobi preconditioner, -%! # which is the diagonal of A -%! -%! N = 100; -%! A = zeros(N,N); -%! for i=1:N-1 # form 1-D Laplacian matrix -%! A(i:i+1,i:i+1) = [2 -1; -1 2]; -%! endfor -%! b = rand(N,1); X = A\b; #X is the true solution -%! maxit = 80; -%! printf('System condition number is %g\n',cond(A)); -%! # No preconditioner: the convergence is very slow! -%! -%! [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],maxit); -%! printf('System condition number estimate is %g\n',eigest(2)/eigest(1)); -%! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); -%! semilogy([0:iter],resvec(:,1),'o-g;NO preconditioning: absolute residual;'); -%! -%! pause(1); -%! # Test Jacobi preconditioner: it will not help much!!! -%! -%! M = diag(diag(A)); # Jacobi preconditioner -%! [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],maxit,M); -%! printf('JACOBI preconditioned system condition number estimate is %g\n',eigest(2)/eigest(1)); -%! hold on; -%! semilogy([0:iter],resvec(:,1),'o-r;JACOBI preconditioner: absolute residual;'); -%! -%! pause(1); -%! # Test nonoverlapping block Jacobi preconditioner: it will help much! -%! -%! M = zeros(N,N);k=4 -%! for i=1:k:N # form 1-D Laplacian matrix -%! M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1); -%! endfor -%! [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],maxit,M); -%! printf('BLOCK JACOBI preconditioned system condition number estimate is %g\n',eigest(2)/eigest(1)); -%! semilogy([0:iter],resvec(:,1),'o-b;BLOCK JACOBI preconditioner: absolute residual;'); -%! hold off; -%!test -%! -%! #solve small diagonal system -%! -%! N = 10; -%! A = diag([1:N]); b = rand(N,1); X = A\b; #X is the true solution -%! [x, flag] = pcg(A,b,[],N+1); -%! assert(norm(x-X)/norm(X),0,1e-10); -%! assert(flag,0); -%! -%!test -%! -%! #solve small indefinite diagonal system -%! #despite A is indefinite, the iteration continues and converges -%! #indefiniteness of A is detected -%! -%! N = 10; -%! A = diag([1:N].*(-ones(1,N).^2)); b = rand(N,1); X = A\b; #X is the true solution -%! [x, flag] = pcg(A,b,[],N+1); -%! assert(norm(x-X)/norm(X),0,1e-10); -%! assert(flag,3); -%! -%!test -%! -%! #solve tridiagonal system, do not converge in default 20 iterations -%! -%! N = 100; -%! A = zeros(N,N); -%! for i=1:N-1 # form 1-D Laplacian matrix -%! A(i:i+1,i:i+1) = [2 -1; -1 2]; -%! endfor -%! b = ones(N,1); X = A\b; #X is the true solution -%! [x, flag, relres, iter, resvec, eigest] = pcg(A,b,1e-12); -%! assert(flag); -%! assert(relres>1.0); -%! assert(iter,20); #should perform max allowable default number of iterations -%! -%!test -%! -%! #solve tridiagonal system with 'prefect' preconditioner -%! #converges in one iteration, so the eigest does not work -%! #and issues a warning -%! -%! N = 100; -%! A = zeros(N,N); -%! for i=1:N-1 # form 1-D Laplacian matrix -%! A(i:i+1,i:i+1) = [2 -1; -1 2]; -%! endfor -%! b = ones(N,1); X = A\b; #X is the true solution -%! [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],[],A,b); -%! assert(norm(x-X)/norm(X),0,1e-6); -%! assert(flag,0); -%! assert(iter,1); #should converge in one iteration -%! assert(isnan(eigest),isnan([NaN NaN])); -%!
--- a/main/sparse/pcr.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,401 +0,0 @@ -function [x, flag, relres, iter, resvec] = pcr(A,b,tol,maxit,M,x0,varargin) -% -%[X, FLAG, RELRES, ITER, RESVEC] = pcr( A, B, TOL, MAXIT, M, X0, ... ) -% -%Solves the linear system of equations A*x = B by means of the Preconditioned -%Conjugate Residuals iterative method. -% -%INPUT ARGUMENTS -%--------------- -% -%A can be either a square (preferably sparse) matrix or a name of a function -%which computes A*x. In principle A should be symmetric and nonsingular; -%if PCR finds A (numerically) singular, you will get a warning message -%and the FLAG output parameter will be set. -% -%B is the right hand side vector. -% -%TOL is the required relative tolerance for the residual error, B-A*X. The -%iteration stops if ||B-A*X|| <= TOL*||B-A*X0||, where ||.|| denotes the -%Euclidean norm. If TOL is empty or is omitted, the function sets TOL = 1e-6 by -%default. -% -%MAXIT is the maximum allowable number of iterations; if [] is supplied for -%MAXIT, or PCR has less arguments, a default value equal to 20 is used. -% -%M is the (left) preconditioning matrix, so that the iteration is -%(theoretically) equivalent to solving by PCR P*x = M\B, with P = M\A. -%Instead of a matrix, the user may pass here a name of a function which returns -%the result of applying the inverse of M to a vector (usually this is the -%preferred way of using the preconditioner), see EXAMPLES for details. -%Usually it is assumed that M is positive definite. Note that a proper choice of -%the preconditioner may dramatically improve the overall performance of the -%method! If [] is supplied for M, or PCR has less arguments, no preconditioning -%is applied. -% -%X0 is the initial guess. If X0 is empty or omitted, the function sets X0 to a -%zero vector by default. -% -%The arguments which follow X0 are treated as parameters, and passed in a proper -%way to any of the functions (A or M) which are passed to PCR. See the EXAMPLES -%for details. -% -%OUTPUT ARGUMENTS -%---------------- -% -%X is the computed approximation to the solution x of A*x=B. -% -%FLAG reports on the convergence. FLAG = 0 means the solution converged -%and the tolerance criterion given by TOL is satisfied. FLAG = 1 means that the -%limit MAXIT for the iteration count was reached. FLAG = 3 reports -%PCR breakdown, see [1] for details. -% -%RELRES is the ratio of the final residual to its initial value, measured in the -%Euclidean norm. -% -%ITER is the actual number of iterations performed. -% -%RESVEC describes the convergence history of the method, so that RESVEC(i) -%contains the Euclidean norms of the residual after the (i-1)-th iteration, -%i = 1,2,...ITER+1. -% -%EXAMPLES -%-------- -% -%Let us consider a trivial problem with a diagonal matrix (we exploit the -%sparsity of A) -% -% N = 10; -% A = diag([1:N]); A = sparse(A); -% b = rand(N,1); -% -%EX. 1. Simplest use of PCR -% -% x = pcr(A,b) -% -%EX. 2. PCR with a function which computes A*x -% -% function y = applyA(x) y = [1:10]'.*x; end -% -% x = pcr('applyA',b) -% -%EX. 3. Preconditioned iteration, with full diagnostics. The preconditioner (quite -%strange, because even the original matrix A is trivial) is defined as a -%function: -% -% function y = applyM(x) -% K = floor(length(x)-2); -% y = x; -% y(1:K) = x(1:K)./[1:K]'; -% end -% -% [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM') -% semilogy([1:iter+1], resvec); -% -%EX. 4. Finally, a preconditioner which depends on a parameter K: -% -% function y = applyM(x, varargin) -% K = varargin{1}; -% y = x; y(1:K) = x(1:K)./[1:K]'; -% end -% -% [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM',[],3) -% -%You can also run -% -% demo('pcr') -% -%from the command line to see simple examples of how the PCR works. -% -%SEE ALSO: sparse, pcg, gmres -% -%REFERENCES -% -% [1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of -% Equations", section 9.5.4; Springer, 1994 - -%% Copyright (C) 2004 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl> -%% -%% This program is free software; you can redistribute it and/or modify -%% it under the terms of the GNU General Public License as published by -%% the Free Software Foundation; either version 2 of the License, or -%% (at your option) any later version. -%% -%% This program is distributed in the hope that it will be useful, -%% but WITHOUT ANY WARRANTY; without even the implied warranty of -%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -%% GNU General Public License for more details. -%% -%% You should have received a copy of the GNU General Public License -%% along with this program; if not, write to the Free Software -%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -%% -%% 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 - -breakdown = false; - -if ((nargin < 6) || isempty(x0)) - x = zeros(size(b)); -else - x = x0; -end - -if (nargin < 5) - M = []; -end - -if ((nargin < 4) || isempty(maxit)) - maxit = 20; -end - -maxit = maxit + 2; - -if ((nargin < 3) || isempty(tol)) - tol = 1e-6; -end - -if (nargin < 2) - usage('x=pcr(A,b)'); -end - -% init -if (isnumeric(A)) % is A a matrix? - r = b-A*x; -else % then A should be a function! - r = b-feval(A,x,varargin{:}); -end - -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; - end -else % then M should be a function! - p = feval(M,r,varargin{:}); -end - -iter = 2; - -b_bot_old = 1; -q_old = p_old = s_old = zeros(size(x)); - - if (isnumeric(A)) % is A a matrix? - q = A*p; - else % then A should be a function! - q = feval(A,p,varargin{:}); - end - -resvec(1) = abs(norm(r)); - -%iteration -while ( (resvec(iter-1) > tol*resvec(1)) && (iter < maxit) ), - - if (isnumeric(M)) % is M a matrix? - if isempty(M) % if M is empty, use no precond - s = q; - else % otherwise, apply the precond - s = M\q; - end - else % then M should be a function! - s = feval(M,q,varargin{:}); - end - b_top = r'*s; - b_bot = q'*s; - - if (b_bot == 0.0) - breakdown = true; - break; - end - lambda = b_top/b_bot; - - x = x + lambda*p; - r = 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{:}); - end - - alpha0 = (t'*s)/b_bot; - alpha1 = (t'*s_old)/b_bot_old; - - p_temp = p; q_temp = q; - p = s - alpha0*p - alpha1*p_old; - q = t - alpha0*q - alpha1*q_old; - - s_old = s; - p_old = p_temp; - q_old = q_temp; - b_bot_old = b_bot; - - - resvec(iter) = abs(norm(r)); - iter = iter + 1; -end - -flag = 0; -relres = resvec(iter-1)./resvec(1); -iter = 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); - end -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); - end -end -if (breakdown) - flag = 3; - if (nargout < 2) - warning('PCR: breakdown occured.\n'); - warning('System matrix singular or preconditioner indefinite?\n'); - end -end - -end -%!demo -%! -%! # Simplest usage of PCR (see also 'help pcr') -%! -%! N = 20; -%! A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution -%! x = pcr(A,b); -%! printf('The solution relative error is %g\n', norm(x-y)/norm(y)); -%! -%! # You shouldn't be afraid if PCR issues some warning messages in this -%! # example: watch out in the second example, why it takes N iterations -%! # of PCR to converge to (a very accurate, by the way) solution -%!demo -%! -%! # Full output from PCR -%! # We use this output to plot the convergence history -%! -%! N = 20; -%! A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution -%! [x, flag, relres, iter, resvec] = pcr(A,b); -%! printf('The solution relative error is %g\n', norm(x-X)/norm(X)); -%! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)'); -%! semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;'); -%!demo -%! -%! # Full output from PCR -%! # We use indefinite matrix based on the Hilbert matrix, with one -%! # strongly negative eigenvalue -%! # Hilbert matrix is extremely ill conditioned, so is ours, -%! # and that's why PCR WILL have problems -%! -%! N = 10; -%! A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution -%! printf('Condition number of A is %g\n', cond(A)); -%! [x, flag, relres, iter, resvec] = pcr(A,b,[],200); -%! if (flag == 3) -%! printf('PCR breakdown. System matrix is [close to] singular\n'); -%! end -%! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); -%! semilogy([0:iter],resvec,'o-g;absolute residual;'); -%!demo -%! -%! # Full output from PCR -%! # We use an indefinite matrix based on the 1-D Laplacian matrix for A, -%! # and here we have cond(A) = O(N^2) -%! # That's the reason we need some preconditioner; here we take -%! # a very simple and not powerful Jacobi preconditioner, -%! # which is the diagonal of A -%! -%! # Note that we use here indefinite preconditioners! -%! -%! N = 100; -%! A = zeros(N,N); -%! for i=1:N-1 # form 1-D Laplacian matrix -%! A(i:i+1,i:i+1) = [2 -1; -1 2]; -%! endfor -%! A = [A, zeros(size(A)); zeros(size(A)), -A]; -%! b = rand(2*N,1); X = A\b; #X is the true solution -%! maxit = 80; -%! printf('System condition number is %g\n',cond(A)); -%! # No preconditioner: the convergence is very slow! -%! -%! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit); -%! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); -%! semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;'); -%! -%! pause(1); -%! # Test Jacobi preconditioner: it will not help much!!! -%! -%! M = diag(diag(A)); # Jacobi preconditioner -%! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M); -%! hold on; -%! semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;'); -%! -%! pause(1); -%! # Test nonoverlapping block Jacobi preconditioner: this one should give -%! # some convergence speedup! -%! -%! M = zeros(N,N);k=4; -%! for i=1:k:N # get k x k diagonal blocks of A -%! M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1); -%! endfor -%! M = [M, zeros(size(M)); zeros(size(M)), -M]; -%! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M); -%! semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;'); -%! hold off; -%!test -%! -%! #solve small indefinite diagonal system -%! -%! N = 10; -%! A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution -%! [x, flag] = pcr(A,b,[],N+1); -%! assert(norm(x-X)/norm(X)<1e-10); -%! assert(flag,0); -%! -%!test -%! -%! #solve tridiagonal system, do not converge in default 20 iterations -%! #should perform max allowable default number of iterations -%! -%! N = 100; -%! A = zeros(N,N); -%! for i=1:N-1 # form 1-D Laplacian matrix -%! A(i:i+1,i:i+1) = [2 -1; -1 2]; -%! endfor -%! b = ones(N,1); X = A\b; #X is the true solution -%! [x, flag, relres, iter, resvec] = pcr(A,b,1e-12); -%! assert(flag,1); -%! assert(relres>0.6); -%! assert(iter,20); -%! -%!test -%! -%! #solve tridiagonal system with 'prefect' preconditioner -%! #converges in one iteration -%! -%! N = 100; -%! A = zeros(N,N); -%! for i=1:N-1 # form 1-D Laplacian matrix -%! A(i:i+1,i:i+1) = [2 -1; -1 2]; -%! endfor -%! b = ones(N,1); X = A\b; #X is the true solution -%! [x, flag, relres, iter] = pcr(A,b,[],[],A,b); -%! assert(norm(x-X)/norm(X)<1e-6); -%! assert(relres<1e-6); -%! assert(flag,0); -%! assert(iter,1); #should converge in one iteration -%!
--- a/main/splines/Makefile Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -sinclude ../../Makeconf - -PCHIP_OBJ = pchip_deriv.o dpchim.o dpchst.o - -all: pchip_deriv.oct - -pchip_deriv.oct: $(PCHIP_OBJ) - $(MKOCTFILE) -v -o $@ $(PCHIP_OBJ) - -clean: ; -$(RM) *.o octave-core core *.oct *~
--- a/main/splines/dpchim.f Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,283 +0,0 @@ -*DECK DPCHIM - SUBROUTINE DPCHIM (N, X, F, D, INCFD, IERR) -C***BEGIN PROLOGUE DPCHIM -C***PURPOSE Set derivatives needed to determine a monotone piecewise -C cubic Hermite interpolant to given data. Boundary values -C are provided which are compatible with monotonicity. The -C interpolant will have an extremum at each point where mono- -C tonicity switches direction. (See DPCHIC if user control -C is desired over boundary or switch conditions.) -C***LIBRARY SLATEC (PCHIP) -C***CATEGORY E1A -C***TYPE DOUBLE PRECISION (PCHIM-S, DPCHIM-D) -C***KEYWORDS CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION, -C PCHIP, PIECEWISE CUBIC INTERPOLATION -C***AUTHOR Fritsch, F. N., (LLNL) -C Lawrence Livermore National Laboratory -C P.O. Box 808 (L-316) -C Livermore, CA 94550 -C FTS 532-4275, (510) 422-4275 -C***DESCRIPTION -C -C DPCHIM: Piecewise Cubic Hermite Interpolation to -C Monotone data. -C -C Sets derivatives needed to determine a monotone piecewise cubic -C Hermite interpolant to the data given in X and F. -C -C Default boundary conditions are provided which are compatible -C with monotonicity. (See DPCHIC if user control of boundary con- -C ditions is desired.) -C -C If the data are only piecewise monotonic, the interpolant will -C have an extremum at each point where monotonicity switches direc- -C tion. (See DPCHIC if user control is desired in such cases.) -C -C To facilitate two-dimensional applications, includes an increment -C between successive values of the F- and D-arrays. -C -C The resulting piecewise cubic Hermite function may be evaluated -C by DPCHFE or DPCHFD. -C -C ---------------------------------------------------------------------- -C -C Calling sequence: -C -C PARAMETER (INCFD = ...) -C INTEGER N, IERR -C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) -C -C CALL DPCHIM (N, X, F, D, INCFD, IERR) -C -C Parameters: -C -C N -- (input) number of data points. (Error return if N.LT.2 .) -C If N=2, simply does linear interpolation. -C -C X -- (input) real*8 array of independent variable values. The -C elements of X must be strictly increasing: -C X(I-1) .LT. X(I), I = 2(1)N. -C (Error return if not.) -C -C F -- (input) real*8 array of dependent variable values to be -C interpolated. F(1+(I-1)*INCFD) is value corresponding to -C X(I). DPCHIM is designed for monotonic data, but it will -C work for any F-array. It will force extrema at points where -C monotonicity switches direction. If some other treatment of -C switch points is desired, DPCHIC should be used instead. -C ----- -C D -- (output) real*8 array of derivative values at the data -C points. If the data are monotonic, these values will -C determine a monotone cubic Hermite function. -C The value corresponding to X(I) is stored in -C D(1+(I-1)*INCFD), I=1(1)N. -C No other entries in D are changed. -C -C INCFD -- (input) increment between successive values in F and D. -C This argument is provided primarily for 2-D applications. -C (Error return if INCFD.LT.1 .) -C -C IERR -- (output) error flag. -C Normal return: -C IERR = 0 (no errors). -C Warning error: -C IERR.GT.0 means that IERR switches in the direction -C of monotonicity were detected. -C "Recoverable" errors: -C IERR = -1 if N.LT.2 . -C IERR = -2 if INCFD.LT.1 . -C IERR = -3 if the X-array is not strictly increasing. -C (The D-array has not been changed in any of these cases.) -C NOTE: The above errors are checked in the order listed, -C and following arguments have **NOT** been validated. -C -C***REFERENCES 1. F. N. Fritsch and J. Butland, A method for construc- -C ting local monotone piecewise cubic interpolants, SIAM -C Journal on Scientific and Statistical Computing 5, 2 -C (June 1984), pp. 300-304. -C 2. F. N. Fritsch and R. E. Carlson, Monotone piecewise -C cubic interpolation, SIAM Journal on Numerical Ana- -C lysis 17, 2 (April 1980), pp. 238-246. -C***ROUTINES CALLED DPCHST, XERMSG -C***REVISION HISTORY (YYMMDD) -C 811103 DATE WRITTEN -C 820201 1. Introduced DPCHST to reduce possible over/under- -C flow problems. -C 2. Rearranged derivative formula for same reason. -C 820602 1. Modified end conditions to be continuous functions -C of data when monotonicity switches in next interval. -C 2. Modified formulas so end conditions are less prone -C of over/underflow problems. -C 820803 Minor cosmetic changes for release 1. -C 870707 Corrected XERROR calls for d.p. name(s). -C 870813 Updated Reference 1. -C 890206 Corrected XERROR calls. -C 890411 Added SAVE statements (Vers. 3.2). -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890703 Corrected category record. (WRB) -C 890831 Modified array declarations. (WRB) -C 891006 Cosmetic changes to prologue. (WRB) -C 891006 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) -C 920429 Revised format and order of references. (WRB,FNF) -C***END PROLOGUE DPCHIM -C Programming notes: -C -C 1. The function DPCHST(ARG1,ARG2) is assumed to return zero if -C either argument is zero, +1 if they are of the same sign, and -C -1 if they are of opposite sign. -C 2. To produce a single precision version, simply: -C a. Change DPCHIM to PCHIM wherever it occurs, -C b. Change DPCHST to PCHST wherever it occurs, -C c. Change all references to the Fortran intrinsics to their -C single precision equivalents, -C d. Change the double precision declarations to real, and -C e. Change the constants ZERO and THREE to single precision. -C -C DECLARE ARGUMENTS. -C - INTEGER N, INCFD, IERR - DOUBLE PRECISION X(*), F(INCFD,*), D(INCFD,*) -C -C DECLARE LOCAL VARIABLES. -C - INTEGER I, NLESS1 - DOUBLE PRECISION DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE, - * H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO - SAVE ZERO, THREE - DOUBLE PRECISION DPCHST - DATA ZERO /0.D0/, THREE/3.D0/ -C -C VALIDITY-CHECK ARGUMENTS. -C -C***FIRST EXECUTABLE STATEMENT DPCHIM - IF ( N.LT.2 ) GO TO 5001 - IF ( INCFD.LT.1 ) GO TO 5002 - DO 1 I = 2, N - IF ( X(I).LE.X(I-1) ) GO TO 5003 - 1 CONTINUE -C -C FUNCTION DEFINITION IS OK, GO ON. -C - IERR = 0 - NLESS1 = N - 1 - H1 = X(2) - X(1) - DEL1 = (F(1,2) - F(1,1))/H1 - DSAVE = DEL1 -C -C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. -C - IF (NLESS1 .GT. 1) GO TO 10 - D(1,1) = DEL1 - D(1,N) = DEL1 - GO TO 5000 -C -C NORMAL CASE (N .GE. 3). -C - 10 CONTINUE - H2 = X(3) - X(2) - DEL2 = (F(1,3) - F(1,2))/H2 -C -C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE -C SHAPE-PRESERVING. -C - HSUM = H1 + H2 - W1 = (H1 + HSUM)/HSUM - W2 = -H1/HSUM - D(1,1) = W1*DEL1 + W2*DEL2 - IF ( DPCHST(D(1,1),DEL1) .LE. ZERO) THEN - D(1,1) = ZERO - ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN -C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. - DMAX = THREE*DEL1 - IF (ABS(D(1,1)) .GT. ABS(DMAX)) D(1,1) = DMAX - ENDIF -C -C LOOP THROUGH INTERIOR POINTS. -C - DO 50 I = 2, NLESS1 - IF (I .EQ. 2) GO TO 40 -C - H1 = H2 - H2 = X(I+1) - X(I) - HSUM = H1 + H2 - DEL1 = DEL2 - DEL2 = (F(1,I+1) - F(1,I))/H2 - 40 CONTINUE -C -C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. -C - D(1,I) = ZERO - IF ( DPCHST(DEL1,DEL2) ) 42, 41, 45 -C -C COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY. -C - 41 CONTINUE - IF (DEL2 .EQ. ZERO) GO TO 50 - IF ( DPCHST(DSAVE,DEL2) .LT. ZERO) IERR = IERR + 1 - DSAVE = DEL2 - GO TO 50 -C - 42 CONTINUE - IERR = IERR + 1 - DSAVE = DEL2 - GO TO 50 -C -C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. -C - 45 CONTINUE - HSUMT3 = HSUM+HSUM+HSUM - W1 = (HSUM + H1)/HSUMT3 - W2 = (HSUM + H2)/HSUMT3 - DMAX = MAX( ABS(DEL1), ABS(DEL2) ) - DMIN = MIN( ABS(DEL1), ABS(DEL2) ) - DRAT1 = DEL1/DMAX - DRAT2 = DEL2/DMAX - D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2) -C - 50 CONTINUE -C -C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE -C SHAPE-PRESERVING. -C - W1 = -H2/HSUM - W2 = (H2 + HSUM)/HSUM - D(1,N) = W1*DEL1 + W2*DEL2 - IF ( DPCHST(D(1,N),DEL2) .LE. ZERO) THEN - D(1,N) = ZERO - ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN -C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. - DMAX = THREE*DEL2 - IF (ABS(D(1,N)) .GT. ABS(DMAX)) D(1,N) = DMAX - ENDIF -C -C NORMAL RETURN. -C - 5000 CONTINUE - RETURN -C -C ERROR RETURNS. -C - 5001 CONTINUE -C N.LT.2 RETURN. - IERR = -1 - CALL XERMSG ('SLATEC', 'DPCHIM', - + 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) - RETURN -C - 5002 CONTINUE -C INCFD.LT.1 RETURN. - IERR = -2 - CALL XERMSG ('SLATEC', 'DPCHIM', 'INCREMENT LESS THAN ONE', IERR, - + 1) - RETURN -C - 5003 CONTINUE -C X-ARRAY NOT STRICTLY INCREASING. - IERR = -3 - CALL XERMSG ('SLATEC', 'DPCHIM', - + 'X-ARRAY NOT STRICTLY INCREASING', IERR, 1) - RETURN -C------------- LAST LINE OF DPCHIM FOLLOWS ----------------------------- - END
--- a/main/splines/dpchst.f Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,59 +0,0 @@ -*DECK DPCHST - DOUBLE PRECISION FUNCTION DPCHST (ARG1, ARG2) -C***BEGIN PROLOGUE DPCHST -C***SUBSIDIARY -C***PURPOSE DPCHIP Sign-Testing Routine -C***LIBRARY SLATEC (PCHIP) -C***TYPE DOUBLE PRECISION (PCHST-S, DPCHST-D) -C***AUTHOR Fritsch, F. N., (LLNL) -C***DESCRIPTION -C -C DPCHST: DPCHIP Sign-Testing Routine. -C -C -C Returns: -C -1. if ARG1 and ARG2 are of opposite sign. -C 0. if either argument is zero. -C +1. if ARG1 and ARG2 are of the same sign. -C -C The object is to do this without multiplying ARG1*ARG2, to avoid -C possible over/underflow problems. -C -C Fortran intrinsics used: SIGN. -C -C***SEE ALSO DPCHCE, DPCHCI, DPCHCS, DPCHIM -C***ROUTINES CALLED (NONE) -C***REVISION HISTORY (YYMMDD) -C 811103 DATE WRITTEN -C 820805 Converted to SLATEC library version. -C 870813 Minor cosmetic changes. -C 890411 Added SAVE statements (Vers. 3.2). -C 890531 Changed all specific intrinsics to generic. (WRB) -C 890531 REVISION DATE from Version 3.2 -C 891214 Prologue converted to Version 4.0 format. (BAB) -C 900328 Added TYPE section. (WRB) -C 910408 Updated AUTHOR and DATE WRITTEN sections in prologue. (WRB) -C 930503 Improved purpose. (FNF) -C***END PROLOGUE DPCHST -C -C**End -C -C DECLARE ARGUMENTS. -C - DOUBLE PRECISION ARG1, ARG2 -C -C DECLARE LOCAL VARIABLES. -C - DOUBLE PRECISION ONE, ZERO - SAVE ZERO, ONE - DATA ZERO /0.D0/, ONE/1.D0/ -C -C PERFORM THE TEST. -C -C***FIRST EXECUTABLE STATEMENT DPCHST - DPCHST = SIGN(ONE,ARG1) * SIGN(ONE,ARG2) - IF ((ARG1.EQ.ZERO) .OR. (ARG2.EQ.ZERO)) DPCHST = ZERO -C - RETURN -C------------- LAST LINE OF DPCHST FOLLOWS ----------------------------- - END
--- a/main/splines/pchip.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ -## Copyright (C) 2001,2002 Kai Habel -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; if not, write to the Free Software -## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{pp} = } pchip (@var{x}, @var{y}) -## @deftypefnx {Function File} {@var{yi} = } pchip (@var{x}, @var{y}, @var{xi}) -## piecewise cubic hermite interpolating polynom. -## @var{x} must be a strictly monotonic vector (either increasing or decreasing). -## @var{y} is a vector of same length as @var{x}, -## or a matrix where the number of columns must match the length -## of @var{x}. In this case the interpolating polynoms are calculated -## for each column. -## In contrast to spline, pchip preserves the monotonicity of (@var{x},@var{y}). -## -## @seealso{ppval, spline, csape} -## @end deftypefn - -## Author: Kai Habel <kai.habel@gmx.de> -## Date: 9. mar 2001 -## -## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom) -## -## 4 conditions: -## S_k(x_k) = y_k; -## S_k(x_k+1) = y_k+1; -## S_k'(x_k) = y_k'; -## S_k'(x_k+1) = y_k+1'; - -function ret = pchip (x, y, xi) - - if (nargin < 2 || nargin > 3) - usage ("pchip (x, y, [xi])"); - endif - - x = x(:); - n = length (x); - - if (columns(y) == n) y = y'; endif - - h = diff(x); - if all(h<0) - x = flipud(x); - h = diff(x); - y = flipud(y); - elseif any(h<=0) - error('x must be strictly monotonic') - endif - - if (rows(y) != n) - error("size of x and y must match"); - endif - - [ry,cy] = size (y); - if (cy > 1) - h = kron (diff (x), ones (1, cy)); - endif - - dy = diff (y) ./ h; - - a = 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; - - 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(:)]; - pp = mkpp (x, coeffs); - - if (nargin == 2) - ret = pp; - else - ret = ppval(pp,xi); - endif - -endfunction
--- a/main/splines/pchip_deriv.cc Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,82 +0,0 @@ -/* Copyright (C) 2002 Kai Habel -** -** This program is free software; you can redistribute it and/or modify -** it under the terms of the GNU General Public License as published by -** the Free Software Foundation; either version 2 of the License, or -** (at your option) any later version. -** -** This program is distributed in the hope that it will be useful, -** but WITHOUT ANY WARRANTY; without even the implied warranty of -** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -** GNU General Public License for more details. -** -** You should have received a copy of the GNU General Public License -** along with this program; if not, write to the Free Software -** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -*/ - -#include <octave/oct.h> -#include <octave/f77-fcn.h> - -// SLATEC/PCHIP function not in libcruft -extern "C" { - extern int F77_FUNC (dpchim, DPCHIM) - (const int &n, double *x, double *f, - double *d, const int &incfd, int *ierr); -} - -DEFUN_DLD(pchip_deriv, args, ,"\ -pchip_deriv(x,y):\n\ -wrapper for SLATEC/PCHIP function DPCHIM to calculate derivates for\n\ -piecewise polynomials. You shouldn't use pchip_deriv, use pchip instead.\n\ -") -{ - octave_value_list retval; - 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(); - - if (nx != nyr) - { - error("number of rows for x and y must match"); - return retval; - } - - 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); - } - - 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); - return retval; - } - for (int r=0; r<nx; r++) - { - dmat(r,c) = dvec(r); - } - } - - retval(0) = dmat; - - } - return retval; - -}
--- a/main/strings/mat2str.m Thu Jun 01 22:42:42 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,88 +0,0 @@ -##USAGE s = mat2str( x, n, PLUS ) -## format real/complex numerical matrix x as string s -## suitable for usage by 'eval' -function -## -##n digits of precision (default n=17) -## n(1) : precision of real parts format -## n(2) : precision of imag parts format -##PLUS 1 : print '+' for pos. real parts (0 by def.) -##NOTE * scalar n sets n(2) = n(1) = n -## * for real x any n(2) is ignored -## * may fail for Octave V2.0.X and complex input -##EXA mat2str( [ -1/3 + i/7; 1/3 - i/7 ], [4 2] ) -## |- [-0.3333+0.14i;0.3333-0.14i] -## mat2str( [ -1/3 +i/7; 1/3 -i/7 ], [4 2] ) -## |- [-0.3333+0i,0+0.14i;0.3333+0i,-0-0.14i] -## mat2str( [1-i -1+i],[],1 ) |- [+1-1i,-1+1i] -##HINT better use commas to seperate row-elements of x -##ASSOC sprintf, int2str -##Copyright (C) 2002 Rolf Fabian <fabian@tu-cottbus.de> -## published under current GNU GENERAL PUBLIC LICENSE 020531 - -function s=mat2str(x,n,PLUS) - -if ( nargin<2||isempty(n) ) - n=17; # default precision -endif - -if ( nargin<3||isempty(PLUS) ) - PLUS=0; # def. PLUS : DO NOT print leading '+' - # for positive real elements -endif - -if ( nargin<1||nargin>3||ischar(x)||isstruct(x)||\ - ischar(n)||isstruct(n)||ischar(PLUS)||isstruct(PLUS) ) - usage ("mat2str( NUMERIC x, NUMERIC n, PLUS 0|1 )"); -endif - -if ( !(COMPLEX=is_complex(x)) ) - if ( !PLUS ) - FMT=sprintf("%%.%dg", n(1)); - else - FMT=sprintf("%%+.%dg",n(1)); - endif -else - if ( length(n)==1 ) - n=[n,n]; - endif - if ( !PLUS ) - FMT=sprintf("%%.%dg%%+.%dgi", n(1),n(2)); - else - FMT=sprintf("%%+.%dg%%+.%dgi",n(1),n(2)); - endif -endif - -[nr,nc] = size(x); - -if ( nr*nc==0 ) # empty .. only print brackets - s = "[]"; - -elseif ( nr*nc==1 ) # scalar x .. don't print brackets - if ( !COMPLEX ) - s = sprintf( FMT, x ); - else - s = sprintf( FMT, real(x), imag(x) ); - endif - -else # non-scalar x .. print brackets - FMT=[FMT,',']; - - if ( !COMPLEX ) - - s = sprintf( FMT, x.' ); - - else - - x = x.'; - s = sprintf( FMT, [ real(x(:))'; imag(x(:))' ] ); - - endif - - s=["[", s]; - s(length(s))= "]"; - IND= find(s == ","); - s( IND(nc:nc:length(IND)) )= ";"; - -endif - -endfunction