*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