Skip to content

Commit

Permalink
chfe replaced with call to chfd with NULL de which can handle
Browse files Browse the repository at this point in the history
  • Loading branch information
@mohawk2
mohawk2 committed Dec 13, 2024
1 parent 5341779 commit eb61697
Showing 1 changed file with 24 additions and 264 deletions.
288 changes: 24 additions & 264 deletions lib/PDL/Primitive-pchip.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,6 @@ void dpchkt(integer, doublereal *, integer *, doublereal *);
void dchfdv(doublereal, doublereal,
doublereal, doublereal, doublereal, doublereal, integer,
doublereal *, doublereal *, doublereal *, integer *, integer *);
void dchfev(doublereal, doublereal,
doublereal, doublereal, doublereal, doublereal, integer,
doublereal *, doublereal *, integer *, integer *);
doublereal dchfie(doublereal, doublereal, doublereal,
doublereal, doublereal, doublereal, doublereal, doublereal);
doublereal dpchid(integer, doublereal *, doublereal *,
Expand Down Expand Up @@ -426,7 +423,7 @@ void dchfdv(doublereal x1, doublereal x2, doublereal f1,
for (i = 0; i < ne; ++i) {
x = xe[i] - x1;
fe[i] = f1 + x * (d1 + x * (c2 + x * c3));
de[i] = d1 + x * (c2t2 + x * c3t3);
if (de) de[i] = d1 + x * (c2t2 + x * c3t3);
/* COUNT EXTRAPOLATION POINTS. */
if (x < xmi) {
++next[0];
Expand All @@ -439,10 +436,6 @@ void dchfdv(doublereal x1, doublereal x2, doublereal f1,
}

/* Programming notes: */
/* 1. To produce a single precision version, simply: */
/* a. Change DPCHFD to PCHFD, and DCHFDV to CHFDV, wherever they */
/* occur, */
/* b. Change the double precision declaration to real, */
/* 2. Most of the coding between the call to DCHFDV and the end of */
/* the IR-loop could be eliminated if it were permissible to */
/* assume that XE is ordered relative to X. */
Expand All @@ -460,8 +453,7 @@ void dpchfd(integer n, doublereal *x, doublereal *f,
doublereal *xe, doublereal *fe, doublereal *de, integer *ierr)
{
/* Local variables */
integer i, j, nj, ir, ierc, next[2], jfirst;
integer located;
integer i, j, nj, ir, ierc, next[2], jfirst, located;
/* VALIDITY-CHECK ARGUMENTS. */
if (ne < 1) {
*ierr = -4;
Expand Down Expand Up @@ -507,270 +499,31 @@ void dpchfd(integer n, doublereal *x, doublereal *f,
}
/* HAVE LOCATED FIRST POINT BEYOND INTERVAL. */
nj = j - jfirst;
/* SKIP EVALUATION IF NO POINTS IN INTERVAL. */
if (nj != 0) {
/* EVALUATE CUBIC AT XE(I), I = JFIRST (1) J-1 . */
/* ---------------------------------------------------------------- */
dchfdv(x[ir - 1], x[ir], f[ir - 1], f[ir],
d[ir - 1], d[ir], nj,
&xe[jfirst], &fe[jfirst], &de[jfirst], next, &ierc);
/* ---------------------------------------------------------------- */
if (ierc < 0) {
*ierr = -5;
xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr);
return;
}
if (next[1] != 0) {
/* IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE */
/* RIGHT OF X(IR). */
if (ir < n-1) {
/* WE SHOULD NEVER HAVE GOTTEN HERE. */
*ierr = -5;
xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr);
return;
}
/* THESE ARE ACTUALLY EXTRAPOLATION POINTS. */
*ierr += next[1];
}
if (next[0] != 0) {
/* IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE */
/* LEFT OF X(IR-1). */
if (ir < 2) {
/* THESE ARE ACTUALLY EXTRAPOLATION POINTS. */
*ierr += next[0];
goto L49;
}
/* XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST */
/* EVALUATION INTERVAL. */
/* FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). */
located = 0;
for (i = jfirst; i < j; ++i) {
if (xe[i] < x[ir - 1]) {
located = 1;
break;
}
}
if (!located) {
/* NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR */
/* IN DCHFDV. */
*ierr = -5;
xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr);
return;
}
/* RESET J. (THIS WILL BE THE NEW JFIRST.) */
j = i;
/* NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. */
for (i = 0; i < ir; ++i) {
if (xe[j] < x[i]) {
break;
}
}
/* NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1). */
/* AT THIS POINT, EITHER XE(J) .LT. X(1) */
/* OR X(I-1) .LE. XE(J) .LT. X(I) . */
/* RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE */
/* CYCLING. */
/* Computing MAX */
ir = max(0,i-1);
}
L49:
jfirst = j;
/* END OF IR-LOOP. */
}
++ir;
if (ir > n-1) break;
}
}

/* ***PURPOSE Evaluate a cubic polynomial given in Hermite form at an */
/* array of points. While designed for use by DPCHFE, it may */
/* be useful directly as an evaluator for a piecewise cubic */
/* Hermite function in applications, such as graphing, where */
/* the interval is known in advance. */
/* DCHFEV: Cubic Hermite Function EValuator */
/* Evaluates the cubic polynomial determined by function values */
/* F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points */
/* XE(J), J=1(1)NE. */
/* ---------------------------------------------------------------------- */
/* Calling sequence: */
/* INTEGER NE, NEXT(2), IERR */
/* DOUBLE PRECISION X1, X2, F1, F2, D1, D2, XE(NE), FE(NE) */
/* CALL DCHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR) */
/* Parameters: */
/* X1,X2 -- (input) endpoints of interval of definition of cubic. */
/* (Error return if X1.EQ.X2 .) */
/* F1,F2 -- (input) values of function at X1 and X2, respectively. */
/* D1,D2 -- (input) values of derivative at X1 and X2, respectively. */
/* NE -- (input) number of evaluation points. (Error return if */
/* NE.LT.1 .) */
/* XE -- (input) real*8 array of points at which the function is to */
/* be evaluated. If any of the XE are outside the interval */
/* [X1,X2], a warning error is returned in NEXT. */
/* FE -- (output) real*8 array of values of the cubic function */
/* defined by X1,X2, F1,F2, D1,D2 at the points XE. */
/* NEXT -- (output) integer array indicating number of extrapolation */
/* points: */
/* NEXT(1) = number of evaluation points to left of interval. */
/* NEXT(2) = number of evaluation points to right of interval. */
/* IERR -- (output) error flag. */
/* Normal return: */
/* IERR = 0 (no errors). */
/* "Recoverable" errors: */
/* IERR = -1 if NE.LT.1 . */
/* IERR = -2 if X1.EQ.X2 . */
/* (The FE-array has not been changed in either case.) */
/* Programming notes: */
/* To produce a single precision version, simply: */
/* a. Change DCHFEV to CHFEV wherever it occurs, */
/* b. Change the double precision declaration to real, and */
/* c. Change the constant ZERO to single precision. */
void dchfev(doublereal x1, doublereal x2, doublereal f1,
doublereal f2, doublereal d1, doublereal d2, integer ne,
doublereal *xe, doublereal *fe, integer *next, integer *ierr)
{
/* Local variables */
doublereal h;
integer i;
doublereal x, c2, c3, xma, xmi, del1, del2, delta;
if (ne < 1) {
*ierr = -1;
xermsg_("SLATEC", "DCHFEV", "NUMBER OF EVALUATION POINTS LESS THAN ONE", *ierr);
return;
}
h = x2 - x1;
if (h == 0.) {
*ierr = -2;
xermsg_("SLATEC", "DCHFEV", "INTERVAL ENDPOINTS EQUAL", *ierr);
return;
}
/* INITIALIZE. */
*ierr = 0;
next[0] = 0;
next[1] = 0;
xmi = min(0.,h);
xma = max(0.,h);
/* COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). */
delta = (f2 - f1) / h;
del1 = (d1 - delta) / h;
del2 = (d2 - delta) / h;
/* (DELTA IS NO LONGER NEEDED.) */
c2 = -(del1 + del1 + del2);
c3 = (del1 + del2) / h;
/* (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) */
/* EVALUATION LOOP. */
for (i = 0; i < ne; ++i) {
x = xe[i] - x1;
fe[i] = f1 + x * (d1 + x * (c2 + x * c3));
/* COUNT EXTRAPOLATION POINTS. */
if (x < xmi) {
++next[0];
}
if (x > xma) {
++next[1];
}
/* (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) */
}
}

/* Programming notes: */
/* 1. To produce a single precision version, simply: */
/* a. Change DPCHFE to PCHFE, and DCHFEV to CHFEV, wherever they */
/* occur, */
/* b. Change the double precision declaration to real, */
/* 2. Most of the coding between the call to DCHFEV and the end of */
/* the IR-loop could be eliminated if it were permissible to */
/* assume that XE is ordered relative to X. */
/* 3. DCHFEV does not assume that X1 is less than X2. thus, it would */
/* be possible to write a version of DPCHFE that assumes a */
/* decreasing X-array by simply running the IR-loop backwards */
/* (and reversing the order of appropriate tests). */
/* 4. The present code has a minor bug, which I have decided is not */
/* worth the effort that would be required to fix it. */
/* If XE contains points in [X(N-1),X(N)], followed by points .LT. */
/* X(N-1), followed by points .GT.X(N), the extrapolation points */
/* will be counted (at least) twice in the total returned in IERR. */
void dpchfe(integer n, doublereal *x, doublereal *f,
doublereal *d, logical *skip, integer ne,
doublereal *xe, doublereal *fe, integer *ierr)
{
/* Local variables */
integer i, j, nj, ir, ierc, next[2];
integer jfirst, located;
/* VALIDITY-CHECK ARGUMENTS. */
if (ne < 1) {
*ierr = -4;
xermsg_("SLATEC", "DPCHFE", "NUMBER OF EVALUATION POINTS LESS THAN ONE", *ierr);
return;
}
if (!*skip) {
if (n < 2) {
*ierr = -1;
xermsg_("SLATEC", "DPCHFE", "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", "DPCHFE", "X-ARRAY NOT STRICTLY INCREASING", *ierr);
return;
}
}
/* FUNCTION DEFINITION IS OK, GO ON. */
}
*ierr = 0;
*skip = TRUE_;
/* LOOP OVER INTERVALS. ( INTERVAL INDEX IS IL = IR-1 . ) */
/* ( INTERVAL IS X(IL).LE.X.LT.X(IR) . ) */
jfirst = 1;
ir = 2;
while (1) {
/* SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS. */
if (jfirst > ne) {
return;
}
/* LOCATE ALL POINTS IN INTERVAL. */
located = 0;
for (j = jfirst; j <= ne; ++j) {
if (xe[j-1] >= x[ir-1]) {
located = 1;
break;
}
}
if (located) {
/* HAVE LOCATED FIRST POINT BEYOND INTERVAL. */
if (ir == n) {
j = ne + 1;
}
} else {
j = ne + 1;
}
nj = j - jfirst;
/* SKIP EVALUATION IF NO POINTS IN INTERVAL. */
if (nj == 0) {
++ir;
if (ir <= n) {
if (ir < n) {
continue;
}
return;
}
/* EVALUATE CUBIC AT XE(I), I = JFIRST (1) J-1 . */
/* ---------------------------------------------------------------- */
dchfev(x[ir - 2], x[ir-1], f[ir - 2], f[ir-1],
d[ir - 2], d[ir - 1], nj,
&xe[jfirst-1], &fe[jfirst-1], next, &ierc);
dchfdv(x[ir - 1], x[ir], f[ir - 1], f[ir], d[ir - 1], d[ir], nj,
&xe[jfirst], &fe[jfirst], de ? &de[jfirst] : NULL, next, &ierc);
/* ---------------------------------------------------------------- */
if (ierc < 0) {
*ierr = -5;
xermsg_("SLATEC", "DPCHFE", "ERROR RETURN FROM DCHFEV -- FATAL", *ierr);
xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr);
return;
}
if (next[1] != 0) {
/* IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE */
/* RIGHT OF X(IR). */
if (ir < n) {
if (ir < n-1) {
/* WE SHOULD NEVER HAVE GOTTEN HERE. */
*ierr = -5;
xermsg_("SLATEC", "DPCHFE", "ERROR RETURN FROM DCHFEV -- FATAL", *ierr);
xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr);
return;
}
/* THESE ARE ACTUALLY EXTRAPOLATION POINTS. */
Expand All @@ -779,7 +532,7 @@ void dpchfe(integer n, doublereal *x, doublereal *f,
if (next[0] != 0) {
/* IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE */
/* LEFT OF X(IR-1). */
if (ir <= 2) {
if (ir < 2) {
/* THESE ARE ACTUALLY EXTRAPOLATION POINTS. */
*ierr += next[0];
goto L49;
Expand All @@ -789,23 +542,23 @@ void dpchfe(integer n, doublereal *x, doublereal *f,
/* FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). */
located = 0;
for (i = jfirst; i < j; ++i) {
if (xe[i-1] < x[ir - 2]) {
if (xe[i] < x[ir - 1]) {
located = 1;
break;
}
}
if (!located) {
/* NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR */
/* IN DCHFEV. */
/* IN DCHFDV. */
*ierr = -5;
xermsg_("SLATEC", "DPCHFE", "ERROR RETURN FROM DCHFEV -- FATAL", *ierr);
xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr);
return;
}
/* RESET J. (THIS WILL BE THE NEW JFIRST.) */
j = i;
/* NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. */
for (i = 1; i < ir; ++i) {
if (xe[j-1] < x[i-1]) {
for (i = 0; i < ir; ++i) {
if (xe[j] < x[i]) {
break;
}
}
Expand All @@ -815,16 +568,23 @@ void dpchfe(integer n, doublereal *x, doublereal *f,
/* RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE */
/* CYCLING. */
/* Computing MAX */
ir = max(1,i - 1);
L49:;
ir = max(0,i-1);
}
L49:
jfirst = j;
/* END OF IR-LOOP. */
++ir;
if (ir > n) break;
if (ir >= n) break;
}
}

void dpchfe(integer n, doublereal *x, doublereal *f,
doublereal *d, logical *skip, integer ne,
doublereal *xe, doublereal *fe, integer *ierr)
{
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 */
Expand Down

0 comments on commit eb61697

Please sign in to comment.