# dpchsp.f

          SUBROUTINE DPCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
    C***BEGIN PROLOGUE  DPCHSP
    C***PURPOSE  Set derivatives needed to determine the Hermite represen-
    C            tation of the cubic spline interpolant to given data, with
    C            specified boundary conditions.
    C***LIBRARY   SLATEC (PCHIP)
    C***CATEGORY  E1A
    C***TYPE      DOUBLE PRECISION (PCHSP-S, DPCHSP-D)
    C***KEYWORDS  CUBIC HERMITE INTERPOLATION, PCHIP,
    C             PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION
    C***AUTHOR  Fritsch, F. N., (LLNL)
    C             Lawrence Livermore National Laboratory
    C             P.O. Box 808  (L-316)
    C             Livermore, CA  94550
    C             FTS 532-4275, (510) 422-4275
    C***DESCRIPTION
    C
    C          DPCHSP:   Piecewise Cubic Hermite Spline
    C
    C     Computes the Hermite representation of the cubic spline inter-
    C     polant to the data given in X and F satisfying the boundary
    C     conditions specified by IC and VC.
    C
    C     To facilitate two-dimensional applications, includes an increment
    C     between successive values of the F- and D-arrays.
    C
    C     The resulting piecewise cubic Hermite function may be evaluated
    C     by DPCHFE or DPCHFD.
    C
    C     NOTE:  This is a modified version of C. de Boor's cubic spline
    C            routine CUBSPL.
    C
    C ----------------------------------------------------------------------
    C
    C  Calling sequence:
    C
    C        PARAMETER  (INCFD = ...)
    C        INTEGER  IC(2), N, NWK, IERR
    C        DOUBLE PRECISION  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
    C
    C        CALL  DPCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
    C
    C   Parameters:
    C
    C     IC -- (input) integer array of length 2 specifying desired
    C           boundary conditions:
    C           IC(1) = IBEG, desired condition at beginning of data.
    C           IC(2) = IEND, desired condition at end of data.
    C
    C           IBEG = 0  to set D(1) so that the third derivative is con-
    C              tinuous at X(2).  This is the "not a knot" condition
    C              provided by de Boor's cubic spline routine CUBSPL.
    C              < This is the default boundary condition. >
    C           IBEG = 1  if first derivative at X(1) is given in VC(1).
    C           IBEG = 2  if second derivative at X(1) is given in VC(1).
    C           IBEG = 3  to use the 3-point difference formula for D(1).
    C                     (Reverts to the default b.c. if N.LT.3 .)
    C           IBEG = 4  to use the 4-point difference formula for D(1).
    C                     (Reverts to the default b.c. if N.LT.4 .)
    C          NOTES:
    C           1. An error return is taken if IBEG is out of range.
    C           2. For the "natural" boundary condition, use IBEG=2 and
    C              VC(1)=0.
    C
    C           IEND may take on the same values as IBEG, but applied to
    C           derivative at X(N).  In case IEND = 1 or 2, the value is
    C           given in VC(2).
    C
    C          NOTES:
    C           1. An error return is taken if IEND is out of range.
    C           2. For the "natural" boundary condition, use IEND=2 and
    C              VC(2)=0.
    C
    C     VC -- (input) real*8 array of length 2 specifying desired boundary
    C           values, as indicated above.
    C           VC(1) need be set only if IC(1) = 1 or 2 .
    C           VC(2) need be set only if IC(2) = 1 or 2 .
    C
    C     N -- (input) number of data points.  (Error return if N.LT.2 .)
    C
    C     X -- (input) real*8 array of independent variable values.  The
    C           elements of X must be strictly increasing:
    C                X(I-1) .LT. X(I),  I = 2(1)N.
    C           (Error return if not.)
    C
    C     F -- (input) real*8 array of dependent variable values to be
    C           interpolated.  F(1+(I-1)*INCFD) is value corresponding to
    C           X(I).
    C
    C     D -- (output) real*8 array of derivative values at the data
    C           points.  These values will determine the cubic spline
    C           interpolant with the requested boundary conditions.
    C           The value corresponding to X(I) is stored in
    C                D(1+(I-1)*INCFD),  I=1(1)N.
    C           No other entries in D are changed.
    C
    C     INCFD -- (input) increment between successive values in F and D.
    C           This argument is provided primarily for 2-D applications.
    C           (Error return if  INCFD.LT.1 .)
    C
    C     WK -- (scratch) real*8 array of working storage.
    C
    C     NWK -- (input) length of work array.
    C           (Error return if NWK.LT.2*N .)
    C
    C     IERR -- (output) error flag.
    C           Normal return:
    C              IERR = 0  (no errors).
    C           "Recoverable" errors:
    C              IERR = -1  if N.LT.2 .
    C              IERR = -2  if INCFD.LT.1 .
    C              IERR = -3  if the X-array is not strictly increasing.
    C              IERR = -4  if IBEG.LT.0 or IBEG.GT.4 .
    C              IERR = -5  if IEND.LT.0 of IEND.GT.4 .
    C              IERR = -6  if both of the above are true.
    C              IERR = -7  if NWK is too small.
    C               NOTE:  The above errors are checked in the order listed,
    C                   and following arguments have **NOT** been validated.
    C             (The D-array has not been changed in any of these cases.)
    C              IERR = -8  in case of trouble solving the linear system
    C                         for the interior derivative values.
    C             (The D-array may have been changed in this case.)
    C             (             Do **NOT** use it!                )
    C
    C***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
    C                 Verlag, New York, 1978, pp. 53-59.
    C***ROUTINES CALLED  DPCHDF, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   820503  DATE WRITTEN
    C   820804  Converted to SLATEC library version.
    C   870707  Corrected XERROR calls for d.p. name(s).
    C   890206  Corrected XERROR calls.
    C   890411  Added SAVE statements (Vers. 3.2).
    C   890703  Corrected category record.  (WRB)
    C   890831  Modified array declarations.  (WRB)
    C   891006  Cosmetic changes to prologue.  (WRB)
    C   891006  REVISION DATE from Version 3.2
    C   891214  Prologue converted to Version 4.0 format.  (BAB)
    C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    C   920429  Revised format and order of references.  (WRB,FNF)
    C***END PROLOGUE  DPCHSP