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