1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
|
*DECK PCHID
REAL FUNCTION PCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERR)
C***BEGIN PROLOGUE PCHID
C***PURPOSE Evaluate the definite integral of a piecewise cubic
C Hermite function over an interval whose endpoints are data
C points.
C***LIBRARY SLATEC (PCHIP)
C***CATEGORY E3, H2A1B2
C***TYPE SINGLE PRECISION (PCHID-S, DPCHID-D)
C***KEYWORDS CUBIC HERMITE INTERPOLATION, NUMERICAL INTEGRATION, PCHIP,
C QUADRATURE
C***AUTHOR Fritsch, F. N., (LLNL)
C Lawrence Livermore National Laboratory
C P.O. Box 808 (L-316)
C Livermore, CA 94550
C FTS 532-4275, (510) 422-4275
C***DESCRIPTION
C
C PCHID: Piecewise Cubic Hermite Integrator, Data limits
C
C Evaluates the definite integral of the cubic Hermite function
C defined by N, X, F, D over the interval [X(IA), X(IB)].
C
C To provide compatibility with PCHIM and PCHIC, includes an
C increment between successive values of the F- and D-arrays.
C
C ----------------------------------------------------------------------
C
C Calling sequence:
C
C PARAMETER (INCFD = ...)
C INTEGER N, IA, IB, IERR
C REAL X(N), F(INCFD,N), D(INCFD,N)
C LOGICAL SKIP
C
C VALUE = PCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERR)
C
C Parameters:
C
C VALUE -- (output) value of the requested integral.
C
C N -- (input) number of data points. (Error return if N.LT.2 .)
C
C X -- (input) real array of independent variable values. The
C elements of X must be strictly increasing:
C X(I-1) .LT. X(I), I = 2(1)N.
C (Error return if not.)
C
C F -- (input) real array of function values. F(1+(I-1)*INCFD) is
C the value corresponding to X(I).
C
C D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is
C the value corresponding to X(I).
C
C INCFD -- (input) increment between successive values in F and D.
C (Error return if INCFD.LT.1 .)
C
C SKIP -- (input/output) logical variable which should be set to
C .TRUE. if the user wishes to skip checks for validity of
C preceding parameters, or to .FALSE. otherwise.
C This will save time in case these checks have already
C been performed (say, in PCHIM or PCHIC).
C SKIP will be set to .TRUE. on return with IERR = 0 or -4.
C
C IA,IB -- (input) indices in X-array for the limits of integration.
C both must be in the range [1,N]. (Error return if not.)
C No restrictions on their relative values.
C
C IERR -- (output) error flag.
C Normal return:
C IERR = 0 (no errors).
C "Recoverable" errors:
C IERR = -1 if N.LT.2 .
C IERR = -2 if INCFD.LT.1 .
C IERR = -3 if the X-array is not strictly increasing.
C IERR = -4 if IA or IB is out of range.
C (VALUE will be zero in any of these cases.)
C NOTE: The above errors are checked in the order listed,
C and following arguments have **NOT** been validated.
C
C***REFERENCES (NONE)
C***ROUTINES CALLED XERMSG
C***REVISION HISTORY (YYMMDD)
C 820723 DATE WRITTEN
C 820804 Converted to SLATEC library version.
C 870813 Minor cosmetic changes.
C 890411 Added SAVE statements (Vers. 3.2).
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890703 Corrected category record. (WRB)
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 930504 Corrected to set VALUE=0 when IERR.ne.0. (FNF)
C***END PROLOGUE PCHID
C
C Programming notes:
C 1. This routine uses a special formula that is valid only for
C integrals whose limits coincide with data values. This is
C mathematically equivalent to, but much more efficient than,
C calls to CHFIE.
C**End
C
C DECLARE ARGUMENTS.
C
INTEGER N, INCFD, IA, IB, IERR
REAL X(*), F(INCFD,*), D(INCFD,*)
LOGICAL SKIP
C
C DECLARE LOCAL VARIABLES.
C
INTEGER I, IUP, LOW
REAL H, HALF, SIX, SUM, VALUE, ZERO
SAVE ZERO, HALF, SIX
C
C INITIALIZE.
C
DATA ZERO /0./, HALF /0.5/, SIX /6./
C***FIRST EXECUTABLE STATEMENT PCHID
VALUE = ZERO
C
C VALIDITY-CHECK ARGUMENTS.
C
IF (SKIP) GO TO 5
C
IF ( N.LT.2 ) GO TO 5001
IF ( INCFD.LT.1 ) GO TO 5002
DO 1 I = 2, N
IF ( X(I).LE.X(I-1) ) GO TO 5003
1 CONTINUE
C
C FUNCTION DEFINITION IS OK, GO ON.
C
5 CONTINUE
SKIP = .TRUE.
IF ((IA.LT.1) .OR. (IA.GT.N)) GO TO 5004
IF ((IB.LT.1) .OR. (IB.GT.N)) GO TO 5004
IERR = 0
C
C COMPUTE INTEGRAL VALUE.
C
IF (IA .NE. IB) THEN
LOW = MIN(IA, IB)
IUP = MAX(IA, IB) - 1
SUM = ZERO
DO 10 I = LOW, IUP
H = X(I+1) - X(I)
SUM = SUM + H*( (F(1,I) + F(1,I+1)) +
* (D(1,I) - D(1,I+1))*(H/SIX) )
10 CONTINUE
VALUE = HALF * SUM
IF (IA .GT. IB) VALUE = -VALUE
ENDIF
C
C NORMAL RETURN.
C
5000 CONTINUE
PCHID = VALUE
RETURN
C
C ERROR RETURNS.
C
5001 CONTINUE
C N.LT.2 RETURN.
IERR = -1
CALL XERMSG ('SLATEC', 'PCHID',
+ 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
GO TO 5000
C
5002 CONTINUE
C INCFD.LT.1 RETURN.
IERR = -2
CALL XERMSG ('SLATEC', 'PCHID', 'INCREMENT LESS THAN ONE', IERR,
+ 1)
GO TO 5000
C
5003 CONTINUE
C X-ARRAY NOT STRICTLY INCREASING.
IERR = -3
CALL XERMSG ('SLATEC', 'PCHID',
+ 'X-ARRAY NOT STRICTLY INCREASING', IERR, 1)
GO TO 5000
C
5004 CONTINUE
C IA OR IB OUT OF RANGE RETURN.
IERR = -4
CALL XERMSG ('SLATEC', 'PCHID', 'IA OR IB OUT OF RANGE', IERR, 1)
GO TO 5000
C------------- LAST LINE OF PCHID FOLLOWS ------------------------------
END
|
