From 87300e7e0132f3b3b4356113be92a4904051469f Mon Sep 17 00:00:00 2001 From: Ed J Date: Mon, 16 Dec 2024 07:39:46 +0000 Subject: [PATCH] make chia into PP kernel --- lib/PDL/Primitive-pchip.c | 241 -------------------------------------- lib/PDL/Primitive.pd | 167 ++++++++++++++++++++------ 2 files changed, 133 insertions(+), 275 deletions(-) diff --git a/lib/PDL/Primitive-pchip.c b/lib/PDL/Primitive-pchip.c index e044aba5b..37e13d004 100644 --- a/lib/PDL/Primitive-pchip.c +++ b/lib/PDL/Primitive-pchip.c @@ -9,18 +9,10 @@ typedef long int logical; #define TRUE_ (1) -#define abs(x) ((x) >= 0 ? (x) : -(x)) #define min(a,b) ((a) <= (b) ? (a) : (b)) #define max(a,b) ((a) >= (b) ? (a) : (b)) /* end f2c.h */ -/* from libf2c */ -double d_sign(doublereal a, doublereal b) -{ - double x = (a >= 0 ? a : - a); - return b >= 0 ? x : -x; -} - #define xermsg_(lib, func, errstr, nerr, ...) \ fprintf(stderr, "%s::%s: %s (err=%ld)\n", lib, func, errstr, nerr) @@ -299,236 +291,3 @@ void dpchfe(integer n, doublereal *x, doublereal *f, { dpchfd(n, x, f, d, skip, ne, xe, fe, NULL, ierr); } - -/* ***PURPOSE Evaluates integral of a single cubic for DPCHIA */ -/* DCHFIE: Cubic Hermite Function Integral Evaluator. */ -/* Called by DPCHIA to evaluate the integral of a single cubic (in */ -/* Hermite form) over an arbitrary interval (A,B). */ -/* ---------------------------------------------------------------------- */ -/* Calling sequence: */ -/* DOUBLE PRECISION X1, X2, F1, F2, D1, D2, A, B */ -/* DOUBLE PRECISION VALUE, DCHFIE */ -/* VALUE = DCHFIE (X1, X2, F1, F2, D1, D2, A, B) */ -/* Parameters: */ -/* VALUE -- (output) value of the requested integral. */ -/* X1,X2 -- (input) endpoints if interval of definition of cubic. */ -/* F1,F2 -- (input) function values at the ends of the interval. */ -/* D1,D2 -- (input) derivative values at the ends of the interval. */ -/* A,B -- (input) endpoints of interval of integration. */ -/* ***SEE ALSO DPCHIA */ -/* Programming notes: */ -/* 1. There is no error return from this routine because zero is */ -/* indeed the mathematically correct answer when X1.EQ.X2 . */ -#define X(ctype, ppsym) \ - static inline ctype chfie_ ## ppsym(ctype x1, ctype x2, ctype f1, ctype f2, \ - ctype d1, ctype d2, ctype a, ctype b) \ - { \ -/* Local variables */ \ - ctype h, ta1, ta2, tb1, tb2, ua1, ua2, ub1, ub2, phia1, \ - phia2, phib1, phib2, psia1, psia2, psib1, psib2, dterm, fterm; \ - if (x1 == x2) { \ - return 0.; \ - } \ - h = x2 - x1; \ - ta1 = (a - x1) / h; \ - ta2 = (x2 - a) / h; \ - tb1 = (b - x1) / h; \ - tb2 = (x2 - b) / h; \ -/* Computing 3rd power */ \ - ua1 = ta1 * (ta1 * ta1); \ - phia1 = ua1 * (2. - ta1); \ - psia1 = ua1 * (3. * ta1 - 4.); \ -/* Computing 3rd power */ \ - ua2 = ta2 * (ta2 * ta2); \ - phia2 = ua2 * (2. - ta2); \ - psia2 = -ua2 * (3. * ta2 - 4.); \ -/* Computing 3rd power */ \ - ub1 = tb1 * (tb1 * tb1); \ - phib1 = ub1 * (2. - tb1); \ - psib1 = ub1 * (3. * tb1 - 4.); \ -/* Computing 3rd power */ \ - ub2 = tb2 * (tb2 * tb2); \ - phib2 = ub2 * (2. - tb2); \ - psib2 = -ub2 * (3. * tb2 - 4.); \ - fterm = f1 * (phia2 - phib2) + f2 * (phib1 - phia1); \ - dterm = (d1 * (psia2 - psib2) + d2 * (psib1 - psia1)) * (h / 6.); \ - return 0.5 * h * (fterm + dterm); \ - } -X(doublereal, D) -#undef X - -/* Programming notes: */ -/* 1. This routine uses a special formula that is valid only for */ -/* integrals whose limits coincide with data values. This is */ -/* mathematically equivalent to, but much more efficient than, */ -/* calls to DCHFIE. */ -doublereal dpchid(integer n, doublereal *x, doublereal *f, doublereal *d, - logical *skip, integer ia, integer ib, integer * - ierr) -{ -/* Local variables */ - doublereal h; - integer i, iup, low; - doublereal sum, value; - value = 0.; -/* VALIDITY-CHECK ARGUMENTS. */ - if (!*skip) { - if (n < 2) { - *ierr = -1; - xermsg_("SLATEC", "DPCHID", "NUMBER OF DATA POINTS LESS THAN TWO", *ierr); - return value; - } - for (i = 1; i < n; ++i) { - if (x[i] <= x[i - 1]) { - *ierr = -3; - xermsg_("SLATEC", "DPCHID", "X-ARRAY NOT STRICTLY INCREASING", *ierr); - return value; - } - } - } -/* FUNCTION DEFINITION IS OK, GO ON. */ - *skip = TRUE_; - if (ia < 0 || ia > n-1 || ib < 0 || ib > n-1) { - *ierr = -4; - xermsg_("SLATEC", "DPCHID", "IA OR IB OUT OF RANGE", *ierr); - return value; - } - *ierr = 0; -/* COMPUTE INTEGRAL VALUE. */ - if (ia != ib) { - low = min(ia,ib); - iup = max(ia,ib); - sum = 0.; - for (i = low; i < iup; ++i) { - h = x[i+1] - x[i]; - sum += h * (f[i] + f[i + 1] + (d[i] - d[i + 1]) * (h / 6.)); - } - value = 0.5 * sum; - if (ia > ib) { - value = -value; - } - } - return value; -} - -/* Programming notes: */ -/* 1. The error flag from DPCHID is tested, because a logic flaw */ -/* could conceivably result in IERD=-4, which should be reported. */ -doublereal dpchia(integer n, doublereal *x, doublereal *f, doublereal *d, - logical *skip, doublereal a, doublereal b, integer *ierr) -{ -/* Local variables */ - integer i, ia, ib, il; - doublereal xa, xb; - integer ir, ierd; - doublereal value; - value = 0.; - if (!*skip) { - if (n < 2) { - *ierr = -1; - xermsg_("SLATEC", "DPCHIA", "NUMBER OF DATA POINTS LESS THAN TWO", *ierr); - return value; - } - for (i = 1; i < n; ++i) { - if (x[i] <= x[i - 1]) { - *ierr = -3; - xermsg_("SLATEC", "DPCHIA", "X-ARRAY NOT STRICTLY INCREASING", *ierr); - return value; - } - } - } -/* FUNCTION DEFINITION IS OK, GO ON. */ - *skip = TRUE_; - *ierr = 0; - if (a < x[0] || a > x[n-1]) { - ++(*ierr); - } - if (b < x[0] || b > x[n-1]) { - *ierr += 2; - } -/* COMPUTE INTEGRAL VALUE. */ - if (a != b) { - xa = min(a,b); - xb = max(a,b); - if (xb <= x[1]) { -/* INTERVAL IS TO LEFT OF X(2), SO USE FIRST CUBIC. */ -/* --------------------------------------- */ - value = chfie_D(x[0], x[1], f[0], f[1], - d[0], d[1], a, b); -/* --------------------------------------- */ - } else if (xa >= x[n - 2]) { -/* INTERVAL IS TO RIGHT OF X(N-1), SO USE LAST CUBIC. */ -/* ------------------------------------------ */ - value = chfie_D(x[n - 2], x[n-1], f[n - 2], f[n-1], d[n-2], d[n-1], a, b); -/* ------------------------------------------ */ - } else { -/* 'NORMAL' CASE -- XA.LT.XB, XA.LT.X(N-1), XB.GT.X(2). */ -/* ......LOCATE IA AND IB SUCH THAT */ -/* X(IA-1).LT.XA.LE.X(IA).LE.X(IB).LE.XB.LE.X(IB+1) */ - ia = 0; - for (i = 0; i < n-1; ++i) { - if (xa > x[i]) { - ia = i + 1; - } - } -/* IA = 1 IMPLIES XA.LT.X(1) . OTHERWISE, */ -/* IA IS LARGEST INDEX SUCH THAT X(IA-1).LT.XA,. */ - ib = n - 1; - for (i = n-1; i >= ia; --i) { - if (xb < x[i]) { - ib = i - 1; - } - } -/* IB = N IMPLIES XB.GT.X(N) . OTHERWISE, */ -/* IB IS SMALLEST INDEX SUCH THAT XB.LT.X(IB+1) . */ -/* ......COMPUTE THE INTEGRAL. */ - if (ib <= ia) { -/* THIS MEANS IB = IA-1 AND */ -/* (A,B) IS A SUBSET OF (X(IB),X(IA)). */ -/* ------------------------------------------- */ - value = chfie_D(x[ib], x[ia], f[ib], - f[ia], d[ib], d[ia], - a, b); -/* ------------------------------------------- */ - } else { -/* FIRST COMPUTE INTEGRAL OVER (X(IA),X(IB)). */ -/* (Case (IB .EQ. IA) is taken care of by initialization */ -/* of VALUE to ZERO.) */ - if (ib > ia-1) { -/* --------------------------------------------- */ - value = dpchid(n, &x[0], &f[0], &d[0], - skip, ia, ib, &ierd); -/* --------------------------------------------- */ - if (ierd < 0) { - *ierr = -4; - xermsg_("SLATEC", "DPCHIA", "TROUBLE IN DPCHID", *ierr); - return value; - } - } -/* THEN ADD ON INTEGRAL OVER (XA,X(IA)). */ - if (xa < x[ia]) { -/* Computing MAX */ - il = max(0,ia - 1); - ir = il + 1; -/* ------------------------------------- */ - value += chfie_D(x[il], x[ir], f[il], f[ir], d[il], d[ir], xa, x[ia]); -/* ------------------------------------- */ - } -/* THEN ADD ON INTEGRAL OVER (X(IB),XB). */ - if (xb > x[ib]) { -/* Computing MIN */ - ir = min(ib + 1,n-1); - il = ir - 1; -/* ------------------------------------- */ - value += chfie_D(x[il], x[ir], f[il], f[ir], d[il], d[ir], x[ib], xb); -/* ------------------------------------- */ - } -/* FINALLY, ADJUST SIGN IF NECESSARY. */ - if (a > b) { - value = -value; - } - } - } - } - return value; -} diff --git a/lib/PDL/Primitive.pd b/lib/PDL/Primitive.pd index 494bc4681..ccbb9da88 100644 --- a/lib/PDL/Primitive.pd +++ b/lib/PDL/Primitive.pd @@ -4330,16 +4330,14 @@ CHFIE => sub {my ($x1, $x2, $f1, $f2, $d1, $d2, $a, $b, $out) = @_; ' /* 1. There is no error return from this routine because zero is */ /* indeed the mathematically correct answer when X1.EQ.X2 . */ do { - $GENERIC() x1='.$x1.', x2='.$x2.', f1='.$f1.', f2='.$f2.', - d1='.$d1.', d2='.$d2.', a='.$a.', b='.$b.'; if ('.$x1.' == '.$x2.') { '.$out.' = 0.; break; } - $GENERIC() h = x2 - x1; - $GENERIC() ta1 = (a - x1) / h; - $GENERIC() ta2 = (x2 - a) / h; - $GENERIC() tb1 = (b - x1) / h; - $GENERIC() tb2 = (x2 - b) / h; + $GENERIC() h = '.$x2.' - '.$x1.'; + $GENERIC() ta1 = ('.$a.' - '.$x1.') / h; + $GENERIC() ta2 = ('.$x2.' - '.$a.') / h; + $GENERIC() tb1 = ('.$b.' - '.$x1.') / h; + $GENERIC() tb2 = ('.$x2.' - '.$b.') / h; /* Computing 3rd power */ $GENERIC() ua1 = ta1 * (ta1 * ta1); $GENERIC() phia1 = ua1 * (2. - ta1); @@ -4356,8 +4354,8 @@ do { $GENERIC() ub2 = tb2 * (tb2 * tb2); $GENERIC() phib2 = ub2 * (2. - tb2); $GENERIC() psib2 = -ub2 * (3. * tb2 - 4.); - $GENERIC() fterm = f1 * (phia2 - phib2) + f2 * (phib1 - phia1); - $GENERIC() dterm = (d1 * (psia2 - psib2) + d2 * (psib1 - psia1)) * (h / 6.); + $GENERIC() fterm = '.$f1.' * (phia2 - phib2) + '.$f2.' * (phib1 - phia1); + $GENERIC() dterm = ('.$d1.' * (psia2 - psib2) + '.$d2.' * (psib1 - psia1)) * (h / 6.); '.$out.' = 0.5 * h * (fterm + dterm); } while(0) '}, @@ -4370,13 +4368,6 @@ PCHID => sub {my ($ia, $ib, $out) = @_; ' /* VALIDITY-CHECK ARGUMENTS. */ do { PDL_Indx ia = '.$ia.', ib = '.$ib.'; - if (!$skip()) { - loop (n=1) %{ - if ($x() > $x(n=>n-1)) continue; - $ierr() = -1; - $CROAK("X-ARRAY NOT STRICTLY INCREASING"); - %} - } /* FUNCTION DEFINITION IS OK, GO ON. */ $skip() = 1; if (ia < 0 || ia > $SIZE(n)-1 || ib < 0 || ib > $SIZE(n)-1) { @@ -4385,14 +4376,14 @@ do { } $ierr() = 0; /* COMPUTE INTEGRAL VALUE. */ - if ($ia() == $ib()) { $ans() = 0; continue; } - PDL_Indx low = PDLMIN($ia(),$ib()), iup = PDLMAX($ia(),$ib()); + if (ia == ib) { '.$out.' = 0; continue; } + PDL_Indx low = PDLMIN(ia,ib), iup = PDLMAX(ia,ib); $GENERIC() sum = 0.; loop (n=low:iup) %{ $GENERIC() h = $x(n=>n+1) - $x(); sum += h * ($f() + $f(n=>n+1) + ($d() - $d(n=>n+1)) * (h / 6.)); %} - '.$out.' = 0.5 * ($ia() > $ib() ? -sum : sum); + '.$out.' = 0.5 * (ia > ib ? -sum : sum); } while(0) '}, ); @@ -5640,19 +5631,121 @@ May be used by itself for Hermite interpolation, or as an evaluator for L or L.' ); -defslatec('pchip_chia', - {D => 'dpchia'}, - 'IntVal n=$SIZE(n); - Mat x (n); - Mat f (n); - Mat d (n); - Logical [io] skip (); - MatScalar la (); - MatScalar lb (); - FuncRet [o] ans (); - Integer [o] ierr (); - ', -'Evaluate the definite integral of a piecewise +pp_def('pchip_chia', + Pars => 'x(n); f(n); d(n); int [io]skip(); la(); lb(); + [o]ans(); longlong [o]ierr()', + GenericTypes => $F, + RedoDimsCode => <<'EOF', +if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO"); +EOF + Code => pp_line_numbers(__LINE__, <<'EOF'), +PDL_Indx ia, ib, il; +$GENERIC() a = $la(), b = $lb(), xa, xb; +PDL_Indx ir, n = $SIZE(n); +$GENERIC() value = 0.; +if (!$skip()) { + loop (n=1) %{ + if ($x() > $x(n=>n-1)) continue; + $ierr() = -1; + $CROAK("X-ARRAY NOT STRICTLY INCREASING"); + %} +} +/* FUNCTION DEFINITION IS OK, GO ON. */ +$skip() = 1; +$ierr() = 0; +if (a < $x(n=>0) || a > $x(n=>n-1)) { + ++($ierr()); +} +if (b < $x(n=>0) || b > $x(n=>n-1)) { + $ierr() += 2; +} +/* COMPUTE INTEGRAL VALUE. */ +if (a != b) { + xa = PDLMIN(a,b); + xb = PDLMAX(a,b); + if (xb <= $x(n=>1)) { +/* INTERVAL IS TO LEFT OF X(2), SO USE FIRST CUBIC. */ +/* --------------------------------------- */ + $CHFIE($x(n=>0), $x(n=>1), $f(n=>0), $f(n=>1), + $d(n=>0), $d(n=>1), a, b, value); +/* --------------------------------------- */ + } else if (xa >= $x(n=>n-2)) { +/* INTERVAL IS TO RIGHT OF X(N-1), SO USE LAST CUBIC. */ +/* ------------------------------------------ */ + $CHFIE($x(n=>n-2), $x(n=>n-1), $f(n=>n-2), $f(n=>n-1), $d(n=>n-2), $d(n=>n-1), a, b, value); +/* ------------------------------------------ */ + } else { +/* 'NORMAL' CASE -- XA.LT.XB, XA.LT.X(N-1), XB.GT.X(2). */ +/* ......LOCATE IA AND IB SUCH THAT */ +/* X(IA-1).LT.XA.LE.X(IA).LE.X(IB).LE.XB.LE.X(IB+1) */ + ia = 0; + loop (n=:-1) %{ + if (xa > $x()) + ia = n + 1; + %} +/* IA = 1 IMPLIES XA.LT.X(1) . OTHERWISE, */ +/* IA IS LARGEST INDEX SUCH THAT X(IA-1).LT.XA,. */ + ib = n - 1; + loop (n=:ia:-1) %{ + if (xb < $x()) + ib = n - 1; + %} +/* IB = N IMPLIES XB.GT.X(N) . OTHERWISE, */ +/* IB IS SMALLEST INDEX SUCH THAT XB.LT.X(IB+1) . */ +/* ......COMPUTE THE INTEGRAL. */ + if (ib <= ia) { +/* THIS MEANS IB = IA-1 AND */ +/* (A,B) IS A SUBSET OF (X(IB),X(IA)). */ +/* ------------------------------------------- */ + $CHFIE($x(n=>ib), $x(n=>ia), $f(n=>ib), + $f(n=>ia), $d(n=>ib), $d(n=>ia), a, b, value); +/* ------------------------------------------- */ + } else { +/* FIRST COMPUTE INTEGRAL OVER (X(IA),X(IB)). */ +/* (Case (IB .EQ. IA) is taken care of by initialization */ +/* of VALUE to ZERO.) */ + if (ib > ia-1) { +/* --------------------------------------------- */ + $PCHID(ia, ib, value); +/* --------------------------------------------- */ + } +/* THEN ADD ON INTEGRAL OVER (XA,X(IA)). */ + if (xa < $x(n=>ia)) { +/* Computing MAX */ + il = PDLMAX(0,ia - 1); + ir = il + 1; +/* ------------------------------------- */ + $GENERIC() chfie_ans = 0; + $CHFIE($x(n=>il), $x(n=>ir), $f(n=>il), $f(n=>ir), $d(n=>il), $d(n=>ir), xa, $x(n=>ia), chfie_ans); + value += chfie_ans; +/* ------------------------------------- */ + } +/* THEN ADD ON INTEGRAL OVER (X(IB),XB). */ + if (xb > $x(n=>ib)) { +/* Computing MIN */ + ir = PDLMIN(ib + 1,n-1); + il = ir - 1; +/* ------------------------------------- */ + $GENERIC() chfie_ans = 0; + $CHFIE($x(n=>il), $x(n=>ir), $f(n=>il), $f(n=>ir), $d(n=>il), $d(n=>ir), $x(n=>ib), xb, chfie_ans); + value += chfie_ans; +/* ------------------------------------- */ + } +/* FINALLY, ADJUST SIGN IF NECESSARY. */ + if (a > b) { + value = -value; + } + } + } +} +$ans() = value; +EOF + Doc => <<'EOF', +=for ref + +Integrate (x,f(x)) over arbitrary limits. + +Evaluate the definite integral of a piecewise cubic Hermite function over an arbitrary interval, given by C<[$la,$lb]>. @@ -5713,8 +5806,7 @@ contains data interval or the intervals do not intersect at all.) which should never happen. =back -', -'Integrate (x,f(x)) over arbitrary limits.' +EOF ); pp_def('pchip_chid', @@ -5725,6 +5817,13 @@ pp_def('pchip_chid', if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO"); EOF Code => <<'EOF', +if (!$skip()) { + loop (n=1) %{ + if ($x() > $x(n=>n-1)) continue; + $ierr() = -1; + $CROAK("X-ARRAY NOT STRICTLY INCREASING"); + %} +} $PCHID($ia(), $ib(), $ans()); EOF GenericTypes => $F,