PROGRAM GENLS C C GENERALIZED LEAST SQUARES PROGRAM: ERRORS IN VARIABLES C GENLS.FOR Version 1.1 C C J. ROSS MACDONALD, DEPT. OF PHYSICS AND ASTRONOMY C UNIVERSITY OF NORTH CAROLINA C CHAPEL HILL, NC 27599-3255 C [(919)-967-5005] C 2/16/92 C C partly based on FORTRAN IV program GNLS written by C Don R. Powell, 1970-1971 C see D.R. Powell and J.R. Macdonald, C Computer Journal, 15(1972)148 C C See also J.R. Macdonald and W.J. Thompson C American Journal of Physics, 1992 C C INPUT DATA FILE: Any name C Its form and contents are described in the C file README.GLS and in Menu 0, which appears C when the program is started C C OUTPUT FILE NAME: GENOUTC C This file contains a record of the commands C and inputs used in an interactive session C C OUTPUT FILE NAME: GENOUTF C This file contains (possible intermediate) C and final fitting outputs C C OUTPUT FILE NAME: GENOUTR C This file contains (possible intermediate) and C final residual rough plots for possible printing C C OUTPUT FILE NAME: GENOUTT C This file contains data and residual-related C tables C C IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL SFUNC CHARACTER*78 ALPHA DIMENSION FX(151),RX(151),RY(151),XKEEP(151) DIMENSION YKEEP(151),XHOLD(151),SDER(10),FY(151), 1 FNY(151),FNX(151),SXI(151),SYI(151),SI(151),DI(151), 3 YN(151),RL(151),AA(10),SAK(10),SAKR(10),RC(10,10), 4 SXKP(151),SYKP(151),ALAM(2),AO(10), 5 SSE(3),WX(151),WY(151) PARAMETER (DM9=1.D-9) C COMMON /IN/ N,NF,NP,MAXD,G,EPS,A(10),NOUT,NS(10),ALPHA COMMON /DATA/ Y(151),X(151),SX(151),SY(151) COMMON /EXTRA/ FAC COMMON /SECDER/ P1,P2,DEL,XX COMMON /EXTRY/ DEL1,DEL2,JJJ,KKK,XPART COMMON /DON/ RX,RY,XKEEP,YPART,XPORT,XHOLD COMMON /WTI/ WXI(151),WYI(151) SAVE C C READ OBSERVATIONS, STANDARD DEVIATIONS, AND SET DEFAULTS C CALL OBSWT(N,ALPHA) C OPEN(9,FILE='GENOUTC') C WRITE(9,929) ALPHA IWE = 0 IDE = 0 IWF = 1 IPD = 0 IPR = 0 IFCT = 0 IFPM = 0 ICONV = 0 IMN5 = 0 ISD = 0 ILV = 0 ITST = 0 DO 120 I = 1,N YKEEP(I) = Y(I) XKEEP(I) = X(I) YN(I) = Y(I) SXKP(I) = SX(I) SYKP(I) = SY(I) 120 CONTINUE C 888 EPS1 = 1.D-8 EPS2 = 1.D-9 EPS3 = 1.D-10 FAC = 1.D-4 G = 1.D0 EPS = 1.D-6 MAXG = 36 MAXD = 10 MAXGO = 36 MAXDO = 10 IF(ICONV.EQ.1) GOTO 322 C NF = 1 NP = 2 NOUT = 0 C ML = 1 AO(1) = 0.1D-0 AO(2) = 1.D0 DO 219 K = 3,10 AO(K) = 1.D-10 A(K) = AO(K) 219 CONTINUE A(1) = AO(1) A(2) = AO(2) C C C ####### MENU NO. 1: DEFAULTS ####### C 317 CONTINUE WRITE(*,321) 321 FORMAT(3X,'@@ MENU #1: DEFAULT CHOICES @@'/ 1 5X,'1: ACCEPT ALL DEFAULTS AND CONTINUE; MIN OUTPUT; ST.-LINE F 1IT'/ 2 5X,'2: CHANGE CONVERGENCE CHOICES'/ 3 5X,'3: CHANGE PRINTING CHOICES'/ 5 5X,'4: CONTINUE TO MENU #2 - SET INPUT VALUES'/) READ(*,*) IMN1 WRITE(9,76) IMN1 76 FORMAT(2X,'MENU 1: DEFAULT CHOICE =',I3) C GOTO(311,312,313,314) IMN1 311 CONTINUE NP = 2 NOUT = 0 IWF = 1 IPR = 0 NF = 1 GOTO 535 C 312 CONTINUE IFCT = 0 C ####### MENU NO. 1A: CONVERGENCE ####### 322 WRITE(*,323) 323 FORMAT(3X,'MENU #1A: CHANGE CONVERGENCE CHOICES'/ 1 5X,'1: RESET MAX # DEMING ITERATIONS, MAXD'/ 2 5X,'2: RESET MAX # GENERAL ITERATIONS, MAXG'/ 4 5X,'3: RESET DEMING PARAMETERS CONV. LIMIT, EPS'/ 6 5X,'4: RESET GEN. PARAMETERS CONV. LIMIT, EPS1'/ 7 5X,'5: TURN ON/OFF LEVENBERG-MARQUARDT CONVERGENCE METHOD'/ 3 5X,'6: RESET CONVERGENCE FACTOR, G'/ 5 5X,'7: RESET RELATIVE INCREMENT FOR DERIVATIVES, FAC'/ 5 5X,'8: RESET DERIVATIVE CONV. LIMIT, EPS2'/ 5 5X,'9: RESET CONV. LIMIT FOR SUM-OF-SQUARES, EPS3'/ 7 5X,'10: RESET ALL DEFAULTS IN 1A TO ORIGINAL VALUES'/ 8 5X,'11: RETURN TO DEFAULT MENU 1'/) READ(*,*) IMN1A WRITE(9,77) IMN1A 77 FORMAT(2X,'MENU 1: DEFAULT CHOICE =',I3) GOTO(331,332,333,334,335,336,337,338,339,340,341) IMN1A 331 CONTINUE WRITE(*,461) 461 FORMAT(3X,'RESET MAX # OF DEMING ITERATIONS (DEF=10)') READ(*,*) MAXD GOTO 322 332 CONTINUE C WRITE(*,462) 462 FORMAT(3X,'RESET MAX # OF GENERAL ITERATIONS (DEF=36)') READ(*,*) MAXG GOTO 322 333 CONTINUE WRITE(*,463) 463 FORMAT(3X,'RESET DEMING PARAMETERS CONVERGENCE LIMIT (DEF=1.D-6)') READ(*,*) EPS GOTO 322 334 CONTINUE WRITE(*,464) 464 FORMAT(3X,'RESET GENERAL PARAMETERS CONVERG. LIMIT (DEF=1.D-8)') READ(*,*) EPS1 GOTO 322 335 CONTINUE WRITE(*,391) 391 FORMAT(3X,'ENTER 1 TO TURN ON; 0 TO TURN OFF') READ(*,*) ILV GOTO 322 336 CONTINUE WRITE(*,465) 465 FORMAT(3X,'RESET CONVERGENCE FACTOR G (DEF=1; 05 FOR THIS CHOICE AND PICK'/ 1 5X,'PROPER NOUT VALUES TO DEFINE FUNCTION:'/ 2 5X,'[1 + A1X +...+ A5(X**5)]/[A6 + A7X +...+ A10(X**5)]'/) ELSEIF(NF.EQ.6) THEN DO 115 I = 1,N IF(YKEEP(I).LE.0.D0) THEN PAUSE 'BAD Y DATA FOR LN TRANSFORMATION' STOP ENDIF Y(I) = DLOG(YKEEP(I)) 115 CONTINUE WRITE(9,184) 184 FORMAT(5X,'LN TRANSFORMATION OF Y DATA') INF = 1 WRITE(*,279) 279 FORMAT(5X,'NOW SELECT NEW FIT FUNCTION, e.g. NF =1',/ 1 5X,'USUALLY UNITY WEIGHTING CHOICE APPROPRIATE HERE'/) GOTO 578 ENDIF GOTO 518 C 532 CONTINUE WRITE(*,622) 622 FORMAT(3X,'MAKE SURE THAT THERE ARE NP INITIAL VALUES SPEC 1IFIED'/5X,'ENTER NEW NP VALUE (< 11)') READ(*,*) NP WRITE(9,84) NP 84 FORMAT(5X,'NEW TOTAL # OF PARAMETERS =',I3) GOTO 518 C 533 CONTINUE WRITE(*,507) 507 FORMAT(3X,'ENTER NEW VALUE OF NOUT') READ(*,*) NOUT WRITE(9,75) NOUT 75 FORMAT(5X,'NEW # OF FIXED PARAMETERS =',I3) IF(NOUT.NE.0) THEN WRITE(*,508) 508 FORMAT(3X,'ENTER INDICES OF NOUT FIXED PARAMETERS') 781 FORMAT(5X,10I4) READ(*,*) (NS(I),I=1,NOUT) WRITE(9,781) (NS(I),I=1,NOUT) ENDIF GOTO 518 C 534 CONTINUE WRITE(*,61) 61 FORMAT(3X,'SET/RESET INITIAL VALUES OF NP PARAMETERS') 739 FORMAT(5X,1P5E11.3) READ(*,*) (A(K),K=1,NP) WRITE(9,739) (A(K),K=1,NP) DO 71 K = 1,NP A(K) = A(K) + 1.D-25 AO(K) = A(K) 71 CONTINUE GOTO 518 C 535 CONTINUE C C IMEND = 0 PRINT 20, ALPHA,N,NP,NOUT,NF 20 FORMAT(/15X,'GENLS PROGRAM V.1.0'// 2X,A78// 1 5X,24HNUMBER OF OBSERVATIONS = I5,3X, 2 22HNUMBER OF PARAMETERS = I4/5X, 3 'NUMBER OF FIXED PARAMETERS =' I4/5X, 4 'FUNCTION CHOICE NUMBER =' I4/) C PRINT 21, EPS,EPS1,G,FAC,EPS2,EPS3 21 FORMAT(5X,'PARAMETER CONVERGENCE TOLERANCE (DEMING) ='1PE10.3/ X 5X, 'PARAMETER CONVERGENCE TOLERANCE (GENERAL) ='E10.3/ X 5X, 'G = CONVERGENCE FACTOR ='E10.3/ X 5X, 'H = RELATIVE INCREMENT FOR NUM. DERIVATIVES ='E10.3/ X 5X, 'DERIVATIVE CONVERGENCE TOLERANCE ='E10.3/ X 5X, 'S CONVERGENCE TOLERANCE ='E10.3/) C 361 IF(NOUT.NE.0) THEN IF(IMEND.GT.1) THEN WRITE(10,30) ELSE PRINT 30 ENDIF 30 FORMAT(/5X,5HNOTE:) DO 40 I = 1,NOUT K = NS(I) IF(IMEND.GT.1) THEN WRITE(10,50) K,A(K) ELSE PRINT 50, K,A(K) ENDIF 40 CONTINUE IF(IMEND.GT.1) WRITE(10,70) PRINT 70 70 FORMAT(/) ENDIF 50 FORMAT(7X, 9HPARAMETER I3,23H IS FIXED AT THE VALUE 1PE16.7) C IF(IMEND.GE.2) GOTO 666 C C INP = 0 C C ICOMP = 0 C SDOLD = 1.D0 C SCOMP = 1.D25 C C RETURN HERE FOR NEW WEIGHTING POSSIBILITIES C 125 CONTINUE DIV = N - NP + NOUT DIVI = 1.D0/DIV IF(IDE.EQ.1) IWE = 0 C C ####### MENU NO. 3: SELECT WEIGHTING ####### C 118 CONTINUE ITD = 0 C IDIRC = 0 ITST = 0 IF(INF.EQ.1) WRITE(*,721) 721 FORMAT(5X,'NAT. LOG. TRANSFORM OF Y DATA'/) C WRITE(*,192) 192 FORMAT(5X,'INITIAL OR CURRENT FITTING VALUES OF THE PARAMETERS') WRITE(*,191) (A(K),K=1,NP) 191 FORMAT(5X,1P5E13.4) PRINT 70 C WRITE(*,121) 121 FORMAT(3X,'MENU #3: SELECT NEW WEIGHTING'/ 1 5X,'0: RESET X-Y DATA TO INPUT VALUES AND RETURN TO THIS MENU'/ 1 5X,'1: RETURN TO DEFAULT MENU #1'// 2 5X,'2: USE LAST CURRENT VALUES FOR FIT'/ 3 5X,'3: ALL SIG(X) CONSTANT; ALL SIG(Y) CONSTANT: SIGX,SIGY'/ 4 5X,'4: RATIO OF VARIABLE X AND Y SIGMAS CONSTANT: IDWMC'/ 5 5X,'5: GENERAL WEIGHTING CHOICES FOR X AND Y'/ 7 5X,'6: EXIT PROGRAM') READ(*,*) IWTC WRITE(9,81) IWTC 81 FORMAT(2X,'MENU 3: WEIGHTING CHOICE =',I3) ILN = 2 IF(IWTC.EQ.6) THEN CLOSE(3) CLOSE(9) CLOSE(10) CLOSE(19) STOP ELSEIF(IWTC.EQ.1) THEN GOTO 317 ELSEIF(IWTC.EQ.0) THEN DO 127 I = 1,N Y(I) = YKEEP(I) X(I) = XKEEP(I) YN(I) = Y(I) SX(I) = SXKP(I) SY(I) = SYKP(I) INF = 0 127 CONTINUE GOTO 118 ELSEIF(IWTC.EQ.2) THEN CONTINUE ELSEIF(IWTC.EQ.5) THEN IPWC = 0 WRITE(*,514) 514 FORMAT(3X,'DETERMINE FORM OF GENERAL ITERATIVE WEIGHTING; LET 1 Q = X OR Y'/ 13X,'THEN FOR EACH POINT: WGTQ = 1/[UQ*UQ + PQ*{|Q|**(2*XIQ)}]'/) C WRITE(*,524) 524 FORMAT(3X,'ENTER ADDITIVE ERROR TERMS UX AND UY') READ(*,*) UX,UY USX = UX*UX USY = UY*UY C WRITE(9,182) UX,UY 182 FORMAT(3X,'UX =',1PE12.3,4X,'UY =',1PE12.3) WRITE(*,124) 124 FORMAT(3X,'ENTER CONSTANT PX; ENTER CONSTANT PY') READ(*,*) PX,PY PX = PX + 1.D-35 PY = PY + 1.D-35 PXS = PX*PX PYS = PY*PY WRITE(*,104) 104 FORMAT(3X,'ENTER |X| POWER EXPONENT, XIX; ENTER |Y| POWER EXPONENT 1, XIY') READ(*,*) XIX,XIY WRITE(*,117) XIX2 = 2.D0*XIX XIY2 = 2.D0*XIY C WRITE(9,82) PX,XIX,PY,XIY 82 FORMAT(2X,'PX =',1PE12.3,2X,'XIX =',1PE12.3,2X,'PY =',1PE12.3, 1 2X,'XIY =',1PE12.3) C DO 128 I = 1,N IF(XIX.NE.0) THEN WX(I) = PXS*(DABS(X(I))**XIX2 + USX + 1.D-20) ELSE WX(I) = PXS + USX ENDIF IF(XIY.NE.0) THEN XX = X(I) YN(I) = SFUNC(ISD) WY(I) = PYS*(DABS(YN(I))**XIY2 + USY + 1.D-20) ELSE WY(I) = PYS + USY ENDIF 128 CONTINUE IF(PX.EQ.PY.AND.XIX.EQ.XIY) THEN IRXY = 1 ELSE IRXY = 0 ENDIF ENDIF C IF(IWTC.EQ.2.OR.IWTC.EQ.5) THEN I23 = 0 GOTO 129 ENDIF C IF(IWTC.EQ.4) THEN KM = 1 ELSE KM = 2 ENDIF C C ####### MENU NO. 4 ####### C WRITE(*,422) 422 FORMAT(3X,'MENU 4: ENTER NUMBER TO SELECT SOLUTION METHO +D'/ 1 5X,'1: STRAIGHT-LINE FIT: PW AND GM ONLY (NF=1,NP<3)'/ 2 5X,'2: DEMING AND GENERALIZED LEAST SQUARES (GLS) METHODS'/ 3 5X,'3: COMBINATION OF CHOICES 1 AND 2 (NF=1,NP<3)'/) READ(*,*) ILN C WRITE(9,83) ILN 83 FORMAT(2X,'MENU 4: SOLUTION CHOICE =',I3) C IF(NF.EQ.1.AND.NP.EQ.2.AND.NOUT.EQ.0.AND.(IWTC.EQ.3.OR.IWTC. 1 EQ.4).AND.ILN.NE.2) THEN I23 = 1 ELSE I23 = 0 ENDIF C IF(IWTC.EQ.3) THEN IF(NF.EQ.1.AND.NP.GT.2.AND.ILN.NE.2) GOTO 125 WRITE(*,122) 122 FORMAT(3X,'ENTER CONSTANT SIGX; ENTER CONSTANT SIGY') READ(*,*) SIGX,SIGY SIGX = SIGX + 1.D-35 SIGY = SIGY + 1.D-35 WRITE(9,93) SIGX,SIGY 93 FORMAT(2X,'MENU 4: SIGX =',1PE14.4,5X,'SIGY ='1PE14.4) ALAM(1) = (SIGY/SIGX)**2 WRITE(*,178) ALAM(1) 178 FORMAT(5X,'LAMBDA =',1PE13.4) WRITE(*,117) 117 FORMAT(/) DO 123 I = 1,N SX(I) = SIGX SY(I) = SIGY X(I) = XKEEP(I) YN(I) = YKEEP(I) 123 CONTINUE ELSEIF(IWTC.EQ.4) THEN C FIRST, TEST WEIGHT CHOICE FOR CONSISTENCY DO 223 I = 1,N SX(I) = SXKP(I) SY(I) = SYKP(I) RL(I) = SY(I)/SX(I) IF(I.GT.1) THEN DIFF = DABS(RL(I) - RL(I-1)) IF(DIFF.GT.1.D-3) THEN PAUSE 'WEIGHT CHOICE IS NOT VALID!!' GOTO 125 ENDIF ENDIF 223 CONTINUE ALAM(1) = RL(2)**2 WRITE(*,85) ALAM(1) WRITE(9,85) ALAM(1) 85 FORMAT(2X,'IDWM LAMBDA =',1PE13.4) ENDIF C 129 CONTINUE C HERE WX IS INVERSE WT AND WXI IS ORDINARY WEIGHT EXPRESSION ******* DO 229 I = 1,N IF(IWTC.EQ.5) THEN SX(I) = DSQRT(WX(I)) SY(I) = DSQRT(WY(I)) ELSE WX(I) = SX(I)*SX(I) WY(I) = SY(I)*SY(I) ENDIF WXI(I) = 1.D0/WX(I) WYI(I) = 1.D0/WY(I) 229 CONTINUE C C PRINT CURRENT VALUES OF X AND Y AND STANDARD DEVIATIONS C IF(IPD.EQ.1) THEN PRINT 130 130 FORMAT(5X,20HINPUT TO THE ROUTINE /) WRITE(*,140) 140 FORMAT(12X,'X',14X,'SIG(X)',10X,'YCALC',11X,'SIG(Y)') C DO 150 I = 1,N 150 PRINT 160, X(I),SX(I),YN(I),SY(I) 160 FORMAT(5X,1PE16.7,E16.7,E16.7,E16.7) ENDIF C C WRITE(*,*) I23,IFPM,ILN IF(I23.EQ.1) THEN CALL PWGM(ALAM,DIVI,WX,WY,WXI,WYI,XKEEP,YKEEP,A,N,KM,IFPM) IF(ILN.EQ.1) GOTO 918 ENDIF C C IRET = 0 CALL MCALC(IPR,IPD,IWF,ITD,ILV,IWE,IPWC,IWTC,WX,WY,XIX,XIX2, 1XIY,XIY2,USX,USY,PXS,PYS,EPS1,EPS2,EPS3,DIVI,IRET,MAXG,SUMSXR, 2SUMSYR,S2,AA,SAK,SAKR,RC,SDER,IDIRC) C IF(IRET.EQ.1) GOTO 518 C C C CALCULATIONS AFTER CONVERGENCE C CALL PRINTOUT(SUMSXR,SUMSYR,S2,A,AA,SAK,SAKR,RC,SDER,1,NP) IMEND = 0 C SB = 0.D0 FMY = 0.D0 FMX = 0.D0 FKY = 0.D0 FKX = 0.D0 FJ = 0.D0 SUMSY = 0.D0 SUMSX = 0.D0 SUMDI = 0.D0 SMURX = 0.D0 SMURY = 0.D0 SMCOR = 0.D0 SMRX = 0.D0 SMRY = 0.D0 C C WX=SIGX**2; WXI=SIGX(-2)= ACTUAL WT DO 980 I = 1,N XX = X(I) CALL DERXE(X,DX,I) DX = -DX DL = WY(I) + WX(I)*DX*DX YN(I) = SFUNC(ISD) DF = Y(I) - YN(I) IF(WYI(I)/WXI(I).LT.1.0D9) THEN RY(I) = DF ELSE RY(I) = 0 ENDIF RX(I) = XKEEP(I) - X(I) SMURX = SMURX + RX(I) SMURY = SMURY + RY(I) FX(I) = XKEEP(I) - RX(I) FY(I) = Y(I) - RY(I) FNY(I) = RY(I)/SY(I) FNX(I) = RX(I)/SX(I) SMRX = SMRX + FNX(I) SMRY = SMRY + FNY(I) SYI(I) = FNY(I)*FNY(I) SXI(I) = FNX(I)*FNX(I) SI(I) = SYI(I) + SXI(I) SCOR = RY(I) -DX*RX(I) DSGS = DSIGN(1.D0,SCOR) DI(I) = DSGS*DSQRT(SI(I)) SMCOR = SMCOR + DI(I) FMY = FMY + FNY(I) FMX = FMX + FNX(I) FKY = FKY + FY(I)*FNY(I) FKX = FKX + FX(I)*FNX(I) FJ = FJ + FNX(I)*FNY(I) SUMSX = SUMSX + SXI(I) SUMSY = SUMSY + SYI(I) SUMDI = SUMDI + DI(I) 980 CONTINUE SB = SUMSX + SUMSY C DIVN = 1.D0/DFLOAT(N) AVURX = SMURX*DIVN AVURY = SMURY*DIVN AVRX = SMRX*DIVN AVRY = SMRY*DIVN AVCOR = SMCOR*DIVN 351 IF(IMN5.EQ.5) THEN WRITE(*,1033) 1033 FORMAT(10X,'MEANS OF UNWEIGHTED X AND Y RESIDUALS') WRITE(*,1035) AVURX,AVURY WRITE(*,1034) WRITE(*,1035) AVRX,AVRY,AVCOR 1034 FORMAT(/10X,'MEANS OF WEIGHTED X,Y, AND COMPOSITE RESIDUALS') 1035 FORMAT(10X,1P3E15.4//) GOTO 305 ENDIF C PIL = 1.D0/DLOG(DFLOAT(N)) IF(SUMSX.LT.1.D-30) THEN MPS = 2 MPF = 2 ELSEIF(SUMSY.LT.1.D-30) THEN MPS = 1 MPF = 1 ELSE MPS = 1 MPF = 3 ENDIF C IF(SUMSX.NE.0) SUMSXI = 1.D0/SUMSX IF(SUMSY.NE.0) SUMSYI = 1.D0/SUMSY SSIN = 1.D0/SB DO 996 MP = MPS,MPF SUMQ = 0.D0 DO 993 I = 1,N IF(MP.EQ.1) THEN QQ = SXI(I)*SUMSXI ELSEIF(MP.EQ.2) THEN QQ = SYI(I)*SUMSYI ELSEIF(MP.EQ.3) THEN QQ = SI(I)*SSIN ENDIF IF(QQ.EQ.0.D0) THEN PPK = 0.D0 ELSE PPK = QQ*DLOG(QQ) ENDIF SUMQ = SUMQ + PPK 993 CONTINUE SSE(MP) = 1.D0 + PIL*SUMQ 996 CONTINUE 305 IF(IMN5.EQ.5) THEN DO 306 MP = MPS,MPF IF(MP.EQ.MPS) WRITE(*,994) IF(MP.EQ.1) WRITE(*,995) SSE(MP) IF(MP.EQ.2) WRITE(*,997) SSE(MP) IF(MP.EQ.3) WRITE(*,998) SSE(MP) 306 CONTINUE ENDIF 994 FORMAT(5X,'NEGATIVE ENTROPY MEASURES FOR WEIGHTED, SQUARED 1 RESIDUALS') 995 FORMAT(8X,'X-PART: SSE =',E13.4/) 997 FORMAT(8X,'Y-PART: SSE =',E13.4/) 998 FORMAT(8X,'(X AND Y)-PARTS: SSE =',E13.4/) C DSX = DIVI*SUMSX DSY = DIVI*SUMSY DSQ = SB*DIVI DDX = DSQRT(DSX) DDY = DSQRT(DSY) DD = DSQRT(DSQ) C WRITE(*,1110) SB 1110 FORMAT(/5X,'WEIGHTED SUM OF SQUARES, S =',1PE17.7) WRITE(*,1141) DDX,DDY 1141 FORMAT(/5X,'STDEV OF X PART =',1PE11.3,5X,'STDEV OF Y PART =' 1,1PE11.3) WRITE(*,1135) DD 1135 FORMAT(/5X,'STANDARD DEVIATION OF FIT, SF =',1PE14.5/) C C C ####### MENU NO. 5: CHOICES AFTER FITTING ####### C 918 WRITE(*,924) IWE = 0 IPWC = 0 ALLV = 0 924 FORMAT(3X,'MENU #5: AFTER-FIT CHOICES'/ 1 5X,'1: RETURN TO MENU #1: DEFAULT CHOICES'/ 1 5X,'2: RETURN TO MENU #2: INPUT CHOICES'/ 1 5X,'3: RETURN TO MENU #3: WEIGHT CHOICES'// 2 5X,'4: PLOT RESIDUALS ON SCREEN AND IN GENOUTR; RTRN TO MENU'/ 2 5X,'5: SHOW RESIDUAL-RELATED QUANTITIES; RETURN TO MENU'/ 3 5X,'6: MAKE OR ADD TO RESIDUALS FILE(GENOUTT); RETURN TO MENU'/ 3 5X,'7: MAKE OR ADD TO OUTPUT FILE (GENOUTF); RETURN TO MENU'/ 4 5X,'8: EXIT FROM PROGRAM'/) READ(*,*) IMN5 C GOTO(931,932,933,934,935,936,937,938) IMN5 931 GOTO 317 932 GOTO 518 933 GOTO 125 C C PLOT RESIDUALS CHOICES BELOW 934 CONTINUE C WRITE(*,987) 987 FORMAT(3X,'PRODUCE VARIOUS RESIDUAL ROUGH PLOTS'/ 1 3X,'USE POSITIVE VALUES FOR SCREEN-ONLY PLOTTING'/ 1 3X,'* USE NEG. VALUES FOR SCREEN AND FILE (GENOUTR) OUTPUT * '/ 1 3X,'DEPENDING ON WEIGHTING CHOICES, NOT ALL PLOT CHOICES BELOW'/ 1 13X,'ARE ALWAYS OPERATIVE'// 1 5X,'1: X-RESIDUALS (RX) VS X'/ 1 5X,'2: Y-RESIDUALS (RY) VS YCALC'/ 2 5X,'3: Y-RESIDUALS VS X-RESIDUALS'/ 2 5X,'4: WEIGHTED X-RESIDUALS (FNX) VS X'/ 3 5X,'5: WEIGHTED Y-RESIDUALS (FNY) VS YCALC'/ 4 5X,'6: WEIGHTED Y-RESIDUALS VS WEIGHTED X-RESIDUALS'/ 5 5X,'7: WEIGHTED COMPOSITE RESIDUALS VS X'/ 6 5X,'8: RETURN TO MENU 5'/) READ(*,*) IRES IRX = 0 IRY = 0 SXX = 0 SYY = 0 DO 303 IJ = 1,4 SXX = SXX + SX(IJ) SYY = SYY + SY(IJ) 303 CONTINUE IF(SXX.LT.DM9) IRX = 1 IF(SYY.LT.DM9) IRY = 1 IF(IWTC.NE.5) THEN IF(SXX.EQ.SYY) THEN IRXY = 1 ELSE IRXY = 0 ENDIF ENDIF IF(IRES.LT.0) THEN OPEN(3,FILE='GENOUTR') WRITE(3,929) ALPHA ENDIF WRITE(*,929) ALPHA 929 FORMAT(/2X,A78) C GOTO(741,742,743,744,745,746,747,748) IABS(IRES) C 741 IF(IRX.EQ.0) CALL PLTR(X,RX,N,IRES) GOTO 934 742 IF(IRY.EQ.0) CALL PLTR(YN,RY,N,IRES) GOTO 934 743 IF(IRX.EQ.0.AND.IRY.EQ.0) CALL PLTR(RX,RY,N,IRES) GOTO 934 744 IF(IRX.EQ.0.AND.IRXY.EQ.0) CALL PLTR(X,FNX,N,IRES) GOTO 934 745 IF(IRY.EQ.0.AND.IRXY.EQ.0) CALL PLTR(YN,FNY,N,IRES) GOTO 934 746 IF(IRX.EQ.0.AND.IRY.EQ.0.AND.IRXY.EQ.0) CALL PLTR(FNX,FNY,N,IRES) GOTO 934 747 CALL PLTR(X,DI,N,IRES) GOTO 934 C 748 GOTO 918 C 935 GOTO 351 C 936 CONTINUE OPEN(19,FILE='GENOUTT') WRITE(19,929) ALPHA WRITE(19,81) IWTC IF(IWTC.EQ.3.OR.IWTC.EQ.4) WRITE(19,83) ILN IF(IWTC.EQ.3) WRITE(19,93) SIGX,SIGY IF(IWTC.EQ.5) THEN WRITE(19,182) UX,UY WRITE(19,82) PX,XIX,PY,XIY ENDIF C WRITE(19,913) WRITE(19,816) DO 947 I = 1,N WRITE(19,856) I,SX(I),XKEEP(I),X(I),DI(I) 947 CONTINUE C WRITE(19,817) DO 948 I = 1,N WRITE(19,856) I,SY(I),YKEEP(I),YN(I),SI(I) 948 CONTINUE C WRITE(19,818) DO 949 I = 1,N WRITE(19,856) I,RX(I),RY(I),FNX(I),FNY(I) 949 CONTINUE C 913 FORMAT(3X,'TABLES OF X-WEIGHTS:SX; X-DATA:x; X-ADJ:X; AND'/ 1 12X,' COMPOSITE WEIGHTED RESIDUALS:DI'/ 2 3X,'TABLES OF Y-WEIGHTS:SY; Y-DATA:y; Y-ADJ:Y; AND SQUARED'/ 2 12X,' WEIGHTED RESIDUALS: SI = DI*DI'/ 4 3X,'TABLES OF X-RESIDUALS:RX; Y-RESIDUALS:RY; WEIGHTED'/ 4 12X,' X-RESIDUALS:FNX; AND WEIGHTED Y-RESIDUALS:FNY'/) 816 FORMAT(5X,'I',8X,'SX',13X,'x',14X,'X',12X,'DI'/) 817 FORMAT(/5X,'I',8X,'SY',13X,'y',14X,'Y',12X,'SI'/) 818 FORMAT(/5X,'I',8X,'RX',13X,'RY',14X,'FNX',12X,'FNY'/) 856 FORMAT(1X,I5,1P4E15.5) C GOTO 918 C 937 CONTINUE OPEN(10,FILE='GENOUTF') WRITE(10,20) ALPHA,N,NP,NOUT,NF IF(IFCT.EQ.0) WRITE(10,21) EPS,EPS1,G,FAC,EPS2,EPS3 IMEND = 2 IFCT = 1 GOTO 361 666 CONTINUE C WRITE(10,81) IWTC IF(IWTC.EQ.3.OR.IWTC.EQ.4) WRITE(10,83) ILN IF(IWTC.EQ.3) WRITE(10,93) SIGX,SIGY IF(IWTC.EQ.5) THEN WRITE(10,182) UX,UY WRITE(10,82) PX,XIX,PY,XIY ENDIF WRITE(10,680) 680 FORMAT(/4X,19H *-*-*-*-*-*-*-*-* ) C IF(I23.EQ.1) THEN IFPM = 1 CALL PWGM(ALAM,DIVI,WX,WY,WXI,WYI,XKEEP,YKEEP,A,N,KM,IFPM) IFPM = 0 ENDIF C IF(IWTC.EQ.5) THEN WRITE(10,751) 751 FORMAT(/10X,'** ITERATIVE RE-WEIGHTING SOLUTION **') IF(IPWC.GT.MAXG) WRITE(10,157) IPWC 157 FORMAT(/5X,'! ITERATIVE RE-WEIGHING DID NOT CON 1VERGE IN',I3,2X,'ITERATIONS !'/) IF(IDIRC.GT.MAXG) THEN WRITE(10,902) IDIRC 902 FORMAT(8X,'!!$$$** N0 CONVERGENCE IN',I3,' ITERATIONS ** 1$$$!!'// 1 5X,'### RESET PARAMETER AND DATA VALUES, THEN DO MENU 1A #5'/ 1 5X,'$$$ AND/OR: GO TO MENU 1A #6, SET G<1, AND RERUN'/ 1 5X,'### AFTER CONVERGENCE, USE NEW PARAMETER VALUES WITH G=1'/ 1 5X,'AND RUN AGAIN AS A CHECK'/) ELSE WRITE(10,905) IDIRC 905 FORMAT(5X,27H ** CONVERGENCE ACHIEVED IN,I3,14H ITERATIONS **/) ENDIF ELSE C TEST WEIGHT CHOICE FOR CONSTANT LAMBDA DO 224 I = 1,N RL(I) = SY(I)/SX(I) IF(I.GT.1) THEN DIFF = DABS(RL(I) - RL(I-1)) IF(DIFF.GT.1.D-3) GOTO 86 ENDIF 224 CONTINUE ALAM(1) = RL(2)**2 WRITE(10,87) ALAM(1) 87 FORMAT(/5X,'LAMBDA =',1PE13.4) C 86 WRITE(10,900) 900 FORMAT(/10X,' *** GENERAL METHOD SOLUTION ***') IF(IDIRC.GT.MAXG) THEN WRITE(10,902) IDIRC ELSE WRITE(10,905) IDIRC ENDIF ENDIF CALL PRINTOUT(SUMSXR,SUMSYR,S2,A,AA,SAK,SAKR,RC,SDER,IMEND,NP) C WRITE(10,1033) WRITE(10,1035) AVURX,AVURY WRITE(10,1034) WRITE(10,1035) AVRX,AVRY,AVCOR CONTINUE C DO 308 MP = MPS,MPF IF(MP.EQ.MPS) WRITE(10,994) IF(MP.EQ.1) WRITE(10,995) SSE(MP) IF(MP.EQ.2) WRITE(10,997) SSE(MP) IF(MP.EQ.3) WRITE(10,998) SSE(MP) 308 CONTINUE C SUMDI = DIVN*SUMDI FM = SB - SUMDI*SUMDI*N D1SQ = FM*DIVI D1 = DSQRT(D1SQ) C IF(SUMSY.GT.1.D-20.AND.SUMSX.GT.1.D-20) THEN IRR = 1 ELSE IRR = 0 ENDIF C IF(MAXG.NE.0) THEN SSUMX = SUMSXR SSUMY = SUMSYR SB = SSUMX + SSUMY ELSE SSUMX = SUMSX SSUMY = SUMSY ENDIF WRITE(10,1110) SB WRITE(10,1141) DDX,DDY WRITE(10,1135) DD C IF(IRR.EQ.1) THEN CALL PEARSN(RX,RY,N,RU,PROB,Z) WRITE(10,1120) RU C CALL PEARSN(FNX,FNY,N,RR,PROB,Z) WRITE(10,1122) RR ENDIF 1120 FORMAT(5X,'CORR. COEFF. OF X-Y UNWEIGHTED RESIDUALS =',E14.4) 1122 FORMAT(5X,'CORR. COEFF. OF X-Y *WEIGHTED* RESIDUALS =',E14.4) C IMEND = 0 WRITE(10,680) WRITE(10,680) GOTO 918 938 CONTINUE CLOSE(10) CLOSE(3) CLOSE(9) CLOSE(19) STOP END