Mercurial > hg > octave-nkf
changeset 5837:55404f3b0da1
[project @ 2006-06-01 19:05:31 by jwe]
author | jwe |
---|---|
date | Thu, 01 Jun 2006 19:05:32 +0000 |
parents | ed69a3b5b3d0 |
children | 376e02b2ce70 |
files | doc/interpreter/sparse.txi libcruft/ChangeLog libcruft/slatec-fn/dpchim.f libcruft/slatec-fn/dpchst.f liboctave/Array.cc liboctave/ChangeLog scripts/ChangeLog 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 |
diffstat | 25 files changed, 3501 insertions(+), 14 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/interpreter/sparse.txi +++ b/doc/interpreter/sparse.txi @@ -427,8 +427,8 @@ sprank, svds, spaugment) @item Iterative techniques: - @dfn{luinc}, (bicg, bicgstab, cholinc, cgs, gmres, lsqr, minres, - pcg, pcr, qmr, symmlq) + @dfn{luinc}, @dfn{pcg}, @dfn{pcr}, (bicg, bicgstab, cholinc, cgs, + gmres, lsqr, minres, qmr, symmlq) @item Miscellaneous: @dfn{spparms}, @dfn{symbfact}, @dfn{spstats}, @@ -744,7 +744,10 @@ @node Iterative Techniques, Real Life Example, Sparse Linear Algebra, Sparse Matrices @section Iterative Techniques applied to sparse matrices -WRITE ME OR DELETE ME IF THERE ARE NO ITERATIVE SOLVERS IN OCTAVE 3.0 +There are three functions currently to document here, these being +@dfn{luinc}, @dfn{pcg} and @dfn{pcr}. + +WRITE ME. @node Real Life Example, Oct-Files, Iterative Techniques, Sparse Matrices @section Real Life Example of the use of Sparse Matrices @@ -1478,9 +1481,13 @@ @item minres @emph{Not implemented} @item pcg -@emph{Not implemented} +Solves the linear system of equations @code{@var{A} * @var{x} = +@var{b}} by means of the Preconditioned Conjugate Gradient iterative +method. @item pcr -@emph{Not implemented} +Solves the linear system of equations @code{@var{A} * @var{x} = +@var{b}} by means of the Preconditioned Conjugate Residual iterative +method. @item qmr @emph{Not implemented} @item symmlq @@ -1548,6 +1555,10 @@ matrix S * nzmax:: Returns the amount of storage allocated to the sparse matrix SM. +* pcg:: Solves linear system of equations by means of the + Preconditioned Conjugate Gradient iterative method. +* pcr:: Solves linear system of equations by means of the + Preconditioned Conjugate Residual iterative method. * spalloc:: Returns an empty sparse matrix of size R-by-C. * sparse:: SPARSE: create a sparse matrix * spatan2:: Compute atan (Y / X) for corresponding sparse matrix @@ -1672,12 +1683,22 @@ @DOCSTRING(nonzeros) -@node nzmax, spalloc, nonzeros, Function Reference +@node nzmax, pcg, nonzeros, Function Reference @subsubsection nzmax @DOCSTRING(nzmax) -@node spalloc, sparse, nzmax, Function Reference +@node pcg, pcr, nzmax, Function Reference +@subsubsection pcg + +@DOCSTRING(pcg) + +@node pcr, spalloc, pcg, Function Reference +@subsubsection pcr + +@DOCSTRING(pcr) + +@node spalloc, sparse, pcr, Function Reference @subsubsection spalloc @DOCSTRING(spalloc)
--- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,3 +1,7 @@ +2006-06-01 David Bateman <dbateman@free.fr> + + * slatec-fn/dpchim.f, slatec-fn/dpchst.f: New files. + 2006-05-22 John W. Eaton <jwe@octave.org> * lapack/dlantr.f, lapack/zlantr.f: New files.
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/dpchim.f @@ -0,0 +1,285 @@ +*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) .LT. 0.) GO TO 42 + IF ( DPCHST(DEL1,DEL2) .EQ. 0.) GO TO 41 + GO TO 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
new file mode 100644 --- /dev/null +++ b/libcruft/slatec-fn/dpchst.f @@ -0,0 +1,59 @@ +*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/liboctave/Array.cc +++ b/liboctave/Array.cc @@ -3104,14 +3104,15 @@ else final_lhs_dims = new_dims; - lhs_dims = new_dims; - - lhs_dims_len = lhs_dims.length (); - - frozen_len = freeze (idx, lhs_dims, true); + lhs_dims_len = new_dims.length (); + + frozen_len = freeze (idx, new_dims, true); if (rhs_is_scalar) { + if (n_idx < orig_lhs_dims_len) + lhs = lhs.reshape (lhs_dims); + lhs.resize_and_fill (new_dims, rfv); if (! final_lhs_dims.any_zero ()) @@ -3147,6 +3148,9 @@ } else { + if (n_idx < orig_lhs_dims_len) + lhs = lhs.reshape (lhs_dims); + lhs.resize_and_fill (new_dims, rfv); if (! final_lhs_dims.any_zero ())
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,7 @@ +2006-05-31 David Bateman <dbateman@free.fr> + + * Array.cc (assignN): Maybe reshape LHS before doing assignment. + 2006-05-23 John W. Eaton <jwe@octave.org> * oct-types.h.in: Include stdint.h or inttypes.h for integer
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,11 @@ +2006-06-01 David Bateman <dbateman@free.fr> + + * general/interpft.m, general/quadl.m, general/polyarea.m, + general/interp1.m, general/gradient.m, general/interp2.m, + general/bicubic.m, miscellaneous/inputname.m, plot/__plt3__.m, + plot/ndgrid.m, plot/plot3.m, polynomial/pchip.m, sparse/pcg.m, + sparse/pcr.m, strings/mat2str.m: New files from Octave Forge. + 2006-05-31 Bill Denney <bill@givebillmoney.com> * miscellaneous/fileparts.m, miscellaneous/fullfile.m: Add seealso.
new file mode 100644 --- /dev/null +++ b/scripts/general/bicubic.m @@ -0,0 +1,197 @@ +## Copyright (C) 2005 Hoxide Ma +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, 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, spline_alpha) + + if (nargin < 1 || nargin > 6) + print_usage (); + endif + + 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 == 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 || nargin == 6) + [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 + print_usage (); + 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;
new file mode 100644 --- /dev/null +++ b/scripts/general/gradient.m @@ -0,0 +1,132 @@ +## Copyright (C) 2000 Kai Habel +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} = } gradient (@var{M}) +## @deftypefnx {Function File} {[@var{x}, @var{y}, @dots{}] = } gradient (@var{M}) +## @deftypefnx {Function File} {[@dots{}] = } gradient (@var{M}, @var{s}) +## @deftypefnx {Function File} {[@dots{}] = } gradient (@var{M}, @var{dx}, @var{dy}, @dots{}) +## +## Calculates the gradient. @code{@var{x} = gradient (@var{M})} +## calculates the one dimensional gradient if @var{M} is a vector. If +## @var{M} is a matrix the gradient is calculated for each row. +## +## @code{[@var{x}, @var{y}] = gradient (@var{M})} calculates the one +## dimensional gradient for each direction if @var{M} if @var{M} is a +## matrix. Additional return arguments can be use for multi-dimensional +## matrices. +## +## Spacing values between two points can be provided by the +## @var{dx}, @var{dy} or @var{h} parameters. If @var{h} is supplied it +## is assumed to be the spacing in all directions. Otherwise, seperate +## values of the spacing can be supplied by the @var{dx}, etc variables. +## A scalar value specifies an equidistant spacing, while 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> +## Modified: David Bateman <dbateman@free.fr> Added NDArray support + +function [varargout] = gradient (M, varargin) + + if ((nargin < 1)) + print_usage () + endif + + transposed = false; + if (isvector (M)) + ## make a column vector + transposed = (size(M,2) == 1); + M = M(:)'; + endif + + nd = ndims (M); + sz = size (M); + if (nargin > 2 && nargin != nd + 1) + print_usage () + endif + + d = cell(1,nd); + if (nargin == 1) + for i=1:nd + 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); + endfor + else + 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); + else + d{i} = varargin{i}; + endif + endfor + + ## Why the hell did matlab decide to swap these two values? + tmp = d{1}; + d{1} = d{2}; + d{2} = tmp; + endif + + 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)); + + 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); + 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)); + endif + varargout{i} = ipermute (Y, [i:nd,1:i-1]); + M = permute (M, [2:nd,1]); + endfor + + ## Why the hell did matlab decide to swap these two values? + tmp = varargout{1}; + varargout{1} = varargout{2}; + varargout{2} = tmp; + + if (transposed) + varargout{1} = varargout{1}.'; + endif +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/interp1.m @@ -0,0 +1,442 @@ +## Copyright (C) 2000 Paul Kienzle +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{yi} =} interp1 (@var{x}, @var{y}, @var{xi}) +## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, @var{method}) +## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, @var{extrap}) +## @deftypefnx {Function File} {@var{pp} =} interp1 (@dots{}, 'pp') +## +## One-dimensional interpolation. Interpolate @var{y}, defined at the +## points @var{x}, at the points @var{xi}. The sample points @var{x} +## must be strictly monotonic. If @var{y} is an array, treat the columns +## of @var{y} seperately. +## +## Method is one of: +## +## @table @asis +## @item 'nearest' +## Return the nearest neighbour. +## @item 'linear' +## Linear interpolation from nearest neighbours +## @item 'pchip' +## Piece-wise cubic hermite interpolating polynomial +## @item 'cubic' +## Cubic interpolation from four nearest neighbours +## @item 'spline' +## Cubic spline interpolation--smooth first and second derivatives +## throughout the curve +## @end table +## +## Appending '*' to the start of the above method forces @code{interp1} +## to assume that @var{x} is uniformly spaced, and only @code{@var{x} +## (1)} and @code{@var{x} (2)} are referenced. This is usually faster, +## and is never slower. The default method is 'linear'. +## +## If @var{extrap} is the string 'extrap', then extrapolate values beyond +## the endpoints. If @var{extrap} is a number, replace values beyond the +## endpoints with that number. If @var{extrap} is missing, assume NaN. +## +## If the string argument 'pp' is specified, then @var{xi} should not be +## supplied and @code{interp1} returns the piece-wise polynomial that +## can later be used with @code{ppval} to evaluate the interpolation. +## There is an equivalence, such that @code{ppval (interp1 (@var{x}, +## @var{y}, @var{method}, 'pp'), @var{xi}) == interp1 (@var{x}, @var{y}, +## @var{xi}, @var{method}, 'extrap')}. +## +## An example of the use of @code{interp1} is +## +## @example +## @group +## 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,'*;;'); +## @end group +## @end example +## +## @seealso{interpft} +## @end deftypefn + +## 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, varargin) + + if ( nargin < 3 || nargin > 6) + print_usage (); + endif + + method = "linear"; + extrap = NaN; + xi = []; + pp = false; + firstnumeric = true; + + if (nargin > 2) + for i = 1:length(varargin) + arg = varargin{i}; + if (ischar(arg)) + arg = tolower (arg); + if (strcmp("extrap",arg)) + extrap = "extrap"; + elseif (strcmp("pp",arg)) + pp = true; + else + method = arg; + endif + else + if (firstnumeric) + xi = arg; + firstnumeric = false; + else + extrap = arg; + endif + endif + endfor + endif + + ## reshape matrices for convenience + x = x(:); + nx = size(x,1); + if (isvector(y) && size (y, 1) == 1) + y = y (:); + endif + ndy = ndims (y); + szy = size(y); + ny = szy(1); + nc = prod (szy(2:end)); + y = reshape (y, ny, nc); + szx = size(xi); + xi = xi(:); + + ## determine sizes + 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 (!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)) + if (szx(1) == 1) + yi = reshape (yi, [szx(2), szy(2:end)]); + else + yi = reshape (yi, [szx(1), szy(2:end)]); + endif + else + yi = reshape (yi, [szx, szy(2:end)]); + endif + return; + endif + xi = xi(range); + endif + endif + + if strcmp(method, "nearest") + if (pp) + 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; + yi(range,:) = y(idx,:); + endif + elseif strcmp(method, "*nearest") + if (pp) + 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))); + yi(range,:) = y(idx,:); + endif + 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)); + else + ## 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 + s = (xi - x(idx))./dx(idx); + yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); + endif + 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)); + else + ## find the interval containing the test point + t = (xi - minx)/dx + 1; + idx = max(1,min(ny,floor(t))); + + ## use the endpoints of the interval to define a line + dy = [y(2:ny,:) - y(1:ny-1,:); y(ny,:) - y(ny-1,:)]; + s = t - idx; + yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); + endif + elseif strcmp(method, "pchip") || strcmp(method, "*pchip") + if (nx == 2 || method(1) == "*") + x = linspace(minx, maxx, ny); + endif + ## Note that pchip's arguments are transposed relative to interp1 + if (pp) + yi = pchip(x.', y.'); + yi.d = szy(2:end); + else + yi(range,:) = pchip(x.', y.', xi.').'; + endif + + 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).'; + nx = ny; + endif + + 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); + + if (pp) + xs = [x(1);x(3:nx-2)]; + yi = mkpp ([x(1);x(3:nx-2);x(nx)], + [a(:), (b(:) + 3.*xs(:,J).*a(:)), ... + (c(:) + 2.*xs(:,J).*b(:) + 3.*xs(:,J)(:).^2.*a(:)), ... + (d(:) + xs(:,J).*c(:) + xs(:,J).^2.*b(:) + ... + xs(:,J).^3.*a(:))], szy(2:end)); + else + yi(range,:) = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ... + + c(idx,:)).*xi(:,J) + d(idx,:); + endif + 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 || 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.d = szy(2:end); + else + yi(range,:) = spline(x.', y.', xi.').'; + endif + else + error(["interp1 doesn't understand method '", method, "'"]); + endif + + 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 + yi = reshape (yi, [szx(1), szy(2:end)]); + endif + else + yi = reshape (yi, [szx, szy(2:end)]); + endif + 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); + +%!assert (ppval(interp1(xp,yp,"nearest","pp"),xi), +%! interp1(xp,yp,xi,"nearest","extrap"),10*eps); +%!assert (ppval(interp1(xp,yp,"linear","pp"),xi), +%! interp1(xp,yp,xi,"linear","extrap"),10*eps); +%!assert (ppval(interp1(xp,yp,"cubic","pp"),xi), +%! interp1(xp,yp,xi,"cubic","extrap"),10*eps); +%!assert (ppval(interp1(xp,yp,"pchip","pp"),xi), +%! interp1(xp,yp,xi,"pchip","extrap"),10*eps); +%!assert (ppval(interp1(xp,yp,"spline","pp"),xi), +%! interp1(xp,yp,xi,"spline","extrap"),10*eps); + +%!assert (ppval(interp1(xp,yp,"*nearest","pp"),xi), +%! interp1(xp,yp,xi,"*nearest","extrap"),10*eps); +%!assert (ppval(interp1(xp,yp,"*linear","pp"),xi), +%! interp1(xp,yp,xi,"*linear","extrap"),10*eps); +%!assert (ppval(interp1(xp,yp,"*cubic","pp"),xi), +%! interp1(xp,yp,xi,"*cubic","extrap"),10*eps); +%!assert (ppval(interp1(xp,yp,"*pchip","pp"),xi), +%! interp1(xp,yp,xi,"*pchip","extrap"),10*eps); +%!assert (ppval(interp1(xp,yp,"*spline","pp"),xi), +%! interp1(xp,yp,xi,"*spline","extrap"),10*eps);
new file mode 100644 --- /dev/null +++ b/scripts/general/interp2.m @@ -0,0 +1,316 @@ +## Copyright (C) 2000 Kai Habel +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, 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 (@dots{}, @var{method}) +## @deftypefnx {Function File} {@var{zi}=} interp2 (@dots{}, @var{method}, @var{extrapval}) +## +## 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{Yy must be +## monotonic. If they are matrices they must have the @code{meshgrid} +## format. +## +## @table @code +## @item interp2 (@var{x}, @var{y}, @var{Z}, @var{xi}, @var{yi}, @dots{}) +## 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. +## +## @item interp2 (@var{z}, @var{xi}, @var{yi}) +## Assumes @code{@var{x} = 1:rows (@var{z})} and @code{@var{y} = +## 1:columns (@var{z})} +## +## @item interp2 (@var{z}, @var{n}) +## Interleaves the Matrix @var{z} n-times. If @var{n} is ommited a value +## of @code{@var{n} = 1} is assumed. +## @end table +## +## The variable @var{method} defines the method to use for the +## interpolation. It can take one of the values +## +## @table @asis +## @item 'nearest' +## Return the nearest neighbour. +## @item 'linear' +## Linear interpolation from nearest neighbours +## @item 'pchip' +## Piece-wise cubic hermite interpolating polynomial +## @item 'cubic' +## Cubic interpolation from four nearest neighbours +## @item 'spline' +## Cubic spline interpolation--smooth first and second derivatives +## throughout the curve (Not implemented yet). +## @end table +## +## If a scalar value @var{extrapval} is defined as the final value, then +## values outside the mesh as set to this value. Note that in this case +## @var{method} must be defined as well. If @var{extrapval} is not +## defined then NaN 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"; + 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 (); + 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 + if (!isscalar(extrapval)) + error("interp2 expected n extrapval"); + 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)); + + elseif (strcmp (method, "spline")) + ## FIXME Implement 2-D (or in fact ND) spline interpolation + error ("interp2, spline interpolation not yet implemented"); + + else + error ("interpolation method not recognized"); + 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); + +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]); +%! assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]); + +%!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); +
new file mode 100644 --- /dev/null +++ b/scripts/general/interpft.m @@ -0,0 +1,116 @@ +## Copyright (C) 2001 Paul Kienzle +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} interpft (@var{x}, @var{n}) +## @deftypefnx {Function File} {} interpft (@var{x}, @var{n}, @var{dim}) +## +## Fourier interpolation. If @var{x} is a vector, then @var{x} is +## resampled with @var{n} points. The data in @var{x} is assumed to be +## equispaced. If @var{x} is an array, then operate along each column of +## the array seperately. If @var{dim} is specified, then interpolate +## along the dimension @var{dim}. +## +## @code{interpft} assumes that the interpolated function is periodic, +## and so assumption are made about the end points of the inetrpolation. +## +## @seealso{interp1} +## @end deftypefn + +## Author: Paul Kienzle +## 2001-02-11 +## * initial version +## 2002-03-17 aadler +## * added code to work on matrices as well +## 2006-05-25 dbateman +## * Make it matlab compatiable, cutting out the 2-D interpolation + +function z = interpft (x, n, dim) + + if (nargin < 2 || nargin > 3) + print_usage (); + endif + + if (nargin == 2) + if (isvector(x) && size(x,1) == 1) + dim = 2; + else + dim = 1; + endif + endif + + if (!isscalar (n)) + error ("interpft: n must be an integer scalar"); + endif + + 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); + m = size (x, 1); + + inc = 1; + while (inc * n < m) + inc++; + endwhile + y = fft (x) / m; + k = floor (m / 2); + sz = size(x); + sz(1) = n * inc - m; + 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)); + 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); + endif + + z = ipermute (z, perm); +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;'); + +%!shared n,y +%! x = [0:10]'; y = sin(x); n = length (x); +%!assert (interpft(y, n), y, eps); +%!assert (interpft(y', n), y', eps); +%!assert (interpft([y,y],n), [y,y], eps); + +%!error (interpft(y,n,0)) +%!error (interpft(y,[n,n]))
new file mode 100644 --- /dev/null +++ b/scripts/general/polyarea.m @@ -0,0 +1,71 @@ +## Copyright (C) 1999 David M. Doolin <doolin@ce.berkeley.edu> +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} polyarea (@var{x}, @var{y}) +## @deftypefnx {Function File} {} polyarea (@var{x}, @var{y}, @var{dim}) +## +## Determines area of a polygon by triangle method. The variables +## @var{x} and @var{y} define the vertex pairs, and must therefore have +## the same shape. Then might be either vectors or arrays. If they are +## arrays then the columns of @var{x} and @var{y} are treated seperately +## and an area returned for each. +## +## If the optional @var{dim} argument is given, then @code{polyarea} +## works along this dimension of the arrays @var{x} and @var{y}. +## +## @end deftypefn + +## 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, 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 + 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 + endif +endfunction + +%!shared x, y +%! x = [1;1;3;3;1]; +%! y = [1;3;3;1;1]; +%!assert (polyarea(x,y), 4, eps) +%!assert (polyarea([x,x],[y,y]), [4,4], eps) +%!assert (polyarea([x,x],[y,y],1), [4,4], eps) +%!assert (polyarea([x,x]',[y,y]',2), [4;4], eps)
new file mode 100644 --- /dev/null +++ b/scripts/general/quadl.m @@ -0,0 +1,173 @@ +## Copyright (C) 1998 Walter Gautschi +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}) +## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}) +## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}) +## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{}) +## +## Numerically evaluate integral using adaptive Lobatto rule. +## @code{quadl (@var{f}, @var{a}, @var{b})} approximates the integral of +## @code{@var{f}(@var{x})} to machine precision. @var{f} is either a +## function handle, inline function or string containing the name of +## the function to evaluate. The function @var{f} must return a vector +## of output values if given a vector of input values. +## +## If defined, @var{tol} defines the relative tolerance to which to +## which to integrate @code{@var{f}(@var{x})}. While if @var{trace} is +## defined, displays the left end point of the current interval, the +## interval length, and the partial integral. +## +## Additional arguments @var{p1}, etc, are passed directly to @var{f}. +## To use default values for @var{tol} and @var{trace}, one may pass +## empty matrices. +## +## Reference: W. Gander and W. Gautschi, 'Adaptive Quadrature - +## Revisited', BIT Vol. 40, No. 1, March 2000, pp. 84--101. +## @url{http://www.inf.ethz.ch/personal/gander/} +## +## @end deftypefn + +## 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=[]; + endif + if (nargin < 5) + trace = []; + endif + if (isempty(tol)) + tol = eps; + endif + if (isempty(trace)) + trace=0; + endif + if (tol < eps) + tol = eps; + endif + + 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; + endif + erri1 = abs(i1-is); + erri2 = abs(i2-is); + R = 1; + if (erri2 != 0) + R = erri1/erri2; + endif + if (R > 0 && R < 1) + 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{:}); +endfunction + +## ADAPTLOBSTP Recursive function used by QUADL. +## +## 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. +## +## Walter Gautschi, 08/03/98 + +function Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin) + 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); + endif + Q = i1; + if (trace) + 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{:}); + endif +endfunction + +function r = need_warning(v) + persistent w = []; + if (nargin == 0) + r = w; + else + w = v; + endif +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/miscellaneous/inputname.m @@ -0,0 +1,31 @@ +## Copyright (C) 2004 Paul Kienzle +## +## This file is part of Octave. +## +## This program is free software and is in the public domain + +## -*- texinfo -*- +## @deftypefn {Function File} {} inputname (@var{n}) +## +## Return the text defining n-th input to the function. +## +## @end deftypefn + +function s = inputname (n) + s = evalin ('caller', sprintf('deblank(argn(%d,:));', n)); +endfunction + +## 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');
new file mode 100644 --- /dev/null +++ b/scripts/plot/__plt3__.m @@ -0,0 +1,65 @@ +## Copyright (C) 1996 John W. Eaton +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} __plt3__ (@var{x}, @var{y}, @var{z}, @var{fmt}) +## @end deftypefn + +## 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 +## Modified to use NaN as seperator for gnuplot, so multiple calls +## aren't needed. +## David Bateman <dbateman@free.fr> +## May 25, 2006 + +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"); + + 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); + unwind_protect_cleanup + __gnuplot_set__ noparametric; + end_unwind_protect +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/ndgrid.m @@ -0,0 +1,74 @@ +## Copyright (C) 2006, Alexander Barth +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{y1}, @var{y2}, ..., @var{y}n]} = ndgrid (@var{x1}, @var{x2}, ..., @var{x}n) +## @deftypefnx {Function File} {[@var{y1}, @var{y2}, ..., @var{y}n]} = ndgrid (@var{x}) +## Given n vectors @var{x1}, ... @var{x}n, ndgrid returns n arrays of dimension n. +## The elements of the ith output argument contains the elements of the vector +## @var{x}i repeated over all dimensions different from the ith dimension. +## Calling ndgrid with only one input argument @var{x} is equivalent of calling ndgrid with +## all n input arguments equal to @var{x}: +## +## [@var{y1}, @var{y2}, ..., @var{y}n] = ndgrid (@var{x}, ..., @var{x}) +## @seealso{meshgrid} +## @end deftypefn + +## Author: Alexander Barth <abarth@marine.usf.edu> + +function varargout = ndgrid(varargin) + + if (nargin == 1) + 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]); + else + error ("ndgrid: wrong number of input arguments"); + endif + + + ## Determine the size of the output arguments + + shape = zeros(1,n); + + for i=1:n + if (~isvector(varargin{i})) + error ("ndgrid: arguments must be vectors"); + endif + + shape(i) = length(varargin{i}); + endfor + + for i=1:n + ## size for reshape + r = ones(1,n); + r(i) = shape(i); + + ## size for repmat + s = shape; + s(i) = 1; + + varargout{i} = repmat(reshape(varargin{i},r) ,s); + endfor + + + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/plot3.m @@ -0,0 +1,199 @@ +## Copyright (C) 1996 John W. Eaton +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} plot (@var{args}) +## +## This function produces three-dimensional plots. Many different +## combinations of arguments are possible. The simplest form is +## +## @example +## plot3 (@var{x}, @var{y}, @var{z}) +## @end example +## +## @noindent +## where the arguments are taken to be the vertices of the points to be +## plotted in three dimensions. If all arguments are vectors of the same +## length, then a single continuous line is drawn. If all arguments are +## matrices, then each column of the matrices is treated as a seperate +## line. No attempt is made to transpose the arguments to make the +## number of rows match. +## +## To save a plot, in one of several image formats such as PostScript +## or PNG, use the @code{print} command. +## +## An optional format argument can be given as +## +## @example +## plot3 (@var{x}, @var{y}, @var{y}, @var{fmt}) +## @end example +## +## If the @var{fmt} argument is supplied, it is interpreted as +## follows. If @var{fmt} is missing, the default gnuplot line style +## is assumed. +## +## @table @samp +## @item - +## Set lines plot style (default). +## +## @item . +## Set dots plot style. +## +## @item @@ +## Set points plot style. +## +## @item -@@ +## Set linespoints plot style. +## +## @item ^ +## Set impulses plot style. +## +## @item L +## Set steps plot style. +## +## @item @var{n} +## Interpreted as the plot color if @var{n} is an integer in the range 1 to +## 6. +## +## @item @var{nm} +## If @var{nm} is a two digit integer and @var{m} is an integer in the +## range 1 to 6, @var{m} is interpreted as the point style. This is only +## valid in combination with the @code{@@} or @code{-@@} specifiers. +## +## @item @var{c} +## If @var{c} is one of @code{"k"}, @code{"r"}, @code{"g"}, @code{"b"}, +## @code{"m"}, @code{"c"}, or @code{"w"}, it is interpreted as the plot +## color (black, red, green, blue, magenta, cyan, or white). +## +## @item ";title;" +## Here @code{"title"} is the label for the key. +## +## @item + +## @itemx * +## @itemx o +## @itemx x +## Used in combination with the points or linespoints styles, set the point +## style. +## @end table +## +## The color line styles have the following meanings on terminals that +## support color. +## +## @example +## Number Gnuplot colors (lines)points style +## 1 red * +## 2 green + +## 3 blue o +## 4 magenta x +## 5 cyan house +## 6 brown there exists +## @end example +## +## The @var{fmt} argument can also be used to assign key titles. +## To do so, include the desired title between semi-colons after the +## formatting sequence described above, e.g. "+3;Key Title;" +## Note that the last semi-colon is required and will generate an error if +## it is left out. +## +## Arguments can also be given in groups of three as +## +## @example +## plot3 (@var{x1}, @var{y1}, @var{y1}, @var{x2}, @var{y2}, @var{y2}, @dots{}) +## @end example +## +## @noindent +## where each set of three arguments are treated as seperate lines or +## sets of lines in three dimensions. +## +## An example of the use of plot3 is +## +## @example +## @group +## z = [0:0.05:5]; +## plot3(cos(2*pi*z), sin(2*pi*z), z, ";helix;"); +## @end group +## @end example +## +## @seealso{plot, semilogx, semilogy, loglog, polar, mesh, contour, __pltopt__ +## bar, stairs, errorbar, replot, xlabel, ylabel, title, print} +## @end deftypefn + +## Author: Paul Kienzle +## (modified from __plt__.m) + +function plot3(varargin) + + 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 = 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 + + 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
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/pchip.m @@ -0,0 +1,154 @@ +## 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 polynomial. Called with two +## arguments, the piece-wise polynomial @var{pp} is returned, that may +## later be used with @code{ppval} to evaluate the polynomial at +## specific points. +## +## The variable @var{x} must be a strictly monotonic vector (either +## increasing or decreasing). While @var{y} can be either a vector or +## array. In the case where @var{y} is a vector, it must have a length +## of @var{n}. If @var{y} is an array, then the size of @var{y} must +## have the form +## @iftex +## @tex +## $$[s_1, s_2, \cdots, s_k, n]$$ +## @end tex +## @end iftex +## @ifinfo +## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]} +## @end ifinfo +## The array is then reshaped internally to a matrix where to leading +## dimension is given by +## @iftex +## @tex +## $$s_1 s_2 \cdots s_k$$ +## @end tex +## @end iftex +## @ifinfo +## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}} +## @end ifinfo +## and each row this matrix is then treated seperately. Note that this +## is exactly the opposite treatment than @code{interp1} and is done +## for compatiability. +## +## Called with a third input argument, @code{pchip} evaluates the +## piece-wise polynomial at the points @var{xi}. There is an equivalence +## between @code{ppval (pchip (@var{x}, @var{y}), @var{xi})} and +## @code{pchip (@var{x}, @var{y}, @var{xi})}. +## +## @seealso{spline, ppval, mkpp, unmkpp} +## @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) + print_usage (); + endif + + x = x(:); + n = length (x); + + ## Check the size and shape of y + ndy = ndims (y); + szy = size (y); + if (ndy == 2 && (szy(1) == 1 || szy(2) == 1)) + if (szy(1) == 1) + y = y'; + else + szy = fliplr (szy); + endif + else + 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)) + error("pchip: x must be strictly monotonic") + endif + + if (rows(y) != n) + error("pchip: 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, szy(1:end-1)); + + if (nargin == 2) + ret = pp; + else + ret = ppval (pp, xi); + endif + +endfunction + +%!demo +%! x = 0:8; +%! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; +%! xi = 0:0.01:8; +%! yspline = spline(x,y,xi); +%! ypchip = pchip(x,y,xi); +%! title("pchip and spline fit to discontinuous function"); +%! plot(xi,yspline,";spline;",... +%! xi,ypchip,"-;pchip;",... +%! x,y,"+;data;"); +%! %------------------------------------------------------------------- +%! % confirm that pchip agreed better to discontinuous data than spline + +%!shared x,y +%! x = 0:8; +%! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; +%!assert (pchip(x,y,x), y); +%!assert (pchip(x,y,x'), y'); +%!assert (pchip(x',y',x'), y'); +%!assert (pchip(x',y',x), y); +%!assert (isempty(pchip(x',y',[]))); +%!assert (isempty(pchip(x,y,[]))); +%!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
--- a/scripts/polynomial/spline.m +++ b/scripts/polynomial/spline.m @@ -23,12 +23,12 @@ ## @deftypefnx {Function File} {@var{yi} = } spline (@var{x}, @var{y}, @var{xi}) ## ## Returns the cubic spline interpolation of @var{y} at the point -## @var{x}. Called with two arguments the piece-wse polynomial @var{pp} +## @var{x}. Called with two arguments the piece-wise polynomial @var{pp} ## that may later be used with @code{ppval} to evaluate the polynomial ## at specific points. ## ## The variable @var{x} must be a vector of length @var{n}, and @var{y} -## can be either a vector of array. In the case where @var{y} is a +## can be either a vector or array. In the case where @var{y} is a ## vector, it can have a length of either @var{n} or @code{@var{n} + 2}. ## If the length of @var{y} is @var{n}, then the 'not-a-knot' end ## condition is used. If the length of @var{y} is @code{@var{n} + 2},
new file mode 100644 --- /dev/null +++ b/scripts/sparse/pcg.m @@ -0,0 +1,482 @@ +## Copyright (C) 2004 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl> +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @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} = +## @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 +## 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 +## message and the @var{flag} output parameter will be set. +## +## @item +## @var{b} is the right hand side vector. +## +## @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} * +## @var{x0})}. If @var{tol} is empty or is omitted, the function sets +## @code{@var{tol} = 1e-6} by default. +## +## @item +## @var{maxit} is the maximum allowable number of iterations; if +## @code{[]} is supplied for @code{maxit}, or @code{pcg} has less +## arguments, a default value equal to 20 is used. +## +## @item +## @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}}. +## 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 +## preferred way of using the preconditioner). If @code{[]} is supplied +## 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 +## function sets @var{x0} to a zero vector by default. +## @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}) +## 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}}. +## +## @item +## @var{flag} reports on the convergence. @code{@var{flag} = 0} means +## the solution converged and the tolerance criterion given by @var{tol} +## is satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit +## for the iteration count was reached. @code{@var{flag} = 3} reports that +## the (preconditioned) matrix was found not positive definite. +## +## @item +## @var{relres} is the ratio of the final residual to its initial value, +## measured in the Euclidean norm. +## +## @item +## @var{iter} is the actual number of iterations performed. +## +## @item +## @var{resvec} describes the convergence history of the method. +## @code{@var{resvec} (i,1)} is the Euclidean norm of the residual, and +## @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{@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 +## particular, if no preconditioning is used, the extimates for the +## 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 +## verifying this assumption. +## @end itemize +## +## Let us consider a trivial problem with a diagonal matrix (we exploit the +## sparsity of A) +## +## @example +## @group +## N = 10; +## A = diag([1:N]); A = sparse(A); +## b = rand(N,1); +## @end group +## @end example +## +## @sc{Example 1:} Simplest use of @code{pcg} +## +## @example +## x = pcg(A,b) +## @end example +## +## @sc{Example 2:} @code{pcg} with a function which computes +## @code{@var{A} * @var{x}} +## +## @example +## @group +## function y = applyA(x) +## y = [1:N]'.*x; +## endfunction +## +## x = pcg('applyA',b) +## @end group +## @end example +## +## @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 +## +## @example +## @group +## function y = applyM(x) +## K = floor(length(x)-2); +## y = x; +## y(1:K) = x(1:K)./[1:K]'; +## endfunction +## +## [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],[],'applyM') +## semilogy([1:iter+1], resvec); +## @end group +## @end example +## +## @sc{Example 4:} Finally, a preconditioner which depends on a +## parameter @var{k}. +## +## @example +## @group +## function y = applyM(x, varargin) +## K = varargin@{1@}; +## y = x; y(1:K) = x(1:K)./[1:K]'; +## endfuntion +## +## [x, flag, relres, iter, resvec, eigest] = ... +## pcg(A,b,[],[],'applyM',[],3) +## @end group +## @end example +## +## @sc{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 +## +## +## @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 +## + +function [x, flag, relres, iter, resvec, eigest] = ... + pcg( A, b, tol, maxit, M, x0, varargin ) + + if ((nargin < 6) || isempty(x0)) + x = zeros(size(b)); + else + x = x0; + endif + + if (nargin < 5) + M = []; + endif + + if ((nargin < 4) || isempty(maxit)) + maxit = min(size(b,1),20); + endif + + maxit = maxit + 2; + + if ((nargin < 3) || isempty(tol)) + tol = 1e-6; + endif + + preconditioned_residual_out = false; + if (nargout > 5) + T = zeros(maxit,maxit); + preconditioned_residual_out = true; + endif + + matrix_positive_definite = true; # assume A is positive definite + + 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{:}); + endif + + 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; + endif + else # then M should be a function! + z = feval(M,r,varargin{:}); + endif + 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{:}); + endif + oldalpha = alpha; # needed only for eigest + 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)) + 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; + endwhile + + 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."); + endif + else + 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 + z = r; + else # otherwise, apply the precond + z = M\r; + endif + else # then M should be a function! + z = feval(M,r,varargin{:}); + endif + resvec(iter-1,2) = sqrt(r'*z); + else + resvec = resvec(:,1); + endif + + 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); + 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 + endif + + if (!matrix_positive_definite) + flag = 3; + if (nargout < 2) + warning("PCG: matrix not positive definite?\n"); + endif + endif +endfunction + +%!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])); +%!
new file mode 100644 --- /dev/null +++ b/scripts/sparse/pcr.m @@ -0,0 +1,437 @@ +## Copyright (C) 2004 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl> +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @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} = +## @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 +## 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 +## message and the @var{flag} output parameter will be set. +## +## @item +## @var{b} is the right hand side vector. +## +## @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} * +## @var{x0})}. If @var{tol} is empty or is omitted, the function sets +## @code{@var{tol} = 1e-6} by default. +## +## @item +## @var{maxit} is the maximum allowable number of iterations; if +## @code{[]} is supplied for @code{maxit}, or @code{pcr} has less +## arguments, a default value equal to 20 is used. +## +## @item +## @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}}. +## 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 +## preferred way of using the preconditioner). If @code{[]} is supplied +## 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 +## function sets @var{x0} to a zero vector by default. +## @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}) +## 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}}. +## +## @item +## @var{flag} reports on the convergence. @code{@var{flag} = 0} means +## the solution converged and the tolerance criterion given by @var{tol} +## is satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit +## for the iteration count was reached. @code{@var{flag} = 3} reports t +## @code{pcr} breakdown, see [1] for details. +## +## @item +## @var{relres} is the ratio of the final residual to its initial value, +## measured in the Euclidean norm. +## +## @item +## @var{iter} is the actual number of iterations performed. +## +## @item +## @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}. +## @end itemize +## +## Let us consider a trivial problem with a diagonal matrix (we exploit the +## sparsity of A) +## +## @example +## @group +## N = 10; +## A = diag([1:N]); A = sparse(A); +## b = rand(N,1); +## @end group +## @end example +## +## @sc{Example 1:} Simplest use of @code{pcr} +## +## @example +## x = pcr(A, b) +## @end example +## +## @sc{Example 2:} @code{pcr} with a function which computes +## @code{@var{A} * @var{x}}. +## +## @example +## @group +## function y = applyA(x) +## y = [1:10]'.*x; +## endfunction +## +## x = pcr('applyA',b) +## @end group +## @end example +## +## @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 +## +## @example +## @group +## function y = applyM(x) +## K = floor(length(x)-2); +## y = x; +## y(1:K) = x(1:K)./[1:K]'; +## endfunction +## +## [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM') +## semilogy([1:iter+1], resvec); +## @end group +## @end example +## +## @sc{Example 4:} Finally, a preconditioner which depends on a +## parameter @var{k}. +## +## @example +## @group +## function y = applyM(x, varargin) +## K = varargin@{1@}; +## y = x; y(1:K) = x(1:K)./[1:K]'; +## endfunction +## +## [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM',[],3) +## @end group +## @end example +## +## @sc{References} +## +## [1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of +## Equations", section 9.5.4; Springer, 1994 +## +## @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 + +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)); + else + x = x0; + endif + + if (nargin < 5) + M = []; + endif + + if ((nargin < 4) || isempty(maxit)) + maxit = 20; + endif + + maxit = maxit + 2; + + if ((nargin < 3) || isempty(tol)) + tol = 1e-6; + endif + + if (nargin < 2) + print_usage(); + endif + + ## 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{:}); + endif + + 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; + endif + else # then M should be a function! + p = feval(M,r,varargin{:}); + endif + + 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{:}); + endif + + resvec(1) = abs(norm(r)); + + ## iteration + while ((resvec(iter-1) > tol*resvec(1)) && (iter < maxit)) + + if (isnumeric(M)) # is M a matrix? + if isempty(M) # if M is empty, use no precond + s = q; + else # otherwise, apply the precond + s = M\q; + endif + else # then M should be a function! + s = feval(M,q,varargin{:}); + endif + b_top = r'*s; + b_bot = q'*s; + + if (b_bot == 0.0) + breakdown = true; + break; + endif + lambda = b_top/b_bot; + + x = 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{:}); + endif + + alpha0 = (t'*s)/b_bot; + alpha1 = (t'*s_old)/b_bot_old; + + p_temp = p; q_temp = q; + p = s - alpha0*p - alpha1*p_old; + q = t - alpha0*q - alpha1*q_old; + + s_old = s; + p_old = p_temp; + q_old = q_temp; + b_bot_old = b_bot; + + + resvec(iter) = abs(norm(r)); + iter = iter + 1; + endwhile + + 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); + 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 + endif + if (breakdown) + flag = 3; + if (nargout < 2) + warning("PCR: breakdown occured.\n"); + warning("System matrix singular or preconditioner indefinite?\n"); + endif + endif + +endfunction + +%!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 +%!
new file mode 100644 --- /dev/null +++ b/scripts/strings/mat2str.m @@ -0,0 +1,113 @@ +## Copyright (C) 2002 Rolf Fabian <fabian@tu-cottbus.de> +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{s} =} mat2str (@var{x}, @var{n}) +## @deftypefnx {Function File} {@var{s} =} mat2str (@dots{}, 'class') +## +## Format real/complex numerial matrices as strings. This function +## returns values that are suitable for the use of the @code{eval} +## function. +## +## The precision of the values is given by @var{n}. If @var{n} is a +## scalar then both real and imaginary parts of the matrix are printed +## to the same precision. Otherwise @code{@var{n} (1)} defines the +## precision of the real part and @code{@var{n} (2)} defines the +## precision of the imaginary part. The default for @var{n} is 17. +## +## If the argument 'class' is given, then the class of @var{x} is +## included in the string in such a way that the eval will result in the +## construction of a matrix of the same class. +## +## @example +## @group +## mat2str( [ -1/3 + i/7; 1/3 - i/7 ], [4 2] ) +## @result{} '[-0.3333+0.14i;0.3333-0.14i]' +## mat2str( [ -1/3 +i/7; 1/3 -i/7 ], [4 2] ) +## @result{} '[-0.3333+0i,0+0.14i;0.3333+0i,-0-0.14i]' +## mat2str( int16([1 -1]), 'class') +## @result{} 'int16([1,-1])' +## @end group +## @end example +## +## @seealso{sprintf, int2str} +## @end deftypefn + +function s = mat2str(x,n,cls) + + if (nargin < 2 || isempty(n)) + n = 17; # default precision + endif + + if (nargin < 3) + if (ischar(n)) + cls = n; + n = 17; + else + cls = ''; + endif + endif + + if (nargin < 1 || nargin > 3 || ischar(x) || isstruct(x) || + ischar(n) || isstruct(n) || isstruct(cls)) + usage ("mat2str"); + endif + + if (ndims (x) > 2) + error ("mat2str: x must be two dimensional"); + endif + + if (!(COMPLEX = is_complex(x))) + FMT = sprintf("%%.%dg", n(1)); + else + if (length(n) == 1 ) + n = [n, n]; + endif + FMT = sprintf("%%.%dg%%+.%dgi", n(1), n(2)); + 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 + + if (strcmp ("class", cls)) + s = [class(x), "(", s, ")"] + endif +endfunction
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,7 @@ +2006-06-01 David Bateman <dbateman@free.fr> + + * DLD-FUNCTIONS/__pchip_deriv__.cc: New file. + 2006-05-31 John W. Eaton <jwe@octave.org> * load-path.h (load_path::set_command_line_path): Make it
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/__pchip_deriv__.cc @@ -0,0 +1,96 @@ +/* + +Copyright (C) 2002 Kai Habel + +This file is part of Octave. + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "oct-obj.h" +#include "utils.h" +#include "f77-fcn.h" + +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, , + "-*- 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\ +@end deftypefn") +{ + octave_value 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 = dmat; + } + + return retval; +}