1*DECK PCHCI
2      SUBROUTINE PCHCI (N, H, SLOPE, D, INCFD)
3C***BEGIN PROLOGUE  PCHCI
4C***SUBSIDIARY
5C***PURPOSE  Set interior derivatives for PCHIC
6C***LIBRARY   SLATEC (PCHIP)
7C***TYPE      SINGLE PRECISION (PCHCI-S, DPCHCI-D)
8C***AUTHOR  Fritsch, F. N., (LLNL)
9C***DESCRIPTION
10C
11C          PCHCI:  PCHIC Initial Derivative Setter.
12C
13C    Called by PCHIC to set derivatives needed to determine a monotone
14C    piecewise cubic Hermite interpolant to the data.
15C
16C    Default boundary conditions are provided which are compatible
17C    with monotonicity.  If the data are only piecewise monotonic, the
18C    interpolant will have an extremum at each point where monotonicity
19C    switches direction.
20C
21C    To facilitate two-dimensional applications, includes an increment
22C    between successive values of the D-array.
23C
24C    The resulting piecewise cubic Hermite function should be identical
25C    (within roundoff error) to that produced by PCHIM.
26C
27C ----------------------------------------------------------------------
28C
29C  Calling sequence:
30C
31C        PARAMETER  (INCFD = ...)
32C        INTEGER  N
33C        REAL  H(N), SLOPE(N), D(INCFD,N)
34C
35C        CALL  PCHCI (N, H, SLOPE, D, INCFD)
36C
37C   Parameters:
38C
39C     N -- (input) number of data points.
40C           If N=2, simply does linear interpolation.
41C
42C     H -- (input) real array of interval lengths.
43C     SLOPE -- (input) real array of data slopes.
44C           If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:
45C                  H(I) =  X(I+1)-X(I),
46C              SLOPE(I) = (Y(I+1)-Y(I))/H(I),  I=1(1)N-1.
47C
48C     D -- (output) real array of derivative values at the data points.
49C           If the data are monotonic, these values will determine a
50C           a monotone cubic Hermite function.
51C           The value corresponding to X(I) is stored in
52C                D(1+(I-1)*INCFD),  I=1(1)N.
53C           No other entries in D are changed.
54C
55C     INCFD -- (input) increment between successive values in D.
56C           This argument is provided primarily for 2-D applications.
57C
58C    -------
59C    WARNING:  This routine does no validity-checking of arguments.
60C    -------
61C
62C  Fortran intrinsics used:  ABS, MAX, MIN.
63C
64C***SEE ALSO  PCHIC
65C***ROUTINES CALLED  PCHST
66C***REVISION HISTORY  (YYMMDD)
67C   820218  DATE WRITTEN
68C   820601  Modified end conditions to be continuous functions of
69C           data when monotonicity switches in next interval.
70C   820602  1. Modified formulas so end conditions are less prone
71C             to over/underflow problems.
72C           2. Minor modification to HSUM calculation.
73C   820805  Converted to SLATEC library version.
74C   890411  Added SAVE statements (Vers. 3.2).
75C   890531  Changed all specific intrinsics to generic.  (WRB)
76C   890831  Modified array declarations.  (WRB)
77C   890831  REVISION DATE from Version 3.2
78C   891214  Prologue converted to Version 4.0 format.  (BAB)
79C   900328  Added TYPE section.  (WRB)
80C   910408  Updated AUTHOR section in prologue.  (WRB)
81C   930503  Improved purpose.  (FNF)
82C***END PROLOGUE  PCHCI
83C
84C  Programming notes:
85C     1. The function  PCHST(ARG1,ARG2)  is assumed to return zero if
86C        either argument is zero, +1 if they are of the same sign, and
87C        -1 if they are of opposite sign.
88C**End
89C
90C  DECLARE ARGUMENTS.
91C
92      INTEGER  N, INCFD
93      REAL  H(*), SLOPE(*), D(INCFD,*)
94C
95C  DECLARE LOCAL VARIABLES.
96C
97      INTEGER  I, NLESS1
98      REAL  DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, HSUM, HSUMT3, THREE,
99     *      W1, W2, ZERO
100      SAVE ZERO, THREE
101      REAL  PCHST
102C
103C  INITIALIZE.
104C
105      DATA  ZERO /0./,  THREE /3./
106C***FIRST EXECUTABLE STATEMENT  PCHCI
107      NLESS1 = N - 1
108      DEL1 = SLOPE(1)
109C
110C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
111C
112      IF (NLESS1 .GT. 1)  GO TO 10
113      D(1,1) = DEL1
114      D(1,N) = DEL1
115      GO TO 5000
116C
117C  NORMAL CASE  (N .GE. 3).
118C
119   10 CONTINUE
120      DEL2 = SLOPE(2)
121C
122C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
123C     SHAPE-PRESERVING.
124C
125      HSUM = H(1) + H(2)
126      W1 = (H(1) + HSUM)/HSUM
127      W2 = -H(1)/HSUM
128      D(1,1) = W1*DEL1 + W2*DEL2
129      IF ( PCHST(D(1,1),DEL1) .LE. ZERO)  THEN
130         D(1,1) = ZERO
131      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
132C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
133         DMAX = THREE*DEL1
134         IF (ABS(D(1,1)) .GT. ABS(DMAX))  D(1,1) = DMAX
135      ENDIF
136C
137C  LOOP THROUGH INTERIOR POINTS.
138C
139      DO 50  I = 2, NLESS1
140         IF (I .EQ. 2)  GO TO 40
141C
142         HSUM = H(I-1) + H(I)
143         DEL1 = DEL2
144         DEL2 = SLOPE(I)
145   40    CONTINUE
146C
147C        SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
148C
149         D(1,I) = ZERO
150         IF ( PCHST(DEL1,DEL2) .LE. ZERO)  GO TO 50
151C
152C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
153C
154         HSUMT3 = HSUM+HSUM+HSUM
155         W1 = (HSUM + H(I-1))/HSUMT3
156         W2 = (HSUM + H(I)  )/HSUMT3
157         DMAX = MAX( ABS(DEL1), ABS(DEL2) )
158         DMIN = MIN( ABS(DEL1), ABS(DEL2) )
159         DRAT1 = DEL1/DMAX
160         DRAT2 = DEL2/DMAX
161         D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2)
162C
163   50 CONTINUE
164C
165C  SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
166C     SHAPE-PRESERVING.
167C
168      W1 = -H(N-1)/HSUM
169      W2 = (H(N-1) + HSUM)/HSUM
170      D(1,N) = W1*DEL1 + W2*DEL2
171      IF ( PCHST(D(1,N),DEL2) .LE. ZERO)  THEN
172         D(1,N) = ZERO
173      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
174C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
175         DMAX = THREE*DEL2
176         IF (ABS(D(1,N)) .GT. ABS(DMAX))  D(1,N) = DMAX
177      ENDIF
178C
179C  NORMAL RETURN.
180C
181 5000 CONTINUE
182      RETURN
183C------------- LAST LINE OF PCHCI FOLLOWS ------------------------------
184      END
185