1*DECK DPOLFT 2 SUBROUTINE DPOLFT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A) 3C***BEGIN PROLOGUE DPOLFT 4C***PURPOSE Fit discrete data in a least squares sense by polynomials 5C in one variable. 6C***LIBRARY SLATEC 7C***CATEGORY K1A1A2 8C***TYPE DOUBLE PRECISION (POLFIT-S, DPOLFT-D) 9C***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT 10C***AUTHOR Shampine, L. F., (SNLA) 11C Davenport, S. M., (SNLA) 12C Huddleston, R. E., (SNLL) 13C***DESCRIPTION 14C 15C Abstract 16C 17C Given a collection of points X(I) and a set of values Y(I) which 18C correspond to some function or measurement at each of the X(I), 19C subroutine DPOLFT computes the weighted least-squares polynomial 20C fits of all degrees up to some degree either specified by the user 21C or determined by the routine. The fits thus obtained are in 22C orthogonal polynomial form. Subroutine DP1VLU may then be 23C called to evaluate the fitted polynomials and any of their 24C derivatives at any point. The subroutine DPCOEF may be used to 25C express the polynomial fits as powers of (X-C) for any specified 26C point C. 27C 28C The parameters for DPOLFT are 29C 30C Input -- All TYPE REAL variables are DOUBLE PRECISION 31C N - the number of data points. The arrays X, Y and W 32C must be dimensioned at least N (N .GE. 1). 33C X - array of values of the independent variable. These 34C values may appear in any order and need not all be 35C distinct. 36C Y - array of corresponding function values. 37C W - array of positive values to be used as weights. If 38C W(1) is negative, DPOLFT will set all the weights 39C to 1.0, which means unweighted least squares error 40C will be minimized. To minimize relative error, the 41C user should set the weights to: W(I) = 1.0/Y(I)**2, 42C I = 1,...,N . 43C MAXDEG - maximum degree to be allowed for polynomial fit. 44C MAXDEG may be any non-negative integer less than N. 45C Note -- MAXDEG cannot be equal to N-1 when a 46C statistical test is to be used for degree selection, 47C i.e., when input value of EPS is negative. 48C EPS - specifies the criterion to be used in determining 49C the degree of fit to be computed. 50C (1) If EPS is input negative, DPOLFT chooses the 51C degree based on a statistical F test of 52C significance. One of three possible 53C significance levels will be used: .01, .05 or 54C .10. If EPS=-1.0 , the routine will 55C automatically select one of these levels based 56C on the number of data points and the maximum 57C degree to be considered. If EPS is input as 58C -.01, -.05, or -.10, a significance level of 59C .01, .05, or .10, respectively, will be used. 60C (2) If EPS is set to 0., DPOLFT computes the 61C polynomials of degrees 0 through MAXDEG . 62C (3) If EPS is input positive, EPS is the RMS 63C error tolerance which must be satisfied by the 64C fitted polynomial. DPOLFT will increase the 65C degree of fit until this criterion is met or 66C until the maximum degree is reached. 67C 68C Output -- All TYPE REAL variables are DOUBLE PRECISION 69C NDEG - degree of the highest degree fit computed. 70C EPS - RMS error of the polynomial of degree NDEG . 71C R - vector of dimension at least NDEG containing values 72C of the fit of degree NDEG at each of the X(I) . 73C Except when the statistical test is used, these 74C values are more accurate than results from subroutine 75C DP1VLU normally are. 76C IERR - error flag with the following possible values. 77C 1 -- indicates normal execution, i.e., either 78C (1) the input value of EPS was negative, and the 79C computed polynomial fit of degree NDEG 80C satisfies the specified F test, or 81C (2) the input value of EPS was 0., and the fits of 82C all degrees up to MAXDEG are complete, or 83C (3) the input value of EPS was positive, and the 84C polynomial of degree NDEG satisfies the RMS 85C error requirement. 86C 2 -- invalid input parameter. At least one of the input 87C parameters has an illegal value and must be corrected 88C before DPOLFT can proceed. Valid input results 89C when the following restrictions are observed 90C N .GE. 1 91C 0 .LE. MAXDEG .LE. N-1 for EPS .GE. 0. 92C 0 .LE. MAXDEG .LE. N-2 for EPS .LT. 0. 93C W(1)=-1.0 or W(I) .GT. 0., I=1,...,N . 94C 3 -- cannot satisfy the RMS error requirement with a 95C polynomial of degree no greater than MAXDEG . Best 96C fit found is of degree MAXDEG . 97C 4 -- cannot satisfy the test for significance using 98C current value of MAXDEG . Statistically, the 99C best fit found is of order NORD . (In this case, 100C NDEG will have one of the values: MAXDEG-2, 101C MAXDEG-1, or MAXDEG). Using a higher value of 102C MAXDEG may result in passing the test. 103C A - work and output array having at least 3N+3MAXDEG+3 104C locations 105C 106C Note - DPOLFT calculates all fits of degrees up to and including 107C NDEG . Any or all of these fits can be evaluated or 108C expressed as powers of (X-C) using DP1VLU and DPCOEF 109C after just one call to DPOLFT . 110C 111C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, 112C Curve fitting by polynomials in one variable, Report 113C SLA-74-0270, Sandia Laboratories, June 1974. 114C***ROUTINES CALLED DP1VLU, XERMSG 115C***REVISION HISTORY (YYMMDD) 116C 740601 DATE WRITTEN 117C 890531 Changed all specific intrinsics to generic. (WRB) 118C 891006 Cosmetic changes to prologue. (WRB) 119C 891006 REVISION DATE from Version 3.2 120C 891214 Prologue converted to Version 4.0 format. (BAB) 121C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 122C 900911 Added variable YP to DOUBLE PRECISION declaration. (WRB) 123C 920501 Reformatted the REFERENCES section. (WRB) 124C 920527 Corrected erroneous statements in DESCRIPTION. (WRB) 125C***END PROLOGUE DPOLFT 126 INTEGER I,IDEGF,IERR,J,JP1,JPAS,K1,K1PJ,K2,K2PJ,K3,K3PI,K4, 127 * K4PI,K5,K5PI,KSIG,M,MAXDEG,MOP1,NDEG,NDER,NFAIL 128 DOUBLE PRECISION TEMD1,TEMD2 129 DOUBLE PRECISION A(*),DEGF,DEN,EPS,ETST,F,FCRIT,R(*),SIG,SIGJ, 130 * SIGJM1,SIGPAS,TEMP,X(*),XM,Y(*),YP,W(*),W1,W11 131 DOUBLE PRECISION CO(4,3) 132 SAVE CO 133 DATA CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2), 134 1 CO(3,2), CO(4,2), CO(1,3), CO(2,3), CO(3,3), 135 2 CO(4,3)/-13.086850D0,-2.4648165D0,-3.3846535D0,-1.2973162D0, 136 3 -3.3381146D0,-1.7812271D0,-3.2578406D0,-1.6589279D0, 137 4 -1.6282703D0,-1.3152745D0,-3.2640179D0,-1.9829776D0/ 138C***FIRST EXECUTABLE STATEMENT DPOLFT 139 M = ABS(N) 140 IF (M .EQ. 0) GO TO 30 141 IF (MAXDEG .LT. 0) GO TO 30 142 A(1) = MAXDEG 143 MOP1 = MAXDEG + 1 144 IF (M .LT. MOP1) GO TO 30 145 IF (EPS .LT. 0.0D0 .AND. M .EQ. MOP1) GO TO 30 146 XM = M 147 ETST = EPS*EPS*XM 148 IF (W(1) .LT. 0.0D0) GO TO 2 149 DO 1 I = 1,M 150 IF (W(I) .LE. 0.0D0) GO TO 30 151 1 CONTINUE 152 GO TO 4 153 2 DO 3 I = 1,M 154 3 W(I) = 1.0D0 155 4 IF (EPS .GE. 0.0D0) GO TO 8 156C 157C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR 158C CHOOSING DEGREE OF POLYNOMIAL FIT 159C 160 IF (EPS .GT. (-.55D0)) GO TO 5 161 IDEGF = M - MAXDEG - 1 162 KSIG = 1 163 IF (IDEGF .LT. 10) KSIG = 2 164 IF (IDEGF .LT. 5) KSIG = 3 165 GO TO 8 166 5 KSIG = 1 167 IF (EPS .LT. (-.03D0)) KSIG = 2 168 IF (EPS .LT. (-.07D0)) KSIG = 3 169C 170C INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING 171C 172 8 K1 = MAXDEG + 1 173 K2 = K1 + MAXDEG 174 K3 = K2 + MAXDEG + 2 175 K4 = K3 + M 176 K5 = K4 + M 177 DO 9 I = 2,K4 178 9 A(I) = 0.0D0 179 W11 = 0.0D0 180 IF (N .LT. 0) GO TO 11 181C 182C UNCONSTRAINED CASE 183C 184 DO 10 I = 1,M 185 K4PI = K4 + I 186 A(K4PI) = 1.0D0 187 10 W11 = W11 + W(I) 188 GO TO 13 189C 190C CONSTRAINED CASE 191C 192 11 DO 12 I = 1,M 193 K4PI = K4 + I 194 12 W11 = W11 + W(I)*A(K4PI)**2 195C 196C COMPUTE FIT OF DEGREE ZERO 197C 198 13 TEMD1 = 0.0D0 199 DO 14 I = 1,M 200 K4PI = K4 + I 201 TEMD1 = TEMD1 + W(I)*Y(I)*A(K4PI) 202 14 CONTINUE 203 TEMD1 = TEMD1/W11 204 A(K2+1) = TEMD1 205 SIGJ = 0.0D0 206 DO 15 I = 1,M 207 K4PI = K4 + I 208 K5PI = K5 + I 209 TEMD2 = TEMD1*A(K4PI) 210 R(I) = TEMD2 211 A(K5PI) = TEMD2 - R(I) 212 15 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2 213 J = 0 214C 215C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION 216C 217 IF (EPS) 24,26,27 218C 219C INCREMENT DEGREE 220C 221 16 J = J + 1 222 JP1 = J + 1 223 K1PJ = K1 + J 224 K2PJ = K2 + J 225 SIGJM1 = SIGJ 226C 227C COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1 228C 229 IF (J .GT. 1) A(K1PJ) = W11/W1 230C 231C COMPUTE NEW A COEFFICIENT 232C 233 TEMD1 = 0.0D0 234 DO 18 I = 1,M 235 K4PI = K4 + I 236 TEMD2 = A(K4PI) 237 TEMD1 = TEMD1 + X(I)*W(I)*TEMD2*TEMD2 238 18 CONTINUE 239 A(JP1) = TEMD1/W11 240C 241C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS 242C 243 W1 = W11 244 W11 = 0.0D0 245 DO 19 I = 1,M 246 K3PI = K3 + I 247 K4PI = K4 + I 248 TEMP = A(K3PI) 249 A(K3PI) = A(K4PI) 250 A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP 251 19 W11 = W11 + W(I)*A(K4PI)**2 252C 253C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE 254C PRECISION 255C 256 TEMD1 = 0.0D0 257 DO 20 I = 1,M 258 K4PI = K4 + I 259 K5PI = K5 + I 260 TEMD2 = W(I)*((Y(I)-R(I))-A(K5PI))*A(K4PI) 261 20 TEMD1 = TEMD1 + TEMD2 262 TEMD1 = TEMD1/W11 263 A(K2PJ+1) = TEMD1 264C 265C UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND 266C ACCUMULATE SUM OF SQUARES OF ERRORS. THE POLYNOMIAL EVALUATIONS ARE 267C COMPUTED AND STORED IN EXTENDED PRECISION. FOR THE I-TH DATA POINT, 268C THE MOST SIGNIFICANT BITS ARE STORED IN R(I) , AND THE LEAST 269C SIGNIFICANT BITS ARE IN A(K5PI) . 270C 271 SIGJ = 0.0D0 272 DO 21 I = 1,M 273 K4PI = K4 + I 274 K5PI = K5 + I 275 TEMD2 = R(I) + A(K5PI) + TEMD1*A(K4PI) 276 R(I) = TEMD2 277 A(K5PI) = TEMD2 - R(I) 278 21 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2 279C 280C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE 281C MAXDEG HAS BEEN REACHED 282C 283 IF (EPS) 23,26,27 284C 285C COMPUTE F STATISTICS (INPUT EPS .LT. 0.) 286C 287 23 IF (SIGJ .EQ. 0.0D0) GO TO 29 288 DEGF = M - J - 1 289 DEN = (CO(4,KSIG)*DEGF + 1.0D0)*DEGF 290 FCRIT = (((CO(3,KSIG)*DEGF) + CO(2,KSIG))*DEGF + CO(1,KSIG))/DEN 291 FCRIT = FCRIT*FCRIT 292 F = (SIGJM1 - SIGJ)*DEGF/SIGJ 293 IF (F .LT. FCRIT) GO TO 25 294C 295C POLYNOMIAL OF DEGREE J SATISFIES F TEST 296C 297 24 SIGPAS = SIGJ 298 JPAS = J 299 NFAIL = 0 300 IF (MAXDEG .EQ. J) GO TO 32 301 GO TO 16 302C 303C POLYNOMIAL OF DEGREE J FAILS F TEST. IF THERE HAVE BEEN THREE 304C SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND. 305C 306 25 NFAIL = NFAIL + 1 307 IF (NFAIL .GE. 3) GO TO 29 308 IF (MAXDEG .EQ. J) GO TO 32 309 GO TO 16 310C 311C RAISE THE DEGREE IF DEGREE MAXDEG HAS NOT YET BEEN REACHED (INPUT 312C EPS = 0.) 313C 314 26 IF (MAXDEG .EQ. J) GO TO 28 315 GO TO 16 316C 317C SEE IF RMS ERROR CRITERION IS SATISFIED (INPUT EPS .GT. 0.) 318C 319 27 IF (SIGJ .LE. ETST) GO TO 28 320 IF (MAXDEG .EQ. J) GO TO 31 321 GO TO 16 322C 323C RETURNS 324C 325 28 IERR = 1 326 NDEG = J 327 SIG = SIGJ 328 GO TO 33 329 29 IERR = 1 330 NDEG = JPAS 331 SIG = SIGPAS 332 GO TO 33 333 30 IERR = 2 334 CALL XERMSG ('SLATEC', 'DPOLFT', 'INVALID INPUT PARAMETER.', 2, 335 + 1) 336 GO TO 37 337 31 IERR = 3 338 NDEG = MAXDEG 339 SIG = SIGJ 340 GO TO 33 341 32 IERR = 4 342 NDEG = JPAS 343 SIG = SIGPAS 344C 345 33 A(K3) = NDEG 346C 347C WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT 348C ALL THE DATA POINTS IF R DOES NOT ALREADY CONTAIN THESE VALUES 349C 350 IF(EPS .GE. 0.0 .OR. NDEG .EQ. MAXDEG) GO TO 36 351 NDER = 0 352 DO 35 I = 1,M 353 CALL DP1VLU (NDEG,NDER,X(I),R(I),YP,A) 354 35 CONTINUE 355 36 EPS = SQRT(SIG/XM) 356 37 RETURN 357 END 358