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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
|
*DECK DPOLFT
SUBROUTINE DPOLFT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A)
C***BEGIN PROLOGUE DPOLFT
C***PURPOSE Fit discrete data in a least squares sense by polynomials
C in one variable.
C***LIBRARY SLATEC
C***CATEGORY K1A1A2
C***TYPE DOUBLE PRECISION (POLFIT-S, DPOLFT-D)
C***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT
C***AUTHOR Shampine, L. F., (SNLA)
C Davenport, S. M., (SNLA)
C Huddleston, R. E., (SNLL)
C***DESCRIPTION
C
C Abstract
C
C Given a collection of points X(I) and a set of values Y(I) which
C correspond to some function or measurement at each of the X(I),
C subroutine DPOLFT computes the weighted least-squares polynomial
C fits of all degrees up to some degree either specified by the user
C or determined by the routine. The fits thus obtained are in
C orthogonal polynomial form. Subroutine DP1VLU may then be
C called to evaluate the fitted polynomials and any of their
C derivatives at any point. The subroutine DPCOEF may be used to
C express the polynomial fits as powers of (X-C) for any specified
C point C.
C
C The parameters for DPOLFT are
C
C Input -- All TYPE REAL variables are DOUBLE PRECISION
C N - the number of data points. The arrays X, Y and W
C must be dimensioned at least N (N .GE. 1).
C X - array of values of the independent variable. These
C values may appear in any order and need not all be
C distinct.
C Y - array of corresponding function values.
C W - array of positive values to be used as weights. If
C W(1) is negative, DPOLFT will set all the weights
C to 1.0, which means unweighted least squares error
C will be minimized. To minimize relative error, the
C user should set the weights to: W(I) = 1.0/Y(I)**2,
C I = 1,...,N .
C MAXDEG - maximum degree to be allowed for polynomial fit.
C MAXDEG may be any non-negative integer less than N.
C Note -- MAXDEG cannot be equal to N-1 when a
C statistical test is to be used for degree selection,
C i.e., when input value of EPS is negative.
C EPS - specifies the criterion to be used in determining
C the degree of fit to be computed.
C (1) If EPS is input negative, DPOLFT chooses the
C degree based on a statistical F test of
C significance. One of three possible
C significance levels will be used: .01, .05 or
C .10. If EPS=-1.0 , the routine will
C automatically select one of these levels based
C on the number of data points and the maximum
C degree to be considered. If EPS is input as
C -.01, -.05, or -.10, a significance level of
C .01, .05, or .10, respectively, will be used.
C (2) If EPS is set to 0., DPOLFT computes the
C polynomials of degrees 0 through MAXDEG .
C (3) If EPS is input positive, EPS is the RMS
C error tolerance which must be satisfied by the
C fitted polynomial. DPOLFT will increase the
C degree of fit until this criterion is met or
C until the maximum degree is reached.
C
C Output -- All TYPE REAL variables are DOUBLE PRECISION
C NDEG - degree of the highest degree fit computed.
C EPS - RMS error of the polynomial of degree NDEG .
C R - vector of dimension at least NDEG containing values
C of the fit of degree NDEG at each of the X(I) .
C Except when the statistical test is used, these
C values are more accurate than results from subroutine
C DP1VLU normally are.
C IERR - error flag with the following possible values.
C 1 -- indicates normal execution, i.e., either
C (1) the input value of EPS was negative, and the
C computed polynomial fit of degree NDEG
C satisfies the specified F test, or
C (2) the input value of EPS was 0., and the fits of
C all degrees up to MAXDEG are complete, or
C (3) the input value of EPS was positive, and the
C polynomial of degree NDEG satisfies the RMS
C error requirement.
C 2 -- invalid input parameter. At least one of the input
C parameters has an illegal value and must be corrected
C before DPOLFT can proceed. Valid input results
C when the following restrictions are observed
C N .GE. 1
C 0 .LE. MAXDEG .LE. N-1 for EPS .GE. 0.
C 0 .LE. MAXDEG .LE. N-2 for EPS .LT. 0.
C W(1)=-1.0 or W(I) .GT. 0., I=1,...,N .
C 3 -- cannot satisfy the RMS error requirement with a
C polynomial of degree no greater than MAXDEG . Best
C fit found is of degree MAXDEG .
C 4 -- cannot satisfy the test for significance using
C current value of MAXDEG . Statistically, the
C best fit found is of order NORD . (In this case,
C NDEG will have one of the values: MAXDEG-2,
C MAXDEG-1, or MAXDEG). Using a higher value of
C MAXDEG may result in passing the test.
C A - work and output array having at least 3N+3MAXDEG+3
C locations
C
C Note - DPOLFT calculates all fits of degrees up to and including
C NDEG . Any or all of these fits can be evaluated or
C expressed as powers of (X-C) using DP1VLU and DPCOEF
C after just one call to DPOLFT .
C
C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston,
C Curve fitting by polynomials in one variable, Report
C SLA-74-0270, Sandia Laboratories, June 1974.
C***ROUTINES CALLED DP1VLU, XERMSG
C***REVISION HISTORY (YYMMDD)
C 740601 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 891006 Cosmetic changes to prologue. (WRB)
C 891006 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 900911 Added variable YP to DOUBLE PRECISION declaration. (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C 920527 Corrected erroneous statements in DESCRIPTION. (WRB)
C***END PROLOGUE DPOLFT
INTEGER I,IDEGF,IERR,J,JP1,JPAS,K1,K1PJ,K2,K2PJ,K3,K3PI,K4,
* K4PI,K5,K5PI,KSIG,M,MAXDEG,MOP1,NDEG,NDER,NFAIL
DOUBLE PRECISION TEMD1,TEMD2
DOUBLE PRECISION A(*),DEGF,DEN,EPS,ETST,F,FCRIT,R(*),SIG,SIGJ,
* SIGJM1,SIGPAS,TEMP,X(*),XM,Y(*),YP,W(*),W1,W11
DOUBLE PRECISION CO(4,3)
SAVE CO
DATA CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2),
1 CO(3,2), CO(4,2), CO(1,3), CO(2,3), CO(3,3),
2 CO(4,3)/-13.086850D0,-2.4648165D0,-3.3846535D0,-1.2973162D0,
3 -3.3381146D0,-1.7812271D0,-3.2578406D0,-1.6589279D0,
4 -1.6282703D0,-1.3152745D0,-3.2640179D0,-1.9829776D0/
C***FIRST EXECUTABLE STATEMENT DPOLFT
M = ABS(N)
IF (M .EQ. 0) GO TO 30
IF (MAXDEG .LT. 0) GO TO 30
A(1) = MAXDEG
MOP1 = MAXDEG + 1
IF (M .LT. MOP1) GO TO 30
IF (EPS .LT. 0.0D0 .AND. M .EQ. MOP1) GO TO 30
XM = M
ETST = EPS*EPS*XM
IF (W(1) .LT. 0.0D0) GO TO 2
DO 1 I = 1,M
IF (W(I) .LE. 0.0D0) GO TO 30
1 CONTINUE
GO TO 4
2 DO 3 I = 1,M
3 W(I) = 1.0D0
4 IF (EPS .GE. 0.0D0) GO TO 8
C
C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR
C CHOOSING DEGREE OF POLYNOMIAL FIT
C
IF (EPS .GT. (-.55D0)) GO TO 5
IDEGF = M - MAXDEG - 1
KSIG = 1
IF (IDEGF .LT. 10) KSIG = 2
IF (IDEGF .LT. 5) KSIG = 3
GO TO 8
5 KSIG = 1
IF (EPS .LT. (-.03D0)) KSIG = 2
IF (EPS .LT. (-.07D0)) KSIG = 3
C
C INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING
C
8 K1 = MAXDEG + 1
K2 = K1 + MAXDEG
K3 = K2 + MAXDEG + 2
K4 = K3 + M
K5 = K4 + M
DO 9 I = 2,K4
9 A(I) = 0.0D0
W11 = 0.0D0
IF (N .LT. 0) GO TO 11
C
C UNCONSTRAINED CASE
C
DO 10 I = 1,M
K4PI = K4 + I
A(K4PI) = 1.0D0
10 W11 = W11 + W(I)
GO TO 13
C
C CONSTRAINED CASE
C
11 DO 12 I = 1,M
K4PI = K4 + I
12 W11 = W11 + W(I)*A(K4PI)**2
C
C COMPUTE FIT OF DEGREE ZERO
C
13 TEMD1 = 0.0D0
DO 14 I = 1,M
K4PI = K4 + I
TEMD1 = TEMD1 + W(I)*Y(I)*A(K4PI)
14 CONTINUE
TEMD1 = TEMD1/W11
A(K2+1) = TEMD1
SIGJ = 0.0D0
DO 15 I = 1,M
K4PI = K4 + I
K5PI = K5 + I
TEMD2 = TEMD1*A(K4PI)
R(I) = TEMD2
A(K5PI) = TEMD2 - R(I)
15 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
J = 0
C
C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION
C
IF (EPS) 24,26,27
C
C INCREMENT DEGREE
C
16 J = J + 1
JP1 = J + 1
K1PJ = K1 + J
K2PJ = K2 + J
SIGJM1 = SIGJ
C
C COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1
C
IF (J .GT. 1) A(K1PJ) = W11/W1
C
C COMPUTE NEW A COEFFICIENT
C
TEMD1 = 0.0D0
DO 18 I = 1,M
K4PI = K4 + I
TEMD2 = A(K4PI)
TEMD1 = TEMD1 + X(I)*W(I)*TEMD2*TEMD2
18 CONTINUE
A(JP1) = TEMD1/W11
C
C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS
C
W1 = W11
W11 = 0.0D0
DO 19 I = 1,M
K3PI = K3 + I
K4PI = K4 + I
TEMP = A(K3PI)
A(K3PI) = A(K4PI)
A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP
19 W11 = W11 + W(I)*A(K4PI)**2
C
C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE
C PRECISION
C
TEMD1 = 0.0D0
DO 20 I = 1,M
K4PI = K4 + I
K5PI = K5 + I
TEMD2 = W(I)*((Y(I)-R(I))-A(K5PI))*A(K4PI)
20 TEMD1 = TEMD1 + TEMD2
TEMD1 = TEMD1/W11
A(K2PJ+1) = TEMD1
C
C UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND
C ACCUMULATE SUM OF SQUARES OF ERRORS. THE POLYNOMIAL EVALUATIONS ARE
C COMPUTED AND STORED IN EXTENDED PRECISION. FOR THE I-TH DATA POINT,
C THE MOST SIGNIFICANT BITS ARE STORED IN R(I) , AND THE LEAST
C SIGNIFICANT BITS ARE IN A(K5PI) .
C
SIGJ = 0.0D0
DO 21 I = 1,M
K4PI = K4 + I
K5PI = K5 + I
TEMD2 = R(I) + A(K5PI) + TEMD1*A(K4PI)
R(I) = TEMD2
A(K5PI) = TEMD2 - R(I)
21 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
C
C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE
C MAXDEG HAS BEEN REACHED
C
IF (EPS) 23,26,27
C
C COMPUTE F STATISTICS (INPUT EPS .LT. 0.)
C
23 IF (SIGJ .EQ. 0.0D0) GO TO 29
DEGF = M - J - 1
DEN = (CO(4,KSIG)*DEGF + 1.0D0)*DEGF
FCRIT = (((CO(3,KSIG)*DEGF) + CO(2,KSIG))*DEGF + CO(1,KSIG))/DEN
FCRIT = FCRIT*FCRIT
F = (SIGJM1 - SIGJ)*DEGF/SIGJ
IF (F .LT. FCRIT) GO TO 25
C
C POLYNOMIAL OF DEGREE J SATISFIES F TEST
C
24 SIGPAS = SIGJ
JPAS = J
NFAIL = 0
IF (MAXDEG .EQ. J) GO TO 32
GO TO 16
C
C POLYNOMIAL OF DEGREE J FAILS F TEST. IF THERE HAVE BEEN THREE
C SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND.
C
25 NFAIL = NFAIL + 1
IF (NFAIL .GE. 3) GO TO 29
IF (MAXDEG .EQ. J) GO TO 32
GO TO 16
C
C RAISE THE DEGREE IF DEGREE MAXDEG HAS NOT YET BEEN REACHED (INPUT
C EPS = 0.)
C
26 IF (MAXDEG .EQ. J) GO TO 28
GO TO 16
C
C SEE IF RMS ERROR CRITERION IS SATISFIED (INPUT EPS .GT. 0.)
C
27 IF (SIGJ .LE. ETST) GO TO 28
IF (MAXDEG .EQ. J) GO TO 31
GO TO 16
C
C RETURNS
C
28 IERR = 1
NDEG = J
SIG = SIGJ
GO TO 33
29 IERR = 1
NDEG = JPAS
SIG = SIGPAS
GO TO 33
30 IERR = 2
CALL XERMSG ('SLATEC', 'DPOLFT', 'INVALID INPUT PARAMETER.', 2,
+ 1)
GO TO 37
31 IERR = 3
NDEG = MAXDEG
SIG = SIGJ
GO TO 33
32 IERR = 4
NDEG = JPAS
SIG = SIGPAS
C
33 A(K3) = NDEG
C
C WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT
C ALL THE DATA POINTS IF R DOES NOT ALREADY CONTAIN THESE VALUES
C
IF(EPS .GE. 0.0 .OR. NDEG .EQ. MAXDEG) GO TO 36
NDER = 0
DO 35 I = 1,M
CALL DP1VLU (NDEG,NDER,X(I),R(I),YP,A)
35 CONTINUE
36 EPS = SQRT(SIG/XM)
37 RETURN
END
|
