Skip to content

Commit

Permalink
make chia 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 b4b7d9f commit fb1f01c
Show file tree
Hide file tree
Showing 2 changed files with 134 additions and 277 deletions.
241 changes: 0 additions & 241 deletions lib/PDL/Primitive-pchip.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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;
}
Loading

0 comments on commit fb1f01c

Please sign in to comment.