Skip to content

Commit

Permalink
make chsp 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 bf7768e commit e3ed1de
Show file tree
Hide file tree
Showing 2 changed files with 195 additions and 238 deletions.
223 changes: 0 additions & 223 deletions lib/PDL/Primitive-pchip.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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. */
Expand Down Expand Up @@ -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 )-------------------- */
}
Loading

0 comments on commit e3ed1de

Please sign in to comment.