Mercurial > octave-nkf

comparison libcruft/slatec-fn/pchim.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
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