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