GCC Code Coverage Report


Directory: ./
File: misc/pchsp.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 111 0.0%
Branches: 0 80 0.0%

Line Branch Exec Source
1 *DECK PCHSP
2 SUBROUTINE PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
3 C***BEGIN PROLOGUE PCHSP
4 C***PURPOSE Set derivatives needed to determine the Hermite represen-
5 C tation of the cubic spline interpolant to given data, with
6 C specified boundary conditions.
7 C***LIBRARY SLATEC (PCHIP)
8 C***CATEGORY E1A
9 C***TYPE SINGLE PRECISION (PCHSP-S, DPCHSP-D)
10 C***KEYWORDS CUBIC HERMITE INTERPOLATION, PCHIP,
11 C PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION
12 C***AUTHOR Fritsch, F. N., (LLNL)
13 C Lawrence Livermore National Laboratory
14 C P.O. Box 808 (L-316)
15 C Livermore, CA 94550
16 C FTS 532-4275, (510) 422-4275
17 C***DESCRIPTION
18 C
19 C PCHSP: Piecewise Cubic Hermite Spline
20 C
21 C Computes the Hermite representation of the cubic spline inter-
22 C polant to the data given in X and F satisfying the boundary
23 C conditions specified by IC and VC.
24 C
25 C To facilitate two-dimensional applications, includes an increment
26 C between successive values of the F- and D-arrays.
27 C
28 C The resulting piecewise cubic Hermite function may be evaluated
29 C by PCHFE or PCHFD.
30 C
31 C NOTE: This is a modified version of C. de Boor's cubic spline
32 C routine CUBSPL.
33 C
34 C ----------------------------------------------------------------------
35 C
36 C Calling sequence:
37 C
38 C PARAMETER (INCFD = ...)
39 C INTEGER IC(2), N, NWK, IERR
40 C REAL VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
41 C
42 C CALL PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
43 C
44 C Parameters:
45 C
46 C IC -- (input) integer array of length 2 specifying desired
47 C boundary conditions:
48 C IC(1) = IBEG, desired condition at beginning of data.
49 C IC(2) = IEND, desired condition at end of data.
50 C
51 C IBEG = 0 to set D(1) so that the third derivative is con-
52 C tinuous at X(2). This is the "not a knot" condition
53 C provided by de Boor's cubic spline routine CUBSPL.
54 C < This is the default boundary condition. >
55 C IBEG = 1 if first derivative at X(1) is given in VC(1).
56 C IBEG = 2 if second derivative at X(1) is given in VC(1).
57 C IBEG = 3 to use the 3-point difference formula for D(1).
58 C (Reverts to the default b.c. if N.LT.3 .)
59 C IBEG = 4 to use the 4-point difference formula for D(1).
60 C (Reverts to the default b.c. if N.LT.4 .)
61 C NOTES:
62 C 1. An error return is taken if IBEG is out of range.
63 C 2. For the "natural" boundary condition, use IBEG=2 and
64 C VC(1)=0.
65 C
66 C IEND may take on the same values as IBEG, but applied to
67 C derivative at X(N). In case IEND = 1 or 2, the value is
68 C given in VC(2).
69 C
70 C NOTES:
71 C 1. An error return is taken if IEND is out of range.
72 C 2. For the "natural" boundary condition, use IEND=2 and
73 C VC(2)=0.
74 C
75 C VC -- (input) real array of length 2 specifying desired boundary
76 C values, as indicated above.
77 C VC(1) need be set only if IC(1) = 1 or 2 .
78 C VC(2) need be set only if IC(2) = 1 or 2 .
79 C
80 C N -- (input) number of data points. (Error return if N.LT.2 .)
81 C
82 C X -- (input) real array of independent variable values. The
83 C elements of X must be strictly increasing:
84 C X(I-1) .LT. X(I), I = 2(1)N.
85 C (Error return if not.)
86 C
87 C F -- (input) real array of dependent variable values to be inter-
88 C polated. F(1+(I-1)*INCFD) is value corresponding to X(I).
89 C
90 C D -- (output) real array of derivative values at the data points.
91 C These values will determine the cubic spline interpolant
92 C with the requested boundary conditions.
93 C The value corresponding to X(I) is stored in
94 C D(1+(I-1)*INCFD), I=1(1)N.
95 C No other entries in D are changed.
96 C
97 C INCFD -- (input) increment between successive values in F and D.
98 C This argument is provided primarily for 2-D applications.
99 C (Error return if INCFD.LT.1 .)
100 C
101 C WK -- (scratch) real array of working storage.
102 C
103 C NWK -- (input) length of work array.
104 C (Error return if NWK.LT.2*N .)
105 C
106 C IERR -- (output) error flag.
107 C Normal return:
108 C IERR = 0 (no errors).
109 C "Recoverable" errors:
110 C IERR = -1 if N.LT.2 .
111 C IERR = -2 if INCFD.LT.1 .
112 C IERR = -3 if the X-array is not strictly increasing.
113 C IERR = -4 if IBEG.LT.0 or IBEG.GT.4 .
114 C IERR = -5 if IEND.LT.0 of IEND.GT.4 .
115 C IERR = -6 if both of the above are true.
116 C IERR = -7 if NWK is too small.
117 C NOTE: The above errors are checked in the order listed,
118 C and following arguments have **NOT** been validated.
119 C (The D-array has not been changed in any of these cases.)
120 C IERR = -8 in case of trouble solving the linear system
121 C for the interior derivative values.
122 C (The D-array may have been changed in this case.)
123 C ( Do **NOT** use it! )
124 C
125 C***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer-
126 C Verlag, New York, 1978, pp. 53-59.
127 C***ROUTINES CALLED PCHDF, XERMSG
128 C***REVISION HISTORY (YYMMDD)
129 C 820503 DATE WRITTEN
130 C 820804 Converted to SLATEC library version.
131 C 870707 Minor cosmetic changes to prologue.
132 C 890411 Added SAVE statements (Vers. 3.2).
133 C 890703 Corrected category record. (WRB)
134 C 890831 Modified array declarations. (WRB)
135 C 890831 REVISION DATE from Version 3.2
136 C 891214 Prologue converted to Version 4.0 format. (BAB)
137 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
138 C 920429 Revised format and order of references. (WRB,FNF)
139 C***END PROLOGUE PCHSP
140 C Programming notes:
141 C
142 C To produce a double precision version, simply:
143 C a. Change PCHSP to DPCHSP wherever it occurs,
144 C b. Change the real declarations to double precision, and
145 C c. Change the constants ZERO, HALF, ... to double precision.
146 C
147 C DECLARE ARGUMENTS.
148 C
149 INTEGER IC(2), N, INCFD, NWK, IERR
150 REAL VC(2), X(*), F(INCFD,*), D(INCFD,*), WK(2,*)
151 C
152 C DECLARE LOCAL VARIABLES.
153 C
154 INTEGER IBEG, IEND, INDEX, J, NM1
155 REAL G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO
156 SAVE ZERO, HALF, ONE, TWO, THREE
157 REAL PCHDF
158 C
159 DATA ZERO /0./, HALF /0.5/, ONE /1./, TWO /2./, THREE /3./
160 C
161 C VALIDITY-CHECK ARGUMENTS.
162 C
163 C***FIRST EXECUTABLE STATEMENT PCHSP
164 IF ( N.LT.2 ) GO TO 5001
165 IF ( INCFD.LT.1 ) GO TO 5002
166 DO 1 J = 2, N
167 IF ( X(J).LE.X(J-1) ) GO TO 5003
168 1 CONTINUE
169 C
170 IBEG = IC(1)
171 IEND = IC(2)
172 IERR = 0
173 IF ( (IBEG.LT.0).OR.(IBEG.GT.4) ) IERR = IERR - 1
174 IF ( (IEND.LT.0).OR.(IEND.GT.4) ) IERR = IERR - 2
175 IF ( IERR.LT.0 ) GO TO 5004
176 C
177 C FUNCTION DEFINITION IS OK -- GO ON.
178 C
179 IF ( NWK .LT. 2*N ) GO TO 5007
180 C
181 C COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
182 C COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
183 DO 5 J=2,N
184 WK(1,J) = X(J) - X(J-1)
185 WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J)
186 5 CONTINUE
187 C
188 C SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
189 C
190 IF ( IBEG.GT.N ) IBEG = 0
191 IF ( IEND.GT.N ) IEND = 0
192 C
193 C SET UP FOR BOUNDARY CONDITIONS.
194 C
195 IF ( (IBEG.EQ.1).OR.(IBEG.EQ.2) ) THEN
196 D(1,1) = VC(1)
197 ELSE IF (IBEG .GT. 2) THEN
198 C PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
199 DO 10 J = 1, IBEG
200 INDEX = IBEG-J+1
201 C INDEX RUNS FROM IBEG DOWN TO 1.
202 XTEMP(J) = X(INDEX)
203 IF (J .LT. IBEG) STEMP(J) = WK(2,INDEX)
204 10 CONTINUE
205 C --------------------------------
206 D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
207 C --------------------------------
208 IF (IERR .NE. 0) GO TO 5009
209 IBEG = 1
210 ENDIF
211 C
212 IF ( (IEND.EQ.1).OR.(IEND.EQ.2) ) THEN
213 D(1,N) = VC(2)
214 ELSE IF (IEND .GT. 2) THEN
215 C PICK UP LAST IEND POINTS.
216 DO 15 J = 1, IEND
217 INDEX = N-IEND+J
218 C INDEX RUNS FROM N+1-IEND UP TO N.
219 XTEMP(J) = X(INDEX)
220 IF (J .LT. IEND) STEMP(J) = WK(2,INDEX+1)
221 15 CONTINUE
222 C --------------------------------
223 D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR)
224 C --------------------------------
225 IF (IERR .NE. 0) GO TO 5009
226 IEND = 1
227 ENDIF
228 C
229 C --------------------( BEGIN CODING FROM CUBSPL )--------------------
230 C
231 C **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
232 C F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
233 C INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
234 C WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
235 C
236 C CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
237 C WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
238 C
239 IF (IBEG .EQ. 0) THEN
240 IF (N .EQ. 2) THEN
241 C NO CONDITION AT LEFT END AND N = 2.
242 WK(2,1) = ONE
243 WK(1,1) = ONE
244 D(1,1) = TWO*WK(2,2)
245 ELSE
246 C NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.
247 WK(2,1) = WK(1,3)
248 WK(1,1) = WK(1,2) + WK(1,3)
249 D(1,1) =((WK(1,2) + TWO*WK(1,1))*WK(2,2)*WK(1,3)
250 * + WK(1,2)**2*WK(2,3)) / WK(1,1)
251 ENDIF
252 ELSE IF (IBEG .EQ. 1) THEN
253 C SLOPE PRESCRIBED AT LEFT END.
254 WK(2,1) = ONE
255 WK(1,1) = ZERO
256 ELSE
257 C SECOND DERIVATIVE PRESCRIBED AT LEFT END.
258 WK(2,1) = TWO
259 WK(1,1) = ONE
260 D(1,1) = THREE*WK(2,2) - HALF*WK(1,2)*D(1,1)
261 ENDIF
262 C
263 C IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
264 C CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
265 C EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
266 C
267 NM1 = N-1
268 IF (NM1 .GT. 1) THEN
269 DO 20 J=2,NM1
270 IF (WK(2,J-1) .EQ. ZERO) GO TO 5008
271 G = -WK(1,J+1)/WK(2,J-1)
272 D(1,J) = G*D(1,J-1)
273 * + THREE*(WK(1,J)*WK(2,J+1) + WK(1,J+1)*WK(2,J))
274 WK(2,J) = G*WK(1,J-1) + TWO*(WK(1,J) + WK(1,J+1))
275 20 CONTINUE
276 ENDIF
277 C
278 C CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
279 C (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
280 C
281 C IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
282 C SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
283 C AT THIS POINT.
284 IF (IEND .EQ. 1) GO TO 30
285 C
286 IF (IEND .EQ. 0) THEN
287 IF (N.EQ.2 .AND. IBEG.EQ.0) THEN
288 C NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.
289 D(1,2) = WK(2,2)
290 GO TO 30
291 ELSE IF ((N.EQ.2) .OR. (N.EQ.3 .AND. IBEG.EQ.0)) THEN
292 C EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
293 C NOT-A-KNOT AT LEFT END POINT).
294 D(1,N) = TWO*WK(2,N)
295 WK(2,N) = ONE
296 IF (WK(2,N-1) .EQ. ZERO) GO TO 5008
297 G = -ONE/WK(2,N-1)
298 ELSE
299 C NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A-
300 C KNOT AT LEFT END POINT.
301 G = WK(1,N-1) + WK(1,N)
302 C DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES).
303 D(1,N) = ((WK(1,N)+TWO*G)*WK(2,N)*WK(1,N-1)
304 * + WK(1,N)**2*(F(1,N-1)-F(1,N-2))/WK(1,N-1))/G
305 IF (WK(2,N-1) .EQ. ZERO) GO TO 5008
306 G = -G/WK(2,N-1)
307 WK(2,N) = WK(1,N-1)
308 ENDIF
309 ELSE
310 C SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.
311 D(1,N) = THREE*WK(2,N) + HALF*WK(1,N)*D(1,N)
312 WK(2,N) = TWO
313 IF (WK(2,N-1) .EQ. ZERO) GO TO 5008
314 G = -ONE/WK(2,N-1)
315 ENDIF
316 C
317 C COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
318 C
319 WK(2,N) = G*WK(1,N-1) + WK(2,N)
320 IF (WK(2,N) .EQ. ZERO) GO TO 5008
321 D(1,N) = (G*D(1,N-1) + D(1,N))/WK(2,N)
322 C
323 C CARRY OUT BACK SUBSTITUTION
324 C
325 30 CONTINUE
326 DO 40 J=NM1,1,-1
327 IF (WK(2,J) .EQ. ZERO) GO TO 5008
328 D(1,J) = (D(1,J) - WK(1,J)*D(1,J+1))/WK(2,J)
329 40 CONTINUE
330 C --------------------( END CODING FROM CUBSPL )--------------------
331 C
332 C NORMAL RETURN.
333 C
334 RETURN
335 C
336 C ERROR RETURNS.
337 C
338 5001 CONTINUE
339 C N.LT.2 RETURN.
340 IERR = -1
341 CALL XERMSG ('SLATEC', 'PCHSP',
342 + 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
343 RETURN
344 C
345 5002 CONTINUE
346 C INCFD.LT.1 RETURN.
347 IERR = -2
348 CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR,
349 + 1)
350 RETURN
351 C
352 5003 CONTINUE
353 C X-ARRAY NOT STRICTLY INCREASING.
354 IERR = -3
355 CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING'
356 + , IERR, 1)
357 RETURN
358 C
359 5004 CONTINUE
360 C IC OUT OF RANGE RETURN.
361 IERR = IERR - 3
362 CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1)
363 RETURN
364 C
365 5007 CONTINUE
366 C NWK TOO SMALL RETURN.
367 IERR = -7
368 CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1)
369 RETURN
370 C
371 5008 CONTINUE
372 C SINGULAR SYSTEM.
373 C *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES ***
374 C *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
375 IERR = -8
376 CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR,
377 + 1)
378 RETURN
379 C
380 5009 CONTINUE
381 C ERROR RETURN FROM PCHDF.
382 C *** THIS CASE SHOULD NEVER OCCUR ***
383 IERR = -9
384 CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR,
385 + 1)
386 RETURN
387 C------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
388 END
389