diff --git a/lib/PDL/Primitive-pchip.c b/lib/PDL/Primitive-pchip.c index 54225896a..e044aba5b 100644 --- a/lib/PDL/Primitive-pchip.c +++ b/lib/PDL/Primitive-pchip.c @@ -24,23 +24,6 @@ double d_sign(doublereal a, doublereal b) #define xermsg_(lib, func, errstr, nerr, ...) \ fprintf(stderr, "%s::%s: %s (err=%ld)\n", lib, func, errstr, nerr) -/* ***PURPOSE DPCHIP Sign-Testing Routine */ -/* Returns: */ -/* -1. if ARG1 and ARG2 are of opposite sign. */ -/* 0. if either argument is zero. */ -/* +1. if ARG1 and ARG2 are of the same sign. */ -/* The object is to do this without multiplying ARG1*ARG2, to avoid */ -/* possible over/underflow problems. */ -/* Fortran intrinsics used: SIGN. */ -/* ***SEE ALSO DPCHCE, DPCHCI, DPCHCS, DPCHIM */ -#define X(ctype, ppsym) \ - static inline ctype pchst_ ## ppsym(ctype arg1, ctype arg2) \ - { \ - return (arg1 == 0. || arg2 == 0.) ? 0. : d_sign(1, arg1) * d_sign(1, arg2); \ - } -X(doublereal, D) -#undef X - void dpchbs(integer n, doublereal *x, doublereal *f, doublereal *d, integer *knotyp, integer *nknots, doublereal *t, doublereal *bcoef, integer *ndim, integer *kord, @@ -549,178 +532,3 @@ doublereal dpchia(integer n, doublereal *x, doublereal *f, doublereal *d, } return value; } - -/* ***PURPOSE Computes divided differences for DPCHCE and DPCHSP */ -/* DPCHDF: DPCHIP Finite Difference Formula */ -/* Uses a divided difference formulation to compute a K-point approx- */ -/* imation to the derivative at X(K) based on the data in X and S. */ -/* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */ -/* derivative approximations. */ -/* ---------------------------------------------------------------------- */ -/* On input: */ -/* K is the order of the desired derivative approximation. */ -/* K must be at least 3 (error return if not). */ -/* X contains the K values of the independent variable. */ -/* X need not be ordered, but the values **MUST** be */ -/* distinct. (Not checked here.) */ -/* S contains the associated slope values: */ -/* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */ -/* (Note that S need only be of length K-1.) */ -/* On return: */ -/* S will be destroyed. */ -/* IERR will be set to -1 if K.LT.2 . */ -/* DPCHDF will be set to the desired derivative approximation if */ -/* IERR=0 or to zero if IERR=-1. */ -/* ---------------------------------------------------------------------- */ -/* ***SEE ALSO DPCHCE, DPCHSP */ -/* ***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- */ -/* Verlag, New York, 1978, pp. 10-16. */ -/* CHECK FOR LEGAL VALUE OF K. */ -#define X(ctype, ppsym) \ - static inline ctype pchdf_ ## ppsym(integer k, ctype *x, ctype *s, integer *ierr) \ - { \ -/* Local variables */ \ - integer i, j; \ - ctype value; \ - if (k < 3) { \ - *ierr = -1; \ - xermsg_("SLATEC", "DPCHDF", "K LESS THAN THREE", *ierr); \ - return 0.; \ - } \ -/* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */ \ - for (j = 2; j < k; ++j) { \ - integer itmp = k - j; \ - for (i = 0; i < itmp; ++i) { \ - s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]); \ - } \ - } \ -/* EVALUATE DERIVATIVE AT X(K). */ \ - value = s[0]; \ - for (i = 1; i < k-1; ++i) { \ - value = s[i] + value * (x[k-1] - x[i]); \ - } \ - *ierr = 0; \ - return value; \ - } -X(doublereal, D) -#undef X - -/* Programming notes: */ -/* 1. The function DPCHST(ARG1,ARG2) is assumed to return zero if */ -/* either argument is zero, +1 if they are of the same sign, and */ -/* -1 if they are of opposite sign. */ -/* 2. To produce a single precision version, simply: */ -/* a. Change DPCHIM to PCHIM wherever it occurs, */ -/* b. Change DPCHST to PCHST wherever it occurs, */ -/* c. Change all references to the Fortran intrinsics to their */ -/* single precision equivalents, */ -/* d. Change the double precision declarations to real, and */ -/* e. Change the constants ZERO and THREE to single precision. */ -void dpchim(integer n, doublereal *x, doublereal *f, - doublereal *d, integer *ierr) -{ -/* System generated locals */ - doublereal dtmp; -/* Local variables */ - integer i; - doublereal h1, h2, w1, w2, del1, del2, dmin, dmax, hsum, drat1, - drat2, dsave; - integer nless1; - doublereal hsumt3; -/* VALIDITY-CHECK ARGUMENTS. */ - if (n < 2) { - *ierr = -1; - xermsg_("SLATEC", "DPCHIM", "NUMBER OF DATA POINTS LESS THAN TWO", *ierr); - return; - } - for (i = 1; i < n; ++i) { - if (x[i] <= x[i - 1]) { - *ierr = -3; - xermsg_("SLATEC", "DPCHIM", "X-ARRAY NOT STRICTLY INCREASING", *ierr); - return; - } - } -/* FUNCTION DEFINITION IS OK, GO ON. */ - *ierr = 0; - nless1 = n - 1; - h1 = x[1] - x[0]; - del1 = (f[1] - f[0]) / h1; - dsave = del1; -/* SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. */ - if (nless1 <= 1) { - d[0] = d[n-1] = del1; - return; - } -/* NORMAL CASE (N .GE. 3). */ - h2 = x[2] - x[1]; - del2 = (f[2] - f[1]) / h2; -/* SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */ -/* SHAPE-PRESERVING. */ - hsum = h1 + h2; - w1 = (h1 + hsum) / hsum; - w2 = -h1 / hsum; - d[0] = w1 * del1 + w2 * del2; - if (pchst_D(d[0], del1) <= 0.) { - d[0] = 0.; - } else if (pchst_D(del1, del2) < 0.) { -/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */ - dmax = 3. * del1; - if (abs(d[0]) > abs(dmax)) { - d[0] = dmax; - } - } -/* LOOP THROUGH INTERIOR POINTS. */ - for (i = 1; i < nless1; ++i) { - if (i != 1) { - h1 = h2; - h2 = x[i+1] - x[i]; - hsum = h1 + h2; - del1 = del2; - del2 = (f[i+1] - f[i]) / h2; - } -/* SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. */ - d[i] = 0.; - dtmp = pchst_D(del1, del2); - if (dtmp <= 0) { - if (dtmp == 0.) { -/* COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY. */ - if (del2 == 0.) { - continue; - } - if (pchst_D(dsave, del2) < 0.) { - ++(*ierr); - } - dsave = del2; - continue; - } - ++(*ierr); - dsave = del2; - continue; - } -/* USE BRODLIE MODIFICATION OF BUTLAND FORMULA. */ - hsumt3 = hsum + hsum + hsum; - w1 = (hsum + h1) / hsumt3; - w2 = (hsum + h2) / hsumt3; -/* Computing MAX */ - dmax = max(abs(del1),abs(del2)); -/* Computing MIN */ - dmin = min(abs(del1),abs(del2)); - drat1 = del1 / dmax; - drat2 = del2 / dmax; - d[i] = dmin / (w1 * drat1 + w2 * drat2); - } -/* SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */ -/* SHAPE-PRESERVING. */ - w1 = -h2 / hsum; - w2 = (h2 + hsum) / hsum; - d[n-1] = w1 * del1 + w2 * del2; - if (pchst_D(d[n-1], del2) <= 0.) { - d[n-1] = 0.; - } else if (pchst_D(del1, del2) < 0.) { -/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */ - dmax = 3. * del2; - if (abs(d[n-1]) > abs(dmax)) { - d[n-1] = dmax; - } - } -} diff --git a/lib/PDL/Primitive.pd b/lib/PDL/Primitive.pd index 0cb247cf5..d233fd811 100644 --- a/lib/PDL/Primitive.pd +++ b/lib/PDL/Primitive.pd @@ -4314,15 +4314,110 @@ PCHST => sub {my ($a, $b) = @_; }, ); -defslatec('pchip_chim', - {D => 'dpchim'}, - 'IntVal n=$SIZE(n); - Mat x (n); - Mat f (n); - Mat [o] d (n); - Integer [o] ierr (); - ', -'Calculate the derivatives needed to determine a monotone piecewise +pp_def('pchip_chim', + Pars => 'x(n); f(n); [o]d(n); longlong [o]ierr();', + RedoDimsCode => <<'EOF', +if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO"); +EOF + Code => pp_line_numbers(__LINE__, <<'EOF'), +/* Local variables */ +$GENERIC() dmax, hsumt3; +PDL_Indx n = $SIZE(n); +/* VALIDITY-CHECK ARGUMENTS. */ +loop (n=1) %{ + if ($x() > $x(n=>n-1)) continue; + $ierr() = -1; + $CROAK("X-ARRAY NOT STRICTLY INCREASING"); +%} +/* FUNCTION DEFINITION IS OK, GO ON. */ +$ierr() = 0; +$GENERIC() h1 = $x(n=>1) - $x(n=>0); +$GENERIC() del1 = ($f(n=>1) - $f(n=>0)) / h1; +$GENERIC() dsave = del1; +/* SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. */ +if (n <= 2) { + $d(n=>0) = $d(n=>1) = del1; + continue; +} +/* NORMAL CASE (N .GE. 3). */ +$GENERIC() h2 = $x(n=>2) - $x(n=>1); +$GENERIC() del2 = ($f(n=>2) - $f(n=>1)) / h2; +/* SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */ +/* SHAPE-PRESERVING. */ +$GENERIC() hsum = h1 + h2; +$GENERIC() w1 = (h1 + hsum) / hsum; +$GENERIC() w2 = -h1 / hsum; +$d(n=>0) = w1 * del1 + w2 * del2; +if ($PCHST($d(n=>0), del1) <= 0.) { + $d(n=>0) = 0.; +} else if ($PCHST(del1, del2) < 0.) { +/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */ + dmax = 3. * del1; + if (PDL_ABS($d(n=>0)) > PDL_ABS(dmax)) { + $d(n=>0) = dmax; + } +} +/* LOOP THROUGH INTERIOR POINTS. */ +loop (n=1:-1) %{ + if (n != 1) { + h1 = h2; + h2 = $x(n=>n+1) - $x(); + hsum = h1 + h2; + del1 = del2; + del2 = ($f(n=>n+1) - $f()) / h2; + } +/* SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. */ + $d() = 0.; + $GENERIC() dtmp = $PCHST(del1, del2); + if (dtmp <= 0) { + if (dtmp == 0.) { +/* COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY. */ + if (del2 == 0.) { + continue; + } + if ($PCHST(dsave, del2) < 0.) { + ++($ierr()); + } + dsave = del2; + continue; + } + ++($ierr()); + dsave = del2; + continue; + } +/* USE BRODLIE MODIFICATION OF BUTLAND FORMULA. */ + hsumt3 = hsum + hsum + hsum; + w1 = (hsum + h1) / hsumt3; + w2 = (hsum + h2) / hsumt3; +/* Computing MAX */ + dmax = PDLMAX(PDL_ABS(del1),PDL_ABS(del2)); +/* Computing MIN */ + $GENERIC() dmin = PDLMIN(PDL_ABS(del1),PDL_ABS(del2)); + $GENERIC() drat1 = del1 / dmax; + $GENERIC() drat2 = del2 / dmax; + $d() = dmin / (w1 * drat1 + w2 * drat2); +%} +/* SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */ +/* SHAPE-PRESERVING. */ +w1 = -h2 / hsum; +w2 = (h2 + hsum) / hsum; +$d(n=>n-1) = w1 * del1 + w2 * del2; +if ($PCHST($d(n=>n-1), del2) <= 0.) { + $d(n=>n-1) = 0.; +} else if ($PCHST(del1, del2) < 0.) { +/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */ + dmax = 3. * del2; + if (PDL_ABS($d(n=>n-1)) > PDL_ABS(dmax)) { + $d(n=>n-1) = dmax; + } +} +EOF + Doc => <<'EOF', +=for ref + +Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation. + +Calculate the derivatives needed to determine a monotone piecewise cubic Hermite interpolant to the given set of points (C<$x,$f>, where C<$x> is strictly increasing). The resulting set of points - C<$x,$f,$d>, referred to @@ -4383,8 +4478,7 @@ Journal on Scientific and Statistical Computing 5, 2 F. N. Fritsch and R. E. Carlson, Monotone piecewise cubic interpolation, SIAM Journal on Numerical Analysis 17, 2 (April 1980), pp. 238-246. -', -'Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation.' +EOF ); pp_def('pchip_chic',