GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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