# pchic.f

          SUBROUTINE PCHIC (IC, VC, SWITCH, N, X, F, D, INCFD, WK, NWK,
         +   IERR)
    C***BEGIN PROLOGUE  PCHIC
    C***PURPOSE  Set derivatives needed to determine a piecewise monotone
    C            piecewise cubic Hermite interpolant to given data.
    C            User control is available over boundary conditions and/or
    C            treatment of points where monotonicity switches direction.
    C***LIBRARY   SLATEC (PCHIP)
    C***CATEGORY  E1A
    C***TYPE      SINGLE PRECISION (PCHIC-S, DPCHIC-D)
    C***KEYWORDS  CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION,
    C             PCHIP, PIECEWISE CUBIC INTERPOLATION,
    C             SHAPE-PRESERVING 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         PCHIC:  Piecewise Cubic Hermite Interpolation Coefficients.
    C
    C     Sets derivatives needed to determine a piecewise monotone piece-
    C     wise cubic interpolant to the data given in X and F satisfying the
    C     boundary conditions specified by IC and VC.
    C
    C     The treatment of points where monotonicity switches direction is
    C     controlled by argument SWITCH.
    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 PCHFE or PCHFD.
    C
    C ----------------------------------------------------------------------
    C
    C  Calling sequence:
    C
    C        PARAMETER  (INCFD = ...)
    C        INTEGER  IC(2), N, NWK, IERR
    C        REAL  VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
    C
    C        CALL  PCHIC (IC, VC, SWITCH, 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  for the default boundary condition (the same as
    C                     used by PCHIM).
    C           If IBEG.NE.0, then its sign indicates whether the boundary
    C                     derivative is to be adjusted, if necessary, to be
    C                     compatible with monotonicity:
    C              IBEG.GT.0  if no adjustment is to be performed.
    C              IBEG.LT.0  if the derivative is to be adjusted for
    C                     monotonicity.
    C
    C           Allowable values for the magnitude of IBEG are:
    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           IBEG = 5  to set D(1) so that the second derivative is con-
    C              tinuous at X(2). (Reverts to the default b.c. if N.LT.4.)
    C              This option is somewhat analogous to the "not a knot"
    C              boundary condition provided by PCHSP.
    C
    C          NOTES (IBEG):
    C           1. An error return is taken if ABS(IBEG).GT.5 .
    C           2. Only in case  IBEG.LE.0  is it guaranteed that the
    C              interpolant will be monotonic in the first interval.
    C              If the returned value of D(1) lies between zero and
    C              3*SLOPE(1), the interpolant will be monotonic.  This
    C              is **NOT** checked if IBEG.GT.0 .
    C           3. If IBEG.LT.0 and D(1) had to be changed to achieve mono-
    C              tonicity, a warning error is returned.
    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 (IEND):
    C           1. An error return is taken if ABS(IEND).GT.5 .
    C           2. Only in case  IEND.LE.0  is it guaranteed that the
    C              interpolant will be monotonic in the last interval.
    C              If the returned value of D(1+(N-1)*INCFD) lies between
    C              zero and 3*SLOPE(N-1), the interpolant will be monotonic.
    C              This is **NOT** checked if IEND.GT.0 .
    C           3. If IEND.LT.0 and D(1+(N-1)*INCFD) had to be changed to
    C              achieve monotonicity, a warning error is returned.
    C
    C     VC -- (input) real 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     SWITCH -- (input) indicates desired treatment of points where
    C           direction of monotonicity switches:
    C           Set SWITCH to zero if interpolant is required to be mono-
    C           tonic in each interval, regardless of monotonicity of data.
    C             NOTES:
    C              1. This will cause D to be set to zero at all switch
    C                 points, thus forcing extrema there.
    C              2. The result of using this option with the default boun-
    C                 dary conditions will be identical to using PCHIM, but
    C                 will generally cost more compute time.
    C                 This option is provided only to facilitate comparison
    C                 of different switch and/or boundary conditions.
    C           Set SWITCH nonzero to use a formula based on the 3-point
    C              difference formula in the vicinity of switch points.
    C           If SWITCH is positive, the interpolant on each interval
    C              containing an extremum is controlled to not deviate from
    C              the data by more than SWITCH*DFLOC, where DFLOC is the
    C              maximum of the change of F on this interval and its two
    C              immediate neighbors.
    C           If SWITCH is negative, no such control is to be imposed.
    C
    C     N -- (input) number of data points.  (Error return if N.LT.2 .)
    C
    C     X -- (input) real 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 array of dependent variable values to be inter-
    C           polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
    C
    C     D -- (output) real array of derivative values at the data points.
    C           These values will determine a monotone cubic Hermite func-
    C           tion on each subinterval on which the data are monotonic,
    C           except possibly adjacent to switches in monotonicity.
    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 array of working storage.  The user may wish
    C           to know that the returned values are:
    C              WK(I)     = H(I)     = X(I+1) - X(I) ;
    C              WK(N-1+I) = SLOPE(I) = (F(1,I+1) - F(1,I)) / H(I)
    C           for  I = 1(1)N-1.
    C
    C     NWK -- (input) length of work array.
    C           (Error return if  NWK.LT.2*(N-1) .)
    C
    C     IERR -- (output) error flag.
    C           Normal return:
    C              IERR = 0  (no errors).
    C           Warning errors:
    C              IERR = 1  if IBEG.LT.0 and D(1) had to be adjusted for
    C                        monotonicity.
    C              IERR = 2  if IEND.LT.0 and D(1+(N-1)*INCFD) had to be
    C                        adjusted for monotonicity.
    C              IERR = 3  if both of the above are true.
    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 ABS(IBEG).GT.5 .
    C              IERR = -5  if ABS(IEND).GT.5 .
    C              IERR = -6  if both of the above are true.
    C              IERR = -7  if NWK.LT.2*(N-1) .
    C             (The D-array has not been changed in any of these cases.)
    C               NOTE:  The above errors are checked in the order listed,
    C                   and following arguments have **NOT** been validated.
    C
    C***REFERENCES  1. F. N. Fritsch, Piecewise Cubic Hermite Interpolation
    C                 Package, Report UCRL-87285, Lawrence Livermore Nation-
    C                 al Laboratory, July 1982.  [Poster presented at the
    C                 SIAM 30th Anniversary Meeting, 19-23 July 1982.]
    C               2. F. N. Fritsch and J. Butland, A method for construc-
    C                 ting local monotone piecewise cubic interpolants, SIAM
    C                 Journal on Scientific and Statistical Computing 5, 2
    C                 (June 1984), pp. 300-304.
    C               3. F. N. Fritsch and R. E. Carlson, Monotone piecewise
    C                 cubic interpolation, SIAM Journal on Numerical Ana-
    C                 lysis 17, 2 (April 1980), pp. 238-246.
    C***ROUTINES CALLED  PCHCE, PCHCI, PCHCS, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   820218  DATE WRITTEN
    C   820804  Converted to SLATEC library version.
    C   870813  Updated Reference 2.
    C   890411  Added SAVE statements (Vers. 3.2).
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    C   890703  Corrected category record.  (WRB)
    C   890831  Modified array declarations.  (WRB)
    C   890831  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  PCHIC