text: *DECK DPCHIC
        SUBROUTINE DPCHIC (IC, VC, SWITCH, N, X, F, D, INCFD, WK, NWK,
       +   IERR)
  C***BEGIN PROLOGUE  DPCHIC
  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      DOUBLE 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         DPCHIC:  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 DPCHFE or DPCHFD.
  C
  C ----------------------------------------------------------------------
  C
  C  Calling sequence:
  C
  C        PARAMETER  (INCFD = ...)
  C        INTEGER  IC(2), N, NWK, IERR
  C        DOUBLE PRECISION  VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N),
  C                          WK(NWK)
  C
  C        CALL DPCHIC (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 DPCHIM).
  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 DPCHSP.
  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*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     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 DPCHIM, 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*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 a monotone cubic
  C           Hermite function on each subinterval on which the data
  C           are monotonic, except possibly adjacent to switches in
  C           monotonicity. 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.  The user may
  C           wish 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 Natio-
  C                 nal 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  DPCHCE, DPCHCI, DPCHCS, XERMSG
  C***REVISION HISTORY  (YYMMDD)
  C   820218  DATE WRITTEN
  C   820804  Converted to SLATEC library version.
  C   870707  Corrected XERROR calls for d.p. name(s).
  C   870813  Updated Reference 2.
  C   890206  Corrected XERROR calls.
  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   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  DPCHIC
  C  Programming notes:
  C
  C     To produce a single precision version, simply:
  C        a. Change DPCHIC to PCHIC wherever it occurs,
  C        b. Change DPCHCE to PCHCE wherever it occurs,
  C        c. Change DPCHCI to PCHCI wherever it occurs,
  C        d. Change DPCHCS to PCHCS wherever it occurs,
  C        e. Change the double precision declarations to real, and
  C        f. Change the constant  ZERO  to single precision.
  C
  C  DECLARE ARGUMENTS.
  C
        INTEGER  IC(2), N, INCFD, NWK, IERR
        DOUBLE PRECISION  VC(2), SWITCH, X(*), F(INCFD,*), D(INCFD,*),
       * WK(NWK)
  C
  C  DECLARE LOCAL VARIABLES.
  C
        INTEGER  I, IBEG, IEND, NLESS1
        DOUBLE PRECISION  ZERO
        SAVE ZERO
        DATA  ZERO /0.D0/
  C
  C  VALIDITY-CHECK ARGUMENTS.
  C
  C***FIRST EXECUTABLE STATEMENT  DPCHIC
        IF ( N.LT.2 )  GO TO 5001
        IF ( INCFD.LT.1 )  GO TO 5002
        DO 1  I = 2, N
           IF ( X(I).LE.X(I-1) )  GO TO 5003
      1 CONTINUE
  C
        IBEG = IC(1)
        IEND = IC(2)
        IERR = 0
        IF (ABS(IBEG) .GT. 5)  IERR = IERR - 1
        IF (ABS(IEND) .GT. 5)  IERR = IERR - 2
        IF (IERR .LT. 0)  GO TO 5004
  C
  C  FUNCTION DEFINITION IS OK -- GO ON.
  C
        NLESS1 = N - 1
        IF ( NWK .LT. 2*NLESS1 )  GO TO 5007
  C
  C  SET UP H AND SLOPE ARRAYS.
  C
        DO 20  I = 1, NLESS1
           WK(I) = X(I+1) - X(I)
           WK(NLESS1+I) = (F(1,I+1) - F(1,I)) / WK(I)
     20 CONTINUE
  C
  C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
  C
        IF (NLESS1 .GT. 1)  GO TO 1000
        D(1,1) = WK(2)
        D(1,N) = WK(2)
        GO TO 3000
  C
  C  NORMAL CASE  (N .GE. 3) .
  C
   1000 CONTINUE
  C
  C  SET INTERIOR DERIVATIVES AND DEFAULT END CONDITIONS.
  C
  C     --------------------------------------
        CALL DPCHCI (N, WK(1), WK(N), D, INCFD)
  C     --------------------------------------
  C
  C  SET DERIVATIVES AT POINTS WHERE MONOTONICITY SWITCHES DIRECTION.
  C
        IF (SWITCH .EQ. ZERO)  GO TO 3000
  C     ----------------------------------------------------
        CALL DPCHCS (SWITCH, N, WK(1), WK(N), D, INCFD, IERR)
  C     ----------------------------------------------------
        IF (IERR .NE. 0)  GO TO 5008
  C
  C  SET END CONDITIONS.
  C
   3000 CONTINUE
        IF ( (IBEG.EQ.0) .AND. (IEND.EQ.0) )  GO TO 5000
  C     -------------------------------------------------------
        CALL DPCHCE (IC, VC, N, X, WK(1), WK(N), D, INCFD, IERR)
  C     -------------------------------------------------------
        IF (IERR .LT. 0)  GO TO 5009
  C
  C  NORMAL RETURN.
  C
   5000 CONTINUE
        RETURN
  C
  C  ERROR RETURNS.
  C
   5001 CONTINUE
  C     N.LT.2 RETURN.
        IERR = -1
        CALL XERMSG ('SLATEC', 'DPCHIC',
       +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
        RETURN
  C
   5002 CONTINUE
  C     INCFD.LT.1 RETURN.
        IERR = -2
        CALL XERMSG ('SLATEC', 'DPCHIC', 'INCREMENT LESS THAN ONE', IERR,
       +   1)
        RETURN
  C
   5003 CONTINUE
  C     X-ARRAY NOT STRICTLY INCREASING.
        IERR = -3
        CALL XERMSG ('SLATEC', 'DPCHIC',
       +   'X-ARRAY NOT STRICTLY INCREASING', IERR, 1)
        RETURN
  C
   5004 CONTINUE
  C     IC OUT OF RANGE RETURN.
        IERR = IERR - 3
        CALL XERMSG ('SLATEC', 'DPCHIC', 'IC OUT OF RANGE', IERR, 1)
        RETURN
  C
   5007 CONTINUE
  C     NWK .LT. 2*(N-1)  RETURN.
        IERR = -7
        CALL XERMSG ('SLATEC', 'DPCHIC', 'WORK ARRAY TOO SMALL', IERR, 1)
        RETURN
  C
   5008 CONTINUE
  C     ERROR RETURN FROM DPCHCS.
        IERR = -8
        CALL XERMSG ('SLATEC', 'DPCHIC', 'ERROR RETURN FROM DPCHCS',
       +   IERR, 1)
        RETURN
  C
   5009 CONTINUE
  C     ERROR RETURN FROM DPCHCE.
  C   *** THIS CASE SHOULD NEVER OCCUR ***
        IERR = -9
        CALL XERMSG ('SLATEC', 'DPCHIC', 'ERROR RETURN FROM DPCHCE',
       +   IERR, 1)
        RETURN
  C------------- LAST LINE OF DPCHIC FOLLOWS -----------------------------
        END
  