Skip to content

Commit

Permalink
make chim into PP kernel
Browse files Browse the repository at this point in the history
  • Loading branch information
@mohawk2
mohawk2 committed Dec 16, 2024
1 parent e3ed1de commit 5d02621
Show file tree
Hide file tree
Showing 2 changed files with 105 additions and 203 deletions.
192 changes: 0 additions & 192 deletions lib/PDL/Primitive-pchip.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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;
}
}
}
116 changes: 105 additions & 11 deletions lib/PDL/Primitive.pd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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',
Expand Down

0 comments on commit 5d02621

Please sign in to comment.