From e3ed1dee16005a6742c28e5d44b1a67ba08650b5 Mon Sep 17 00:00:00 2001 From: Ed J Date: Mon, 16 Dec 2024 01:19:13 +0000 Subject: [PATCH] make chsp into PP kernel --- lib/PDL/Primitive-pchip.c | 223 -------------------------------------- lib/PDL/Primitive.pd | 210 ++++++++++++++++++++++++++++++++--- 2 files changed, 195 insertions(+), 238 deletions(-) diff --git a/lib/PDL/Primitive-pchip.c b/lib/PDL/Primitive-pchip.c index 47e9cdeeb..54225896a 100644 --- a/lib/PDL/Primitive-pchip.c +++ b/lib/PDL/Primitive-pchip.c @@ -24,10 +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) -doublereal d1mach() { - return 2.3e-16; /* for float, 1.2e-7 */ -} - /* ***PURPOSE DPCHIP Sign-Testing Routine */ /* Returns: */ /* -1. if ARG1 and ARG2 are of opposite sign. */ @@ -728,222 +724,3 @@ void dpchim(integer n, doublereal *x, doublereal *f, } } } - -/* SINGULAR SYSTEM. */ -/* *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES *** */ -/* *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). *** */ -#define dpchsp_singular \ - *ierr = -8; \ - xermsg_("SLATEC", "DPCHSP", "SINGULAR LINEAR SYSTEM", *ierr); \ - return; -void dpchsp(integer *ic, doublereal *vc, integer n, - doublereal *x, doublereal *f, doublereal *d, - doublereal *wk, integer nwk, integer *ierr) -{ -/* System generated locals */ - doublereal dtmp; -/* Local variables */ - doublereal g; - integer j, nm1, ibeg, iend; - doublereal stemp[3], xtemp[4]; - wk -= 3; -/* VALIDITY-CHECK ARGUMENTS. */ - if (n < 2) { - *ierr = -1; - xermsg_("SLATEC", "DPCHSP", "NUMBER OF DATA POINTS LESS THAN TWO", *ierr); - return; - } - for (j = 1; j < n; ++j) { - if (x[j] <= x[j - 1]) { - *ierr = -3; - xermsg_("SLATEC", "DPCHSP", "X-ARRAY NOT STRICTLY INCREASING", *ierr); - return; - } - } - ibeg = ic[0]; - iend = ic[1]; - *ierr = 0; - if (ibeg < 0 || ibeg > 4) { - --(*ierr); - } - if (iend < 0 || iend > 4) { - *ierr += -2; - } - if (*ierr < 0) { - *ierr += -3; - xermsg_("SLATEC", "DPCHSP", "IC OUT OF RANGE", *ierr); - return; - } -/* FUNCTION DEFINITION IS OK -- GO ON. */ - if (nwk < n << 1) { - *ierr = -7; - xermsg_("SLATEC", "DPCHSP", "WORK ARRAY TOO SMALL", *ierr); - return; - } -/* COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, */ -/* COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). */ - for (j = 1; j < n; ++j) { - wk[(j << 1) + 3] = x[j] - x[j - 1]; - wk[(j << 1) + 4] = (f[j] - f[j - 1]) / wk[(j << 1) + 3]; - } -/* SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. */ - if (ibeg > n) { - ibeg = 0; - } - if (iend > n) { - iend = 0; - } -/* SET UP FOR BOUNDARY CONDITIONS. */ - if (ibeg == 1 || ibeg == 2) { - d[0] = vc[0]; - } else if (ibeg > 2) { -/* PICK UP FIRST IBEG POINTS, IN REVERSE ORDER. */ - for (j = 0; j < ibeg; ++j) { - integer index = ibeg - j; -/* INDEX RUNS FROM IBEG DOWN TO 1. */ - xtemp[j] = x[index+1]; - if (j < ibeg-1) { - stemp[j] = wk[(index << 1) + 6]; - } - } -/* -------------------------------- */ - d[0] = pchdf_D(ibeg, xtemp, stemp, ierr); -/* -------------------------------- */ - if (*ierr != 0) { - *ierr = -9; - xermsg_("SLATEC", "DPCHSP", "ERROR RETURN FROM DPCHDF", *ierr); - return; - } - ibeg = 1; - } - if (iend == 1 || iend == 2) { - d[n-1] = vc[1]; - } else if (iend > 2) { -/* PICK UP LAST IEND POINTS. */ - for (j = 0; j < iend; ++j) { - integer index = n - iend + j; -/* INDEX RUNS FROM N+1-IEND UP TO N. */ - xtemp[j] = x[index]; - if (j < iend-1) { - stemp[j] = wk[((index + 2) << 1) + 2]; - } - } -/* -------------------------------- */ - d[n-1] = pchdf_D(iend, xtemp, stemp, ierr); -/* -------------------------------- */ - if (*ierr != 0) { - *ierr = -9; - xermsg_("SLATEC", "DPCHSP", "ERROR RETURN FROM DPCHDF", *ierr); - return; - } - iend = 1; - } -/* --------------------( BEGIN CODING FROM CUBSPL )-------------------- */ -/* **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF */ -/* F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- */ -/* INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. */ -/* WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. */ -/* CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM */ -/* WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) */ - if (ibeg == 0) { - if (n == 2) { -/* NO CONDITION AT LEFT END AND N = 2. */ - wk[4] = 1.; - wk[3] = 1.; - d[0] = 2. * wk[6]; - } else { -/* NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2. */ - wk[4] = wk[7]; - wk[3] = wk[5] + wk[7]; -/* Computing 2nd power */ - d[0] = ((wk[5] + 2. * wk[3]) * wk[6] * wk[7] + wk[5] * - wk[5] * wk[8]) / wk[3]; - } - } else if (ibeg == 1) { -/* SLOPE PRESCRIBED AT LEFT END. */ - wk[4] = 1.; - wk[3] = 0.; - } else { -/* SECOND DERIVATIVE PRESCRIBED AT LEFT END. */ - wk[4] = 2.; - wk[3] = 1.; - d[0] = 3. * wk[6] - 0.5 * wk[5] * d[0]; - } -/* IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND */ -/* CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH */ -/* EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). */ - nm1 = n - 1; - if (nm1 > 1) { - for (j = 1; j < nm1; ++j) { - if (wk[((j) << 1) + 2] == 0.) { - dpchsp_singular; - } - g = -wk[((j + 2) << 1) + 1] / wk[(j << 1) + 2]; - d[j] = g * d[j - 1] + 3. * - (wk[(j << 1) + 3] * wk[((j + 2) << 1) + 2] + - wk[(j << 1) + 5] * wk[(j << 1) + 4]); - wk[(j << 1) + 4] = g * wk[(j << 1) + 1] + 2. * - (wk[(j << 1) + 3] + wk[((j + 2) << 1) + 1]); - } - } -/* CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM */ -/* (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) */ -/* IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- */ -/* SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT */ -/* AT THIS POINT. */ - if (iend != 1) { - if (iend == 0 && n == 2 && ibeg == 0) { -/* NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2. */ - d[1] = wk[6]; - } else { - if (iend == 0) { - if (n == 2 || n == 3 && ibeg == 0) { -/* EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT* */ -/* NOT-A-KNOT AT LEFT END POINT). */ - d[n-1] = 2. * wk[(n << 1) + 2]; - wk[(n << 1) + 2] = 1.; - if (wk[((n - 1) << 1) + 2] == 0.) { - dpchsp_singular; - } - g = -1. / wk[((n - 1) << 1) + 2]; - } else { -/* NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A- */ -/* KNOT AT LEFT END POINT. */ - g = wk[((n - 1) << 1) + 1] + wk[(n << 1) + 1]; -/* DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES). */ -/* Computing 2nd power */ - dtmp = wk[(n << 1) + 1]; - d[n-1] = ((wk[(n << 1) + 1] + 2. * g) * wk[(n << 1) + 2] * wk[((n - 1) << 1) + 1] + dtmp * dtmp * (f[n - 2] - f[n - 3]) / wk[((n - 1) << 1) + 1]) / g; - if (wk[((n - 1) << 1) + 2] == 0.) { - dpchsp_singular; - } - g = -g / wk[((n - 1) << 1) + 2]; - wk[(n << 1) + 2] = wk[((n - 1) << 1) + 1]; - } - } else { -/* SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT. */ - d[n-1] = 3. * wk[(n << 1) + 2] + 0.5 * wk[(n << 1) + 1] * d[n-1]; - wk[(n << 1) + 2] = 2.; - if (wk[((n - 1) << 1) + 2] == 0.) { - dpchsp_singular; - } - g = -1. / wk[((n - 1) << 1) + 2]; - } -/* COMPLETE FORWARD PASS OF GAUSS ELIMINATION. */ - wk[(n << 1) + 2] = g * wk[((n - 1) << 1) + 1] + wk[(n << 1) + 2]; - if (wk[(n << 1) + 2] == 0.) { - dpchsp_singular; - } - d[n-1] = (g * d[n - 2] + d[n-1]) - / wk[(n << 1) + 2]; - } - } -/* CARRY OUT BACK SUBSTITUTION */ - for (j = nm1-1; j >= 0; --j) { - if (wk[(j << 1) + 4] == 0.) { - dpchsp_singular; - } - d[j] = (d[j] - wk[(j << 1) + 3] * d[j+1]) / wk[(j << 1) + 4]; - } -/* --------------------( END CODING FROM CUBSPL )-------------------- */ -} diff --git a/lib/PDL/Primitive.pd b/lib/PDL/Primitive.pd index 36a23b08f..0cb247cf5 100644 --- a/lib/PDL/Primitive.pd +++ b/lib/PDL/Primitive.pd @@ -5009,19 +5009,200 @@ Analysis 17, 2 (April 1980), pp. 238-246. EOF ); -defslatec('pchip_chsp', - {D => 'dpchsp'}, - 'Integer ic (two=2); - Mat vc (two=2); - IntVal n=$SIZE(n); - Mat x (n); - Mat f (n); - Mat [o] d (n); - Mat [t] wk (nwk=CALC(2*$SIZE(n))); - IntVal nwk=$SIZE(nwk); - Integer [o] ierr (); - ', -'Computes the Hermite representation of the cubic spline interpolant +pp_def('pchip_chsp', + Pars => 'longlong ic(two=2); vc(two=2); x(n); f(n); + [o]d(n); longlong [o]ierr(); + [t]dx(n); [t]dy_dx(n); + ', + GenericTypes => $F, + RedoDimsCode => <<'EOF', +if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO"); +EOF + Code => pp_line_numbers(__LINE__, <<'EOF'), +/* SINGULAR SYSTEM. */ +/* *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES *** */ +/* *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). *** */ +#define dpchsp_singular(x, ind) \ + $ierr() = -8; $CROAK("SINGULAR LINEAR SYSTEM(" #x ") at %td", ind); +/* Local variables */ +$GENERIC() stemp[3], xtemp[4]; +PDL_Indx n = $SIZE(n), nm1 = n - 1; +/* VALIDITY-CHECK ARGUMENTS. */ +loop (n=1) %{ + if ($x() > $x(n=>n-1)) continue; + $ierr() = -1; + $CROAK("X-ARRAY NOT STRICTLY INCREASING"); +%} +PDL_Indx ibeg = $ic(two=>0), iend = $ic(two=>1), j; +$ierr() = 0; +if (PDL_ABS(ibeg) > 5) + --($ierr()); +if (PDL_ABS(iend) > 5) + $ierr() += -2; +if ($ierr() < 0) { + $ierr() += -3; + $CROAK("IC OUT OF RANGE"); +} +/* FUNCTION DEFINITION IS OK -- GO ON. */ +/* COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, */ +/* COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). */ +loop (n=1) %{ + $dx() = $x() - $x(n=>n-1); + $dy_dx() = ($f() - $f(n=>n-1)) / $dx(); +%} +/* SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. */ +if (ibeg > n) { + ibeg = 0; +} +if (iend > n) { + iend = 0; +} +/* SET UP FOR BOUNDARY CONDITIONS. */ +if (ibeg == 1 || ibeg == 2) { + $d(n=>0) = $vc(two=>0); +} else if (ibeg > 2) { +/* PICK UP FIRST IBEG POINTS, IN REVERSE ORDER. */ + for (j = 0; j < ibeg; ++j) { + PDL_Indx index = ibeg - j + 1; +/* INDEX RUNS FROM IBEG DOWN TO 1. */ + xtemp[j] = $x(n=>index); + if (j < ibeg-1) + stemp[j] = $dy_dx(n=>index); + } +/* -------------------------------- */ + $PCHDF(ibeg, xtemp, stemp, $d(n=>0)); +/* -------------------------------- */ + ibeg = 1; +} +if (iend == 1 || iend == 2) { + $d(n=>n-1) = $vc(two=>1); +} else if (iend > 2) { +/* PICK UP LAST IEND POINTS. */ + for (j = 0; j < iend; ++j) { + PDL_Indx index = n - iend + j; +/* INDEX RUNS FROM N+1-IEND UP TO N. */ + xtemp[j] = $x(n=>index); + if (j < iend-1) + stemp[j] = $dy_dx(n=>index+1); + } +/* -------------------------------- */ + $PCHDF(iend, xtemp, stemp, $d(n=>n-1)); +/* -------------------------------- */ + iend = 1; +} +/* --------------------( BEGIN CODING FROM CUBSPL )-------------------- */ +/* **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF */ +/* F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- */ +/* INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. */ +/* WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. */ +/* CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM */ +/* WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) */ +if (ibeg == 0) { + if (n == 2) { +/* NO CONDITION AT LEFT END AND N = 2. */ + $dy_dx(n=>0) = 1.; + $dx(n=>0) = 1.; + $d(n=>0) = 2. * $dy_dx(n=>1); + } else { +/* NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2. */ + $dy_dx(n=>0) = $dx(n=>2); + $dx(n=>0) = $dx(n=>1) + $dx(n=>2); +/* Computing 2nd power */ + $d(n=>0) = (($dx(n=>1) + 2. * $dx(n=>0)) * $dy_dx(n=>1) * $dx(n=>2) + $dx(n=>1) * + $dx(n=>1) * $dy_dx(n=>2)) / $dx(n=>0); + } +} else if (ibeg == 1) { +/* SLOPE PRESCRIBED AT LEFT END. */ + $dy_dx(n=>0) = 1.; + $dx(n=>0) = 0.; +} else { +/* SECOND DERIVATIVE PRESCRIBED AT LEFT END. */ + $dy_dx(n=>0) = 2.; + $dx(n=>0) = 1.; + $d(n=>0) = 3. * $dy_dx(n=>1) - 0.5 * $dx(n=>1) * $d(n=>0); +} +/* IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND */ +/* CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH */ +/* EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). */ +if (n > 2) { + loop (n=1:-1) %{ + if ($dy_dx(n=>n-1) == 0.) { + dpchsp_singular(1, n-1); + } + $GENERIC() g = -$dx(n=>n+1) / $dy_dx(n=>n-1); + $d() = g * $d(n=>n-1) + 3. * ($dx() * $dy_dx(n=>n+1) + + $dx(n=>n+1) * $dy_dx()); + $dy_dx() = g * $dx(n=>n-1) + 2. * ($dx() + $dx(n=>n+1)); + %} +} +/* CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM */ +/* (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) */ +/* IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- */ +/* SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT */ +/* AT THIS POINT. */ +if (iend != 1) { + if (iend == 0 && n == 2 && ibeg == 0) { +/* NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2. */ + $d(n=>1) = $dy_dx(n=>1); + } else { + $GENERIC() g; + if (iend == 0) { + if (n == 2 || n == 3 && ibeg == 0) { +/* EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT* */ +/* NOT-A-KNOT AT LEFT END POINT). */ + $d(n=>n-1) = 2. * $dy_dx(n=>n-1); + $dy_dx(n=>n-1) = 1.; + if ($dy_dx(n=>n-2) == 0.) { + dpchsp_singular(2, n-2); + } + g = -1. / $dy_dx(n=>n-2); + } else { +/* NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A- */ +/* KNOT AT LEFT END POINT. */ + g = $dx(n=>n-2) + $dx(n=>n-1); +/* DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES). */ +/* Computing 2nd power */ + $GENERIC() dtmp = $dx(n=>n-1); + $d(n=>n-1) = (($dx(n=>n-1) + 2. * g) * $dy_dx(n=>n-1) * $dx(n=>n-2) + dtmp * dtmp * ($f(n=>n-2) - $f(n=>n-3)) / $dx(n=>n-2)) / g; + if ($dy_dx(n=>n-2) == 0.) { + dpchsp_singular(3, n-2); + } + g /= -$dy_dx(n=>n-2); + $dy_dx(n=>n-1) = $dx(n=>n-2); + } + } else { +/* SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT. */ + $d(n=>n-1) = 3. * $dy_dx(n=>n-1) + 0.5 * $dx(n=>n-1) * $d(n=>n-1); + $dy_dx(n=>n-1) = 2.; + if ($dy_dx(n=>n-2) == 0.) { + dpchsp_singular(4, n-2); + } + g = -1. / $dy_dx(n=>n-2); + } +/* COMPLETE FORWARD PASS OF GAUSS ELIMINATION. */ + $dy_dx(n=>n-1) = g * $dx(n=>n-2) + $dy_dx(n=>n-1); + if ($dy_dx(n=>n-1) == 0.) { + dpchsp_singular(5, n-1); + } + $d(n=>n-1) = (g * $d(n=>n-2) + $d(n=>n-1)) / $dy_dx(n=>n-1); + } +} +/* CARRY OUT BACK SUBSTITUTION */ +loop (n=nm1-1::-1) %{ + if ($dy_dx() == 0.) { + dpchsp_singular(6, n); + } + $d() = ($d() - $dx() * $d(n=>n+1)) / $dy_dx(); +%} +/* --------------------( END CODING FROM CUBSPL )-------------------- */ +#undef dpchsp_singular +EOF + Doc => <<'EOF', +=for ref + +Calculate the derivatives of (x,f(x)) using cubic spline interpolation. + +Computes the Hermite representation of the cubic spline interpolant to the data given in (C<$x,$f>), with the specified boundary conditions. Control over the boundary conditions is given by the C<$ic> and C<$vc> ndarrays. @@ -5113,8 +5294,7 @@ for the interior derivative values. References: Carl de Boor, A Practical Guide to Splines, Springer-Verlag, New York, 1978, pp. 53-59. -', -'Calculate the derivatives of (x,f(x)) using cubic spline interpolation.' +EOF ); defslatec('pchip_chfd',