1*DECK DPCHIM
2      SUBROUTINE DPCHIM (N, X, F, D, INCFD, IERR)
3C***BEGIN PROLOGUE  DPCHIM
4C***PURPOSE  Set derivatives needed to determine a monotone piecewise
5C            cubic Hermite interpolant to given data.  Boundary values
6C            are provided which are compatible with monotonicity.  The
7C            interpolant will have an extremum at each point where mono-
8C            tonicity switches direction.  (See DPCHIC if user control
9C            is desired over boundary or switch conditions.)
10C***LIBRARY   SLATEC (PCHIP)
11C***CATEGORY  E1A
12C***TYPE      DOUBLE PRECISION (PCHIM-S, DPCHIM-D)
13C***KEYWORDS  CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION,
14C             PCHIP, PIECEWISE CUBIC INTERPOLATION
15C***AUTHOR  Fritsch, F. N., (LLNL)
16C             Lawrence Livermore National Laboratory
17C             P.O. Box 808  (L-316)
18C             Livermore, CA  94550
19C             FTS 532-4275, (510) 422-4275
20C***DESCRIPTION
21C
22C          DPCHIM:  Piecewise Cubic Hermite Interpolation to
23C                  Monotone data.
24C
25C     Sets derivatives needed to determine a monotone piecewise cubic
26C     Hermite interpolant to the data given in X and F.
27C
28C     Default boundary conditions are provided which are compatible
29C     with monotonicity.  (See DPCHIC if user control of boundary con-
30C     ditions is desired.)
31C
32C     If the data are only piecewise monotonic, the interpolant will
33C     have an extremum at each point where monotonicity switches direc-
34C     tion.  (See DPCHIC if user control is desired in such cases.)
35C
36C     To facilitate two-dimensional applications, includes an increment
37C     between successive values of the F- and D-arrays.
38C
39C     The resulting piecewise cubic Hermite function may be evaluated
40C     by DPCHFE or DPCHFD.
41C
42C ----------------------------------------------------------------------
43C
44C  Calling sequence:
45C
46C        PARAMETER  (INCFD = ...)
47C        INTEGER  N, IERR
48C        DOUBLE PRECISION  X(N), F(INCFD,N), D(INCFD,N)
49C
50C        CALL  DPCHIM (N, X, F, D, INCFD, IERR)
51C
52C   Parameters:
53C
54C     N -- (input) number of data points.  (Error return if N.LT.2 .)
55C           If N=2, simply does linear interpolation.
56C
57C     X -- (input) real*8 array of independent variable values.  The
58C           elements of X must be strictly increasing:
59C                X(I-1) .LT. X(I),  I = 2(1)N.
60C           (Error return if not.)
61C
62C     F -- (input) real*8 array of dependent variable values to be
63C           interpolated.  F(1+(I-1)*INCFD) is value corresponding to
64C           X(I).  DPCHIM is designed for monotonic data, but it will
65C           work for any F-array.  It will force extrema at points where
66C           monotonicity switches direction.  If some other treatment of
67C           switch points is desired, DPCHIC should be used instead.
68C                                     -----
69C     D -- (output) real*8 array of derivative values at the data
70C           points.  If the data are monotonic, these values will
71C           determine a monotone cubic Hermite function.
72C           The value corresponding to X(I) is stored in
73C                D(1+(I-1)*INCFD),  I=1(1)N.
74C           No other entries in D are changed.
75C
76C     INCFD -- (input) increment between successive values in F and D.
77C           This argument is provided primarily for 2-D applications.
78C           (Error return if  INCFD.LT.1 .)
79C
80C     IERR -- (output) error flag.
81C           Normal return:
82C              IERR = 0  (no errors).
83C           Warning error:
84C              IERR.GT.0  means that IERR switches in the direction
85C                 of monotonicity were detected.
86C           "Recoverable" errors:
87C              IERR = -1  if N.LT.2 .
88C              IERR = -2  if INCFD.LT.1 .
89C              IERR = -3  if the X-array is not strictly increasing.
90C             (The D-array has not been changed in any of these cases.)
91C               NOTE:  The above errors are checked in the order listed,
92C                   and following arguments have **NOT** been validated.
93C
94C***REFERENCES  1. F. N. Fritsch and J. Butland, A method for construc-
95C                 ting local monotone piecewise cubic interpolants, SIAM
96C                 Journal on Scientific and Statistical Computing 5, 2
97C                 (June 1984), pp. 300-304.
98C               2. F. N. Fritsch and R. E. Carlson, Monotone piecewise
99C                 cubic interpolation, SIAM Journal on Numerical Ana-
100C                 lysis 17, 2 (April 1980), pp. 238-246.
101C***ROUTINES CALLED  DPCHST, XERMSG
102C***REVISION HISTORY  (YYMMDD)
103C   811103  DATE WRITTEN
104C   820201  1. Introduced  DPCHST  to reduce possible over/under-
105C             flow problems.
106C           2. Rearranged derivative formula for same reason.
107C   820602  1. Modified end conditions to be continuous functions
108C             of data when monotonicity switches in next interval.
109C           2. Modified formulas so end conditions are less prone
110C             of over/underflow problems.
111C   820803  Minor cosmetic changes for release 1.
112C   870707  Corrected XERROR calls for d.p. name(s).
113C   870813  Updated Reference 1.
114C   890206  Corrected XERROR calls.
115C   890411  Added SAVE statements (Vers. 3.2).
116C   890531  Changed all specific intrinsics to generic.  (WRB)
117C   890703  Corrected category record.  (WRB)
118C   890831  Modified array declarations.  (WRB)
119C   891006  Cosmetic changes to prologue.  (WRB)
120C   891006  REVISION DATE from Version 3.2
121C   891214  Prologue converted to Version 4.0 format.  (BAB)
122C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
123C   920429  Revised format and order of references.  (WRB,FNF)
124C***END PROLOGUE  DPCHIM
125C  Programming notes:
126C
127C     1. The function  DPCHST(ARG1,ARG2)  is assumed to return zero if
128C        either argument is zero, +1 if they are of the same sign, and
129C        -1 if they are of opposite sign.
130C     2. To produce a single precision version, simply:
131C        a. Change DPCHIM to PCHIM wherever it occurs,
132C        b. Change DPCHST to PCHST wherever it occurs,
133C        c. Change all references to the Fortran intrinsics to their
134C           single precision equivalents,
135C        d. Change the double precision declarations to real, and
136C        e. Change the constants ZERO and THREE to single precision.
137C
138C  DECLARE ARGUMENTS.
139C
140      INTEGER  N, INCFD, IERR
141      DOUBLE PRECISION  X(*), F(INCFD,*), D(INCFD,*)
142C
143C  DECLARE LOCAL VARIABLES.
144C
145      INTEGER  I, NLESS1
146      DOUBLE PRECISION  DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE,
147     *      H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO
148      SAVE ZERO, THREE
149      DOUBLE PRECISION  DPCHST
150      DATA  ZERO /0.D0/, THREE/3.D0/
151C
152C  VALIDITY-CHECK ARGUMENTS.
153C
154C***FIRST EXECUTABLE STATEMENT  DPCHIM
155      IF ( N.LT.2 )  GO TO 5001
156      IF ( INCFD.LT.1 )  GO TO 5002
157      DO 1  I = 2, N
158         IF ( X(I).LE.X(I-1) )  GO TO 5003
159    1 CONTINUE
160C
161C  FUNCTION DEFINITION IS OK, GO ON.
162C
163      IERR = 0
164      NLESS1 = N - 1
165      H1 = X(2) - X(1)
166      DEL1 = (F(1,2) - F(1,1))/H1
167      DSAVE = DEL1
168C
169C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
170C
171      IF (NLESS1 .GT. 1)  GO TO 10
172      D(1,1) = DEL1
173      D(1,N) = DEL1
174      GO TO 5000
175C
176C  NORMAL CASE  (N .GE. 3).
177C
178   10 CONTINUE
179      H2 = X(3) - X(2)
180      DEL2 = (F(1,3) - F(1,2))/H2
181C
182C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
183C     SHAPE-PRESERVING.
184C
185      HSUM = H1 + H2
186      W1 = (H1 + HSUM)/HSUM
187      W2 = -H1/HSUM
188      D(1,1) = W1*DEL1 + W2*DEL2
189      IF ( DPCHST(D(1,1),DEL1) .LE. ZERO)  THEN
190         D(1,1) = ZERO
191      ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO)  THEN
192C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
193         DMAX = THREE*DEL1
194         IF (ABS(D(1,1)) .GT. ABS(DMAX))  D(1,1) = DMAX
195      ENDIF
196C
197C  LOOP THROUGH INTERIOR POINTS.
198C
199      DO 50  I = 2, NLESS1
200         IF (I .EQ. 2)  GO TO 40
201C
202         H1 = H2
203         H2 = X(I+1) - X(I)
204         HSUM = H1 + H2
205         DEL1 = DEL2
206         DEL2 = (F(1,I+1) - F(1,I))/H2
207   40    CONTINUE
208C
209C        SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
210C
211         D(1,I) = ZERO
212         IF ( DPCHST(DEL1,DEL2) .LT. 0.)  GO TO 42
213         IF ( DPCHST(DEL1,DEL2) .EQ. 0.)  GO TO 41
214         GO TO 45
215C
216C        COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY.
217C
218   41    CONTINUE
219         IF (DEL2 .EQ. ZERO)  GO TO 50
220         IF ( DPCHST(DSAVE,DEL2) .LT. ZERO)  IERR = IERR + 1
221         DSAVE = DEL2
222         GO TO 50
223C
224   42    CONTINUE
225         IERR = IERR + 1
226         DSAVE = DEL2
227         GO TO 50
228C
229C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
230C
231   45    CONTINUE
232         HSUMT3 = HSUM+HSUM+HSUM
233         W1 = (HSUM + H1)/HSUMT3
234         W2 = (HSUM + H2)/HSUMT3
235         DMAX = MAX( ABS(DEL1), ABS(DEL2) )
236         DMIN = MIN( ABS(DEL1), ABS(DEL2) )
237         DRAT1 = DEL1/DMAX
238         DRAT2 = DEL2/DMAX
239         D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2)
240C
241   50 CONTINUE
242C
243C  SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
244C     SHAPE-PRESERVING.
245C
246      W1 = -H2/HSUM
247      W2 = (H2 + HSUM)/HSUM
248      D(1,N) = W1*DEL1 + W2*DEL2
249      IF ( DPCHST(D(1,N),DEL2) .LE. ZERO)  THEN
250         D(1,N) = ZERO
251      ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO)  THEN
252C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
253         DMAX = THREE*DEL2
254         IF (ABS(D(1,N)) .GT. ABS(DMAX))  D(1,N) = DMAX
255      ENDIF
256C
257C  NORMAL RETURN.
258C
259 5000 CONTINUE
260      RETURN
261C
262C  ERROR RETURNS.
263C
264 5001 CONTINUE
265C     N.LT.2 RETURN.
266      IERR = -1
267      CALL XERMSG ('SLATEC', 'DPCHIM',
268     +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
269      RETURN
270C
271 5002 CONTINUE
272C     INCFD.LT.1 RETURN.
273      IERR = -2
274      CALL XERMSG ('SLATEC', 'DPCHIM', 'INCREMENT LESS THAN ONE', IERR,
275     +   1)
276      RETURN
277C
278 5003 CONTINUE
279C     X-ARRAY NOT STRICTLY INCREASING.
280      IERR = -3
281      CALL XERMSG ('SLATEC', 'DPCHIM',
282     +   'X-ARRAY NOT STRICTLY INCREASING', IERR, 1)
283      RETURN
284C------------- LAST LINE OF DPCHIM FOLLOWS -----------------------------
285      END
286