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