1*DECK DPCHCI
2      SUBROUTINE DPCHCI (N, H, SLOPE, D, INCFD)
3C***BEGIN PROLOGUE  DPCHCI
4C***SUBSIDIARY
5C***PURPOSE  Set interior derivatives for DPCHIC
6C***LIBRARY   SLATEC (PCHIP)
7C***TYPE      DOUBLE PRECISION (PCHCI-S, DPCHCI-D)
8C***AUTHOR  Fritsch, F. N., (LLNL)
9C***DESCRIPTION
10C
11C          DPCHCI:  DPCHIC Initial Derivative Setter.
12C
13C    Called by DPCHIC 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 DPCHIM.
26C
27C ----------------------------------------------------------------------
28C
29C  Calling sequence:
30C
31C        PARAMETER  (INCFD = ...)
32C        INTEGER  N
33C        DOUBLE PRECISION  H(N), SLOPE(N), D(INCFD,N)
34C
35C        CALL  DPCHCI (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*8 array of interval lengths.
43C     SLOPE -- (input) real*8 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*8 array of derivative values at 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  DPCHIC
65C***ROUTINES CALLED  DPCHST
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   870813  Minor cosmetic changes.
75C   890411  Added SAVE statements (Vers. 3.2).
76C   890531  Changed all specific intrinsics to generic.  (WRB)
77C   890831  Modified array declarations.  (WRB)
78C   890831  REVISION DATE from Version 3.2
79C   891214  Prologue converted to Version 4.0 format.  (BAB)
80C   900328  Added TYPE section.  (WRB)
81C   910408  Updated AUTHOR section in prologue.  (WRB)
82C   930503  Improved purpose.  (FNF)
83C***END PROLOGUE  DPCHCI
84C
85C  Programming notes:
86C     1. The function  DPCHST(ARG1,ARG2)  is assumed to return zero if
87C        either argument is zero, +1 if they are of the same sign, and
88C        -1 if they are of opposite sign.
89C**End
90C
91C  DECLARE ARGUMENTS.
92C
93      INTEGER  N, INCFD
94      DOUBLE PRECISION  H(*), SLOPE(*), D(INCFD,*)
95C
96C  DECLARE LOCAL VARIABLES.
97C
98      INTEGER  I, NLESS1
99      DOUBLE PRECISION  DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, HSUM,
100     *      HSUMT3, THREE, W1, W2, ZERO
101      SAVE ZERO, THREE
102      DOUBLE PRECISION  DPCHST
103C
104C  INITIALIZE.
105C
106      DATA  ZERO /0.D0/, THREE/3.D0/
107C***FIRST EXECUTABLE STATEMENT  DPCHCI
108      NLESS1 = N - 1
109      DEL1 = SLOPE(1)
110C
111C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
112C
113      IF (NLESS1 .GT. 1)  GO TO 10
114      D(1,1) = DEL1
115      D(1,N) = DEL1
116      GO TO 5000
117C
118C  NORMAL CASE  (N .GE. 3).
119C
120   10 CONTINUE
121      DEL2 = SLOPE(2)
122C
123C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
124C     SHAPE-PRESERVING.
125C
126      HSUM = H(1) + H(2)
127      W1 = (H(1) + HSUM)/HSUM
128      W2 = -H(1)/HSUM
129      D(1,1) = W1*DEL1 + W2*DEL2
130      IF ( DPCHST(D(1,1),DEL1) .LE. ZERO)  THEN
131         D(1,1) = ZERO
132      ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO)  THEN
133C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
134         DMAX = THREE*DEL1
135         IF (ABS(D(1,1)) .GT. ABS(DMAX))  D(1,1) = DMAX
136      ENDIF
137C
138C  LOOP THROUGH INTERIOR POINTS.
139C
140      DO 50  I = 2, NLESS1
141         IF (I .EQ. 2)  GO TO 40
142C
143         HSUM = H(I-1) + H(I)
144         DEL1 = DEL2
145         DEL2 = SLOPE(I)
146   40    CONTINUE
147C
148C        SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
149C
150         D(1,I) = ZERO
151         IF ( DPCHST(DEL1,DEL2) .LE. ZERO)  GO TO 50
152C
153C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
154C
155         HSUMT3 = HSUM+HSUM+HSUM
156         W1 = (HSUM + H(I-1))/HSUMT3
157         W2 = (HSUM + H(I)  )/HSUMT3
158         DMAX = MAX( ABS(DEL1), ABS(DEL2) )
159         DMIN = MIN( ABS(DEL1), ABS(DEL2) )
160         DRAT1 = DEL1/DMAX
161         DRAT2 = DEL2/DMAX
162         D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2)
163C
164   50 CONTINUE
165C
166C  SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
167C     SHAPE-PRESERVING.
168C
169      W1 = -H(N-1)/HSUM
170      W2 = (H(N-1) + HSUM)/HSUM
171      D(1,N) = W1*DEL1 + W2*DEL2
172      IF ( DPCHST(D(1,N),DEL2) .LE. ZERO)  THEN
173         D(1,N) = ZERO
174      ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO)  THEN
175C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
176         DMAX = THREE*DEL2
177         IF (ABS(D(1,N)) .GT. ABS(DMAX))  D(1,N) = DMAX
178      ENDIF
179C
180C  NORMAL RETURN.
181C
182 5000 CONTINUE
183      RETURN
184C------------- LAST LINE OF DPCHCI FOLLOWS -----------------------------
185      END
186