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