SUBROUTINE MMRQT(NAME,NN,NMM,NXM,NYM,NM,NX,NY,NP,NS,MW,MD,DP0,DPP, 1 IP,DT,DV,DN,NSTEP,DSTEP,DW,AC,IC,FFIT,FPRT) ************************************************************************ * * * Interactive MicroSoft Fortran 5.1 version -- creates * * printout of the results (88) and fitted values+residuals (99) * * * * Multiple non-linear weighted least squares parameter estimation * * using the MARQUARDT algorithm to find minimum chi-square. Sign of * * parameters is not allowed to change (see label 247). There is no * * upper limit for the number of either the independent or dependent * * variables or parameters, or measured points. The execution time * * is sensitive to the number of parameters x the number of data po- * * ints x the number of functions (The use of analytical derivatives * * reduces computation time almost to half of that using numerical * * derivatives.) * * * * Input data: NAME is an identifier string of length 40 * * NN= number of data blocks * * NMM= max. number of data points * * NXM= max. number of independent variables * * NYM= max. number of dependent variables * * NM= number of data points * * NX= number of independent variables * * NY= number of dependent variables * * NP= number of parameters * * NS= number of param.+const. * * DT(NN*NMM*NXM) contains (measured) INDEPENDENT varbles * * DV(NN*NMM*NYM) contains (measured) DEPENDENT variables * * DN(NN*NMM*NYM) contains the uncertainty values (MW=5) * * (rows: X(1),X(2),....,X(NX) ; Y(1),Y(2),......,Y(NY) * ******** ******** * WARNING !!! The DIMENSION statement in the (main) program * * calling MMRQT s h o u l d specify the actual * * values of NM,NX and NY, if NX>1 or NY>1. If both NX * * and NY are equal 1, the NM dimension might be greater * * then the actual NM value. ( Never NX or NY !) * ******** ******** * MW is a flag: * * if MW=0 : equally weighted Y(I)-s * * if MW=1 : relative weights; DW(I,I)=1/Y(I)**2 * * DW(I,J)=0 * * if MW=2 : weigts assigned before calling MMRQT * * if MW=3 : weights calculated from user defined * * function named SUBROUTINE USERW * * if MW=4 : Poisson weights; DW(I,I)=1/Y(I) * * DW(I,J)=0 * * (default=0, if not specified) * * if MW=5 : uncertainty weights; DW(I,I)=1/U(I)**2 * * DW(I,J)=0 * * MD is another flag: * * if MD=0 : derivatives evaluated NUMERICALLY * * if MD=1 : derivatives evaluated ANALYTICALLY * * in the subroutine DERJACOBI (written * * by the user and linked to MMRQT) * * (default=0, if not specified) * * DW(NY,NY) is the weighting matrix (see above, MW=2) * * DP0(NS) contains the initial guess of the * * parameters. (Signs remain during iteration !!!) * * DPP(NS) contains the in. guesses and const. * * IP(NS) contains the assignment of parameters, based * * on the 0's and 1's specifying fixed or * * fitting parameters. * * DSTEP= convergence criterion in parameter space. * * (default= 1E-9 if not specified) * * NSTEP= maximum number of iterative steps if no con- * * vergence below DSTEP. (default= 50) * * AC is a real dummy array to transmit max. 10 eventual * * fixed parameters to evaluate the model. * * IC is an integer to transmit eventual fixed integer * * values to evaluate the model function. * * RESULTS will be displayed within this subroutine after iteration. * * Two output files will be created. The one with the filename con- * * tained in FPRT (logical unit 88) with the parameters and statis- * * tical analysis results, the other with the filename contained in * * FFIT (logical unit 99) listing the fitted Yi values and calcula- * * ted residuals for each input point to facilitate graphical * * evaluation of the fit. * * * * * * * * * * * * * SUBROUTINE MODEL calculates the actual Y value(s) from the * * actual X value(s). ( See " CALL MODEL ".) * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DT(13,226,1),DV(13,226,1),DN(13,226,1), 1 DW(9,9),DA(15,15),DC(15,15),DU(15,15),DG(15,15),DP(15), 2 DP0(15),DPP(15),IP(15),DD(15),DBW(15),CCY(13,226,1), 3 RES(13,226,1),DC0(15),D0C(15),DX(15),DY(15),DB(15),NM(1), 4 NX(1),NY(1),AC(10) CHARACTER NAME*40,CHR*1,FFIT*20,FPRT*20 NO=0 DO 1112 MN=1,NN NO=NO+NY(MN)*NM(MN) 1112 CONTINUE IF(DSTEP.EQ.0.)DSTEP=1.D-9 IF(NSTEP.EQ.0)NSTEP=50 IF(MW.EQ.0)THEN DO 202 I=1,NYM DO 202 J=1,NYM IF(I.EQ.J)THEN DW(I,J)=1. ELSE DW(I,J)=0. ENDIF 202 CONTINUE ENDIF IF (MW.EQ.1)GOTO 2203 IF (MW.EQ.3)GOTO 2203 IF (MW.EQ.4)GOTO 2203 IF (MW.EQ.5)GOTO 2203 GOTO 203 2203 DSW=0. DWS=1. *** DWS (^) is for scaling rel. or Poisson weights. *************** DO 204 MN=1,NN DO 204 M=1,NM(MN) IF(MW.EQ.1)CALL RELW(NN,MN,NMM,NYM,M,NY(MN),DV,DW,DWS) IF(MW.EQ.5)CALL UNCERW(NN,MN,NMM,NYM,M,NY(MN),DN,DW,DWS) IF(MW.EQ.4)CALL POISSONW(NN,MN,NMM,NYM,M,NY(MN),DV,DW,DWS) IF(MW.EQ.3)THEN IER=0 CALL USERW(NN,MN,NMM,NXM,NYM,M,NX(MN),NY(MN),NP,NS,DT,DV, 1 DW,DWS,DP,AC,IC) IF(IER.GT.0)RETURN ENDIF DO 205 I=1,NY(MN) 205 DSW=DSW+DW(I,I) 204 CONTINUE DWS=DBLE(NO)/DSW 203 NEI=0 NEF=0 NES=0 DO 615 I=1,NP 615 DP(I)=DP0(I) DPM=.01 DO 1990 MN=1,NN WRITE(*,1117)MN 1990 WRITE(*,701)NM(MN) 701 FORMAT(I10,' data points') C *** Starting the iteration.********************* DF=DSSQRS(NN,NMM,NXM,NYM,NM,NX,NY,NP,NS,MW,DT,DV,DW,DWS,DP, 1 DPP,IP,AC,IC) IF(DF.GT.9.D9)THEN WRITE(*,933)DF ELSE WRITE(*,300)DF ENDIF 300 FORMAT(1H0,'Starting sum of squared residuals=',G20.10) 933 FORMAT(1H0,'Starting sum of squared residuals=',E12.6) DO 206 K=1,NP 206 WRITE(*,301)K,DP(K) 301 FORMAT(1H ,16X,2HP(,I2,2H)=,G20.10) C *** Iteration loop begins ********************* DO 207 IT=1,NSTEP DO 208 K=1,NP 208 DU(K,1)=DP(K) DFR=DF DO 209 K=1,NP DB(K)=0. DO 209 L=1,K DC(K,L)=0. 209 CONTINUE DO 210 MN=1,NN DO 210 M=1,NM(MN) DO 211 I=1,NX(MN) 211 DX(I)=DT(MN,M,I) IF(MW.EQ.1)CALL RELW(NN,MN,NMM,NYM,M,NY(MN),DV,DW,DWS) IF(MW.EQ.5)CALL UNCERW(NN,MN,NMM,NYM,M,NY(MN),DN,DW,DWS) IF(MW.EQ.4)CALL POISSONW(NN,MN,NMM,NYM,M,NY(MN),DV,DW,DWS) IF(MW.EQ.3)THEN IER=0 CALL USERW(NN,MN,NMM,NXM,NYM,M,NX(MN),NY(MN),NP,NS,DT,DV, 1 DW,DWS,DP,AC,IC) IF(IER.GT.0)RETURN ENDIF C *** Calculating the Jacobian DGij = dYi/dPj *********** IF(MD.NE.0) GOTO 220 CALL NUMJACOBI(MN,NX(MN),NY(MN),NP,NS,DX,DY,DG,DP,DPP,IP,AC,IC) GOTO 222 220 CALL DERJACOBI(MN,NX(MN),NY(MN),NP,NS,DX,DY,DG,DP,DPP,IP,AC,IC) C *** Calculating second derivatives of chi-square *********** 222 DO 230 K=1,NP DO 231 L=1,K DAA=0. DO 232 I=1,NY(MN) DO 232 J=1,NY(MN) 232 DAA=DAA+DW(I,J)*DG(I,L)*DG(J,K)*DP(L)*DP(K) DC(K,L)=DC(K,L)+DAA 231 CONTINUE DAA=0. DO 233 I=1,NY(MN) DO 233 J=1,NY(MN) 233 DAA=DAA+DW(I,J)*DG(J,K)*(DV(MN,M,I)-DY(I))*DP(K) DB(K)=DB(K)+DAA 230 CONTINUE 210 CONTINUE C *** Normalization of correlation matrix. *************** DTR=0. DO 234 I=1,NP DC0(I)=DC(I,I) DTR=DTR+DC(I,I) 234 CONTINUE DTR=DTR/DBLE(NP)/1000. DO 235 I=1,NP IF(DC0(I).LE.DTR)THEN DC0(I)=1. ELSE DC0(I)=DSQRT(DC0(I)) ENDIF 235 CONTINUE DO 236 I=1,NP DO 236 J=1,I 236 DC(I,J)=DC(I,J)/DC0(I)/DC0(J) C *** Marquardt modification ********************** 253 DO 237 I=1,NP DO 238 J=1,I-1 238 DA(I,J)=DC(I,J) DA(I,I)=DC(I,I)+DPM 237 CONTINUE C *** Inversion of DA matrix ********************** IER=0 N=NP+1 CALL INVPDS(NP,DA,IER) IF(IER.EQ.1)GOTO 801 C *** DA inverted. Solution for estimation next *** DO 245 I=1,NP DS=0. DO 246 J=1,NP 246 DS=DS+DA(I,J)/DC0(J)*DB(J) DD(I)=DS/DC0(I) 245 CONTINUE DSL=0. DXI=1. DO 247 I=1,NP C *** Next statement does not allow sign change of the parameters *** IF(DXI*DD(I).LE.-.95)DXI=-.95/DD(I) DSL=DSL+DD(I)*DD(I) 247 CONTINUE DSL=DSL*DXI*DXI C *** Solution for parameters calculated. DSL=squared step-length. *** DO 248 I=1,NP 248 DP(I)=DU(I,1)*(1.+DXI*DD(I)) ******** start rescaling USER weights ********************** IF (MW.EQ.1)GOTO 3203 IF (MW.EQ.3)GOTO 3203 IF (MW.EQ.4)GOTO 3203 IF (MW.EQ.5)GOTO 3203 GOTO 6203 3203 DSW=0. DWS=1. *** DWS (^) is for rescaling USER weights. *************** DO 6204 MN=1,NN DO 6204 M=1,NM(MN) IF(MW.EQ.1)CALL RELW(NN,MN,NMM,NYM,M,NY(MN),DV,DW,DWS) IF(MW.EQ.5)CALL UNCERW(NN,MN,NMM,NYM,M,NY(MN),DN,DW,DWS) IF(MW.EQ.4)CALL POISSONW(NN,MN,NMM,NYM,M,NY(MN),DV,DW,DWS) IF(MW.EQ.3)THEN IER=0 CALL USERW(NN,MN,NMM,NXM,NYM,M,NX(MN),NY(MN),NP,NS,DT,DV, 1 DW,DWS,DP,AC,IC) IF(IER.GT.0)RETURN ENDIF DO 6205 I=1,NY(MN) 6205 DSW=DSW+DW(I,I) 6204 CONTINUE DWS=DBLE(NO)/DSW 6203 CONTINUE DF=DSSQRS(NN,NMM,NXM,NYM,NM,NX,NY,NP,NS,MW,DT,DV,DW,DWS,DP, 1 DPP,IP,AC,IC) IF(DF.GT.9.D9)THEN WRITE(*,934)IT,DSL,DPM,DF ELSE WRITE(*,302)IT,DSL,DPM,DF ENDIF 302 FORMAT(1X,I3,2X,'Step=',1PE7.1,2X,'Marq.pm=',1PE6.0,2X,'Sum of squ 1d.res.=',0PF21.10) 934 FORMAT(1X,I3,2X,'Step=',1PE7.1,2X,'Marq.pm=',1PE6.0,2X,'Sum of squ 1d.res.=',E12.6) IF(DF.GE.DFR)GOTO 250 DO 249 K=1,NP 249 WRITE(*,301)K,DP(K) 250 IF(DSL.LE.DSTEP)THEN NEI=0 GOTO 252 ENDIF C *** Setting the Marquardt-parameter ************** IF(DF.LE.DFR)GOTO 251 801 DPM=DPM*10. GOTO 253 251 DPM=DPM/10. IF(DPM.LT..000001)DPM=.000001 207 CONTINUE C *** End of the iteration loop ********************* NEI=1 C *** Flag indicating not sufficient convergence **** 252 IF(DF.LT.DFR)GOTO 254 DF=DFR DO 260 I=1,NP 260 DP(I)=DU(I,1) 254 NF=NO-NP DNF=DBLE(NF) DSE=DSQRT(DF/DNF) C *** Calculating variances & correlation matrix of the parameters ***** DO 261 I=1,NP DO 261 J=1,I 261 DA(I,J)=DC(I,J) CALL INVPDS(NP,DA,IER) IF(IER.EQ.1)THEN NES=1 IER=0 C *** Flag indicating singular Hessian matrix. (Too many params.) GOTO 262 ENDIF DO 263 I=1,NP DB(I)=DSQRT(DF/DNF*DA(I,I)/DC0(I)/DC0(I)) D0C(I)=DSQRT(DA(I,I)) 263 CONTINUE DO 264 I=1,NP DO 264 J=1,NP 264 DC(I,J)=DINT(1000.*DA(I,J)/D0C(I)/D0C(J)+.5)/1000 262 DO 265 I=1,NP DO 265 J=1,I DA(I,J)=DC(I,J) DA(J,I)=DC(I,J) 265 CONTINUE C *** Principal component analysis of the correlation matrix ******* CALL EIGSLT(NP,DA,DU) DO 266 I=1,N DO 266 J=2,N 266 DU(I,J)=DNINT(10000.*DU(I,J))/10000. IF(DU(1,N).LE..000001)NEF=1 C *** Flag for too few parameters, or not suitable model ********** C C *** End of fitting & estimation ************************************** C *** Printout of results ********************************************** OPEN(88,FILE=FPRT) Cc NNW=0 Cc416 WRITE(*,400) Cc400 FORMAT(1H0,//////,15X,'Do you want to Display results on the scree Cc 1n (D)',//,20X,'or Skip them and put into PRT.OUT (S) ?' Cc 2,///,15X,'Press a key & , or S & ! ',$) Cc READ(*,401)CHR Cc401 FORMAT(A1) Cc IF(CHR.EQ.'S')NNW=1 Cc IF(CHR.EQ.'s')NNW=1 WRITE (88,310)NAME 777 WRITE(*,310)NAME 310 FORMAT( 1H1, //,10X,'LEAST SQUARES PARAMETER ESTIMATION - Marquard *t method',///,15X,'DATA from file:',A40,///) WRITE(88,3101)NN,NP WRITE(*,3101)NN,NP DO 1113 MN=1,NN WRITE(88,311)MN,NY(MN),NX(MN),NM(MN) 1113 WRITE(*,311)MN,NY(MN),NX(MN),NM(MN) WRITE(88,3102)DSTEP WRITE(*,3102)DSTEP 3101 FORMAT( 15X,'Number of BLOCKS.................',I2,/,15X, 1'Number of PARAMETERS.............',I2,/) 311 FORMAT( 15X, 'Block Number.....................',I2,/,15X, 1 'Number of FUNCTIONS..............',I2,/,15X, 2 'Number of independent VARIABLES..',I2,/,15X, 3 'Number of DATA POINTS............',I3,/) 3102 FORMAT(15X,'Convergency limit................',1PE8.1) IF(MW.EQ.0)WRITE(88,312) IF(MW.EQ.0)WRITE(*,312) IF(MW.EQ.1)WRITE(88,313) IF(MW.EQ.1)WRITE(*,313) IF(MW.EQ.4)WRITE(88,3313) IF(MW.EQ.4)WRITE(*,3313) IF(MW.EQ.5)WRITE(88,3513) IF(MW.EQ.5)WRITE(*,3513) IF(MW.EQ.2)WRITE(88,314) IF(MW.EQ.2)WRITE(*,314) IF(MW.EQ.3)WRITE(88,314) IF(MW.EQ.3)WRITE(*,314) 3312 FORMAT( /,15X,'Starting parameters:',//,2X,15(1PE10.3)) 312 FORMAT( 17X,'Equal weights') 313 FORMAT( 17X,'Relative weights') 3313 FORMAT( 17X,'Poisson weights') 3513 FORMAT( 17X,'Uncertainty weights') 314 FORMAT( 17X,'User defined weights') WRITE(88,3312)(DP0(IKI),IKI=1,NP) WRITE(*,3312)(DP0(IKI),IKI=1,NP) Cc IF(NNW.EQ.0) CALL SCROLL(IP) WRITE(88,315) 315 FORMAT( ////,30X,'D A T A P O I N T S',//,2X,75('-')) DO 1115 MN=1,NN WRITE(88,1117)MN DO 267 I=1,NX(MN) 267 WRITE(88,317)I 317 FORMAT( 9X,'X',I2,$) DO 268 I=1,NY(MN) 268 WRITE(88,318)I 318 FORMAT( 9X,'Y',I2,$) IF (MW.NE.5) GOTO 3181 DO 2681 I=1,NY(MN) 2681 WRITE(88,3182)I 3182 FORMAT( 9X,'U',I2,$) 3181 WRITE(88,319) 319 FORMAT( /,2X,75('-')) DO 270 IM=1,30 WRITE(88,320)IM 320 FORMAT( 2X,I3,1X,$) IF (MW.NE.5) GOTO 3211 WRITE(88,3212)(DT(MN,IM,I),I=1,NX(MN)),(DV(MN,IM,I), 1 I=1,NY(MN)),(DN(MN,IM,I),I=1,NY(MN)) GOTO 3213 3211 WRITE(88,321)(DT(MN,IM,I),I=1,NX(MN)),(DV(MN,IM,I),I=1,NY(MN)) 321 FORMAT( 0P,10(F10.4)) 3212 FORMAT( 0P,15(F10.4)) 3213 IF(MOD(IM,5).EQ.0)WRITE(88,336) IF(IM.EQ.NM(MN))GOTO 277 270 CONTINUE 271 WRITE(88,319) KM=31 272 WRITE(88,322) 322 FORMAT(1H1,//,2X,75('-')) DO 273 I=1,NX(MN) 273 WRITE(88,317)I DO 274 I=1,NY(MN) 274 WRITE(88,318)I WRITE(88,319) DO 275 IM=KM,KM+49 WRITE(88,320)IM IF (MW.NE.5) GOTO 3214 WRITE(88,3212)(DT(MN,IM,I),I=1,NX(MN)),(DV(MN,IM,I), 1 I=1,NY(MN)),(DN(MN,IM,I),I=1,NY(MN)) GOTO 3215 3214 WRITE(88,321)(DT(MN,IM,I),I=1,NX(MN)),(DV(MN,IM,I),I=1,NY(MN)) 3215 IF(MOD(IM,5).EQ.0)WRITE(88,336) IF(IM.EQ.NM(MN))GOTO 277 275 CONTINUE KM=KM+50 GOTO 272 277 WRITE(88,733) 733 FORMAT( 2X,75('-')) 1115 CONTINUE DO 289 I=1,NP 289 DBW(I)=DABS(DB(I)*DP(I)) WRITE(88,323)NAME WRITE(*,323)NAME 323 FORMAT(1H1,///,20X,A40,////,10X,'Estimates of the PARAMETERS',5X, 1'of their STANDARD DEVIATIONS',/) DO 278 I=1,NP WRITE(88,324)I,DP(I),DBW(I) 278 WRITE(*,324)I,DP(I),DBW(I) 324 FORMAT( /,6X,I2,3X,F25.10,5X,1PE12.3) Cc IF(NNW.EQ.0)CALL SCROLL(IP) WRITE(88,30) WRITE(*,30) 30 FORMAT( ////,25X,'95 % confidence intervals',//,15X,'half-width', 18X,'lower limit - upper limit',/) TST=T95CNF(DNF) DO 31 I=1,NP WRITE(88,32)I,DBW(I)*TST,DP(I)-DBW(I)*TST,DP(I)+DBW(I)*TST 31 WRITE(*,32)I,DBW(I)*TST,DP(I)-DBW(I)*TST,DP(I)+DBW(I)*TST 32 FORMAT(8X,'P(',I2,')',2X,1PE12.3,5X,E12.3,2X,'-',E12.3,/) WRITE(88,325) WRITE(*,325) 325 FORMAT( ////, 5X,'Principal component analysis of correlation matr 1ix',///,5X,'EIGENVALUE',$) DO 279 I=1,NP WRITE(88,326)I 279 WRITE(*,326)I 326 FORMAT( 4X,'P(',I2,')',1X,$) WRITE(88,336) WRITE(*,336) DO 280 I=1,NP WRITE(88,327)DU(1,I) WRITE(*,327)DU(1,I) WRITE(88,1327)(DU(J,I),J=2,NP+1) WRITE(*,1327)(DU(J,I),J=2,NP+1) 280 CONTINUE 327 FORMAT( /,8X,F5.3,5X,$) 1327 FORMAT( F5.3,7(4X,F5.3)) WRITE(88,328) WRITE(*,328) IF(NEI.EQ.1)WRITE(88,329) IF(NEI.EQ.1)WRITE(*,329) IF(NEF.EQ.1)WRITE(88,330) IF(NEF.EQ.1)WRITE(*,330) IF(NES.EQ.1)WRITE(88,331) IF(NES.EQ.1)WRITE(*,331) 328 FORMAT( ///) 329 FORMAT( 10X,'Desired convergence not obtained') 330 FORMAT( 10X,'Not correctly defined model.(Maybe too few parameters 1)') 331 FORMAT( 10X,'Singular Hessian matrix.(Maybe too many parameters)') Cc IF(NNW.EQ.0)CALL SCROLL(IP) WRITE(88,332) WRITE(*,332) 332 FORMAT(1H1,///,10X,'Fit of the model at the actual data points', 1//) DSS=0. DVS=0. DO 281 MN=1,NN WRITE(88,1117)MN WRITE(*,1117)MN 1117 FORMAT(//,31X,'BLOCK No. ',I2,/) DO 281 I=1,NY(MN) WRITE(88,333)I WRITE(*,333)I 333 FORMAT(/,30X,I2,'. FUNCTION',/) C***** printout of the fit omitted - results in FIT.OUT ******* DS2=0. DE1=0. DE2=0. DSN = 0. DO 282 M=1,NM(MN) DO 283 J=1,NX(MN) 283 DX(J)=DT(MN,M,J) CALL MODEL(MN,NX(MN),NY(MN),NP,NS,DX,DY,DP,DPP,IP,AC,IC) CCY(MN,M,I)=DY(I) RES(MN,M,I)=DV(MN,M,I)-DY(I) IF(M.GT.1) DSN = DSN + (RES(MN,M-1,I)-RES(MN,M,I))**2 DS2=DS2+RES(MN,M,I)**2 DSH = DSN/DS2 DE1=DE1+DV(MN,M,I) DE2=DE2+DV(MN,M,I)*DV(MN,M,I) IF(DV(MN,M,I).EQ.0.)DV(MN,M,I)=.0000000000001 335 FORMAT(2X,I3,2X,2(G18.10E2),2X,1PE15.7,6X,2PF6.1) 336 FORMAT(1X) 337 FORMAT(1H1,//) 282 CONTINUE DVV=(DE2-DE1*DE1/DBLE(NM(MN)))/DBLE(NM(MN)-1) DRV=DS2/DBLE(NM(MN)-NP-1) DYM=DE1/DBLE(NM(MN)) DGFIT=100*(1-DS2/DVV/DBLE(NM(MN)-1)) DSS=DSS+DS2 DVS=DVS+DVV*DBLE(NM(MN)-1) WRITE(88,338)DYM,DVV,DRV,DGFIT,DSH WRITE(*,338)DYM,DVV,DRV,DGFIT,DSH 338 FORMAT( /,20X,'MEAN=',1PE10.4,8X,'VARIANCE=',1PE10.3,//,16X, 1'RESIDUAL=',1PE10.3,13X,'FIT=',0PF7.3,'%',//,18X, 2'D-STAT=',1PE10.3,//) Cc IF(NNW.EQ.0)CALL SCROLL(IP) 281 CONTINUE C OPEN(99,FILE=FFIT) WRITE(99,3315)NAME 3315 FORMAT(A40) WRITE(99,*)'Estimated parameters:' WRITE(99,3316)(DP(IKI),IKI=1,NP) 3316 FORMAT(15(1PE10.3)) WRITE(99,*)'Calculated values and residuals at the data points:' DO 3270 MN=1,NN WRITE(99,1117)MN DO 3267 I=1,NX(MN) 3267 WRITE(99,3317)I 3317 FORMAT( 10X,'X',I2,$) DO 3268 I=1,NY(MN) 3268 WRITE(99,3318)I 3318 FORMAT( 10X,'Y',I2,$) DO 3269 I=1,NY(MN) IF(I.LT.NY(MN)) WRITE(99,3319)I 3269 CONTINUE WRITE(99,4119)I-1 3319 FORMAT( 8X,'RES',I2,$) 4119 FORMAT( 8X,'RES',I2) Cc WRITE(99,401)' ' DO 3270 M=1,NM(MN) WRITE(99,3320)(DT(MN,M,I),I=1,NX(MN)),(CCY(MN,M,I),I=1,NY(MN)), 1 (RES(MN,M,I),I=1,NY(MN)) 3320 FORMAT(20(1PE12.4)) 3270 CONTINUE CLOSE(99) C DGFIT=100*(1-DSS/DVS) DWWE=DSQRT(DF/DNF) IF(DWS.EQ.0.)DWS=1.D-040 WRITE(88,339)IT,DF,DF/DWS,DF/DWS/NF,NF,DWWE,TST,DGFIT WRITE(*,339)IT,DF,DF/DWS,DF/DWS/NF,NF,DWWE,TST,DGFIT 339 FORMAT( /,20X,'Number of iterations............', 1 I3,/,20X,'Sum of squared residuals........', 2 1PE12.6,/,20X,'Chi2 (if proper weights used)...', 3 0PF7.3,/,20X,'Reduced chi2 (if pr. weights)...', 4 F7.3,/,20X,'Number of degrees of freedom....', 5 I3,/,20X,'Estimated weighted error........', 6 1PE9.3,/,20X,'Critical t (95 % confidence)....', 7 0PF5.3,/,20X,'Goodness of overall fit.........',F7.3,'%', 8 //,13X,'Calculated values of the model function and re 9siduals',/,23X,'are listed in the file FIT.OUT') Cc IF(NNW.EQ.0)CALL SCROLL(IP) WRITE(88,340) WRITE(*,340) 340 FORMAT( ///,15X,'Correlation matrix of the parameters',///,13X,$) DO 285 I=1,NP WRITE(88,226)I 285 WRITE(*,226)I 226 FORMAT( 2X,'P(',I2,')',3X,$) WRITE(88,341) WRITE(*,341) 341 FORMAT( /) DO 286 I=1,NP WRITE(88,342)I WRITE(*,342)I 342 FORMAT(5X,'P(',I2,')',1X,$) WRITE(88,343) (DC(I,J),J=1,I) WRITE(*,343) (DC(I,J),J=1,I) 343 FORMAT( 8( 4X,F5.3)) WRITE(88,336) WRITE(*,336) 286 CONTINUE WRITE(88,344) 344 FORMAT( /////,10X,'Covariance matrix of the parameters',///,8X,$) DO 287 I=1,NP 287 WRITE(88,350)I 350 FORMAT( 8X,'P(',I2,')',1X,$) WRITE(88,341) DO 288 I=1,NP WRITE(88,352)I 352 FORMAT( 3X,'P(',I2,')',3X,$) DO 299 J=1,I DPO=DC(I,J)*DBW(I)*DBW(J) WRITE(88,353)DPO 353 FORMAT( 3X,1PE10.4,$) 299 CONTINUE WRITE(88,888) 288 CONTINUE Cc IF(NNW.EQ.1)RETURN Cc290 WRITE(*,402) 888 FORMAT( /) Cc402 FORMAT(1H0,//, 6X,'To Display the results AGAIN on the screen, pre Cc 1ss D & !',//, 6X,'If you don`t want to see them,just pres Cc 2s ! ',$) Cc READ(*,401)CHR Cc IF(CHR.EQ.' ')RETURN Cc IF(CHR.EQ.'D')NNW=0 Cc IF(CHR.EQ.'d')NNW=0 Cc GOTO 777 END C SUBROUTINE NUMJACOBI(MN,NX,NY,NP,NS,DX,DY,DG,DP,DPP,IP,AC,IC) ************************************************************************ * Calculates the Jacobian matrix of the derivatives dYi/dPj with * * a numerical method using finite differences, which needs twice * * the evaluation of Yi to have one derivative dYi/dPj = DG(i,j). * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DX(NX),DG(15,15),DP(NS),DY(NY),D(15),AC(10) DO 20 J=1,NP DE=.001*DABS(DP(J))+.000001 DP(J)=DP(J)+DE CALL MODEL(MN,NX,NY,NP,NS,DX,DY,DP,DPP,IP,AC,IC) DO 21 I=1,NY 21 DG(I,J)=DY(I)/DE DP(J)=DP(J)-DE D(J)=DE 20 CONTINUE CALL MODEL(MN,NX,NY,NP,NS,DX,DY,DP,DPP,IP,AC,IC) DO 22 I=1,NY DO 22 J=1,NP 22 DG(I,J)=DG(I,J)-DY(I)/D(J) C *** DG is the actual Jacobian. ************************ RETURN END C SUBROUTINE INVPDS(NP,DA1,IER) ************************************************************************ * Inversion of positive definite symmetrical matrix DA(N*N) * * Inverse:DA ; IER=0 ,if O.K. ; IER=1 ,if not positive definit * ************************************************************************ REAL*8 DA(15,15),DA1(15,15),DD N=NP+1 DO 990 I=2,N DO 990 J=2,N 990 DA(I,J)=DA1(I-1,J-1) DO 101 K=N,2,-1 IF (DA(2,2).LE.0.)GOTO 109 DA(1,N)=1/DA(2,2) DO 102 I=3,N DD=DA(I,2)*DA(1,N) IF(I.GT.K)THEN DA(1,I-1)=DD ELSE DA(1,I-1)=-DD ENDIF DO 103 J=3,I 103 DA(I-1,J-1)=DA(I,J)+DA(I,2)*DA(1,J-1) 102 CONTINUE DO 104 I=2,N 104 DA(N,I)=DA(1,I) 101 CONTINUE DO 105 I=2,N-1 DO 105 J=I+1,N DA(I,J)=DA(J,I) 105 CONTINUE IER=0 GOTO 110 109 IER=1 110 DO 991 I=1,N-1 DO 991 J=1,N-1 991 DA1(I,J)=DA(I+1,J+1) RETURN END C SUBROUTINE EIGSLT(NP,DA1,DU) ************************************************************************ * Calculates the eigenvalues of DA(N*N) symmetric matrix using lower * * triangle - JACOBI method. Eigenvalues & eigenvectors in DU(N+1*N). * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DU(15,15),DA(15,15),DA1(15,15) N=NP+1 DO 992 I=2,N DO 992 J=2,N 992 DA(I,J)=DA1(I-1,J-1) DO 111 I=2,N DO 111 J=2,N IF (I.EQ.J)THEN DU(I,J)=1. ELSE DU(I,J)=0. ENDIF 111 CONTINUE DV=0. DO 112 I=3,N DO 112 J=2,I-1 112 DV=DV+DABS(DA(I,J)) IF(DV.EQ.0.)GOTO 696 DN1=DBLE(N) DV0=DV/DN1/DN1*.000005 DV1=0. 130 DV=DV/DN1 131 DO 113 I0=3,N DO 113 J0=2,I0-1 IF(DABS(DA(I0,J0)).LE.DV)GOTO 113 DV1=1. IF(DA(J0,J0).EQ.DA(I0,I0))THEN DT=1. GOTO 114 ELSEIF(DA(J0,J0).GT.DA(I0,I0))THEN DV2=1. ELSE DV2=-1. ENDIF DV3=DABS(DA(J0,J0)-DA(I0,I0)) DV4=DSQRT((DA(J0,J0)-DA(I0,I0))**2+4.*DA(I0,J0)**2) DT=2.*DA(I0,J0)*DV2/(DV3+DV4) 114 DC=1/DSQRT(1+DT**2) DS=DT*DC DC1=DC**2 DS1=DS**2 DT1=DT**2 DV5=DA(I0,I0) DA(I0,I0)=DC1*(DV5-2.*DT*DA(I0,J0)+DT1*DA(J0,J0)) DA(J0,J0)=DC1*(DA(J0,J0)+2.*DT*DA(I0,J0)+DT1*DV5) DA(I0,J0)=0. DO 115 J=2,J0-1 DV5=-DS*DA(J0,J)+DC*DA(I0,J) DA(J0,J)=DC*DA(J0,J)+DS*DA(I0,J) DA(I0,J)=DV5 115 CONTINUE DO 116 I=J0+1,I0-1 DV5=-DS*DA(I,J0)+DC*DA(I0,I) DA(I,J0)=DC*DA(I,J0)+DS*DA(I0,I) DA(I0,I)=DV5 116 CONTINUE DO 117 I=I0+1,N DV5=-DS*DA(I,J0)+DC*DA(I,I0) DA(I,J0)=DC*DA(I,J0)+DS*DA(I,I0) DA(I,I0)=DV5 117 CONTINUE DO 118 I=2,N DV5=DC*DU(I,I0)-DS*DU(I,J0) DU(I,J0)=DS*DU(I,I0)+DC*DU(I,J0) DU(I,I0)=DV5 118 CONTINUE 113 CONTINUE IF(DV1.EQ.1.)THEN DV1=0. GOTO 131 ELSEIF(DV.GE.DV0)THEN GOTO 130 ENDIF C Sorting the eigen's DO 119 I=2,N 119 DU(1,I)=DA(I,I) DO 120 I=3,N L=I 132 IF(DU(1,L-1).GE.DU(1,L))GOTO 120 DO 121 J=1,N DT=DU(J,L-1) DU(J,L-1)=DU(J,L) DU(J,L)=DT 121 CONTINUE IF(L.GT.2)THEN L=L-1 GOTO 132 ENDIF 120 CONTINUE 696 DO 697 I=1,N-1 DO 697 J=1,N-1 697 DU(I,J)=DU(I,J+1) RETURN END C FUNCTION T95CNF(DNF) ************************************************************************ * Calculates the critical t value for 95% confidence. * * (t: realisation of a Student distribution random variable) * * NF= Number of degrees of Freedom. * ************************************************************************ REAL*8 DT,DNF NF=INT(DNF) IF(NF.GT.5)GOTO 140 DT=6.145+(5.-DNF)*(.289+(4.-DNF)*(.0575+(3.-DNF)*(.022+(2.-DNF)* 1.020552))) GOTO 150 140 IF(NF.GT.30)GOTO 141 DT=7.6278-.2316*DNF+.00421*DNF*DNF-.186*DSIN(.2214*DNF)-.0116* 1DSIN(.4428*DNF)+.0186*DSIN(.4428*DNF) GOTO 150 141 IF(NF.LT.40)THEN DT=4.4067+.0296*(30.-DNF) GOTO 150 ELSEIF(NF.LT.60)THEN DT=4.1108+.021*(40.-DNF) GOTO 150 ELSE DT=3.6888+.0116*(60.-DNF) ENDIF 150 DT=DEXP(DT) IF(NF.EQ.7.AND.NF.EQ.8.AND.NF.EQ.13.AND.NF.EQ.15.AND.NF.EQ.24) 1DT=DT+2 DT=DINT(DT+1960.5)/1000. T95CNF=REAL(DT) RETURN END C SUBROUTINE RELW(NN,MN,NMM,NYM,M,NY,DV,DW,DWS) ************************************************************************ * Calculates relative weights for the dependent variables at * * each measured point. * * Input variables: NY= number of dependent varibles * * M= index of the actual measured point * * DV(MN,M,I)= value of the ith dependent variable * * DW(MN,I,I)= actual weight of the ith dep.var. * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DV(NN,NMM,NYM),DW(NYM,NYM) DO 1600 I=1,NY DY1=DABS(DV(MN,M,I)) DO 1600 J=1,NY IF(DY1.LT.1.E-40)DY1=1.E-40 IF(I.EQ.J)THEN DW(I,J)=DWS/DY1/DY1 ELSE DW(I,J)=0. ENDIF 1600 CONTINUE RETURN END C SUBROUTINE UNCERW(NN,MN,NMM,NYM,M,NY,DN,DW,DWS) ************************************************************************ * Calculates relative weights for the dependent variables at * * each measured point. * * Input variables: NY= number of dependent varibles * * M= index of the actual measured point * * DN(MN,M,I)= uncert. value of the ith dep.var. * * DW(I,I)= actual weight of the ith dependent var. * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DN(NN,NMM,NYM),DW(NYM,NYM) DO 1600 I=1,NY DN1=DABS(DN(MN,M,I)) DO 1600 J=1,NY IF(DN1.LT.1.E-40)DN1=1.E-40 IF(I.EQ.J)THEN DW(I,J)=DWS/DN1/DN1 ELSE DW(I,J)=0. ENDIF 1600 CONTINUE RETURN END C SUBROUTINE POISSONW(NN,MN,NMM,NYM,M,NY,DV,DW,DWS) ************************************************************************ * Calculates POISSON weights for the dependent variables at * * each measured point. ( i. e. standard deviation = SQRT(DY(I)) ) * * Input variables: NY= number of dependent varibles * * M= index of the actual measured point * * DV(MN,M,I)= value of the ith dependent variable * * DW(MN,I,I)= actual weight of the ith dep.var. * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DV(NN,NMM,NYM),DW(NYM,NYM) DO 160 I=1,NY DY1=DABS(DV(MN,M,I)) DO 160 J=1,NY IF(DY1.LT.1.E-40)DY1=1.E-40 IF(I.EQ.J)THEN DW(I,J)=DWS/DY1 ELSE DW(I,J)=0. ENDIF 160 CONTINUE RETURN END C FUNCTION DSSQRS(NN,NMM,NXM,NYM,NM,NX,NY,NP,NS,MW,DT,DV,DW,DWS,DP, 1 DPP,IP,AC,IC) ************************************************************************ * Calculates the sum of weighted squared differences between * * measured and calculated values of the independent variables. * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DT(13,226,1),DV(13,226,1),DW(NYM,NYM),DP(15), 1 DPP(NS),IP(NS),DX(15),DY(15),NM(NN),NX(NN),NY(NN),AC(10) DSSQRS=0. DO 170 MN=1,NN DO 170 M=1,NM(MN) DO 171 I=1,NX(MN) 171 DX(I)=DT(MN,M,I) CALL MODEL(MN,NX(MN),NY(MN),NP,NS,DX,DY,DP,DPP,IP,AC,IC) IF(MW.EQ.1)CALL RELW(NN,MN,NMM,NYM,M,NY(MN),DV,DW,DWS) IF(MW.EQ.5)CALL UNCERW(NN,MN,NMM,NYM,M,NY(MN),DN,DW,DWS) IF(MW.EQ.4)CALL POISSONW(NN,MN,NMM,NYM,M,NY(MN),DV,DW,DWS) IF(MW.EQ.3)THEN IER=0 CALL USERW(NN,MN,NMM,NXM,NYM,M,NX(MN),NY(MN),NP,NS,DT,DV, 1 DW,DWS,DP,AC,IC) IF(IER.GT.0)RETURN ENDIF DO 172 I=1,NY(MN) DO 172 J=1,NY(MN) 172 DSSQRS=DSSQRS+DW(I,J)*(DV(MN,M,I)-DY(I))*(DV(MN,M,J)-DY(J)) 170 CONTINUE RETURN END C SUBROUTINE SCROLL(IP) INTEGER IP IP=IP+1 WRITE(*,600) 600 FORMAT(/) PAUSE'Press to SCROLL !' RETURN END