PCHSP Subroutine

subroutine PCHSP(IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)

BEGIN PROLOGUE PCHSP PURPOSE Set derivatives needed to determine the Hermite represen- tation of the cubic spline interpolant to given data, with specified boundary conditions. LIBRARY SLATEC (PCHIP) CATEGORY E1A TYPE SINGLE PRECISION (PCHSP-S, DPCHSP-D) KEYWORDS CUBIC HERMITE INTERPOLATION, PCHIP, PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION AUTHOR Fritsch, F. N., (LLNL) Lawrence Livermore National Laboratory P.O. Box 808 (L-316) Livermore, CA 94550 FTS 532-4275, (510) 422-4275 DESCRIPTION

  PCHSP:   Piecewise Cubic Hermite Spline

Computes the Hermite representation of the cubic spline inter- polant to the data given in X and F satisfying the boundary conditions specified by IC and VC.

To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays.

The resulting piecewise cubic Hermite function may be evaluated by PCHFE or PCHFD.

NOTE: This is a modified version of C. de Boor's cubic spline routine CUBSPL.

Calling sequence:

PARAMETER  (INCFD = ...)
INTEGER  IC(2), N, NWK, IERR
REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)

CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)

Parameters:

IC -- (input) integer array of length 2 specifying desired boundary conditions: IC(1) = IBEG, desired condition at beginning of data. IC(2) = IEND, desired condition at end of data.

   IBEG = 0  to set D(1) so that the third derivative is con-
      tinuous at X(2).  This is the "not a knot" condition
      provided by de Boor's cubic spline routine CUBSPL.
      < This is the default boundary condition. >
   IBEG = 1  if first derivative at X(1) is given in VC(1).
   IBEG = 2  if second derivative at X(1) is given in VC(1).
   IBEG = 3  to use the 3-point difference formula for D(1).
             (Reverts to the default b.c. if N.LT.3 .)
   IBEG = 4  to use the 4-point difference formula for D(1).
             (Reverts to the default b.c. if N.LT.4 .)
  NOTES:
   1. An error return is taken if IBEG is out of range.
   2. For the "natural" boundary condition, use IBEG=2 and
      VC(1)=0.

   IEND may take on the same values as IBEG, but applied to
   derivative at X(N).  In case IEND = 1 or 2, the value is
   given in VC(2).

  NOTES:
   1. An error return is taken if IEND is out of range.
   2. For the "natural" boundary condition, use IEND=2 and
      VC(2)=0.

VC -- (input) real array of length 2 specifying desired boundary values, as indicated above. VC(1) need be set only if IC(1) = 1 or 2 . VC(2) need be set only if IC(2) = 1 or 2 .

N -- (input) number of data points. (Error return if N.LT.2 .)

X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.)

F -- (input) real array of dependent variable values to be inter- polated. F(1+(I-1)*INCFD) is value corresponding to X(I).

D -- (output) real array of derivative values at the data points. These values will determine the cubic spline interpolant with the requested boundary conditions. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed.

INCFD -- (input) increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD.LT.1 .)

WK -- (scratch) real array of working storage.

NWK -- (input) length of work array. (Error return if NWK.LT.2*N .)

IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if IBEG.LT.0 or IBEG.GT.4 . IERR = -5 if IEND.LT.0 of IEND.GT.4 . IERR = -6 if both of the above are true. IERR = -7 if NWK is too small. NOTE: The above errors are checked in the order listed, and following arguments have NOT been validated. (The D-array has not been changed in any of these cases.) IERR = -8 in case of trouble solving the linear system for the interior derivative values. (The D-array may have been changed in this case.) ( Do NOT use it! )

REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- Verlag, New York, 1978, pp. 53-59. ROUTINES CALLED PCHDF, XERMSG REVISION HISTORY (YYMMDD) 820503 DATE WRITTEN 820804 Converted to SLATEC library version. 870707 Minor cosmetic changes to prologue. 890411 Added SAVE statements (Vers. 3.2). 890703 Corrected category record. (WRB) 890831 Modified array declarations. (WRB) 890831 REVISION DATE from Version 3.2 891214 Prologue converted to Version 4.0 format. (BAB) 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 920429 Revised format and order of references. (WRB,FNF) END PROLOGUE PCHSP Programming notes:

To produce a double precision version, simply: a. Change PCHSP to DPCHSP wherever it occurs, b. Change the real declarations to double precision, and c. Change the constants ZERO, HALF, ... to double precision.

DECLARE ARGUMENTS.

**FIRST EXECUTABLE STATEMENT PCHSP

Arguments