SUBROUTINE MMRQT(NAME,NM,NX,NY,NP,MW,MD,DP0,DT,DV,NSTEP,DSTEP,DW, 1 AC,IC) ************************************************************************ * * * Interactive MicroSoft Fortran 5.1 version -- * * creates PRT.OUT and FIT.OUT files * * * * 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-2). 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 * * NM= number of data points * * NX= number of independent variables * * NY= number of dependent variables * * NP= number of parameters * * DT(NM*NX) contains (measured) INDEPENDENT variables * * DV(NM*NY) contains (measured) DEPENDENT variables * * (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) * * 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) * * DP0(NP) contains the initial guess of the * * parameters. (Signs remain during iteration !!!) * * 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= 20) * * 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 the subroutine after iteration. * * PRINTOUT of the results is possible using Ctrl-PrtScreen option. * * An output file named MRQ.OUT (logical unit 99) will be created * * containing the identification and the estimated parameters in its * * heading and Xi's, calculated Yi's and calculated 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(NM,NX),DV(NM,NY),DW(NY,NY),DA(15,15), 1DC(15,15),DU(15,15),DG(15,15),DP(15),DP0(NP),DD(15),DBW(15), 2 CCY(1024,10),RES(1024,10),DC0(15),D0C(15),DX(15),DY(15),DB(15), 3 AC(10) CHARACTER NAME*40,CHR*1 IF(DSTEP.EQ.0.)DSTEP=1.D-9 IF(NSTEP.EQ.0)NSTEP=50 IF(MW.EQ.0)THEN DO 202 I=1,NY DO 202 J=1,NY 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 GOTO 203 2203 DSW=0. DWS=1. *** DWS (^) is for scaling rel. or Poisson weights. *************** DO 204 M=1,NM IF(MW.EQ.1)CALL RELW(NM,M,NY,DV,DW,DWS) IF(MW.EQ.4)CALL POISSONW(NM,M,NY,DV,DW,DWS) IF(MW.EQ.3)CALL USERW(M,NM,NX,NY,NP,DT,DV,DW,DWS,DP,AC,I) DO 205 I=1,NY 205 DSW=DSW+DW(I,I) 204 CONTINUE DWS=DBLE(NM*NY)/DSW 203 NEI=0 NEF=0 NES=0 DO 615 I=1,NP 615 DP(I)=DP0(I) DPM=.01 WRITE(*,701)NM 701 FORMAT(I10,' data points') 721 FORMAT( 5X,'X:',10F12.4) 722 FORMAT( 5X,'Y:',10F12.4) C *** Starting the iteration.********************* DF=DSSQRS(NM,NX,NY,NP,MW,DT,DV,DW,DWS,DP,AC,IC) IF(DF.GT.9.D9)WRITE(*,933)DF WRITE(*,300)DF 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(,I1,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 M=1,NM DO 211 I=1,NX 211 DX(I)=DT(M,I) IF(MW.EQ.1)CALL RELW(NM,M,NY,DV,DW,DWS) IF(MW.EQ.4)CALL POISSONW(NM,M,NY,DV,DW,DWS) IF(MW.EQ.3)THEN IER=0 CALL USERW(M,NM,NX,NY,NP,DT,DV,DW,DWS,DP,AC,I) IF(IER.GT.0)RETURN ENDIF C *** Calculating the Jacobian DGij = dYi/dPj *********** IF(MD.NE.0) GOTO 220 CALL NUMJACOBI(NX,NY,NP,DX,DY,DG,DP,AC,IC) GOTO 222 220 CALL DERJACOBI(NX,NY,NP,DX,DY,DG,DP,AC,IC) C *** Calculating second derivatives of chi-square *********** C 222 write(*,*)(dg(1,j),j=1,4) 222 DO 230 K=1,NP DO 231 L=1,K DAA=0. DO 232 I=1,NY DO 232 J=1,NY 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 DO 233 J=1,NY 233 DAA=DAA+DW(I,J)*DG(J,K)*(DV(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.3)GOTO 3203 GOTO 6203 3203 DSW=0. DWS=1. *** DWS (^) is for rescaling USER weights. *************** DO 6204 M=1,NM IF(MW.EQ.1)CALL RELW(NM,M,NY,DV,DW,DWS) IF(MW.EQ.4)CALL POISSONW(NM,M,NY,DV,DW,DWS) IF(MW.EQ.3)CALL USERW(M,NM,NX,NY,NP,DT,DV,DW,DWS,DP,AC,I) DO 6205 I=1,NY 6205 DSW=DSW+DW(I,I) 6204 CONTINUE DWS=DBLE(NM*NY)/DSW 6203 CONTINUE DF=DSSQRS(NM,NX,NY,NP,MW,DT,DV,DW,DWS,DP,AC,IC) IF(DF.GT.9.D9)WRITE(*,934)IT,DSL,DPM,DF WRITE(*,302)IT,DSL,DPM,DF 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=NM*NY-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 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='PRT.OUT') NNW=0 416 WRITE(*,400) 400 FORMAT(1H0,//////,15X,'Do you want to Display results on the scree 1n (D)',//,20X,'or Skip them and put into PRT.OUT (S) ?' 2,///,15X,'Press D & , or S & ! ',\) READ(*,401)CHR 401 FORMAT(A1) IF(CHR.EQ.'S')NNW=1 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,311)NP,NY,NX,NM,DSTEP WRITE(*,311)NP,NY,NX,NM,DSTEP 311 FORMAT( 15X,'Number of PARAMETERS.............',I2,/,15X, 1 'Number of FUNCTIONS..............',I2,/,15X, 2 'Number of independent VARIABLES..',I2,/,15X, 3 'Number of DATA POINTS............',I3,/,15X, 4 '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.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') 314 FORMAT( 17X,'User defined weights') WRITE(88,3312)(DP0(IKI),IKI=1,NP) WRITE(*,3312)(DP0(IKI),IKI=1,NP) 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('-')) WRITE(88,316) 316 FORMAT( 7X,\) DO 267 I=1,NX 267 WRITE(88,317)I 317 FORMAT( 4X,'X',I1,4X,\) DO 268 I=1,NY 268 WRITE(88,318)I 318 FORMAT( 4X,'Y',I1,4X,\) WRITE(88,319) 319 FORMAT( /,2X,75('-')) DO 270 IM=1,30 WRITE(88,320)IM 320 FORMAT( 2X,I3,1X,\) WRITE(88,321)(DT(IM,I),I=1,NX),(DV(IM,I),I=1,NY) 321 FORMAT( 0P10(F10.4)) IF(MOD(IM,5).EQ.0)WRITE(88,336) IF(IM.EQ.NM)GOTO 277 270 CONTINUE 271 WRITE(88,319) KM=31 272 WRITE(88,322) 322 FORMAT(1H1,//,2X,75('-')) WRITE(88,316) DO 273 I=1,NX 273 WRITE(88,317)I DO 274 I=1,NY 274 WRITE(88,318)I WRITE(88,319) DO 275 IM=KM,KM+49 WRITE(88,320)IM WRITE(88,321)(DT(IM,I),I=1,NX),(DV(IM,I),I=1,NY) IF(MOD(IM,5).EQ.0)WRITE(88,336) IF(IM.EQ.NM)GOTO 277 275 CONTINUE KM=KM+50 GOTO 272 277 WRITE(88,733) 733 FORMAT( 2X,75('-')) 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) 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(',I1,')',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(',I1,')',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)') IF(NNW.EQ.0)CALL SCROLL(IP) WRITE(88,332) 332 FORMAT(1H1,///,10X,'Fit of the model at the actual data points', 1//) DSS=0. DVS=0. DO 281 I=1,NY IF(NY.GT.1)WRITE(*,333)I 333 FORMAT(/,30X,I1,'. FUNCTION',//) C***** printout of the fit omitted - results in FIT.OUT ******* CC WRITE(*,334) CC334 FORMAT( 10X,'Y observed',8X,'Y calculated',9X,'DIFFERENCE',5X, CC 1 '% DIFFERENCE',/) CC WRITE(*,*)' ' DS2=0. DE1=0. DE2=0. DO 282 M=1,NM DO 283 J=1,NX 283 DX(J)=DT(M,J) CALL MODEL(NX,NY,NP,DX,DY,DP,AC,IC) CCY(M,I)=DY(I) RES(M,I)=DV(M,I)-DY(I) DS2=DS2+(DV(M,I)-DY(I))**2 DE1=DE1+DV(M,I) DE2=DE2+DV(M,I)*DV(M,I) IF(DV(M,I).EQ.0.)DV(M,I)=DV(M,I)+.0000000000001 CC WRITE(*,335)M,DV(M,I),DY(I),DV(M,I)-DY(I), CC 1 1.*(DV(M,I)-DY(I))/DV(M,I) CC IF(MOD(M,5).EQ.0)WRITE(*,336) CC IF(MOD(M,15).EQ.0.AND.NNW.EQ.0)CALL SCROLL(IP) CC IF(MOD(M,45).EQ.0)WRITE(*,337) 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))/DBLE(NM-1) DRV=DS2/DBLE(NM-NP-1) DYM=DE1/DBLE(NM) DGFIT=100*(1-DS2/DVV/DBLE(NM-1)) DSS=DSS+DS2 DVS=DVS+DVV*DBLE(NM-1) WRITE(88,338)DYM,DVV,DRV,DGFIT WRITE(*,338)DYM,DVV,DRV,DGFIT 338 FORMAT( /,20X,'MEAN=',1PE10.4,8X,'VARIANCE=',1PE10.3,//,16X, 1'RESIDUAL=',1PE10.3,13X,'FIT=',0PF7.3,'%',//) IF(NNW.EQ.0)CALL SCROLL(IP) CC WRITE(*,337) 281 CONTINUE C OPEN(99,FILE='FIT.OUT') 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 3267 I=1,NX 3267 WRITE(99,3317)I 3317 FORMAT( 10X,'X',I1,\) DO 3268 I=1,NY 3268 WRITE(99,3318)I 3318 FORMAT( 10X,'Y',I1,\) DO 3269 I=1,NY 3269 WRITE(99,3319)I 3319 FORMAT( 8X,'RES',I1,\) WRITE(99,401)' ' DO 3270 M=1,NM WRITE(99,3320)(DT(M,I),I=1,NX),(CCY(M,I),I=1,NY),(RES(M,I),I=1,NY) 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') 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(',I1,')',3X,\) WRITE(88,341) WRITE(*,341) 341 FORMAT( /) DO 286 I=1,NP WRITE(88,342)I WRITE(*,342)I 342 FORMAT(5X,'P(',I1,')',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 IF(NNW.EQ.0)GOTO 290 WRITE(88,344) WRITE(*,344) 344 FORMAT( /////,10X,'Covariance matrix of the parameters',///,8X,\) DO 287 I=1,NP WRITE(88,350)I 287 WRITE(*,350)I 350 FORMAT( 8X,'P(',I1,')',1X,\) WRITE(88,341) WRITE(*,341) DO 288 I=1,NP WRITE(88,352)I WRITE(*,352)I 352 FORMAT( 3X,'P(',I1,')',3X,\) DO 299 J=1,I DPP=DC(I,J)*DBW(I)*DBW(J) WRITE(88,353)DPP WRITE(*,353)DPP 353 FORMAT( 3X,1PE10.4,\) 299 CONTINUE WRITE(88,888) WRITE(*,888) 288 CONTINUE IF(NNW.EQ.1)RETURN 290 WRITE(*,402) 888 FORMAT( /) 402 FORMAT(1H0,//, 6X,'To Display the results AGAIN on the screen, pre 1ss D & !',//, 6X,'If you don`t want to see them,just pres 2s ! ',\) READ(*,401)CHR IF(CHR.EQ.' ')RETURN IF(CHR.EQ.'D')NNW=0 IF(CHR.EQ.'d')NNW=0 GOTO 777 END C SUBROUTINE NUMJACOBI(NX,NY,NP,DX,DY,DG,DP,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(NP),DY(NY),D(15),AC(10) DO 20 J=1,NP DE=.001*DABS(DP(J))+.000001 DP(J)=DP(J)+DE CALL MODEL(NX,NY,NP,DX,DY,DP,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(NX,NY,NP,DX,DY,DP,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 DN=DBLE(N) DV0=DV/DN/DN*.000005 DV1=0. 130 DV=DV/DN 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(NM,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(M,I)= value of the ith dependent variable * * DW(I,I)= actual weight of the ith dep.var. * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DV(NM,NY),DW(NY,NY) DO 160 I=1,NY DY1=DABS(DV(M,I)) IF(DY1.LT.1.E-40)DY1=1.E-40 DW(I,I)=DWS/DY1/DY1 160 CONTINUE RETURN END C SUBROUTINE POISSONW(NM,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(M,I)= value of the ith dependent variable * * DW(I,I)= actual weight of the ith dep.var. * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DV(NM,NY),DW(NY,NY) DO 160 I=1,NY DW(I,I)=DWS/DV(M,I) 160 CONTINUE RETURN END C FUNCTION DSSQRS(NM,NX,NY,NP,MW,DT,DV,DW,DWS,DP,AC,IC) ************************************************************************ * Calculates the sum of weighted squared differences between * * measured and calculated values of the independent variables. * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DIMENSION DT(NM,NX),DV(NM,NY),DW(NY,NY),DP(NP),DX(15),DY(15), 1 AC(10) DSSQRS=0. DO 170 M=1,NM DO 171 I=1,NX 171 DX(I)=DT(M,I) CALL MODEL(NX,NY,NP,DX,DY,DP,AC,IC) IF(MW.EQ.1) CALL RELW(NM,M,NY,DV,DW,DWS) IF(MW.EQ.4) CALL POISSONW(NM,M,NY,DV,DW,DWS) IF(MW.EQ.3) CALL USERW(M,NM,NX,NY,NP,DT,DV,DW,DWS,DP,AC,I) DO 172 I=1,NY DO 172 J=1,NY 172 DSSQRS=DSSQRS+DW(I,J)*(DV(M,I)-DY(I))*(DV(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