1*DECK PCHIM
2      SUBROUTINE PCHIM (N, X, F, D, INCFD, IERR)
3C***BEGIN PROLOGUE  PCHIM
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 PCHIC if user control is
9C            desired over boundary or switch conditions.)
10C***LIBRARY   SLATEC (PCHIP)
11C***CATEGORY  E1A
12C***TYPE      SINGLE 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          PCHIM:  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 PCHIC 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 PCHIC 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 PCHFE or PCHFD.
41C
42C ----------------------------------------------------------------------
43C
44C  Calling sequence:
45C
46C        PARAMETER  (INCFD = ...)
47C        INTEGER  N, IERR
48C        REAL  X(N), F(INCFD,N), D(INCFD,N)
49C
50C        CALL  PCHIM (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 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 array of dependent variable values to be inter-
63C           polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
64C           PCHIM is designed for monotonic data, but it will work for
65C           any F-array.  It will force extrema at points where mono-
66C           tonicity switches direction.  If some other treatment of
67C           switch points is desired, PCHIC should be used instead.
68C                                     -----
69C     D -- (output) real array of derivative values at the data points.
70C           If the data are monotonic, these values will determine a
71C           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  PCHST, XERMSG
102C***REVISION HISTORY  (YYMMDD)
103C   811103  DATE WRITTEN
104C   820201  1. Introduced  PCHST  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   870813  Updated Reference 1.
113C   890411  Added SAVE statements (Vers. 3.2).
114C   890531  Changed all specific intrinsics to generic.  (WRB)
115C   890703  Corrected category record.  (WRB)
116C   890831  Modified array declarations.  (WRB)
117C   890831  REVISION DATE from Version 3.2
118C   891214  Prologue converted to Version 4.0 format.  (BAB)
119C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
120C   920429  Revised format and order of references.  (WRB,FNF)
121C***END PROLOGUE  PCHIM
122C  Programming notes:
123C
124C     1. The function  PCHST(ARG1,ARG2)  is assumed to return zero if
125C        either argument is zero, +1 if they are of the same sign, and
126C        -1 if they are of opposite sign.
127C     2. To produce a double precision version, simply:
128C        a. Change PCHIM to DPCHIM wherever it occurs,
129C        b. Change PCHST to DPCHST wherever it occurs,
130C        c. Change all references to the Fortran intrinsics to their
131C           double precision equivalents,
132C        d. Change the real declarations to double precision, and
133C        e. Change the constants ZERO and THREE to double precision.
134C
135C  DECLARE ARGUMENTS.
136C
137      INTEGER  N, INCFD, IERR
138      REAL  X(*), F(INCFD,*), D(INCFD,*)
139C
140C  DECLARE LOCAL VARIABLES.
141C
142      INTEGER  I, NLESS1
143      REAL  DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE,
144     *      H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO
145      SAVE ZERO, THREE
146      REAL  PCHST
147      DATA  ZERO /0./,  THREE /3./
148C
149C  VALIDITY-CHECK ARGUMENTS.
150C
151C***FIRST EXECUTABLE STATEMENT  PCHIM
152      IF ( N.LT.2 )  GO TO 5001
153      IF ( INCFD.LT.1 )  GO TO 5002
154      DO 1  I = 2, N
155         IF ( X(I).LE.X(I-1) )  GO TO 5003
156    1 CONTINUE
157C
158C  FUNCTION DEFINITION IS OK, GO ON.
159C
160      IERR = 0
161      NLESS1 = N - 1
162      H1 = X(2) - X(1)
163      DEL1 = (F(1,2) - F(1,1))/H1
164      DSAVE = DEL1
165C
166C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
167C
168      IF (NLESS1 .GT. 1)  GO TO 10
169      D(1,1) = DEL1
170      D(1,N) = DEL1
171      GO TO 5000
172C
173C  NORMAL CASE  (N .GE. 3).
174C
175   10 CONTINUE
176      H2 = X(3) - X(2)
177      DEL2 = (F(1,3) - F(1,2))/H2
178C
179C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
180C     SHAPE-PRESERVING.
181C
182      HSUM = H1 + H2
183      W1 = (H1 + HSUM)/HSUM
184      W2 = -H1/HSUM
185      D(1,1) = W1*DEL1 + W2*DEL2
186      IF ( PCHST(D(1,1),DEL1) .LE. ZERO)  THEN
187         D(1,1) = ZERO
188      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
189C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
190         DMAX = THREE*DEL1
191         IF (ABS(D(1,1)) .GT. ABS(DMAX))  D(1,1) = DMAX
192      ENDIF
193C
194C  LOOP THROUGH INTERIOR POINTS.
195C
196      DO 50  I = 2, NLESS1
197         IF (I .EQ. 2)  GO TO 40
198C
199         H1 = H2
200         H2 = X(I+1) - X(I)
201         HSUM = H1 + H2
202         DEL1 = DEL2
203         DEL2 = (F(1,I+1) - F(1,I))/H2
204   40    CONTINUE
205C
206C        SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
207C
208         D(1,I) = ZERO
209         IF ( PCHST(DEL1,DEL2) )  42, 41, 45
210C
211C        COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY.
212C
213   41    CONTINUE
214         IF (DEL2 .EQ. ZERO)  GO TO 50
215         IF ( PCHST(DSAVE,DEL2) .LT. ZERO)  IERR = IERR + 1
216         DSAVE = DEL2
217         GO TO 50
218C
219   42    CONTINUE
220         IERR = IERR + 1
221         DSAVE = DEL2
222         GO TO 50
223C
224C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
225C
226   45    CONTINUE
227         HSUMT3 = HSUM+HSUM+HSUM
228         W1 = (HSUM + H1)/HSUMT3
229         W2 = (HSUM + H2)/HSUMT3
230         DMAX = MAX( ABS(DEL1), ABS(DEL2) )
231         DMIN = MIN( ABS(DEL1), ABS(DEL2) )
232         DRAT1 = DEL1/DMAX
233         DRAT2 = DEL2/DMAX
234         D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2)
235C
236   50 CONTINUE
237C
238C  SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
239C     SHAPE-PRESERVING.
240C
241      W1 = -H2/HSUM
242      W2 = (H2 + HSUM)/HSUM
243      D(1,N) = W1*DEL1 + W2*DEL2
244      IF ( PCHST(D(1,N),DEL2) .LE. ZERO)  THEN
245         D(1,N) = ZERO
246      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
247C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
248         DMAX = THREE*DEL2
249         IF (ABS(D(1,N)) .GT. ABS(DMAX))  D(1,N) = DMAX
250      ENDIF
251C
252C  NORMAL RETURN.
253C
254 5000 CONTINUE
255      RETURN
256C
257C  ERROR RETURNS.
258C
259 5001 CONTINUE
260C     N.LT.2 RETURN.
261      IERR = -1
262      CALL XERMSG ('SLATEC', 'PCHIM',
263     +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
264      RETURN
265C
266 5002 CONTINUE
267C     INCFD.LT.1 RETURN.
268      IERR = -2
269      CALL XERMSG ('SLATEC', 'PCHIM', 'INCREMENT LESS THAN ONE', IERR,
270     +   1)
271      RETURN
272C
273 5003 CONTINUE
274C     X-ARRAY NOT STRICTLY INCREASING.
275      IERR = -3
276      CALL XERMSG ('SLATEC', 'PCHIM', 'X-ARRAY NOT STRICTLY INCREASING'
277     +   , IERR, 1)
278      RETURN
279C------------- LAST LINE OF PCHIM FOLLOWS ------------------------------
280      END
281