2
SUBROUTINE DPCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
3
C***BEGIN PROLOGUE DPCHFE
4
C***PURPOSE Evaluate a piecewise cubic Hermite function at an array of
5
C points. May be used by itself for Hermite interpolation,
6
C or as an evaluator for DPCHIM or DPCHIC.
7
C***LIBRARY SLATEC (PCHIP)
9
C***TYPE DOUBLE PRECISION (PCHFE-S, DPCHFE-D)
10
C***KEYWORDS CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
11
C PIECEWISE CUBIC EVALUATION
12
C***AUTHOR Fritsch, F. N., (LLNL)
13
C Lawrence Livermore National Laboratory
14
C P.O. Box 808 (L-316)
16
C FTS 532-4275, (510) 422-4275
19
C DPCHFE: Piecewise Cubic Hermite Function Evaluator
21
C Evaluates the cubic Hermite function defined by N, X, F, D at
22
C the points XE(J), J=1(1)NE.
24
C To provide compatibility with DPCHIM and DPCHIC, includes an
25
C increment between successive values of the F- and D-arrays.
27
C ----------------------------------------------------------------------
31
C PARAMETER (INCFD = ...)
33
C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
36
C CALL DPCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
40
C N -- (input) number of data points. (Error return if N.LT.2 .)
42
C X -- (input) real*8 array of independent variable values. The
43
C elements of X must be strictly increasing:
44
C X(I-1) .LT. X(I), I = 2(1)N.
45
C (Error return if not.)
47
C F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is
48
C the value corresponding to X(I).
50
C D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD)
51
C is the value corresponding to X(I).
53
C INCFD -- (input) increment between successive values in F and D.
54
C (Error return if INCFD.LT.1 .)
56
C SKIP -- (input/output) logical variable which should be set to
57
C .TRUE. if the user wishes to skip checks for validity of
58
C preceding parameters, or to .FALSE. otherwise.
59
C This will save time in case these checks have already
60
C been performed (say, in DPCHIM or DPCHIC).
61
C SKIP will be set to .TRUE. on normal return.
63
C NE -- (input) number of evaluation points. (Error return if
66
C XE -- (input) real*8 array of points at which the function is to
70
C 1. The evaluation will be most efficient if the elements
71
C of XE are increasing relative to X;
72
C that is, XE(J) .GE. X(I)
73
C implies XE(K) .GE. X(I), all K.GE.J .
74
C 2. If any of the XE are outside the interval [X(1),X(N)],
75
C values are extrapolated from the nearest extreme cubic,
76
C and a warning error is returned.
78
C FE -- (output) real*8 array of values of the cubic Hermite
79
C function defined by N, X, F, D at the points XE.
81
C IERR -- (output) error flag.
83
C IERR = 0 (no errors).
85
C IERR.GT.0 means that extrapolation was performed at
87
C "Recoverable" errors:
88
C IERR = -1 if N.LT.2 .
89
C IERR = -2 if INCFD.LT.1 .
90
C IERR = -3 if the X-array is not strictly increasing.
91
C IERR = -4 if NE.LT.1 .
92
C (The FE-array has not been changed in any of these cases.)
93
C NOTE: The above errors are checked in the order listed,
94
C and following arguments have **NOT** been validated.
97
C***ROUTINES CALLED DCHFEV, XERMSG
98
C***REVISION HISTORY (YYMMDD)
100
C 820803 Minor cosmetic changes for release 1.
101
C 870707 Corrected XERROR calls for d.p. name(s).
102
C 890206 Corrected XERROR calls.
103
C 890531 Changed all specific intrinsics to generic. (WRB)
104
C 890831 Modified array declarations. (WRB)
105
C 891006 Cosmetic changes to prologue. (WRB)
106
C 891006 REVISION DATE from Version 3.2
107
C 891214 Prologue converted to Version 4.0 format. (BAB)
108
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
109
C***END PROLOGUE DPCHFE
112
C 1. To produce a single precision version, simply:
113
C a. Change DPCHFE to PCHFE, and DCHFEV to CHFEV, wherever they
115
C b. Change the double precision declaration to real,
117
C 2. Most of the coding between the call to DCHFEV and the end of
118
C the IR-loop could be eliminated if it were permissible to
119
C assume that XE is ordered relative to X.
121
C 3. DCHFEV does not assume that X1 is less than X2. thus, it would
122
C be possible to write a version of DPCHFE that assumes a
123
C decreasing X-array by simply running the IR-loop backwards
124
C (and reversing the order of appropriate tests).
126
C 4. The present code has a minor bug, which I have decided is not
127
C worth the effort that would be required to fix it.
128
C If XE contains points in [X(N-1),X(N)], followed by points .LT.
129
C X(N-1), followed by points .GT.X(N), the extrapolation points
130
C will be counted (at least) twice in the total returned in IERR.
134
INTEGER N, INCFD, NE, IERR
135
DOUBLE PRECISION X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*)
138
C DECLARE LOCAL VARIABLES.
140
INTEGER I, IERC, IR, J, JFIRST, NEXT(2), NJ
142
C VALIDITY-CHECK ARGUMENTS.
144
C***FIRST EXECUTABLE STATEMENT DPCHFE
147
IF ( N.LT.2 ) GO TO 5001
148
IF ( INCFD.LT.1 ) GO TO 5002
150
IF ( X(I).LE.X(I-1) ) GO TO 5003
153
C FUNCTION DEFINITION IS OK, GO ON.
156
IF ( NE.LT.1 ) GO TO 5004
160
C LOOP OVER INTERVALS. ( INTERVAL INDEX IS IL = IR-1 . )
161
C ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
166
C SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
168
IF (JFIRST .GT. NE) GO TO 5000
170
C LOCATE ALL POINTS IN INTERVAL.
173
IF (XE(J) .GE. X(IR)) GO TO 30
178
C HAVE LOCATED FIRST POINT BEYOND INTERVAL.
181
IF (IR .EQ. N) J = NE + 1
186
C SKIP EVALUATION IF NO POINTS IN INTERVAL.
188
IF (NJ .EQ. 0) GO TO 50
190
C EVALUATE CUBIC AT XE(I), I = JFIRST (1) J-1 .
192
C ----------------------------------------------------------------
193
CALL DCHFEV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR)
194
* ,NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
195
C ----------------------------------------------------------------
196
IF (IERC .LT. 0) GO TO 5005
198
IF (NEXT(2) .EQ. 0) GO TO 42
199
C IF (NEXT(2) .GT. 0) THEN
200
C IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
203
IF (IR .LT. N) GO TO 41
204
C IF (IR .EQ. N) THEN
205
C THESE ARE ACTUALLY EXTRAPOLATION POINTS.
206
IERR = IERR + NEXT(2)
210
C WE SHOULD NEVER HAVE GOTTEN HERE.
216
IF (NEXT(1) .EQ. 0) GO TO 49
217
C IF (NEXT(1) .GT. 0) THEN
218
C IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
221
IF (IR .GT. 2) GO TO 43
222
C IF (IR .EQ. 2) THEN
223
C THESE ARE ACTUALLY EXTRAPOLATION POINTS.
224
IERR = IERR + NEXT(1)
228
C XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
229
C EVALUATION INTERVAL.
231
C FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
232
DO 44 I = JFIRST, J-1
233
IF (XE(I) .LT. X(IR-1)) GO TO 45
235
C NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
240
C RESET J. (THIS WILL BE THE NEW JFIRST.)
243
C NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
245
IF (XE(J) .LT. X(I)) GO TO 47
247
C NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
250
C AT THIS POINT, EITHER XE(J) .LT. X(1)
251
C OR X(I-1) .LE. XE(J) .LT. X(I) .
252
C RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
265
IF (IR .LE. N) GO TO 10
277
CALL XERMSG ('SLATEC', 'DPCHFE',
278
+ 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
284
CALL XERMSG ('SLATEC', 'DPCHFE', 'INCREMENT LESS THAN ONE', IERR,
289
C X-ARRAY NOT STRICTLY INCREASING.
291
CALL XERMSG ('SLATEC', 'DPCHFE',
292
+ 'X-ARRAY NOT STRICTLY INCREASING', IERR, 1)
298
CALL XERMSG ('SLATEC', 'DPCHFE',
299
+ 'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
303
C ERROR RETURN FROM DCHFEV.
304
C *** THIS CASE SHOULD NEVER OCCUR ***
306
CALL XERMSG ('SLATEC', 'DPCHFE',
307
+ 'ERROR RETURN FROM DCHFEV -- FATAL', IERR, 2)
309
C------------- LAST LINE OF DPCHFE FOLLOWS -----------------------------
added
removed
Lib/Slatec/slatec/dpchfe.f
