//U11701AA JOB (11701,123-45-6789),TIME=(1,0),CLASS=2, // MSGCLASS=Y,NOTIFY=* /*JOBPARM ROOM=H,FORMS=9021 // EXEC WATFIV,REGION.WATFIV=500K $JOB ,LIBLIST C C CRICFJC.CNTL CRICF WITH DSHUFF AND NEW GAUSF C C PROGRAM CRICF (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT) C C CRICF 2.3 A.N.S.I. STANDARD FORTRAN JUNE 1992 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C DOYLE E. HILL AND H. OLIN SPIVEY, DEPARTMENT OF BIOCHEMISTRY C OKLAHOMA STATE UNIVERSITY C STILLWATER,OKLAHOMA 74074 C C CRICF PERFORMS LEAST SQUARES FITS TO DATA OF MODELS WHICH C MAY BE THE SOLUTIONS OF SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS. C SEVERAL METHODS ARE AVAILABLE FOR EVALUATING THE GOODNESS-OF-FIT C AND THE ERRORS IN THE FITTED VALUES OF THE PARAMETERS. C C REFERENCE.... J. P. CHANDLER, D. E. HILL, AND H. O. SPIVEY, C COMPUTERS AND BIOMEDICAL RESEARCH 5 (1972) 515-534. C C CRICF IS DISTRIBUTED BY C H. O. SPIVEY, DEPARTMENT OF BIOCHEMISTRY, C OKLAHOMA STATE UNIVERSITY, STILLWATER, OKLAHOMA 74074 C C FOR INSTRUCTIONS ON USAGE OF CRICF, READ THE WRITE-UP. C C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THIS MAIN PROGRAM READS IN DATA AND CONTROLS SUPPORT PLANE C COMPUTATIONS AND MONTE CARLO RE-FITS. C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION FIT,FITSV DOUBLE PRECISION DSEEDA,DSEEDB,DSHUFF,GAUSF C C THE DIMENSIONS OF THE MATRICES ARE C Y(NPTS), YY(NPTS), YSIG(NPTS), FIT((PTS), FITSV(NPTS), X(NPTS), C B(NPAR), BMIN(NPAR), BMAX(NPAR), BERR(NPAR), BERSV(NPAR), C BSAVE(NPAR), ERALF(NPAR,NPAR), SIG(NPAR), LOGB(NPAR), MSK(NPAR), C SPAR(4), KTITL(80), C YDE(MAX(NEQFU,NEQDE)), BD(NPAR), CUSER(NCUSR), DFDB(NPAR) C DIMENSION Y(100),YY(100),YSIG(100),FIT(100),FITSV(100) DIMENSION B(30),BMIN(30),BMAX(30),BERR(30),BERSV(30),BSAVE(30) DIMENSION SIG(30) DIMENSION SPAR(4),KTITL(80) C DIMENSION ERALF(30,30) C COMMON /CRICK/ X(100),DFRAC,ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,KALC,NMAX,NPRNT,KW, * LOGB(30),MSK(30),NPAR,NFUPR COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30), * JSET,NCUSR COMMON /CDSHUF/ TABLE(128) C C THE STEP SIZE FOR INTEGRATING OVER THE FIRST INTERVAL, FOR EACH C JSET (EACH SET OF DATA), IS SAVED IN DHSAV. C COMMON /CDHSAV/ DHSAV(10) COMMON /MINK/ NSW C ZSQRT(ARG)=SQRT(ARG) ZLOGT(ARG)=ALOG(ARG)/2.302585 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INITIALIZE. C KR ... LOGICAL UNIT NUMBER OF THE READER KR=5 C KW ... LOGICAL UNIT NUMBER OF THE PRINTER KW=6 C MXPAR ... MAXIMUM NUMBER OF PARAMETERS, C AND THE DIMENSION OF B, BMIN, ETC. MXPAR=30 C MXPTS ... MAXIMUM NUMBER OF DATA POINTS, C AND THE DIMENSION OF Y, YY, YSIG, ETC. MXPTS=100 C MXNCU ... MAXIMUM NUMBER OF USER C CONSTANTS, AND THE DIMENSION OF CUSER MXNCU=60 C DFRAC ... RELATIVE STEP SIZE FOR DERIV DFRAC=.002 C MXSUP ... MAXIMUM NUMBER OF ITERATIONS IN C SUPPORT PLANE COMPUTATIONS MXSUP=3 C TOL ... TOLERANCE IN SUPPORT PLANE C COMPUTATIONS TOL=.1 C HUGE ... A VERY LARGE REAL NUMBER C (DEFAULT VALUE FOR XMAX AND -XMIN) HUGE=1.0E30 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C READ IN THE DATA. C READ(KR,1000)KTITL 1000 FORMAT(80A1) WRITE(KW,1010)KTITL 1010 FORMAT('1 CRICF PROBLEM......'///10X,80A1) READ(KR,1020)NSETS 1020 FORMAT(I5) C CHECK FOR DEFAULT. IF(NSETS)1030,1030,1040 1030 NSETS=1 1040 WRITE(KW,1050)NSETS 1050 FORMAT(/////1X,I5,' SET(S) OF DATA') NPTS=0 YSDEF=1.0 DO 1140 JSET=1,NSETS DHSAV(JSET)=0.0 READ(KR,1020)MPTS WRITE(KW,1060)MPTS 1060 FORMAT(/' SET OF',I4,' DATA POINTS....'/' ') IF(MPTS)1070,1070,1080 1070 STOP 1080 JPTS(JSET)=MPTS DO 1120 J=1,MPTS NPTS=NPTS+1 READ(KR,1090)X(NPTS),Y(NPTS),YS 1090 FORMAT(3E10.5) IF(YS)1110,1110,1100 1100 YSDEF=YS 1110 YS=YSDEF YSIG(NPTS)=YS 1120 WRITE(KW,1130)NPTS,X(NPTS),Y(NPTS),YS 1130 FORMAT(6X,'J =',I4,6X,'X =',1PE15.7,6X,'Y =',E15.7, * 6X,'YSIG =',E11.3) 1140 CONTINUE C C READ IN THE OTHER INPUT QUANTITIES, SET DEFAULT VALUES, AND PRINT. C READ(KR,1150)NPRNT,MAXIT,NMAX,NCUSR,ER,RER,(SPAR(J),J=1,4) 1150 FORMAT(4I5,6E10.5) IF(MAXIT)1160,1160,1170 1160 MAXIT=100 1170 IF(NMAX)1180,1180,1190 1180 NMAX=5000 1190 IF(ER)1220,1200,1210 1200 IF(RER)1220,1220,1230 1210 IF(RER)1220,1230,1230 1220 ER=0.0 RER=.0001 1230 IF(SPAR(1))1240,1240,1250 1240 SPAR(1)=10.0 1250 IF(SPAR(2))1260,1260,1270 1260 SPAR(2)=.1 1270 WRITE(KW,1280)NPRNT,MAXIT,NMAX,NCUSR,ER,RER,(SPAR(J),J=1,4) 1280 FORMAT(////'1'////' NPRNT =',I2,5X,'MAXIT =',I4,5X, * 'NMAX =',I6,5X,'NCUSR =',I4//' ER =',1PE13.5,5X, * 'RER =',E13.5,5X,'SPAR =',4E15.5) IF(NCUSR)1390,1330,1290 1290 IF(NCUSR-MXNCU)1300,1300,1390 1300 READ(KR,1310)(CUSER(J),J=1,NCUSR) 1310 FORMAT(8D10.5) WRITE(KW,1320)(J,CUSER(J),J=1,NCUSR) 1320 FORMAT(/////' USER CONSTANTS....'// * (6X,'CUSER(',I2,') =',1PG15.7)) 1330 READ(KR,1340)NPAR,NEQFU,NEQDE,KALC,MCSIM,KVSUP,SPFAC,NFUPR 1340 FORMAT(6I5,E10.5,I5) IF(KALC-2)1350,1360,1350 1350 KALC=1 1360 WRITE(KW,1370)NPAR,NEQFU,NEQDE,KALC,MCSIM,KVSUP,SPFAC,NFUPR 1370 FORMAT(////' NPAR =',I4,5X,'NEQFU =',I4,5X,'NEQDE =', * I4,5X,'KALC =',I2//5X,'MCSIM =',I5,5X,'KVSUP =',I3, * 5X,'SPFAC =',1PG13.5,5X,'NFUPR =',I4) IF(NPAR)1390,1390,1380 1380 IF(NPAR-MXPAR)1410,1410,1390 1390 WRITE(KW,1400) 1400 FORMAT(////' ONE OR MORE INPUT VALUES ARE ILLEGAL.'//' ') STOP C C THE QUANTITIES READ IN ABOVE ARE FIXED FOR EACH RUN, AT PRESENT. C NOW READ IN THE QUANTITIES FOR THE NEXT FIT IN THIS RUN. C 1410 CONTINUE READ(KR,1420)(MSK(J),J=1,NPAR) 1420 FORMAT(8I10) IF(MSK(1)-1)1422,1422,1421 1421 STOP 1422 WRITE(KW,1430)(MSK(J),J=1,NPAR) 1430 FORMAT(////' MSK(J).....',I9,4I12/(9X,5I12)) READ(KR,1420)(LOGB(J),J=1,NPAR) NLOGB=0 DO 1470 J=1,NPAR IF(MSK(J))1440,1450,1440 1440 LOGB(J)=0 1450 IF(LOGB(J))1460,1470,1460 1460 NLOGB=NLOGB+1 1470 CONTINUE WRITE(KW,1480)(LOGB(J),J=1,NPAR) 1480 FORMAT(/' LOGB(J)....',I9,8I12/(9X,5I12)) READ(KR,1490)(BSAVE(J),J=1,NPAR) 1490 FORMAT(8E10.5) READ(KR,1490)(BMIN(J),J=1,NPAR) READ(KR,1490)(BMAX(J),J=1,NPAR) C C DEFAULT BMIN AND BMAX. DO 1510 J=1,NPAR IF(BMAX(J)-BMIN(J))1500,1500,1510 1500 BMAX(J)=HUGE BMIN(J)=-HUGE 1510 CONTINUE WRITE(KW,1520)(BSAVE(J),J=1,NPAR) 1520 FORMAT(/' B(J)....... ',1PE12.4,4E12.4/ * (10X,5E12.4)) WRITE(KW,1530)(BMIN(J),J=1,NPAR) 1530 FORMAT(/' BMIN(J).... ',1PE12.4,4E12.4/ * (10X,5E12.4)) WRITE(KW,1540)(BMAX(J),J=1,NPAR) 1540 FORMAT(/' BMAX(J).... ',1PE12.4,4E12.4/ * (10X,5E12.4)) C C MAP B, BMAX, AND BMIN INTO LOGS (BASE 10) C IF REQUESTED. C NLOGB ... NUMBER OF LOGB(J) .NE. ZERO NLOGB=0 DO 1680 J=1,NPAR IF(LOGB(J))1550,1640,1550 1550 IF(BMIN(J))1560,1560,1570 1560 BMIN(J)=1.0/HUGE 1570 BMIN(J)=ZLOGT(BMIN(J)) 1580 IF(BMAX(J))1590,1590,1600 1590 BMAX(J)=HUGE 1600 BMAX(J)=ZLOGT(BMAX(J)) 1610 IF(BSAVE(J))1620,1620,1630 1620 BSAVE(J)=1.0 GO TO 1640 1630 BSAVE(J)=ZLOGT(BSAVE(J)) 1640 IF(BSAVE(J)-BMAX(J))1660,1660,1650 1650 BSAVE(J)=BMAX(J) 1660 IF(BSAVE(J)-BMIN(J))1670,1680,1680 1670 BSAVE(J)=BMIN(J) 1680 BD(J)=BSAVE(J) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C FIT THE DATA. C WRITE(KW,1690) 1690 FORMAT(////'1'///' FIT THE DATA..........'//' ') NFTYP=0 CALL SOLSUB (Y,YSIG,FITSV,BSAVE,BMIN,BMAX,BERSV,SPAR,CH, * ERSCL,NPTS,NFTYP,MAXIT) DO 1700 J=1,NPAR 1700 B(J)=BSAVE(J) IF(ERSCL-1.0)1740,1740,1710 1710 DO 1720 J=1,NPTS 1720 YSIG(J)=ERSCL*YSIG(J) WRITE(KW,1730)ERSCL 1730 FORMAT(/////' ALL YSIG(J) SCALED UP BY THE FACTOR', * ' ERSCL =',1PG13.5) C C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C DO A SUPPORT PLANE COMPUTATION IF KVSUP .GT. ZERO. C 1740 IF(KVSUP)2020,2020,1750 1750 NFTYP=1 C C LOOP OVER THE PARAMETERS FOR WHICH SUPPORT C PLANES ARE TO BE COMPUTED. DO 2010 KMASK=KVSUP,NPAR WRITE(KW,1760)KMASK,SPFAC 1760 FORMAT(////'1SUPPORT PLANE COMPUTATION FOR PARAMETER', * I4//' THE SUPPORT PLANES WILL LIE APPROXIMATELY', * ' SPFAC =',1PG13.5/5X, * 'STANDARD ERRORS FROM THE BEST FIT VALUES.') C C HOLD PARAMETER B(KMASK) FIXED. MSAV=MSK(KMASK) MSK(KMASK)=1 C SET B(KMASK) TO A FIRST GUESS. C DELTA=BERSV(KMASK)*SPFAC BD(KMASK)=BSAVE(KMASK)+DELTA C C RE-FIT, HOLDING B(KMASK) FIXED. C LOOP OVER THE TWO HALF-INTERVALS. DO 1990 JJ=1,2 C ITERATE TO FIND THIS SUPPORT PLANE. DO 1900 JSUP=1,MXSUP B(KMASK)=BD(KMASK) BHOLD=BD(KMASK) IF(LOGB(KMASK))1770,1790,1770 1770 IF(MSAV)1790,1780,1790 1780 BD(KMASK)=10.0**BHOLD C C FIND THE MINIMUM OF PH IN THIS PLANE. C 1790 CALL SOLSUB (Y,YSIG,FIT,B,BMIN,BMAX,BERR,SPAR,PH, * ERSCL,NPTS,NFTYP,MAXIT) BD(KMASK)=BHOLD DIFFR=PH-CH IF(DIFFR)1800,1820,1830 1800 WRITE(KW,1810)DIFFR 1810 FORMAT(////' THE SUM OF SQUARES IS LESS THAN AT', * 'THE SUPPOSED "GLOBAL MINIMUM". DIFFR =', * 1PG13.5//' ') 1820 FACTR=2. GO TO 1870 1830 DD=ZSQRT(DIFFR) FACTR=SPFAC/DD BSP=BD(KMASK) ERB=FACTR*(BSP-BSAVE(KMASK)) WRITE(KW,1840)KMASK,BSP,DIFFR,FACTR,ERB 1840 FORMAT(//' SUPPORT PLANE FOR PARAMETER NUMBER', * I3//6X,'BD(KMASK) =',1PG15.7,5X,'DIFFER =', * G15.7,5X,'FACTR =',G13.5,6X,'ERB =',G13.5) C C TEST FOR CONVERGENCE. C IF(ABS(DD-SPFAC) .LE. TOL) GO TO DIFF=DD-SPFAC IF(DIFF)1850,1860,1860 1850 DIFF=-DIFF 1860 IF(DIFF-TOL)1920,1920,1870 C C MOVE THE PLANE TO A NEW POINT AND C TRY AGAIN. 1870 DO 1890 J=1,NPAR IF(MSK(J))1890,1880,1890 1880 B(J)=B(J)+(FACTR-1.0)*(B(J)-BSAVE(J)) 1890 CONTINUE BSP=BD(KMASK) BD(KMASK)=BSP+(FACTR-1.0)*(BSP-BSAVE(KMASK)) 1900 CONTINUE C WRITE(KW,1910)KMASK,MXSUP 1910 FORMAT(///' THIS SUPPORT PLANE COMPUTATION FOR', * ' PARAMETER',I4/5X,'FAILED TO CONVERGE IN MXSUP =', * I4,' ITERATIONS') GO TO 1960 C CONVERGED. PRINT THE RESULTS. 1920 BSP=BSAVE(KMASK) BFID=BSP+ERB WRITE(KW,1930)KMASK,ERB,BFID 1930 FORMAT(///' RESULTS OF THE SUPPORT PLANE COMPUTATION', * ' FOR PARAMETER',I4//6X, * 'LENGTH OF THE FIDUCIAL HALF-INTERVAL =',1PG13.5// * 6X,'ENDPOINT OF THE HALF-INTERVAL =',G15.7) IF(LOGB(KMASK))1940,1960,1940 1940 BFID=10.0**BFID WRITE(KW,1950)BFID 1950 FORMAT(//6X,'EXPONENTIATED ENDPOINT =',1PG15.7) C C REFLECT THROUGH THE GLOBAL MINIMUM AND C FIND THE PLANE ON THE OTHER SIDE. 1960 DO 1980 J=1,NPAR IF(MSK(J))1980,1970,1980 1970 B(J)=BSAVE(J)-(B(J)-BSAVE(J)) 1980 CONTINUE BSP=BD(KMASK) 1990 BD(KMASK)=BSAVE(KMASK)-(BSP-BSAVE(KMASK)) C C RESET THE PARAMETERS AND THE MASK. DO 2000 J=1,NPAR 2000 B(J)=BSAVE(J) 2010 MSK(KMASK)=MSAV C C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C PERFORM MONTE CARLO RE-FITS OF SIMULATED DATA, IF MCSIM .GT. ZERO. C 2020 IF(MCSIM)1410,1410,2030 2030 NFTYP=2 DSEEDA=123456789.0D0 DSEEDB=DSEEDA INITDS=1 WRITE(KW,2040) 2040 FORMAT(///'1 MONTE CARLO RE-FITS OF SIMULATED SETS', * ' OF DATA....') C C NW ... NUMBER OF RE-FITS WITH A VALUE OF C CHI-SQUARE LARGER THAN -CH- NW=0 C INITIALIZE THE ERROR MATRIX. DO 2050 J=1,NPAR DO 2050 JH=1,NPAR 2050 ERALF(J,JH)=0.0 C LOOP OVER THE NUMBER OF SIMULATIONS. DO 2130 KI=1,MCSIM WRITE(KW,2060)KI 2060 FORMAT(//////' THE FOLLOWING FIT WILL BE MONTE CARLO', * ' SIMULATION NUMBER',I5) C C LOOP OVER THE DATA POINTS. DO 2080 J=1,NPTS YY(J)=0.0 IF(FITSV(J)-Y(J) .EQ. 0.0) GO TO 2080 C C COMPUTE THE RANDOM -DATA POINT-. C 2070 YY(J)=FITSV(J)+YSIG(J)*GAUSF(DSEEDA,DSEEDB,INITDS) C C REJECT NEGATIVE DATA POINTS. C (THE USER MAY WISH TO CHANGE THIS....) C IF(YY(J))2070,2070,2080 2080 CONTINUE C RESET THE PARAMETERS. DO 2090 J=1,NPAR 2090 B(J)=BSAVE(J) C FIT THE SIMULATED DATA. C CALL SOLSUB (YY,YSIG,FIT,B,BMIN,BMAX,BERR,SPAR,PH, * ERSCL,NPTS,NFTYP,MAXIT) C IF(PH-CH)2110,2110,2100 2100 NW=NW+1 C ACCUMULATE THE SUMS FOR THE ERROR MATRIX. 2110 DO 2120 K=1,NPAR DO 2120 L=1,NPAR 2120 ERALF(K,L)=ERALF(K,L)+(B(K)-BSAVE(K))*(B(L)-BSAVE(L)) C 2130 CONTINUE C END OF THE SIMULATIONS. C COMPUTE THE ERROR MATRIX. ANES=MCSIM DO 2140 K=1,NPAR DO 2140 L=1,NPAR 2140 ERALF(K,L)=ERALF(K,L)/ANES C C COMPUTE THE STANDARD DEVIATIONS. DO 2150 K=1,NPAR 2150 SIG(K)=ZSQRT(ERALF(K,K)) C C COMPUTE THE CORRELATIONS. DO 2160 K=1,NPAR DO 2160 L=1,NPAR 2160 ERALF(K,L)=ERALF(K,L)/(SIG(K)*SIG(L)) C C COMPUTE THE CONFIDENCE LEVEL. CL=NW CL=CL/ANES C PRINT ALL RESULTS. WRITE(KW,2170) MCSIM 2170 FORMAT('1 RESULTS OF MONTE CARLO SIMULATIONS......'///// * ' NUMBER OF SIMULATIONS =',I5) WRITE(KW,2180)(SIG(K),K=1,NPAR) 2180 FORMAT(///' STANDARD ERRORS OF THE PARAMETERS....'// * (1X,1PE13.5,4E13.5)) WRITE(KW,2190) 2190 FORMAT(////' ENDPOINTS OF THE MONTE CARLO', * ' CONFIDENCE INTERVALS....'//' ') DO 2200 J=1,NPAR BFA=BSAVE(J)-SIG(J) BFB=BSAVE(J)+SIG(J) 2200 WRITE(KW,2210)J,BFA,BFB 2210 FORMAT(10X,I10,2E20.7) IF(NLOGB)2260,2260,2220 2220 WRITE(KW,2230) 2230 FORMAT(////' EXPONENTIATED ENDPOINTS.....'/' ') DO 2250 J=1,NPAR IF(LOGB(J))2240,2250,2240 2240 BFB=10.0**(BSAVE(J)+SIG(J)) BFA=10.0**(BSAVE(J)-SIG(J)) WRITE(KW,2210)J,BFA,BFB 2250 CONTINUE 2260 WRITE(KW,2270)CL 2270 FORMAT(///' APPROXIMATE CONFIDENCE LEVEL =',1PG15.7///// * ' CORRELATION MATRIX....') DO 2280 LN=1,NPAR 2280 WRITE(KW,2290)(ERALF(LN,LH),LH=1,NPAR) 2290 FORMAT(//(1X,8E13.5)) GO TO 1410 C C END CRICF MAIN PROGRAM. C END SUBROUTINE SOLSUB (YY,YSIG,FITSL,B,BMIN,BMAX,BERR,SPAR,PH, * ERSCL,NPTS,NFTYP,MAXIT) C C SOLSUB SUPERVISES ONE FIT (SOLSUB CALLS SOLVE), COMPUTES THE C VARIOUS STATISTICS, AND PRINTS RESULTS. C C NFTYP=0 TO FIT THE DATA, =1 TO COMPUTE SUPPORT PLANES, =2 FOR C MONTE CARLO RE-FITS OF SIMULATED DATA. C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION FITSL C C THE DIMENSIONS OF MATRICES ARE.... C YY(NPTS), YSIG(NPTS), FITSL(NPTS), B(NPAR), BMIN(NPAR), C BMAX(NPAR), BERR(NPAR), SPAR(4), C BSOLV(NACTV), BMNSL(NACTV), BMXSL(NACTV), BERSL(NACTV), C A(NACTV*(NACTV+1)/2), C(NACTV,5), ROOT(NACTV), C VECT(NACTV,NACTV), RESID(NPTS) C C NACTV IS THE NUMBER OF NON-MASKED B(J). C DIMENSION YY(1),YSIG(1),FITSL(1),B(1),BMIN(1),BMAX(1),BERR(1), * SPAR(1) DIMENSION BSOLV(30),BMNSL(30),BMXSL(30),BERSL(30) DIMENSION A(465),C(30,5),ROOT(30),VECT(30,30) DIMENSION RESID(100) C COMMON /CRICK/ X(100),DFRAC,ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,KALC,NMAX,NPRNT,KW, * LOGB(30),MSK(30),NPAR,NFUPR COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C ZSQRT(ARG)=SQRT(ARG) C C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C NJX ... THE FIRST DIMENSION OF THE ARRAYS C C, ROOT, AND VECT (NJX.GE.NACTV) NJX=30 C CLTHA, CLTHB ... THRESHOLDS FOR -POOR- C AND -VERY POOR- CONFIDENCE LEVELS CLTHA=.01 CLTHB=.001 RTEN=10.0 C C PACK THE PARAMETERS. SOLVE SEES ONLY PACKED PARAMETERS (TO SAVE C STORAGE IF MANY MSK(J) ARE NONZERO). THE PARAMETERS ARE UNPACKED C AND EXPONENTIATED (IF NECESSARY) IN FOFX. DIFFUN SEES ONLY THE C ORIGINAL UNPACKED, EXPONENTIATED PARAMETERS. C KK=0 DO 3010 J=1,NPAR IF(MSK(J))3010,3000,3010 3000 KK=KK+1 BSOLV(KK)=B(J) BMNSL(KK)=BMIN(J) BMXSL(KK)=BMAX(J) 3010 CONTINUE NACTV=KK C WRITE(KW,3020)NPAR,NACTV,(J,B(J),BMIN(J),BMAX(J),MSK(J),LOGB(J), * J=1,NPAR) 3020 FORMAT(//////' SUBROUTINE SOLSUB. INPUT QUANTITIES....'// * 5X,I3,' PARAMETERS,',I3,' ACTIVE'//' PARAMETER',7X, * 'STARTING',9X,'LOWER',11X,'UPPER',14X,'MSK',7X, * 'LOGB'/' NUMBER',10X,'VALUE',2(11X,'LIMIT')// * (1X,I6,7X,3E16.7,2I10)) IF(NACTV)3040,3060,3030 3030 IF(NACTV-NJX)3090,3090,3040 3040 WRITE(KW,3050)NACTV,NJX 3050 FORMAT(////' ILLEGAL VALUE OF NACTV.',5X,'NACTV =',I7, * 5X,'NJX =',I7////' ') STOP C C THERE ARE NO ACTIVE (UNMASKED) PARAMETERS. CALL FUNC TO DO C THIS ZERO-PARAMETER -FIT-. C 3060 CALL FUNC (NACTV,BSOLV,NPTS,FITSL) PH=0.0 DO 3070 J=1,NPTS 3070 PH=PH+((YY(J)-FITSL(J))/YSIG(J))**2 DO 3080 J=1,NACTV 3080 BERSL(J)=0.0 GO TO 3230 C C CALL SOLVE TO DO THE LEAST SQUARES FIT. C USE THE ARRAY VECT FOR THE CORRELATION MATRIX, TO SAVE STORAGE. C 3090 ITER=0 3100 CALL SOLVE (NACTV,NPTS,ITER,ICON,KERR,BSOLV,BMNSL,BMXSL,YY,FITSL, * PH,SPAR,VECT,BERSL,KW,MAXIT,YSIG,NPRNT,NJX) IF(KERR)3110,3130,3110 3110 WRITE(KW,3120) 3120 FORMAT(///' SINGULAR MATRIX.... RANK = 0 .'//' ') GO TO 3230 3130 IF(ICON)3150,3230,3140 3140 IF(ITER-MAXIT)3100,3100,3210 3150 IF(ICON-(-2))3180,3180,3160 3160 WRITE(KW,3170) 3170 FORMAT(////' THE NUMBER OF PARAMETERS TO BE FITTED', * ' EXCEEDS THE NUMBER OF DATA POINTS.') GO TO 3200 3180 WRITE(KW,3190) 3190 FORMAT(////' THE NUMBER OF ACTIVE (NON-MASKED)', * ' PARAMETERS IS .LE. ZERO.'/////' ') 3200 STOP 3210 WRITE(KW,3220)MAXIT 3220 FORMAT(///' MORE THAN MAXIT =',I6,' ITERATIONS REQUIRED.') C C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE FIT IS COMPLETED, FOR BETTER OR FOR WORSE. C PRINT THE RESULTS. C 3230 WRITE(KW,3240)PH 3240 FORMAT(///' SOLSUB FINAL RESULTS..........'// * 6X,'FINAL SUB OF SQUARES =',1PG15.7/' ') IF(NFTYP)3250,3250,3340 3250 DO 3260 J=1,NACTV 3260 WRITE(KW,3270)J,BSOLV(J),BERSL(J) 3270 FORMAT(' J =',I4,5X,'BSOLV(J) =',1PE15.7,5X, * 'ERROR(J) =',E12.5) C C COMPUTE THE RESIDUALS AND NPFIT. NPFIT=0 DO 3290 J=1,NPTS RESID(J)=YY(J)-FITSL(J) IF(RESID(J))3280,3290,3280 3280 NPFIT=NPFIT+1 3290 CONTINUE NDF=NPFIT-NACTV DF=NDF C SCALE UP THE ERRORS IF CHI-SQUARE IS C GREATER THAN THE EXPECTED VALUE. ERSCL=1.0 IF(PH-DF)3340,3340,3300 3300 ERSCL=ZSQRT(PH/DF) DO 3310 J=1,NACTV 3310 BERSL(J)=BERSL(J)*ERSCL 3320 WRITE(KW,3330)(BERSL(J),J=1,NACTV) 3330 FORMAT(//' ERRORS SCALED UP BY THE FACTOR', * ' SQRT(CHI-SQUARE/NO.D.F.)....'//(5X,1PE18.6,4E18.6)) C C UNPACK THE PARAMETERS AND PRINT. 3340 WRITE(KW,3350) 3350 FORMAT(/////' UNPACKED, EXPONENTIATED RESULTS FROM', * ' SOLSUB....'///83X,'ONE-STD.-ERROR'/8X,'PARAMETER', * 7X,'MSK',4X,'LOGB',10X,'BMIN',11X,'BMAX',22X, * 'ENDPOINTS'/' ') J=0 DO 3410 K=1,NPAR IF(MSK(K))3360,3380,3360 3360 WRITE(KW,3370)B(K),MSK(K),LOGB(K) 3370 FORMAT(4X,1PE15.7,2I7,5X,2E15.5,5X,2E15.5) GO TO 3410 3380 J=J+1 B(K)=BSOLV(J) BERR(K)=BERSL(J) BB=B(K) BBPL=BB+BERR(K) BBMI=BB-BERR(K) BBMIN=BMIN(K) BBMAX=BMAX(K) IF(LOGB(K))3390,3400,3390 3390 BB=RTEN**BB BBMAX=RTEN**BBMAX BBMIN=RTEN**BBMIN BBPL=RTEN**BBPL BBMI=RTEN**BBMI 3400 WRITE(KW,3370)BB,MSK(K),LOGB(K),BBMIN,BBMAX,BBMI,BBPL 3410 CONTINUE C PRINT THE DATA, FIT, AND RESIDUALS. IF(NFTYP)3415,3415,3430 3415 WRITE(KW,3420) 3420 FORMAT('1',4X,'ABSCISSA',8X,'DATA',10X,'FIT',9X, * 'RESIDUAL',6X,'PER CENT',6X,'RESIDUAL/YSIG(J)'/ * 19X,'ORDINATE',6X,'ORDINATE',20X,'RESIDUAL') 3430 SUMP=0.0 SUMD=0.0 I=0 DO 3540 JSET=1,NSETS WRITE(KW,3440) 3440 FORMAT(' ') ITOP=JPTS(JSET) DO 3530 IPT=1,ITOP I=I+1 DIF=YY(I)-FITSL(I) RDSIG=DIF/YSIG(I) SUMD=SUMD+DIF*DIF IF(YY(I))3460,3450,3460 3450 PCT=0.0 IF(FITSL(I))3470,3480,3470 3460 PCT=100.0*DIF/YY(I) GO TO 3480 3470 PCT=100.0 3480 IF(PCT)3490,3500,3500 3490 PCT=-PCT 3500 SUMP=SUMP+PCT IF(NFTYP)3510,3510,3530 3510 WRITE(KW,3520)X(I),YY(I),FITSL(I),DIF,PCT,RDSIG 3520 FORMAT(1X,1PE14.5,4E14.5,3X,E14.5) 3530 CONTINUE 3540 CONTINUE C XN=NPTS SEE=ZSQRT(SUMD/XN) AVP=SUMP/XN IF(NFTYP)3550,3550,3630 3550 WRITE(KW,3560)SEE,AVP 3560 FORMAT(///' STANDARD DEVIATION OF FIT =',1PG15.7/// * ' ABSOLUTE AVERAGE PERCENT DEVIATION =',G15.7) C C PRINT THE CORRELATION MATRIX. WRITE(KW,3570) 3570 FORMAT(////' CORRELATION MATRIX....'/' ') DO 3580 J=1,NACTV 3580 WRITE(KW,3590)(VECT(J,JJ),JJ=1,J) 3590 FORMAT(/(1X,1PG13.4,4G13.4)) C COMPUTE THE EIGENVALUES AND EIGENVECTORS C OF THE CORRELATION MATRIX. KKK=1 DO 3600 I=1,NACTV ROOT(I)=0.0 DO 3600 J=1,I A(KKK)=VECT(I,J) 3600 KKK=KKK+1 CALL GIVENS (NACTV,NACTV,NJX,A,C,ROOT,VECT) WRITE(KW,3610)(ROOT(J),J=1,NACTV) 3610 FORMAT(////' EIGENVALUES OF THE CORRELATION MATRIX....'// * (/5X,1PE14.5,4E14.5)) C C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C COMPUTE AND PRINT SOME STATISTICS. C WRITE(KW,3620) 3620 FORMAT('1 STATISTICS....'////' ') C 3630 SECHI=ZSQRT(2.*DF) SPH=PH CHIPR=CHSPR(SPH,NDF) SIGFT=PH/DF SIGFT=ZSQRT(SIGFT) IF(NFTYP)3640,3640,3890 3640 WRITE(KW,3650)PH,NDF,DF,SECHI,CHIPR 3650 FORMAT(///6X,'CHISQUARE =',1PG15.7//6X, * 'NUMBER OF DEGRERES OF FREEDOM =',I4//6X, * 'EXPECTED VALUE OF CHI=SQUARE =',G12.5//6X, * 'STANDARD ERROR OF THE EXPECTED VALUE OF CHI-SQUARE =', * ' SQRT(2*N.D.F.) =',G13.5//6X, * 'CHI-SQUARE PROBABILITY =',G15.7) WRITE(KW,3660)SIGFT 3660 FORMAT(//6X,'SQRT(CHI-SQUARE PER D.F.) =',1PG13.5) C C COMPUTE THE SIGNS AND RUNS CONFIDENCE C LEVELS. WRITE(KW,3670) 3670 FORMAT(////' SIGNS AND RUNS TESTS .....') CLSUB=1.0 NA=0 NB=0 NRA=0 NRB=0 NZERO=0 JH=0 DO 3750 JSET=1,NSETS JL=JH+1 JH=JL+JPTS(JSET)-1 CALL NRUNS (RESID,JL,JH,NMI,NZER,NPL,NRMI,NRPL) PROB1=RUNPR(NPL,NMI,NRPL,NRMI) PROB2=SIGNP(NPL,NMI) IF(PROB1-CLSUB)3680,3690,3690 3680 CLSUB=PROB1 3690 IF(PROB2-CLSUB)3700,3710,3710 3700 CLSUB=PROB2 3710 WRITE(KW,3720)JSET 3720 FORMAT(//6X,'FOR DATA SET NUMBER',I3,',') WRITE(KW,3730)NPL,NMI,NZER,PROB2 3730 FORMAT(//11X,'NUMBER OF POSITIVE RESIDUALS =',I8// * 11X,'NUMBER OF NEGATIVE RESIDUALS =',I8// * 11X,'NUMBER OF ZERO RESIDUALS =',I12// * 11X,'SIGNS CONFIDENCE LEVEL =',1PG17.5) WRITE(KW,3740)NRPL,NRMI,PROB1 3740 FORMAT(//11X,'NUMBER OF RUNS OF POSITIVE RESIDUALS =', * I4//11X,'NUMBER OF RUNS OF NEGATIVE RESIDUALS =',I4// * 11X,'RUNS CONFIDENCE LEVEL =',1PG18.5//' ') 3750 CONTINUE C PRINT A MESSAGE IF THE FIT IS POOR. C IF(CLSUB-SPH)3760,3770,3770 3760 SPH=CLSUB 3770 IF(SPH-CLTHB)3780,3780,3820 3780 WRITE(KW,3790) 3790 FORMAT('1'//////' IN AT LEAST ONE RESPECT THE QUALITY', * ' OF THE FIT IS TERRIBLE.'/// * ' PERHAPS MOTHER NATURE IS TRYING TO TELL YOU', * ' SOMETHING........'//////' ') WRITE(KW,3800) 3800 FORMAT(/22X,'*',39X,'*' * /12X,2('*',9X),'*',19X,'*',9X,'*' * /' *',9X,'*',8X,'****',7X,2('*',9X),'*',8X,'****' * /' *',8X,2('****',6X),'****',7X,'*',8X,'****', * 6X,'****'/' **** ****',16X,2('****',6X),'****' * /' ****',36X,'****') WRITE(KW,3810) 3810 FORMAT(///////////////) GO TO 3850 3820 IF(SPH-CLTHA)3830,3830,3850 3830 WRITE(KW,3840) 3840 FORMAT(////' IN AT LEAST ONE RESPECT THE QUALITY OF THE', * ' FIT IS POOR.') C C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C RE-SCALE YSIG(J) AND PH IF NECESSARY. C 3850 IF(ERSCL-1.0)3890,3890,3860 3860 WRITE(KW,3870)ERSCL,DF 3870 FORMAT(////' THE YSIG(J) ARE BEING SCALED UP BY THE', * ' FACTOR',1PG13.5// * ' CHI-SQUARE IS BEING RE-SCALED TO ITS EXPECTED', * ' VALUE,',G15.7, * ' IN ORDER TO MAKE ANY ERRORS COMPUTED BELOW LARGER.' * //' ') DO 3880 J=1,NPTS 3880 YSIG(J)=YSIG(J)*ERSCL PH=DF C CALL FUNC TO SET THINGS TO FINAL VALUES. C 3890 NPRSV=NPRNT IF(NFTYP)3900,3900,3930 3900 IF(NFUPR)3930,3930,3910 3910 NPRNT=2 WRITE(KW,3920) 3920 FORMAT('1') 3930 CALL FUNC (NACTV,BSOLV,NPTS,FITSL) NPRNT=NPRSV C RETURN C END SOLSUB. END SUBROUTINE SOLVE (NACTV,NPTS,ITER,ICON,KERR,B,BMIN,BMAX,Y,FIT, * PH,SPAR,QSAVE,BERR,KW,MAXIT,YSIG,NPRNT,LQSAV) C C SOLVE 2.4 A.N.S.I. STANDARD FORTRAN AUGUST 1973 C C SOLVE PERFORMS A NONLINEAR LEAST SQUARES FIT OF A GIVEN FUNCTION C TO A GIVEN SET OF DATA, USING MARQUARDT-S METHOD. C REFERENCE.... D. W. MARQUARDT, J. S.I.A.M. JUNE 1963, P. 431 C C NACTV = THE NUMBER OF PARAMETERS, B(J), TO BE DETERMINED C NPTS = THE NUMBER OF OBSERVATIONS OF Y AVAILABLE C ITER = ITERATION COUNT - INITIALLY SHOULD BE ZERO C ICON = THE NUMBER OF B VALUES NOT SATISFYING CONVERGENCE REQUIREMENT C C SOLVE RETURNS TO THE CALLING PROGRAM AFTER EACH ITERATION. C ON RETURN TO THE MAIN PROGRAM C ICON = 0 INDICATES FINAL SOLUTION OF PROBLEM C ICON = -1 INDICATES MORE UNKNOWNS THAN FUNCTIONS C ICON = -2 INDICATES ILLEGAL VALUE OF NACTV OR NPTS C INTEGER NSUB C C KERR IS AN ERROR FLAG WHICH INDICATES THAT THE MATRIX OF NORMAL C EQUATIONS IS SINGULAR. C C B = THE VECTOR OF PARAMETERS C BMIN(I) = MINIMUM VALUE FOR B(I) TO BE ALLOWED (INPUT) C BMAX(I) = MAXIMUM VALUE FOR B(I) TO BE ALLOWED (INPUT) C (THE PARAMETERS STAY WITHIN THESE BOUNDS (EXCEPT FOR A C POSSIBLE SMALL VIOLATION IN SUBROUTINE DERIV), BUT IF ONE C OF THE BOUNDS IS REACHED, THE BEST CONSTRAINED FIT MAY NOT C BE FOUND.) C Y = VECTOR OF THE NPTS DATA ORDINATES (OBSERVATIONS) C FIT = VECTOR OF THE NPTS COMPUTED FITTED VALUES C PH = VALUE OF FUNCTION TO BE MINIMIZED C (THIS IS THE SUM OF ((Y(J)-FIT(J))/YSIG(J))**2) C C SPAR IS A VECTOR WHICH SUPPLIES THE SUBROUTINE WITH THE C FOLLOWING FOUR PARAMETERS C C FNU = THE NU FACTOR GIVEN IN THE SOURCE PAPER C FLA = THE LAMBDA FACTOR CITED IN SOURCE PAPER C EPS1 = CONVERGENCE CRITERIA FOR THE RESIDUAL SUM OF SQUARES C EPS2 = CONVERGENCE CRITERIA FOR THE PARAMETERS C DOUBLE PRECISION FIT,FTEMP C C ON COMPLETION, QSAVE CONTAINS THE CORRELATION MATRIX. C C Q AND QSAVE ARE SYMMETRIC MATRICES, SO ONLY ONE TRIANGLE OF EACH C IS REALLY NEEDED. THE FULL MATRICES ARE USED FOR CLARITY, AND C BECAUSE MATNV IS USED (DUE TO TROUBLE WITH THE ORIGINAL MATRIX C ROUTINES). C C DIMENSIONS OF MATRICES ARE C B(NACTV),BMIN(NACTV),BMAX(NACTV),Y(NPTS),FIT(NPTS),SPAR(4), C QSAVE(NACTV,NACTV),BERR(NACTV),YSIG(NPTS),FTEMP(NPTS), C X(NACTV),GRAD(NACTV),XSAVE(NACTV),P(NPTS,NACTV),Q(NACTV,NACTV) C C NOTE.... P AND Q ARE PRESENTLY DIMENSIONED FOR ONLY TEN (10) C ACTIVE PARAMETERS. C DIMENSION B(1),BMIN(1),BMAX(1),BERR(1),Y(1),YSIG(1),FIT(1), * SPAR(4),QSAVE(LQSAV,NACTV) DIMENSION BTEMP(30),X(30),GRAD(30) DIMENSION FTEMP(100),P(100,30),Q(30,30) C ZSQRT(ARG)=SQRT(ARG) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C LDIMP IS THE DIMENSION OF THE ARRAY FTEMP, AND THE FIRST C DIMENSION OF THE ARRAY P. C THE ARRAY P IS ORDINARILY THE LARGEST SINGLE ARRAY. C LDIMP=100 C C LDIMQ IS THE SECOND DIMENSION OF THE ARRAY P, AND IT IS THE FIRST C AND SECOND DIMENSIONS OF THE ARRAY Q. C LDIMQ=30 C C MXSUB ... MAXIMUM NUMBER OF SUBITERATIONS MXSUB=8 C IF NFLAT=1, CONVERGENCE IS ACHIEVED IF THE C TRIAL SUM OF SQUARES IS EQUAL TO THE C VALUE FROM THE PREVIOUS ITERATION. NFLAT=1 C IF NFLSB=1, CONVERGENCE IS ACHIEVED IF THE C TRIAL SUM OF SQUARES IS EQUAL TO THE C VALUE FROM THE PREVIOUS SUBITERATION. NFLSB=1 NFLSB=0 C FCUT ... FACTOR FOR CUTTING STEP SIZE FCUT=2. C FMULT ... FACTOR BY WHICH FCUT IS C INCREASED AT EACH SUBITERATION FMULT=1.7 C CCRIT ... CRITICAL VALUE OF THE COSINE C OF MARQUARDT-S ANGLE GAMMA CCRIT=0.866 C NSUB ... SUBITERATION COUNTER NSUB=0 KERR=0 IF(NACTV)5040,5040,5010 5010 IF(NACTV-LDIMQ)5020,5020,5040 5020 IF(NPTS)5040,5040,5030 5030 IF(NPTS-LDIMP)5050,5050,5040 5040 ICON=-2 GO TO 5160 5050 IF(ITER)5090,5090,5060 5060 IF(ITER-MAXIT)5170,5070,5070 5070 WRITE(KW,5080)MAXIT 5080 FORMAT(//////' MORE THAN',I5, * ' ITERATIONS ATTEMPTED IN SOLVE.') ITER=ITER+1 GO TO 6000 C C THIS IS THE FIRST CALL TO SOLVE IN THIS FIT. C INITIALIZE. C 5090 FNU=SPAR(1) FLA=SPAR(2) IF(SPAR(3))5100,5110,5110 5100 SPAR(3)=0.00001 5110 EPS2=SPAR(3) IF(SPAR(4))5120,5130,5130 5120 SPAR(4)=0.00001 5130 EPS1=SPAR(4) L=1 DO 5140 J=1,NPTS FIT(J)=0.0 5140 FTEMP(J)=0.0 C C THE INITIALIZATION IS COMPLETED. C IF(NPTS-NACTV)5150,5180,5180 5150 ICON=-1 5160 ITER=ITER+1 RETURN C THIS IS NOT THE FIRST CALL TO SOLVE. 5170 L=2 C 5180 DO 5870 I1=L,2 GO TO (5190,5280),I1 C C COMPUTE THE INITIAL VALUE OF PH, THE SUM OF SQUARES. C THE ROUTE THROUGH THIS SECTION IS TAKEN ONLY ON THE FIRST CALL. C 5190 DO 5200 J=1,NACTV 5200 BTEMP(J)=B(J) IF(NPRNT)5230,5210,5210 5210 WRITE(KW,5220)(BTEMP(J),J=1,NACTV) 5220 FORMAT(/6X,'BSOLV(J)....',4X,1PE15.7,4E15.7/ * (22X,5E15.7)) 5230 CALL FUNC (NACTV,B,NPTS,FIT) PH=0.0 DO 5250 J=1,NPTS IF(YSIG(J))5240,5240,5250 5240 ICRY=1 GO TO 6120 5250 PH=PH+((Y(J)-FIT(J))/YSIG(J))**2 IF(NPRNT)5870,5260,5260 5260 WRITE(KW,5270)PH 5270 FORMAT(/' SOLVE.... INITIAL VALUE OF PHI =',1PG15.7) GO TO 5870 C C CALL DERIV TO COMPUTE P, THE MATRIX OF PARTIAL DERIVATIVES OF C FIT(K) WITH RESPECT TO B(J). P MAY BE APPROXIMATED USING FINITE C DIFFERENCES. C 5280 CALL DERIV (NACTV,B,NPTS,FIT,P,LDIMP) C C COMPUTE THE NORMAL EQUATIONS. C DO 5330 JA=1,NACTV SUM=0.0 DO 5300 J=1,NPTS IF(YSIG(J))5290,5290,5300 5290 ICRY=2 GO TO 6120 5300 SUM=SUM+(Y(J)-FIT(J))*P(J,JA)/YSIG(J)**2 GRAD(JA)=SUM DO 5330 JB=1,JA SUM=0.0 DO 5320 J=1,NPTS IF(YSIG(J))5310,5310,5320 5310 ICRY=3 GO TO 6120 5320 SUM=SUM+P(J,JA)*P(J,JB)/YSIG(J)**2 5330 QSAVE(JA,JB)=SUM C C THE COEFFICIENTS OF THE NORMAL EQUATIONS ARE NOW IN QSAVE. C SCALE THE MATRIX OF COEFFICIENTS SO THAT THE DIAGONAL ELEMENTS C ARE EQUAL TO UNITY. THE GRADIENT VECTOR IS SCALED ACCORDINGLY. C THE GRADIENT VECTOR IS SAVED FOR USE IN COMPUTING COSIN, THE C COSINE OF MARQUARDT-S ANGLE GAMMA (THE ANGLE BETWEEN THE STEP C AND THE NEGATIVE OF THE GRADIENT VECTOR). C C COMPUTE THE SCALE FACTORS. DO 5360 J=1,NACTV IF(QSAVE(J,J))5340,5340,5350 5340 BERR(J)=1.0 GO TO 5360 5350 BERR(J)=ZSQRT(QSAVE(J,J)) 5360 CONTINUE C IF(NPRNT)5410,5370,5390 5370 WRITE(KW,5380) 5380 FORMAT(/' ') 5390 WRITE(KW,5400)(GRAD(J),J=1,NACTV) 5400 FORMAT(/6X,'SCALED GRADIENT....',1PE14.6,4E14.6/ * (26X,5E14.6)) C C SCALE QSAVE AND THE GRADIENT. 5410 DO 5450 J=1,NACTV BERRJ=BERR(J) IF(BERRJ)5420,5420,5430 5420 ICRY=4 GO TO 6120 5430 GRAD(J)=GRAD(J)/BERRJ DO 5450 JJ=1,J DENOM=BERRJ*BERR(JJ) IF(DENOM)5420,5420,5440 5440 QSAVE(J,JJ)=QSAVE(J,JJ)/DENOM 5450 QSAVE(JJ,J)=QSAVE(J,JJ) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C START THE PROCEDURE TO OBTAIN A BETTER B VECTOR. C NSUB=0 IF(ITER)5470,5470,5460 5460 FLA=FLA/FNU C C MOVE THE MATRIX OF COEFFICIENTS INTO THE Q ARRAY. C THIS PRESERVES QSAVE IN CASE FLA MUST BE INCREASED OR THAT C CONVERGENCE IS ACHIEVED AND THE ERROR MATRIX MUST BE RETURNED. C ADD FLA (MARQUARDT-S LAMBDA) TO THE DIAGONAL ELEMENTS. C 5470 DO 5490 J=1,NACTV X(J)=GRAD(J) DO 5480 JJ=1,NACTV 5480 Q(J,JJ)=QSAVE(J,JJ) 5490 Q(J,J)=Q(J,J)+FLA C C SOLVE FOR THE INCREMENTS IN THE PARAMETERS. C M=1 CALL MATNV (Q,NACTV,X,M,DET,LDIMQ) MRANK=NACTV-M C C IF THE COEFFICIENT MATRIX WAS OF RANK ZERO, QUIT. C IF IT WAS RANK DEFICIENT BUT NOT OF RANK ZERO, PRINT A MESSAGE. C IF(M-NACTV)5510,5500,5500 5500 KERR=M GO TO 5160 5510 IF(M)5520,5540,5520 5520 WRITE(KW,5530)MRANK,NACTV 5530 FORMAT(/' RANK-DEFICIENT NORMAL EQUATIONS IN SOLVE.'/ * 6X,'RANK =',I3,6X,'ORDER OF MATRIX =',I3) C C COMPUTE THE COSINE OF GAMMA. 5540 SA=0.0 SB=0.0 SC=0.0 DO 5550 J=1,NACTV SA=SA+X(J)*GRAD(J) SB=SB+X(J)**2 5550 SC=SC+GRAD(J)**2 IF(SB*SC)5560,5560,5570 5560 ICRY=5 GO TO 6120 5570 COSIN=SA/ZSQRT(SB*SC) C C INITIALIZE STFAC, THE CUTSTEP FACTOR. STFAC WILL BE DECREASED C IF THE SUM OF SQUARES HAS NOT DECREASED AND COSIN IS .GT. CCRIT C (IN MARQUARDT-S NOTATION, GAMMA IS LESS THAN GAMMA(ZERO)). C STFAC=1.0 C C INCREMENT THE PARAMETER VECTOR BTEMP AND CHECK FOR C CONSTRAINT VIOLATIONS. C 5580 DO 5620 J=1,NACTV IF(BERR(J))5590,5590,5600 5590 ICRY=6 GO TO 6120 5600 BTEMP(J)=B(J)+STFAC*X(J)/BERR(J) IF(BTEMP(J)-BMIN(J))5630,5610,5610 5610 IF(BTEMP(J)-BMAX(J))5620,5620,5630 5620 CONTINUE GO TO 5650 5630 NSUB=NSUB+1 WRITE(KW,5640)J,J,BTEMP(J) 5640 FORMAT(//' PARAMETER NUMBER',I3, * ' HAS VIOLATED A CONSTRAINT.',6X,' B(',I2,') =', * 1PG15.7) IF(NSUB-MXSUB)5850,5850,5810 C C CHECK FOR A REDUCTION IN THE WEIGHTED SUM OF SQUARES. C 5650 WRITE(KW,5220)(BTEMP(J),J=1,NACTV) CALL FUNC(NACTV,BTEMP,NPTS,FTEMP) SSQ=0.0 DO 5670 J=1,NPTS IF(YSIG(J))5660,5660,5670 5660 ICRY=7 GO TO 6120 5670 SSQ=SSQ+((Y(J)-FTEMP(J))/YSIG(J))**2 IF(NPRNT)5700,5680,5680 5680 WRITE(KW,5690)ITER,PH,SSQ,COSIN,FLA,STFAC 5690 FORMAT(/' SOLVE.... ITERATION',I4,5X,'OLD PHI =', * 1PG15.7,5X,'NEW PHI =',G15.7,5X,'COS(GAMMA) =',G10.3/ * 47X,'LAMBDA =',G11.3,5X,'K =',G11.3) 5700 IF(SSQ-PH)5870,5710,5740 5710 IF(NFLAT)5800,5800,5720 5720 WRITE(KW,5730) 5730 FORMAT(/' SOLVE CONVERGED UNDER THE "NFLAT" OPTION.') GO TO 5790 5740 IF(NSUB)5800,5800,5750 5750 IF(NFLSB)5800,5800,5760 5760 IF(SSQ-SSUB)5800,5770,5800 5770 WRITE(KW,5780) 5780 FORMAT(/' SOLVE CONVERGED UNDER THE "NFLSB" OPTION.') 5790 ICON=0 GO TO 6010 5800 SSUB=SSQ NSUB=NSUB+1 IF(NSUB-MXSUB)5830,5830,5810 5810 WRITE(KW,5820) MXSUB 5820 FORMAT(//' MORE THAN MXSUB =',I4, * ' SUBITERATIONS REQUIRED.'//' ') ICON=0 GO TO 6010 C C IF THERE IS NO REDUCTION, AND THE ANGLE IS LESS THAN THE CRITICAL C ANGLE, USE CUT-STEPS. IF THE ANGLE EXCEEDS THE CRITICAL ANGLE, C INCREASE FLA AND SOLVE THE NORMAL EQUATIONS AGAIN. C 5830 IF(SSQ-PH)5870,5840,5840 5840 IF(COSIN-CCRIT)5860,5860,5850 C C THE ANGLE IS .LT. THE CRITICAL ANGLE. 5850 STFAC=STFAC/FCUT FCUT=FCUT*FMULT GO TO 5580 C THE ANGLE IS .GE. THE CRITICAL ANGLE. 5860 FLA=FLA*FNU GO TO 5470 5870 CONTINUE C END OF DO XXX I1=L,2 . C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C A NEW PARAMETER VECTOR, BTEMP, HAS BEEN FOUND WHICH DECREASES C THE SUM OF SQUARES. TEST FOR CONVERGENCE. C DO 5880 J=1,NPTS 5880 FIT(J)=FTEMP(J) C C FIRST TEST FOR NO LARGE CHANGES IN THE PARAMETERS. C ICON=0 DO 5940 J=1,NACTV PHL=BTEMP(J)-B(J) IF(PHL)5890,5900,5900 5890 PHL=-PHL 5900 B(J)=BTEMP(J) DENOM=B(J) IF(DENOM)5910,5920,5920 5910 DENOM=-DENOM 5920 DENOM=DENOM+1.0E-20 IF(PHL-EPS2*DENOM)5940,5940,5930 5930 ICON=ICON+1 5940 CONTINUE C C IF EPS2 IS GREATER THAN ZERO, THEN ICON WILL CONTAIN THE NUMBER C OF PARAMETERS WHICH FAIL TO MEET THE CONVERGENCE CRITERION. C PHL=PH-SSQ IF(PHL)5950,5960,5960 5950 PHL=-PHL 5960 PH=SSQ IF(EPS1)5990,5990,5970 5970 IF(PHL-EPS1*PH)5990,5990,5980 5980 ICON=ICON+1 5990 IF(ICON)6010,6010,6000 C C IF EPS1 AND EPS2 ARE BOTH NONZERO, THEN ICON WILL BE THE SUM OF C THE NUMBER OF PARAMETERS NOT MEETING THE CONVERGENCE REQUIREMENT, C PLUS ONE IF THE SUM-OF-SQUARES REQUIREMENT IS NOT MET. C C IF EITHER CONVERGENCE TEST FAILS AND THE NUMBER OF ITERATIONS IS C LESS THAN MAXIT, WE RETURN TO THE CALLING PROGRAM. C OTHERWISE, COMPUTE THE ERROR MATRIX. C 6000 IF(ITER-MAXIT)5160,6010,6010 C 6010 M=0 CALL MATNV (QSAVE,NACTV,X,M,DET,LQSAV) C C DE-SCALE THE INVERSE TO FORM THE ERROR MATRIX. C DO 6040 J=1,NACTV BERRJ=BERR(J) DO 6030 JJ=1,J DENOM=BERRJ*BERR(JJ) IF(DENOM)6020,6020,6030 6020 ICRY=8 GO TO 6120 6030 QSAVE(J,JJ)=QSAVE(J,JJ)/DENOM 6040 CONTINUE C C QSAVE NOW CONTAINS THE ERROR MATRIX. COMPUTE THE ERRORS. C DO 6080 J=1,NACTV IF(QSAVE(J,J))6050,6050,6070 6050 WRITE(KW,6060)J,J,QSAVE(J,J) 6060 FORMAT(/' NEGATIVE OR ZERO MEAN SQUARE ERROR.'/5X, * 'QSAVE(',I2,') =',1PG15.7/5X, * 'THEREFORE THE ERROR MATRIX IS MEANINGLESS.'/' ') BERR(J)=1.0 GO TO 6080 6070 BERR(J)=ZSQRT(QSAVE(J,J)) 6080 CONTINUE C COMPUTE THE CORRELATION MATRIX AND RETURN C IT IN QSAVE. DO 6110 J=1,NACTV BERRJ=BERR(J) DO 6110 JJ=1,J DENOM=BERRJ*BERR(JJ) IF(DENOM)6090,6090,6100 6090 ICRY=9 GO TO 6120 6100 QSAVE(J,JJ)=QSAVE(J,JJ)/DENOM 6110 QSAVE(JJ,J)=QSAVE(J,JJ) GO TO 5160 6120 WRITE(KW,6130)ICRY 6130 FORMAT(//////' ILLEGAL VALUE OF DIVISOR AT CHECKPOINT', * I2,' IN SUBROUTINE SOLVE'//////' ') STOP C END SOLVE. END SUBROUTINE FUNC (NACTV,BSOLV,NPTS,FITSL) C C FUNC IS CALLED BY SOLVE TO COMPUTE THE FITTED VALUES, FITSL(J). C FUNC CALLS FOFX ONCE FOR EACH VALUE TO BE COMPUTED. C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION FITSL C DIMENSION BSOLV(1),FITSL(1) C COMMON /CRICK/ X(100),DFRAC,ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,KALC,NMAX,NPRNT,KW, * LOGB(30),MSK(30),NPAR,NFUPR COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C IF(NPRNT-2)120,100,100 100 WRITE(KW,110) 110 FORMAT(//' INTEGRATED SOLUTION OF THE SYSTEM OF', * ' DIFFERENTIAL EQUATIONS....'/' ') C 120 NSTPS=0 JPT=0 C LOOP OVER THE SETS OF DATA. DO 130 JS=1,NSETS JSET=JS MPTS=JPTS(JS) C LOOP OVER THE POINTS IN THE JS-TH SET. DO 130 J=1,MPTS KPTIN=J JPT=JPT+1 C CALL FOFX TO COMPUTE FITJ AT THIS POINT. C CALL FOFX (X(JPT),BSOLV,NACTV,KPTIN,NEQFU,NUSED) FITSL(JPT)=FITJ NSTPS=NSTPS+NUSED 130 CONTINUE C IF(NPRNT)160,160,140 140 WRITE(KW,150)NSTPS 150 FORMAT(/' FUNC.... TOTAL NUMBER OF FUNCTION', * ' EVALUATIONS IN GEM =',I6) 160 RETURN C END FUNC. END SUBROUTINE DERIV (NACTV,BSOLV,NPTS,FITSL,P,LP) C C COMPUTES APPROXIMATIONS TO THE PARTIAL DERIVATIVES OF FIT(J) WITH C RESPECT TO B(L). THE FINITE DIFFERENCES ARE EITHER FIRST ORDER C APPROXIMATIONS (KALC=1) OR SECOND ORDER (KALC=2). THE LATTER C TAKE ABOUT TWICE AS LONG TO COMPUTE. C C IF KALC=1, THE VALUES OF FITSL(J) MUST HAVE BEEN SET (BY A CALL TO C FUNC) BEFORE CALLING DERIV. C FITSL IS NOT ALTERED BY DERIV. C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION FITSL,FFIT,DEL,DEFIT,DPF C DIMENSION BSOLV(1),FITSL(1),P(LP,NACTV) DIMENSION DEFIT(100) C COMMON /CRICK/ X(100),DFRAC,ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,KALC,NMAX,NPRNT,KW, * LOGB(30),MSK(30),NPAR,NFUPR COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C IF(NPRNT)120,120,100 100 WRITE(KW,110) 110 FORMAT(///' DERIV.... FINITE DIFFERENCE', * ' APPROXIMATIONS OF PARTIAL DERIVATIVES OF', * ' FIT(J) WITH RESPECT TO B(L)....') C C LOOP OVER THE ACTIVE PARAMETERS. 120 DO 250 L=1,NACTV DEL=DFRAC*BSOLV(L) IF(DEL)130,140,150 130 DEL=-DEL GO TO 150 140 DEL=DFRAC 150 BSAVE=BSOLV(L) C COMPUTE -KALC- SETS OF FITJ-S FOR THE C FINITE DIFFERENCE APPROXIMATIONS. DO 240 JJ=1,KALC C AVOID MIXED MODE ARITHMETIC.... DPF=3-2*JJ FFIT=BSAVE BSOLV(L)=FFIT+DPF*DEL NSTPS=0 JPT=0 C LOOP OVER THE SETS OF DATA. DO 210 JS=1,NSETS JSET=JS MPTS=JPTS(JS) C LOOP OVER THE POINTS IN THE JS-TH SET. DO 200 J=1,MPTS KPTIN=J JPT=JPT+1 C CALL FOFX TO COMPUTE FITJ AT THIS POINT. C CALL FOFX (X(JPT),BSOLV,NACTV,KPTIN,NEQDE,NUSED) IF(KALC-1)160,160,170 C C KALC=1. COMPUTE THE FINITE DIFFERENCE C APPROXIMATION (FIRST ORDER). C 160 P(JPT,L)=(FITJ-FITSL(JPT))/DEL GO TO 200 C 170 IF(JJ-1)180,180,190 C C KALC=2, JJ=1. SAVE FITJ. 180 DEFIT(JPT)=FITJ GO TO 200 C KALC=2, JJ=2. COMPUTE THE SYMMETRIC C FINITE DIFFERENCE (A SECOND ORDER C APPROXIMATION TO THE FIRST PARTIAL C DERIVATIVES). C 190 P(JPT,L)=(DEFIT(JPT)-FITJ)/(DEL+DEL) 200 NSTPS=NSTPS+NUSED 210 CONTINUE IF(NPRNT)240,240,220 220 WRITE(KW,230)L,NSTPS 230 FORMAT(6X,'L =',I4,6X,I6,' FUNCTION EVALUATIONS', * ' IN GEM') 240 CONTINUE 250 BSOLV(L)=BSAVE C RETURN C END DERIV. END SUBROUTINE GEM(NEQ,KN,INT,X,XF,H,HMAX,Y,ERR,RERR,FUNCT,NMAX,NUSED) C C GEM 3.0 A.N.S.I. STANDARD FORTRAN FEBRUARY 1973 C J. P. CHANDLER, COMPUTER SCIENCE DEPT., OKLAHOMA STATE UNIVERSITY C C INTEGRATES A SYSTEM OF SIMULTANEOUS FIRST-ORDER ORDINARY C DIFFERENTIAL EQUATIONS USING A FOURTH-ORDER KUTTA METHOD WITH C AUTOMATIC STEP SIZE CONTROL. C PRESENTLY SET UP FOR THE GENERALIZED NEWTON-S THREE-EIGHTHS RULE. C FOR THE THREE-EIGHTHS RULE THE STEP SIZE CONTROL DOES NOT BREAK DOWN C IN THE CASE WHERE A DERIVATIVE IS A FUNCTION ONLY OF X, AS IT DOES C IN MERSON-S METHOD. C SEE M.I.T. LINCOLN LABS REPORT 6D-355 (10/1/1957) BY L. D. EARNEST. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C USAGE........ C CALL GEM ONCE FOR EACH INTERVAL OF INTEGRATION. C C INPUT QUANTITIES..... NEQ,KN,INT,X,XF,H,HMAX,Y,ERR,RERR,FUNCT,NMAX C OUTPUT QUANTITIES.... X,H,Y,NUSED C INTEGER MM C C NEQ -- NUMBER OF EQUATIONS IN THE SYSTEM C KN -- SET KN ZERO FOR THE FIRST INTERVAL IN EACH C INTEGRATION, AND SET IT NONZERO THEREAFTER. C INT -- SET INT NONZERO IF ALL DERIVATIVES ARE FUNCTIONS C ONLY OF X. C X -- THE INDEPENDENT VARIABLE. SET X TO THE INITIAL C VALUE BEFORE THE FIRST CALL TO GEM. IT WILL C BE RETURNED EQUAL TO THE FINAL VALUE BY GEM. C XF -- THE DESIRED FINAL VALUE OF X FOR THIS STEP C H -- THE STEP SIZE. SET H BEFORE THE FIRST STEP OF EACH C INTEGRATION. IT WILL BE READJUSTED BY GEM. C HMAX -- THE MAXIMUM VALUE OF H TO BE USED C Y(J) -- THE J-TH DEPENDENT VARIABLE. SET Y TO THE INITIAL C VALUES BEFORE THE FIRST INTERVAL OF INTEGRATION. C ERR(J) -- ABSOLUTE CONTROL ON THE ERROR IN Y(J) PER STEP C RERR(J) -- RELATIVE CONTROL ON THE ERROR IN Y(J) PER STEP. C A STEP IS TAKEN IF THE ABSOLUTE ERROR IS .LE. C ERR OR IF THE RELATIVE ERROR IS .LE. RERR. C FUNCT -- THE NAME OF THE SUBROUTINE THAT EVALUATES THE C DERIVATIVES. CALL FUNCT(X,Y,DYDX) MUST RETURN C THE VECTOR OF DERIVATIVES DYDX(J) OF Y(J) WITH C RESPECT TO X. THE NAME OF EACH FUNCT C SUBROUTINE MUST APPEAR IN AN EXTERNAL STATEMENT C IN THE CALLING PROGRAM. C NMAX -- MAXIMUM NUMBER OF CALLS TO BE MADE TO FUNCT PER STEP C NUSED -- RETURNS THE NUMBER OF CALLS ACTUALLY MADE. IF THE C INTEGRATION IS UNSUCCESSFUL, NUSED WILL BE C RETURNED .LT. ZERO. INTEGER INK C C EXAMPLE OF USAGE.... D2Y/DX2=-Y , Y(0)=0 , DY/DX(0)=1 C WRITE AS TWO FIRST ORDER EQUATIONS... DY(1)/DX=Y(2) , DY(2)/DX=-Y(1) C DIMENSION Y(2),ERR(2),RERR(2) C EXTERNAL TD C X=0.0 $ Y(1)=0.0 $ Y(2)=1.0 C ERR(1)=ERR(2)=RERR(1)=RERR(2)=1.0E-5 C KN=0 $ H=0.5 C DO 1 J=1,50 C CALL GEM(2,KN,0,X,0.5*J,H,1.0,Y,ERR,RERR,TD,1000,NUSED) C IF(NUSED.LT.0) STOP C KN=1 $ D=SIN(X)-Y(1) C 1 PRINT 2,X,H,Y(1),D,NUSED C 2 FORMAT(' X =',1PE13.5,5X,'H =',E13.5,5X,'Y =', C * E15.7/5X,'ERROR =',E13.5,5X,'NUSED =',I5) C END C SUBROUTINE TD (NEQ,X,Y,DYDX) C DIMENSION Y(2),DYDX(2) C DYDX(1)=Y(2) $ DYDX(2)=-Y(1) $ RETURN C END C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE FOLLOWING STATEMENTS CONVERT GEM TO DOUBLE PRECISION. C IF DOUBLE PRECISION IS NOT USED, SUBROUTINE GEM IS THEN WRITTEN C ENTIRELY IN A.N.S.I. STANDARD BASIC FORTRAN. C IN EITHER CASE, GEM CONTAINS NO MIXED MODE ARITHMETIC. C GEM CALLS NO LIBRARY FUNCTIONS. C DOUBLE PRECISION X,XF,H,HMAX,Y,ERR,RERR,YT,FT,F DOUBLE PRECISION XX,FF,HT,HHMAX,HH,XT,ABHH,XL,RATFC,TT,ABERR, * ABRER,XSAVE,A,B,C,D,ALPHA,BETA,GAMMA,DELTA,RUNIT,RTWO,RTHRE, * FRAC,RATIO,UPFAC,DWNFC,SCALE,U,V,EPS,SFOUR,R,FAC,S,T,CST,VMR DOUBLE PRECISION REIGH,RNINE C DIMENSION Y(1),ERR(1),RERR(1) DIMENSION YT(40),FT(40),F(40,4) C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C YES. SET ALL CONSTANTS. C NEQMX ... MAXIMUM VALUE OF NEQ. NEQMX IS C ALSO THE DIMENSION OF YT AND FT, AND C THE FIRST DIMENSION OF F. NEQMX=40 C RATIO ... INITIAL VALUE FOR RATFC RATIO=16. C FRAC ... USED IN CHECKING FOR FINAL STEP FRAC=.99 C UPFAC ... FACTOR FOR INCREASING RATFC UPFAC=1.2 C DWNFC ... FACTOR FOR DECREASING RATFC DWNFC=.98 C SCALE ... FACTOR FOR CHANGING STEP SIZE SCALE=2. C RATIO=AMAX1(RATIO,SCALE**4) SFOUR=SCALE**4 IF(RATIO-SFOUR)20,30,30 20 RATIO=SFOUR 30 RUNIT=1.0 RTWO=2. RTHRE=3. C SET THE CONSTANTS FOR CASE I WITH U=1/3 , C V=2/3 , AND EPS=3/2 . REIGH=8. RNINE=9. U=RUNIT/RTHRE V=RTWO/RTHRE EPS=RTHRE/RTWO R=RUNIT VMR=-U S=-RUNIT T=RUNIT CST=RUNIT A=RUNIT/REIGH B=RTHRE/REIGH C=B D=A ALPHA=-B BETA=RNINE/REIGH GAMMA=-BETA DELTA=GAMMA C ACTIVATE THE FOLLOWING STATEMENTS TO C OBTAIN THE GENERAL CASE I. C DOUBLE PRECISION RFOUR,RFIVE,RSIX,RTWEL,DIDLE C RFOUR=4. C RFIVE=5. C RSIX=6. C RTWEL=12. C R=V*(U-V)/(RTWO*U*(RTWO*U-RUNIT)) C VMR=V-R C FAC=RSIX*U*V-RFOUR*U-RFOUR*V+RTHRE C S=(RUNIT-U)*(U-RFOUR*V**2+RFIVE*V-RTWO)/(RTWO*U*(V-U)*FAC) C T=(RUNIT-RTWO*U)*(RUNIT-U)*(RUNIT-V)/(V*(V-U)*FAC) C CST=RUNIT-S-T C A=(RSIX*U*V-RTWO*U-RTWO*V+RUNIT)/(RTWEL*U*V) C B=(RTWO*V-RUNIT)/(RTWEL*U*(V-U)*(RUNIT-U)) C C=(RTWO*U-RUNIT)/(RTWEL*V*(U-V)*(RUNIT-V)) C D=FAC/(RTWEL*(RUNIT-U)*(RUNIT-V)) C DIDLE=EPS*(RTWO*U-RUNIT)*(RTWO*V-RUNIT)/RTWO C ALPHA=DIDLE/(U*V) C BETA=-DIDLE/(U*(V-U)*(RUNIT-U)) C GAMMA=DIDLE/(V*(V-U)*(RUNIT-V)) C DELTA=-EPS*FAC/(RTWO*(RUNIT-U)*(RUNIT-V)) C C IQUIT ... =1 IF THIS IS THE FINAL STEP 40 IQUIT=0 C MOVE THE NON-SUBSCRIPTED INPUT VARIABLES. XX=X FF=XF HT=H IF(XX-FF)50,770,50 50 MM=NEQ C CHECK VALIDITY OF INPUT VALUES. IF(MM)80,80,60 60 IF(MM .GT. NEQMX) STOP 70 N=NMAX IF(N)80,80,90 80 NDONE=-1 GO TO 780 C HHMAX=ABS(HMAX) 90 HHMAX=HMAX IF(HHMAX)100,110,110 100 HHMAX=-HHMAX C HH=SIGN(AMIN1(ABS(H),HHMAX),XF-X) 110 HH=HT IF(HH)120,130,130 120 HH=-HH 130 IF(HH-HHMAX)150,150,140 140 HH=HHMAX 150 IF(FF-XX)160,170,170 160 HH=-HH C CHECK THAT HH IS NOT NEGLIGIBLE. 170 XT=XX+HH IF(XT-XX)180,230,180 C SEE WHETHER WE CAN FINISH IN ONE STEP. C IF(ABS(HH)-0.99*ABS(XF-X)) 180 ABHH=HH IF(ABHH)190,200,200 190 ABHH=-ABHH 200 XL=FRAC*(FF-XX) IF(XL)210,220,220 210 XL=-XL 220 IF(ABHH-XL)240,230,230 230 HH=FF-XX C IQUIT=1 240 ITCH=0 C H=AMIN1(ABS(HH),HHMAX) HT=ABHH IF(HT-HHMAX)260,260,250 250 HT=HHMAX C CHECK WHETHER THIS STEP IS A CONTINUATION. 260 NDONE=0 IF(KN)270,280,270 270 IF(XX-XSAVE)280,300,280 280 RATFC=RATIO CALL FUNCT (MM,XX,Y,FT) NDONE=1 DO 290 K=1,MM 290 F(K,1)=FT(K) C START OF EACH STEP.... C CHECK THAT HH IS NOT NEGLIGIBLE. 300 XT=XX+HH IF(XT-XX)310,760,310 310 DO 320 K=1,MM 320 YT(K)=Y(K)+U*HH*F(K,1) CALL FUNCT (MM,XX+U*HH,YT,FT) DO 330 K=1,MM F(K,2)=FT(K) 330 YT(K)=Y(K)+HH*(VMR*F(K,1)+R*F(K,2)) CALL FUNCT (MM,XX+V*HH,YT,FT) DO 340 K=1,MM F(K,3)=FT(K) 340 YT(K)=Y(K)+HH*(CST*F(K,1)+S*F(K,2)+T*F(K,3)) CALL FUNCT (MM,XX+HH,YT,FT) DO 350 K=1,MM F(K,4)=FT(K) 350 YT(K)=Y(K)+HH*(A*F(K,1)+B*F(K,2)+C*F(K,3)+D*F(K,4)) NDONE=NDONE+3 IF(INT)370,360,370 360 CALL FUNCT (MM,XX+HH,YT,FT) NDONE=NDONE+1 C YT NOW CONTAINS THE Y VALUES AT THE END OF C THE STEP. ESTIMATE AND CHECK THE C TRUNCATION ERROR FOR EACH K. INK=1 C IF WE CAN INCREASE THE STEP SIZE. 370 INK=1 DO 530 K=1,MM TT=HH*(ALPHA*F(K,1)+BETA*F(K,2)+GAMMA*F(K,3)+DELTA*F(K,4)+ X EPS*FT(K)) C TT=ABS(TT) IF(TT)380,390,390 380 TT=-TT C IF(TT-ABS(ERR(K))) 390 ABERR=ERR(K) IF(ABERR)400,410,410 400 ABERR=-ABERR C IF(TT-ABS(YT(K)*RERR(K))) 410 ABRER=YT(K)*RERR(K) IF(ABRER)420,430,430 420 ABRER=-ABRER 430 IF(TT-ABERR)500,500,440 440 IF(TT-ABRER)500,500,450 C DECREASE THE STEP SIZE. 450 HH=HH/SCALE C H=ABS(HH) HT=HH IF(HH)460,470,470 460 HT=-HH 470 IQUIT=0 C WAS THE STEP SIZE INCREASED LAST STEP.... IF(ITCH)490,490,480 C YES. RATFC IS TOO SMALL. INCREASE IT. 480 RATFC=UPFAC*RATFC C THE STEP SIZE IS BEING DECREASED (ITCH=-1). 490 ITCH=-1 GO TO 750 C CHECK WHETHER STEP SIZE CAN BE INCREASED. 500 IF(RATFC*TT-ABERR)530,530,510 510 IF(RATFC*TT-ABRER)530,530,520 C C NO IT CANNOT. COMPONENT K FAILS. 520 INK=0 530 CONTINUE C END OF TRUNCATION ERROR TESTS. C THE STEP WAS SUCCESSFUL. MOVE YT AND FT. DO 540 K=1,MM Y(K)=YT(K) 540 F(K,1)=FT(K) C WAS THIS THE FINAL STEP.... IF(IQUIT)550,550,770 C NO. INCREMENT X. 550 XX=XX+HH C ARE THE ERRORS SMALL ENOUGH TO INCREASE C THE STEP SIZE.... IF(INK)580,580,560 C YES. WAS IT DECREASED ON THE PREVIOUS C STEP.... 560 IF(ITCH)570,610,610 C YES. RATFC IS TOO SMALL. INCREASE IT. 570 RATFC=UPFAC*RATFC GO TO 600 C ALL IS WELL. DECREASE RATFC SLIGHTLY. C RATFC=AMAX1(DWNFC*RATFC,RATIO) 580 RATFC=DWNFC*RATFC IF(RATFC-RATIO)590,600,600 590 RATFC=RATIO C THE STEP SIZE IS NOT BEING CHANGED. 600 ITCH=0 GO TO 690 C THE STEP SIZE IS BEING INCREASED (ITCH=1). 610 ITCH=+1 C INCREASE THE STEP SIZE. C HH=SIGN(AMIN1(ABS(SCALE*HH),HHMAX),HH) ABHH=SCALE*HH IF(ABHH)620,630,630 620 ABHH=-ABHH 630 IF(ABHH-HHMAX)650,650,640 640 ABHH=HHMAX 650 IF(HH)660,670,670 660 ABHH=-ABHH 670 HH=ABHH C H=ABS(HH) HT=HH IF(HH)680,690,690 680 HT=-HH C SEE WHETHER WE CAN FINISH IN ONE STEP. C IF(ABS(HH)-0.99*ABS(XF-X)) 690 ABHH=HH IF(ABHH)700,710,710 700 ABHH=-ABHH 710 XL=FRAC*(FF-XX) IF(XL)720,730,730 720 XL=-XL 730 IF(ABHH-XL)750,740,740 C YES WE CAN. SET THE IQUIT SWITCH. 740 HH=FF-XX IQUIT=1 GO TO 300 C CHECK NDONE, AND START THE NEXT STEP. 750 IF(NDONE-N)300,300,760 C TOO MANY STEPS -- INTEGRATION UNSUCCESSFUL. 760 NDONE=-NDONE GO TO 780 C THE INTEGRATION WAS SUCCESSFUL. 770 XX=FF XSAVE=XX C RETURN MODIFIED VALUES OF X, H, AND NUSED. 780 X=XX H=HT NUSED=NDONE RETURN C END GEM. END SUBROUTINE FOFX (XX,BSOLV,NACTV,KPTIN,NEQ,NUSED) C C FOFX 1.2 FEBRUARY 1974 C C FOFX COMPUTES AND RETURNS THE VALUE OF THE FITTED FUNCTION, FITJ, C AT XX. C THIS VERSION IS USED WITH MODELS THAT ARE THE SOLUTIONS OF C SYSTEMS OF DIFFERENTIAL EQUATIONS. THE EQUATIONS ARE INTEGRATED C FROM THE PREVIOUS XX TO THE CURRENT XX BY SUBROUTINE GEM. C GEM CALLS SUBROUTINE DIFFUN TO COMPUTE THE VALUES OF THE DERIVATIVES. C FOFX CALLS SUBROUTINE DINIT TO INITIALIZE XZERO AND YDE(J) FOR C THE INTEGRATION TO THE FIRST XX OF EACH DATA SET. C THIS VERSION OF FOFX IS USED WITH VERSIONS 1.2 AND LATER OF DINIT, C AND PERMITS THE INITIAL CONDITIONS TO BE FIXED OR TO BE ADJUSTABLE C PARAMETERS TO BE FITTED. C C JSET IS THE INDEX OF THE SET OF DATA OF WHICH XX IS PART. C KPTIN IS THE INDEX OF THE XX DATA POINT IN THAT SET. C EXTERNAL DIFFUN C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION ERR,RERR,DXX,DXOLD,DH,DHMAX,DHSAV,DYDT C DIMENSION BSOLV(1) DIMENSION ERR(40),RERR(40),DYDT(40) C C THE STEP SIZE FOR INTEGRATING OVER THE FIRST INTERVAL, FOR EACH C JSET (EACH SET OF DATA), IS SAVED IN DHSAV. C COMMON /CDHSAV/ DHSAV(10) C COMMON /CRICK/ X(100),DFRAC,ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,KALC,NMAX,NPRNT,KW, * LOGB(30),MSK(30),NPAR,NFUPR COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C DXX=XX C C UNPACK THE PARAMETER VECTOR, EXPONENTIATING IT IF NECESSARY. C J=0 DO 130 K=1,NPAR IF(MSK(K))130,110,130 110 J=J+1 BD(K)=BSOLV(J) IF(LOGB(K))120,130,120 120 BD(K)=10.0**BSOLV(J) 130 CONTINUE IF(KPTIN-1)140,140,160 C C THIS IS THE FIRST CALL TO FOFX FOR THIS SET OF DATA (THAT IS, THE C CALL TO FOFX FOR THE FIRST POINT IN THIS SET). C 140 KN=0 CALL DINIT (NPAR,BD,JSET,CUSER,NCUSR,NEQ,XZERO,YDE) DXOLD=XZERO JSTEP=0 IF(XZERO-XX)160,150,160 150 DH=0.0 NUSED=1 C FOFX 63 IS DELETED AND FOFX 81 THRU 84 FOLLOW FOFX 62 IF(NPRNT-2)280,220,220 220 WRITE(KW,230) 230 FORMAT(/' ') GO TO 280 C 160 IF(KN)250,170,250 C C THIS IS THE FIRST INTEGRATION IN THIS SET OF DATA. C INITIALIZE. C 170 IF(DHSAV(JSET))190,180,190 180 DH=DXX-DXOLD GO TO 200 C USE THE FIRST STEPSIZE FROM THE PREVIOUS C INTEGRATION FOR THIS SET OF DATA. 190 DH=DHSAV(JSET) 200 DHMAX=1.0E30 DO 210 J=1,NEQ ERR(J)=ER 210 RERR(J)=RER C FOFX 81 THRU 84 ARE NOW MOVED TO BELOW FOFX 62 C C CALL GEM TO INTEGRATE THE SYSTEM OF DIFFERENTIAL EQUATIONS TO XX, C THE ABSCISSA OF THE NEXT DATA POINT. C 250 CALL GEM (NEQ,KN,0,DXOLD,DXX,DH,DHMAX,YDE,ERR,RERR,DIFFUN,NMAX, * NUSED) IF(DHSAV(JSET))270,260,270 260 DHSAV(JSET)=DH 270 KN=1 JSTEP=JSTEP+1 C CALL DIFFUN TO SET FITJ FOR THE FINAL C VALUES OF YDE. C 280 CALL DIFFUN (NEQ,DXOLD,YDE,DYDT) IF(NPRNT-2)290,300,300 290 IF(NUSED)300,340,340 300 NYPR=NEQ IF(NFUPR-NYPR)302,302,301 301 NYPR=NFUPR 302 WRITE(KW,310)JSTEP,XX,DH,FITJ,NUSED,(YDE(J),J=1,NYPR) 310 FORMAT(' STEP',I5,5X,'XT =',1PE15.7,5X,'HH =', * E13.5,5X,'FITJ =',E15.7/5X,'NUSED =',I6,5X, * 'YDE(J)....'/(5X,5E15.7)) IF(NUSED)320,320,340 320 WRITE(KW,330)DXOLD,DXX,DH,NMAX,NUSED 330 FORMAT(///' GEM FAILURE.',31PE20.10,2E20.10,2I8) STOP 340 RETURN C END FOFX. END SUBROUTINE DINIT (NPAR,BD,JSET,CUSER,NCUSR,NEQ,XZERO,YDE) C C DINIT 1.2 FEBRUARY 1974 C C INITIALIZES XZERO AND YDE(J) FOR THE INTEGRATION OF THE SYSTEM OF C DIFFERENTIAL EQUATIONS. C THIS SUBROUTINE IS PROBLEM-DEPENDENT. C C THIS VERSION OF DINIT ALLOWS THE INITIAL VALUES TO BE FIXED C (SET FROM CUSER(J)) OR TO BE ADJUSTABLE PARAMETERS TO BE FITTED C (SET FROM BD(K)). C DOUBLE PRECISION YDE,CUSER,BD DIMENSION CUSER(1),YDE(1),BD(1) C XZERO=0.0 DO 10 J=1,NEQ 10 YDE(J)=0.0 RETURN C END DINIT. END SUBROUTINE MATNV (A,N,B,M,DET,MA) C C CO-OP 61 NBSB MATINV - MATRIX INVERSION WITH SOLUTION OF LINEAR EQNS. C C INVERTS -A-, THE N BY N MATRIX OF COEFFICIENTS, SOLVES M SYSTEMS OF C N LINEAR EQUATIONS, AND EVALUATES THE DETRMINANT OF -A-. C ON ENTRY A CONTAINS THE MATRIX OF COEFFICIENTS AND B CONTAINS THE M C VECTORS OF CONSTANTS. ON EXIT A CONTAINS THE INVERSE OF THE MATRIX C OF COEFFICIENTS AND B CONTAINS THE M SOLUTION VECTORS. C ON EXIT M CONTAINS (N MINUS THE RANK OF -A-). C MA IS THE FIRST DIMENSION OF THE ARRAYS IN WHICH THE MATRICES -A- C AND -B- ARE STORED. C C USES THE GAUSS-JORDAN METHOD WITH COMPLETE (DOUBLE) PIVOTING. C C ORIGINALLY A SHARE PROGRAM BY B. GARBOW (ARGONNE NATIONAL LAB). C THIS VERSION TREATS THE LARGEST NONSINGULAR SUBMATRIX OF -A-.... C IF ANY PIVOT IS .LE. PVBAD (.GE.ZERO), THE REMAINDER OF -A- AND C B IS ZEROED OUT, -A- IS UNSCRAMBLED, AND MATINV RETURNS. C NOTE.... IF PVBAD IS .GT. ZERO, THE MATRIX -A- MUST BE SCALED.... C DIMENSION A(MA,N),B(MA,1),INDEX(50,2) C DIMENSION A(15,15),B(15,1),INDEX(15,2) C C INITIALIZATION C C NINDX IS THE FIRST DIMENSION OF -INDEX-. NINDX=15 NINDX=50 PVBAD=0.001 PVBAD=0.0 ZERO=0.0 UNITY=1.0 NVAR=N IF(NVAR-NINDX)1010,1010,1000 1000 IRANK=0 DETRM=ZERO GO TO 1310 1010 NB=M DETRM=UNITY DO 1020 J=1,NVAR DO 1020 JL=1,2 1020 INDEX(J,JL)=0 C C SEARCH FOR PIVOT ELEMENT C I=0 IRANK=0 1030 AMAX=ZERO DO 1100 J=1,NVAR IF(INDEX(J,1))1100,1040,1100 1040 DO 1090 K=1,NVAR IF(INDEX(K,1))1090,1050,1090 1050 T=A(J,K) IF(T)1060,1070,1070 1060 T=-T 1070 IF(T-AMAX)1090,1090,1080 1080 IROW=J KOLUM=K AMAX=T 1090 CONTINUE 1100 CONTINUE IF(AMAX-PVBAD)1320,1320,1110 1110 INDEX(KOLUM,1)=IROW C C INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL C IF(IROW-KOLUM)1120,1160,1120 1120 DETRM=-DETRM DO 1130 L=1,NVAR SWAP=A(IROW,L) A(IROW,L)=A(KOLUM,L) 1130 A(KOLUM,L)=SWAP IF(NB)1160,1160,1140 1140 DO 1150 L=1, NB SWAP=B(IROW,L) B(IROW,L)=B(KOLUM,L) 1150 B(KOLUM,L)=SWAP 1160 I=I+1 INDEX(I,2)=KOLUM PIVOT=A(KOLUM,KOLUM) DETRM=PIVOT*DETRM IRANK=IRANK+1 C C DIVIDE PIVOT ROW BY PIVOT ELEMENT C A(KOLUM,KOLUM)=UNITY DO 1170 L=1,NVAR 1170 A(KOLUM,L)=A(KOLUM,L)/PIVOT IF(NB)1200,1200,1180 1180 DO 1190 L=1,NB 1190 B(KOLUM,L)=B(KOLUM,L)/PIVOT C C REDUCE NON-PIVOT ROWS C 1200 DO 1250 L1=1,NVAR IF(L1-KOLUM)1210,1250,1210 1210 T=A(L1,KOLUM) A(L1,KOLUM)=ZERO DO 1220 L=1,NVAR 1220 A(L1,L)=A(L1,L)-A(KOLUM,L)*T IF(NB)1250,1250,1230 1230 DO 1240 L=1,NB 1240 B(L1,L)=B(L1,L)-B(KOLUM,L)*T 1250 CONTINUE IF(I-NVAR) 1030,1260,1260 C C INTERCHANGE COLUMNS C 1260 KOLUM=INDEX(I,2) IROW=INDEX(KOLUM,1) IF(IROW-KOLUM) 1270,1290,1270 1270 DO 1280 K=1,NVAR SWAP=A(K,IROW) A(K,IROW)=A(K,KOLUM) 1280 A(K,KOLUM)=SWAP 1290 I=I-1 1300 IF(I)1310,1310,1260 1310 M=NVAR-IRANK DET=DETRM RETURN C C ZERO OUT THE REST OF A AND B. C 1320 DO 1370 J=1,NVAR DO 1330 K=1,NVAR IF(INDEX(K,2)-J)1330,1370,1330 1330 CONTINUE DO 1340 JL=1,NVAR A(JL,J)=ZERO 1340 A(J,JL)=ZERO IF(M)1370,1370,1350 1350 DO 1360 JL=1,M 1360 B(J,JL)=ZERO 1370 CONTINUE GO TO 1300 END SUBROUTINE GIVENS ( NX, NROOTX, NJX, A, B, ROOT, VECT ) C 62.3 GIVENS -EIGENVALUES AND EIGENVECTORS BY THE GIVENS METHOD. C BY FRANKLIN PROSSER, INDIANA UNIVERSITY. C SEPTEMBER, 1967 C CALCULATES EIGENVALUES AND EIGENVECTORS OF REAL SYMMETRIC MATRIX C STORED IN PACKED UPPER TRIANGULAR FORM. C C THANKS ARE DUE TO F. E. HARRIS (STANFORD UNIVERSITY) AND H. H. C MICHELS (UNITED AIRCRAFT RESEARCH LABORATORIES) FOR EXCELLENT C WORK ON NUMERICAL DIFFICULTIES WITH EARLIER VERSIONS OF THIS C PROGRAM. C C THE PARAMETERS FOR THE ROUTINE ARE... C NX ORDER OF MATRIX C NROOTX NUMBER OF ROOTS WANTED. THE NROOTX SMALLEST (MOST C NEGATIVE) ROOTS WILL BE CALCULATED. IF NO VECTORS C ARE WANTED, MAKE THIS NUMBER NEGATIVE. C NJX ROW DIMENSION OF VECT ARRAY. SEE 'VECT' BELOW. C NJX MUST BE NOT LESS THAN NX. C A MATRIX STORED BY COLUMNS IN PACKED UPPER TRIANGULAR C FORM, I.E. OCCUPYING NX*(NX+1)/2 CONSECUTIVE C LOCATIONS. C B SCRATCH ARRAY USED BY GIVENS. MUST BE AT LEAST C NX*5 CELLS. C ROOT ARRAY TO HOLD THE EIGENVALUES. MUST BE AT LEAST C NROOTX CELLS LONG. THE NROOTX SMALLEST ROOTS ARE C ORDERED LARGEST FIRST IN THIS ARRAY. INTEGER N C VECT EIGENVECTOR ARRAY. EACH COLUMN WILL HOLD AN C EIGENVECTOR FOR THE CORRESPONDING ROOT. MUST BE C DIMENSIONED WITH 'NJX' ROWS AND AT LEAST 'NROOTX' C COLUMNS, UNLESS NO VECTORS C ARE REQUESTED (NEGATIVE NROOTX). IN THIS LATTER C CASE, THE ARGUMENT VECT IS JUST A DUMMY, AND THE C STORAGE IS NOT USED. C THE EIGENVECTORS ARE NORMALIZED TO UNITY. C C THE ARRAYS A AND B ARE DESTROYED BY THE COMPUTATION. THE RESULTS C APPEAR IN ROOT AND VECT. C FOR PROPER FUNCTIONING OF THIS ROUTINE, THE RESULT OF A FLOATING C POINT UNDERFLOW SHOULD BE A ZERO. C TO CONVERT THIS ROUTINE TO DOUBLE PRECISION (E.G. ON IBM 360 C MACHINES), BE SURE THAT ALL REAL VARIABLES AND FUNCTION C REFERENCES ARE PROPERLY MADE DOUBLE PRECISION. C THE VALUE OF 'ETA' (SEE BELOW) SHOULD ALSO BE CHANGED, TO REFLECT C THE INCREASED PRECISION. INTEGER NROOT C C THE ORIGINAL REFERENCE TO THE GIVENS TECHNIQUE IS IN OAK RIDGE C REPORT NUMBER ORNL 1574 (PHYSICS), BY WALLACE GIVENS. C THE METHOD AS PRESENTED IN THIS PROGRAM CONSISTS OF FOUR STEPS, C ALL MODIFICATIONS OF THE ORIGINAL METHOD... C FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE C HOUSEHOLDER TECHNIQUE (J. H. WILKINSON, COMP. J. 3, 23 (1960)). C THE ROOTS ARE THEN LOCATED BY THE STURM SEQUENCE METHOD (J. M. C ORTEGA (SEE REFERENCE BELOW). THE VECTORS OF THE TRIDIAGONAL C FORM ARE THEN EVALUATED (J. H. WILKINSON, COMP. J. 1, 90 (1958)), C AND LAST THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE C ORIGINAL ARRAY (FIRST REFERENCE). C VECTORS FOR DEGENERATE (OR NEAR-DEGENERATE) ROOTS ARE FORCED C TO BE ORTHOGONAL, USING A METHOD SUGGESTED BY B. GARBOW, ARGONNE C NATIONAL LABS (PRIVATE COMMUNICATION, 1964). THE GRAM-SCHMIDT C PROCESS IS USED FOR THE ORTHOGONALIZATION. C C AN EXCELLENT PRESENTATION OF THE GIVENS TECHNIQUE IS FOUND IN C J. M. ORTEGA'S ARTICLE IN 'MATHEMATICS FOR DIGITAL COMPUTERS,' C VOLUME 2, ED. BY RALSTON AND WILF, WILEY (1967), PAGE 94. C DIMENSION B(NX,5),A(100),ROOT(10),VECT(NJX,NROOTX) C ZABS(ARG)=ABS(ARG) ZSIGN(ARGA,ARGB)=SIGN(ARGA,ARGB) ZMAX1(ARGA,ARGB)=AMAX1(ARGA,ARGB) ZMIN1(ARGA,ARGB)=AMIN1(ARGA,ARGB) ZMOD(ARGA,ARGB)=AMOD(ARGA,ARGB) ZSQRT(ARG)=SQRT(ARG) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C ** USERS PLEASE NOTE... C ** THE FOLLOWING TWO PARAMETERS, ETA AND THETA, SHOULD BE ADJUSTED C ** BY THE USER FOR HIS PARTICULAR MACHINE. C ** ETA IS AN INDICATION OF THE PRECISION OF THE FLOATING POINT C ** REPRESENTATION ON THE COMPUTER BEING USED (ROUGHLY 10**(-M), C ** WHERE M IS THE NUMBER OF DECIMALS OF PRECISION ). C ** THETA IS AN INDICATION OF THE RANGE OF NUMBERS THAT CAN BE C ** EXPRESSED IN THE FLOATING POINT REPRESENTATION (ROUGHLY THE C ** LARGEST NUMBER). C ** SOME RECOMMENDED VALUES FOLLOW. C ** FOR IBM 7094, UNIVAC 1108, ETC. (27-BIT BINARY FRACTION, 8-BIT C ** BINARY EXPONENT), ETA=1.0E-8, THETA=1.0E37. C ** FOR CONTROL DATA 3600 (36-BIT BINARY FRACTION, 11-BIT BINARY C ** EXPONENT), ETA=1.0E-11, THETA=1.0E307. C ** FOR CONTROL DATA 6600 (48-BIT BINARY FRACTION, 11-BIT BINARY C ** EXPONENT), ETA=1.0E-14, THETA=1.0E307. C ** FOR IBM 360/50 AND 360/65 DOUBLE PRECISION (56-BIT HEXADECIMAL C ** FRACTION, 7-BIT HEXADECIMAL EXPONENT), ETA=1.E-16, THETA=1.E75. C ** ETA=1.0E-7 THETA = 1.0E75 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DEL1 = ETA/100.0 DELTA = ETA**2*100.0 SMALL = ETA**2/100.0 DELBIG = THETA*DELTA/1000.0 THETA1 = 1000.0/THETA C TOLER IS A FACTOR USED TO DETERMINE IF TWO ROOTS ARE CLOSE C ENOUGH TO BE CONSIDERED DEGENERATE FOR PURPOSES OF ORTHOGONALI- C ZING THEIR VECTORS. FOR THE MATRIX NORMED TO UNITY, IF THE C DIFFERENCE BETWEEN TWO ROOTS IS LESS THAN TOLER, THEN C ORTHOGONALIZATION WILL OCCUR. TOLER = ETA*100.0 C C INITIAL VALUE FOR PSEUDORANDOM NUMBER GENERATOR... (2**23)-3 RPOWER = 8388608. RPOW1 = RPOWER/2. RAND1 = RPOWER - 3. C N = NX NROOT = IABS(NROOTX) IF (NROOT.EQ.0) GO TO 1640 IF (N-1) 1640,1000,1010 1000 ROOT(1) = A(1) IF (NROOTX.GT.0) VECT(1,1) = 1.0 GO TO 1640 1010 CONTINUE C NSIZE NUMBER OF ELEMENTS IN THE PACKED ARRAY NSIZE = (N*(N+1))/2 NM1 = N-1 NM2 = N-2 C C SCALE MATRIX TO EUCLIDEAN NORM OF 1. SCALE FACTOR IS ANORM. FACTOR = 0.0 DO 1020 I=1,NSIZE 1020 FACTOR= ZMAX1 ( FACTOR, ZABS( A(I) ) ) IF (FACTOR.NE.0.0) GO TO 1050 C NULL MATRIX. FIX UP ROOTS AND VECTORS, THEN EXIT. DO 1040 I=1,NROOT IF (NROOTX.LT.0) GO TO 1040 DO 1030 J=1,N 1030 VECT(J,I) = 0.0 VECT(I,I) = 1.0 1040 ROOT(I) = 0.0 GO TO 1640 C 1050 ANORM = 0.0 J = 1 K = 1 DO 1070 I=1,NSIZE IF (I.NE.J) GO TO 1060 ANORM = ANORM + (A(I)/FACTOR)**2/2. K = K+1 J = J+K GO TO 1070 1060 ANORM = ANORM + (A(I)/FACTOR)**2 1070 CONTINUE ANORM = ZSQRT( ANORM * 2. ) * FACTOR DO 1080 I=1,NSIZE 1080 A(I) = A(I)/ANORM ALIMIT = 1.0 C C TRIDIA SECTION. C TRIDIAGONALIZATION OF SYMMETRIC MATRIX ID = 0 IA = 1 IF (NM2.EQ.0) GO TO 1190 DO 1180 J=1,NM2 C J COUNTS ROW OF A-MATRIX TO BE DIAGONALIZED C IA START OF NON-CODIAGONAL ELEMENTS IN THE ROW C ID INDEX OF CODIAGONAL ELEMENT ON ROW BEING CODIAGONALIZED. IA = IA+J+2 ID = ID + J + 1 JP2 = J+2 C SUM SQUARES OF NON-CODIAGONAL ELEMENTS IN ROW J II = IA SUM = 0.0 DO 1090 I=JP2,N SUM=SUM+A(II)**2 1090 II = II + I TEMP = A(ID) IF (SUM.GT.SMALL) GO TO 1100 C NO TRANSFORMATION NECESSARY IF ALL THE NON-CODIAGONAL C ELEMENTS ARE TINY. B(J,1) = TEMP A(ID) = 0.0 GO TO 1180 C NOW COMPLETE THE SUM OF OFF-DIAGONAL SQUARES 1100 SUM = ZSQRT ( SUM + TEMP ** 2 ) C NEW CODIAGONAL ELEMENT B(J,1) =-ZSIGN(SUM,TEMP) C FIRST NON-ZERO ELEMENT OF THIS W-VECTOR B(J+1,2) = ZSQRT (( 1.0 + ZABS ( TEMP ) / SUM ) / 2.0 ) C FORM REST OF THE W-VECTOR ELEMENTS TEMP =ZSIGN(0.5/(B(J+1,2)*SUM),TEMP) II = IA DO 1110 I=JP2,N B(I,2) = A(II)*TEMP 1110 II = II + I C FORM P-VECTOR AND SCALAR. P-VECTOR = A-MATRIX*W-VECTOR. C SCALAR = W-VECTOR*P-VECTOR. AK = 0.0 C IC LOCATION OF NEXT DIAGONAL ELEMENT IC = ID + 1 J1 = J + 1 DO 1140 I=J1,N JJ = IC TEMP = 0.0 DO 1130 II=J1,N C I RUNS OVER THE NON-ZERO P-ELEMENTS C II RUNS OVER ELEMENTS OF W-VECTOR TEMP = TEMP + B(II,2)*A(JJ) C CHANGE INCREMENTING MODE AT THE DIAGONAL ELEMENTS. IF (II.LT.I) GO TO 1120 JJ = JJ + II GO TO 1130 1120 JJ = JJ + 1 1130 CONTINUE C BUILD UP THE K-SCALAR (AK) AK = AK + TEMP*B(I,2) B(I,1) = TEMP C MOVE IC TO TOP OF NEXT A-MATRIX 'ROW' 1140 IC = IC + I C FORM THE Q-VECTOR DO 1150 I=J1,N 1150 B(I,1) = B(I,1) - AK*B(I,2) C TRANSFORM THE REST OF THE A-MATRIX C JJ START-1 OF THE REST OF THE A-MATRIX JJ = ID C MOVE W-VECTOR INTO THE OLD A-MATRIX LOCATIONS TO SAVE SPACE C I RUNS OVER THE SIGNIFICANT ELEMENTS OF THE W-VECTOR DO 1170 I=J1,N A(JJ) = B(I,2) DO 1160 II=J1,I JJ = JJ + 1 1160 A(JJ) = A(JJ) - 2.0*(B(I,1)*B(II,2) + B(I,2)*B(II,1)) 1170 JJ = JJ + J 1180 CONTINUE C MOVE LAST CODIAGONAL ELEMENT OUT INTO ITS PROPER PLACE 1190 CONTINUE B(NM1,1) = A(NSIZE-1) A(NSIZE-1) = 0.0 C C STURM SECTION. C STURM SEQUENCE ITERATION TO OBTAIN ROOTS OF TRIDIAGONAL FORM. C MOVE DIAGONAL ELEMENTS INTO SECOND N ELEMENTS OF B-VECTOR. C THIS IS A MORE CONVENIENT INDEXING POSITION. C ALSO, PUT SQUARE OF CODIAGONAL ELEMENTS IN THIRD N ELEMENTS. JUMP=1 DO 1200 J=1,N B(J,2)=A(JUMP) B(J,3) = B(J,1)**2 1200 JUMP = JUMP+J+1 DO 1210 I=1,NROOT 1210 ROOT(I) = +ALIMIT ROOTL = -ALIMIT C ISOLATE THE ROOTS. THE NROOT LOWEST ROOTS ARE FOUND, LOWEST FIRST DO 1290 I=1,NROOT C FIND CURRENT 'BEST' UPPER BOUND ROOTX = +ALIMIT DO 1220 J=I,NROOT 1220 ROOTX = ZMIN1 ( ROOTX, ROOT(J) ) ROOT(I) = ROOTX C GET IMPROVED TRIAL ROOT 1230 TRIAL = (ROOTL + ROOT(I))*0.5 IF (TRIAL.EQ.ROOTL.OR.TRIAL.EQ.ROOT(I)) GO TO 1290 C FORM STURM SEQUENCE RATIOS, USING ORTEGA'S ALGORITHM (MODIFIED). C NOMTCH IS THE NUMBER OF ROOTS LESS THAN THE TRIAL VALUE. CONTINUE NOMTCH = N J = 1 1240 F0 = B(J,2) - TRIAL 1250 CONTINUE IF ( ZABS ( F0 ).LT.THETA1 ) GO TO 1260 IF (F0.GE.0.0) NOMTCH = NOMTCH - 1 J = J + 1 IF (J.GT.N) GO TO 1270 C SINCE MATRIX IS NORMED TO UNITY, MAGNITUDE OF B(J,3) IS LESS THAN C ONE, AND SINCE F0 IS GREATER THAN THETA1, OVERFLOW IS NOT POSSIBLE C AT THE DIVISION STEP. F0 = B(J,2) - TRIAL - B(J-1,3)/F0 GO TO 1250 1260 J = J + 2 NOMTCH = NOMTCH - 1 IF (J.LE.N) GO TO 1240 1270 CONTINUE C FIX NEW BOUNDS ON ROOTS IF (NOMTCH.GE.I) GO TO 1280 ROOTL = TRIAL GO TO 1230 1280 ROOT(I) = TRIAL NOM = MIN0(NROOT,NOMTCH) ROOT(NOM) = TRIAL GO TO 1230 1290 CONTINUE C C TRIVEC SECTION. C EIGENVECTORS OF CODIAGONAL FORM CONTINUE C QUIT NOW IF NO VECTORS WERE REQUESTED. IF (NROOTX.LT.0) GO TO 1620 C INITIALIZE VECTOR ARRAY. DO 1300 I=1,N DO 1300 J=1,NROOT VECT ( N+1 , J ) = 0.0 VECT ( N+2 , J ) = 0.0 1300 VECT(I,J) = 1.0 DO 1560 I=1,NROOT AROOT = ROOT(I) C ORTHOGONALIZE IF ROOTS ARE CLOSE. IF (I.EQ.1) GO TO 1310 C THE ABSOLUTE VALUE IN THE NEXT TEST IS TO ASSURE THAT THE TRIVEC C SECTION IS INDEPENDENT OF THE ORDER OF THE EIGENVALUES. IF ( ZABS ( ROOT ( I-1 ) - AROOT ) .LT. TOLER ) GO TO 1320 1310 IA = -1 1320 IA = IA + 1 ELIM1 = A(1) - AROOT ELIM2 = B(1,1) JUMP = 1 DO 1350 J=1,NM1 JUMP = JUMP+J+1 C GET THE CORRECT PIVOT EQUATION FOR THIS STEP. IF ( ZABS ( ELIM1 ) .LE. ZABS ( B (J,1 ) )) GO TO 1330 C FIRST (ELIM1) EQUATION IS THE PIVOT THIS TIME. CASE 1. B(J,2) = ELIM1 B(J,3) = ELIM2 B(J,4) = 0.0 TEMP = B(J,1)/ELIM1 ELIM1 = A(JUMP) - AROOT - TEMP*ELIM2 ELIM2 = B(J+1,1) GO TO 1340 C SECOND EQUATION IS THE PIVOT THIS TIME. CASE 2. 1330 B(J,2) = B(J,1) B(J,3) = A(JUMP) - AROOT B(J,4) = B(J+1,1) TEMP = 1.0 IF ( ZABS ( B (J,1) ) .GT. THETA1 ) TEMP = ELIM1 / B(J,1) ELIM1 = ELIM2 - TEMP*B(J,3) ELIM2 = -TEMP*B(J+1,1) C SAVE FACTOR FOR THE SECOND ITERATION. 1340 B(J,5) = TEMP 1350 CONTINUE B(N,2) = ELIM1 B(N,3) = 0.0 B(N,4) = 0.0 B(NM1,4) = 0.0 ITER = 1 IF (IA.NE.0) GO TO 1450 C BACK SUBSTITUTE TO GET THIS VECTOR. 1360 L = N + 1 DO 1400 J=1,N L = L - 1 1370 CONTINUE ELIM1=VECT(L,I)-VECT(L+1,I)*B(L,3)-VECT(L+2,I)*B(L,4) C IF OVERFLOW IS CONCEIVABLE, SCALE THE VECTOR DOWN. C THIS APPROACH IS USED TO AVOID MACHINE-DEPENDENT AND SYSTEM- C DEPENDENT CALLS TO OVERFLOW ROUTINES. IF ( ZABS( ELIM1 ) .GT. DELBIG ) GO TO 1380 TEMP = B(L,2) IF ( ZABS( B(L,2) ) .LT.DELTA ) TEMP = DELTA VECT(L,I) = ELIM1/TEMP GO TO 1400 C VECTOR IS TOO BIG. SCALE IT DOWN. 1380 DO 1390 K=1,N 1390 VECT(K,I) = VECT(K,I)/DELBIG GO TO 1370 1400 CONTINUE GO TO (1410,1470), ITER C SECOND ITERATION. (BOTH ITERATIONS FOR REPEATED-ROOT VECTORS). 1410 ITER = ITER + 1 1420 ELIM1 = VECT(1,I) DO 1440 J=1,NM1 IF (B(J,2).EQ.B(J,1)) GO TO 1430 C CASE ONE. VECT(J,I) = ELIM1 ELIM1 = VECT(J+1,I) - ELIM1*B(J,5) GO TO 1440 C CASE TWO. 1430 VECT(J,I) = VECT(J+1,I) ELIM1 = ELIM1 - VECT(J+1,I)*TEMP 1440 CONTINUE VECT(N,I) = ELIM1 GO TO 1360 C PRODUCE A RANDOM VECTOR 1450 CONTINUE DO 1460 J=1,N C GENERATE PSEUDORANDOM NUMBERS WITH UNIFORM DISTRIBUTION IN (-1,1). C THIS RANDOM NUMBER SCHEME IS OF THE FORM... C RAND1 = AMOD((2**12+3)*RAND1,2**23) C IT HAS A PERIOD OF 2**21 NUMBERS. RAND1 = ZMOD(4099.*RAND1,RPOWER) 1460 VECT(J,I) = RAND1/RPOW1 - 1.0 GO TO 1360 C C ORTHOGONALIZE THIS REPEATED-ROOT VECTOR TO OTHERS WITH THIS ROOT. 1470 IF (IA.EQ.0) GO TO 1510 DO 1500 J1=1,IA K = I - J1 TEMP = 0.0 DO 1480 J=1,N 1480 TEMP = TEMP + VECT(J,I)*VECT(J,K) DO 1490 J=1,N 1490 VECT(J,I) = VECT(J,I) - TEMP*VECT(J,K) 1500 CONTINUE 1510 GO TO (1420,1520), ITER C NORMALIZE THE VECTOR 1520 ELIM1 = 0.0 DO 1530 J=1,N 1530 ELIM1 = ZMAX1 ( ZABS(VECT(J,I)),ELIM1 ) TEMP=0.0 DO 1540 J=1,N ELIM2=VECT(J,I)/ELIM1 TEMP=TEMP+ELIM2**2 1540 CONTINUE TEMP = 1.0 /(ZSQRT ( TEMP ) * ELIM1 ) DO 1550 J=1,N VECT(J,I) = VECT(J,I)*TEMP IF ( ZABS ( VECT(J,I) ) .LT. DEL1 ) VECT (J,I ) = 0.0 1550 CONTINUE 1560 CONTINUE C C SIMVEC SECTION. C ROTATE CODIAGONAL VECTORS INTO VECTORS OF ORIGINAL ARRAY C LOOP OVER ALL THE TRANSFORMATION VECTORS IF (NM2.EQ.0) GO TO 1620 JUMP = NSIZE - (N+1) IM = NM1 DO 1610 I=1,NM2 J1 = JUMP C MOVE A TRANSFORMATION VECTOR OUT INTO BETTER INDEXING POSITION. DO 1570 J=IM,N B(J,2) = A(J1) 1570 J1 = J1 + J C MODIFY ALL REQUESTED VECTORS. DO 1600 K=1,NROOT TEMP = 0.0 C FORM SCALAR PRODUCT OF TRANSFORMATION VECTOR WITH EIGENVECTOR DO 1580 J=IM,N 1580 TEMP = TEMP + B(J,2)*VECT(J,K) TEMP = TEMP + TEMP DO 1590 J=IM,N 1590 VECT(J,K) = VECT(J,K) - TEMP*B(J,2) 1600 CONTINUE JUMP = JUMP - IM 1610 IM = IM - 1 1620 CONTINUE C RESTORE ROOTS TO THEIR PROPER SIZE. DO 1630 I=1,NROOT 1630 ROOT(I) = ROOT(I)*ANORM 1640 RETURN END DOUBLE PRECISION FUNCTION GAUSF (DSEEDA,DSEEDB,INITDS) C C GAUSF 1.3 MAY 1992 C C GENERATES RANDOM NORMAL DEVIATES... C RANDOM NUMBERS FROM A GAUSSIAN DISTRIBUTION WITH ZERO C MEAN AND UNIT VARIANCE. C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C G. MARSAGLIA AND T. A. BRAY, S.I.A.M. REVIEW 6 (1964) 260. C C FOR A RANDOM NUMBER FROM THE GAUSSIAN DISTRIBUTION WITH MEAN C VALUE EQUAL TO 'AMEAN' AND STANDARD DEVIATION EQUAL TO C 'SIGMA', USE ... C C RNG=AMEAN+SIGMA*GAUSF(DSEEDA,DSEEDB,INITDS) C DOUBLE PRECISION DSEEDA,DSEEDB, RJUMP,A,B,X,Y,EX,ABX, * G,VA,VB,SUMSQ,FAC, DSHUFF,DABS,DSQRT,DLOG,DEXP C INTEGER INITDS C RJUMP=DSHUFF(DSEEDA,DSEEDB,INITDS) INITDS=0 IF(RJUMP.LT.0.1362D0) GO TO 20 C 10 A=(RJUMP-0.1362D0)/0.8638D0 B=DSHUFF(DSEEDA,DSEEDB,INITDS) GAUSF=2.0D0*(A+B+DSHUFF(DSEEDA,DSEEDB,INITDS)-1.5D0) GO TO 150 C 20 IF(RJUMP.LT.0.0255D0) GO TO 40 C 30 A=(RJUMP-0.0255D0)/0.1107D0 GAUSF=1.5D0*(A+DSHUFF(DSEEDA,DSEEDB,INITDS)-1.0D0) GO TO 150 C 40 IF(RJUMP.LT.0.0026997961D0) GO TO 110 C 50 X=6.0D0*DSHUFF(DSEEDA,DSEEDB,INITDS)-3.0D0 Y=0.358D0*DSHUFF(DSEEDA,DSEEDB,INITDS) EX=17.49731196D0*DEXP(-0.5D0*X**2) ABX=DABS(X) IF(ABX.GE.1.0D0) GO TO 70 C 60 G=EX-4.73570326D0*(3.0D0-X**2)-2.15787544D0*(1.5D0-ABX) GO TO 90 C 70 G=EX-2.36785163D0*(3.0D0-ABX)**2 IF(ABX.LT.1.5D0) G=G-2.15787544D0*(1.5D0-ABX) 90 IF(Y.GT.G) GO TO 50 C 100 GAUSF=X GO TO 150 C 110 VA=2.0D0*DSHUFF(DSEEDA,DSEEDB,INITDS)-1.0D0 VB=2.0D0*DSHUFF(DSEEDA,DSEEDB,INITDS)-1.0D0 SUMSQ=VA**2+VB**2 IF(SUMSQ.GT.1.0D0 .OR. SUMSQ.LE.0.0D0) GO TO 110 C 130 FAC=DSQRT((9.0D0-2.0D0*DLOG(SUMSQ))/SUMSQ) GAUSF=VA*FAC IF(DABS(GAUSF).GE.3.0D0) GO TO 150 GAUSF=VB*FAC IF(DABS(GAUSF).LT.3.0D0) GO TO 110 C 150 RETURN C C END GAUSF C END DOUBLE PRECISION FUNCTION DSHUFF (DSEEDA,DSEEDB,INITDS) C C DSHUFF 1.2 MAY 1992 C C DSHUFF IS A PORTABLE DOUBLE PRECISION UNIFORM (0,1) C PSEUDORANDOM NUMBER GENERATOR. C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C DSHUFF USES A SHUFFLING METHOD. THE METHOD HAS A VERY LONG C PERIOD (AT LEAST AS LARGE AS DMODA*DMODB, AND PERHAPS MUCH C LARGER) AND SHOULD HAVE GOOD RANDOMNESS PROPERTIES. C C DSEEDA -- ONE CURRENT "SEED" FOR THE GENERATOR. C ON THE FIRST ENTRY, DSEEDA SHOULD BE SET C TO A DOUBLE PRECISION INTEGER VALUE C BETWEEN 1.0D0 AND 2147483646.0D0, C INCLUSIVE, FOR EXAMPLE C DSEEDA = 123456789.0D0 . C AFTER THAT, THE USER MUST NOT CHANGE THE C VALUE OF DSEEDA. C (SUBROUTINE DSHUFF WILL ALTER DSEEDA C AS NEEDED.) C C DSEEDB -- THE OTHER CURRENT SEED FOR THE GENERATOR. C ON THE FIRST ENTRY, DSEEDB SHOULD BE SET C TO A DOUBLE PRECISION INTEGER VALUE C BETWEEN 1.0D0 AND 2147483710.0D0, C INCLUSIVE. AFTER THAT, THE USER MUST NOT C CHANGE THE VALUE OF DSEEDB. C THE SAME INITIAL VALUE MAY BE USED FOR C DSEEDA AND DSEEDB, IF DESIRED. C C INITDS -- =1 ON THE FIRST CALL TO SUBROUTINE DSHUFF, C =0 ON EACH SUBSEQUENT CALL C C DOUBLE PRECISION DSEEDA,DSEEDB, TABLE, * DMODA,DMODB,DMULT,DLTABL C INTEGER INITDS, LTABLE,J C C SAVE TABLE C COMMON /CDSHUF/ TABLE(128) C C DMODA = 2**31 - 1 , A PRIME. DMODB IS ALSO A PRIME. C DMODA=2147483647.0D0 DMODB=2147483711.0D0 C C DMULT = 16807 = 7**5 . C DMULT=16807.0D0 C C LTABLE IS THE DIMENSION OF THE ARRAY TABLE(*). C LTABLE=128 DLTABL=128.0D0 C C IF INITDS=1, FILL THE TABLE. C IF(INITDS.EQ.1) THEN C DO 10 J=1,LTABLE DSEEDA=DMOD(DMULT*DSEEDA,DMODA) DSEEDB=DMOD(DMULT*DSEEDB,DMODB) TABLE(J)=(DSEEDA+DSEEDB/DMODB)/DMODA 10 CONTINUE C INITDS=0 C ENDIF C C GENERATE A NEW UNIFORM (0,1) PSEUDORANDOM DEVIATE. C 20 DSEEDA=DMOD(DMULT*DSEEDA,DMODA) DSEEDB=DMOD(DMULT*DSEEDB,DMODB) C C SHUFFLE BY PICKING A NUMBER PSEUDORANDOMLY FROM THE TABLE, C RETURNING IT AS OUTPUT, AND REPLACING IT. C J=DMOD(DSEEDB,DLTABL)+1.0D0 IF(J.LT.1 .OR. J.GT.LTABLE) GO TO 20 DSHUFF=TABLE(J) TABLE(J)=(DSEEDA+DSEEDB/DMODB)/DMODA IF(DSHUFF.LE.0.0D0 .OR. DSHUFF.GE.1.0D0) GO TO 20 C RETURN C C END DSHUFF C END FUNCTION CHSPR(CHISQ,NU) C C CHSPR 1.1 JUNE 1992 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C CHSPR IS THE CHI-SQUARE PROBABILITY (CONFIDENCE LEVEL) OF A C GIVEN VALUE OF CHI-SQUARE (CHISQ) FOR NU DEGREES OF FREEDOM. C THAT IS, IT IS THE INTEGRAL FROM CHISQ TO INFINITY OF THE C CHI-SQUARE PROBABILITY DISTRIBUTION FOR NU DEGREES OF C FREEDOM. C OR, IF YOU LIKE, IT IS THE A PRIORI PROBABILITY OF GETTING C A WORSE (LARGER) VALUE OF CHI-SQUARE THAN THE VALUE YOU C ACTUALLY GOT, ASSUMING YOUR HYPOTHESIS WAS CORRECT. C C CHSPR IS ALSO RELATED TO THE FOLLOWING FUNCTIONS.... C C CHSPR(2.*A,2*L) IS THE CUMULATIVE SUM (FROM ZERO TO L-1) OF C THE POISSON DISTRIBUTION WITH MEAN VALUE EQUAL TO -A-. C C CHSPR(2.*X,2*L)=GAMMA(L,X)/GAMMA(L) , C THE INCOMPLETE GAMMA FUNCTION. C C SEE EQUATIONS 26.4.4 AND 26.4.5 (P. 941) OF THE C HANDBOOK OF MATHEMATICAL FUNCTIONS C (NBS APPLIED MATH. SERIES, NO. 55) C DOUBLE PRECISION IS USED ONLY INTERNALLY. C DOUBLE PRECISION TWOLN,CHIS,CHISL,TRMLN,FACT,RZXSM,ROOTP DOUBLE PRECISION CONLX,SAMD,ZERO,UNITY,TWO,Q,UPLIM DOUBLE PRECISION DLOG,DEXP,ZLOG,ZEXP,DATAN,PI C ZLOG(Q)=DLOG(Q) ZEXP(Q)=DEXP(Q) C C ERFCP(X) IS THE COMPLEMENTARY ERROR FUNCTION OF ABSF(X), C X REAL. HANDBOOK OF MATHEMATICAL FUNCTIONS, P.299 . C ERFCP(X,ETA)=EXP (-X**2)*ETA*(.254829592+ETA*(-0.284496736+ * ETA*(1.421413741+ETA*(-1.453152027+ETA*1.061405429)))) C TT(X)=1.0/(1.0+.3275911*ABS (X)) C KW=6 ZERO=0.0D0 UNITY=1.0D0 TWO=2.0D0 UPLIM=1.0001D0 TWOLN=DLOG(2.0) PI=4.0D0*DATAN(1.0D0) ROOTP=DSQRT(2.0D0*PI) C CHSQQ=CHISQ CHIS=CHSQQ N=NU CONLX=-UNITY IF(CHIS)160,10,20 10 CONLX=UNITY GO TO 190 20 IF(N)160,160,30 30 RZXSM=ZERO IF(N-2)110,120,40 40 NUODD=N-2*(N/2) CHISL=ZLOG(CHIS) IF(NUODD)60,60,50 50 TRMLN=CHISL/TWO FACT=UNITY GO TO 70 60 TRMLN=CHISL-TWOLN FACT=TWO C C COMPUTE SQRTF(2*PI)*Z*SUM IN SUCH A WAY AS TO AVOID OVERFLOW C OR HARMFUL UNDERFLOW. C 70 LIMIT=(N-1)/2 DO 90 J=1,LIMIT RZXSM=RZXSM+ZEXP(TRMLN-CHIS/TWO) IF(J-LIMIT)80,100,100 80 FACT=FACT+TWO 90 TRMLN=TRMLN+CHISL-ZLOG(FACT) C 100 IF(NUODD)120,120,110 110 ARG=SQRT(0.5*CHSQQ) SAMD=ERFCP(ARG,TT(ARG)) CONLX=SAMD+(RZXSM/ROOTP)*TWO GO TO 130 120 CONLX=RZXSM+DEXP(-CHIS/TWO) 130 CL=CONLX IF(CONLX)150,140,140 140 IF(CONLX-UPLIM)190,190,160 C 150 NERR=2 CONLX=ZERO GO TO 170 160 NERR=1 CONLX=UNITY 170 WRITE(KW,180)NERR,CHSQQ,N,CL 180 FORMAT(/' ERROR IN CHSPR.'/5X,'NERR =',I5,5X, * 'CHISQ =',1PG15.7,5X,'NU =',I10,5X,'CHSPR =', * G15.7) C 190 CHSPR=CONLX RETURN END SUBROUTINE NRUNS (R,JL,JH,NMI,NZER,NPL,NRMI,NRPL) C C NRUNS 2.0 A.N.S.I. STANDARD BASIC FORTRAN MARCH 1973 C C J. P. CHANDLER, COMPUTER SCIENCE DEPT., OKLAHOMA STATE UNIVERSITY C C GIVEN A LIST OF REAL NUMBERS, R(JL),...,R(JH), SUBROUTINE NRUNS C COUNTS THE NUMBER OF NEGATIVE, ZERO, AND POSITIVE NUMBERS, AND THE C NUMBERS OF RUNS OF NEGATIVE AND POSITIVE NUMBERS. C IN COUNTING THE NUMBERS OF RUNS, ZERO VALUES OF R(J) ARE IGNORED. C DIMENSION R(1) C INITIALIZE. NM=0 NZ=0 NP=0 NRM=0 NRP=0 IF(JH-JL)150,10,10 10 IF(R(JL))20,30,40 20 NM=1 NPREV=-1 NRM=1 GO TO 50 30 NZ=1 NPREV=0 GO TO 50 40 NP=1 NPREV=1 NRP=1 50 IF(JH-JL)150,150,60 C PROCESS THE R(J). 60 JLP=JL+1 DO 140 J=JLP,JH IF(R(J))70,100,110 70 IF(NPREV)90,80,80 80 NRM=NRM+1 90 NM=NM+1 NPREV=-1 GO TO 140 100 NZ=NZ+1 GO TO 140 110 IF(NPREV)120,120,130 120 NRP=NRP+1 130 NP=NP+1 NPREV=1 140 CONTINUE C RETURN THE OUTPUT QUANTITIES. 150 NMI=NM NZER=NZ NPL=NP NRMI=NRM NRPL=NRP RETURN C END NRUNS. END FUNCTION SIGNP(NA,NB) C C SIGNP 1.1 A.N.S.I. STANDARD BASIC FORTRAN MARCH 1973 C C DOYLE HILL, DEPT. OF BIOCHEMISTRY, OKLAHOMA STATE UNIVERSITY C C COMPUTES THE CONFIDENCE LEVEL FOR THE SIGNS TEST ON A SET OF C RESIDUALS. THIS IS THE PROBABILITY OF OBTAINING AS UNBALANCED A C A DISTRIBUTION OF SIGNS AS WAS OBTAINED, IF THE SIGNS WERE RANDOM. C ZEXP(ARG)=EXP(ARG) ZLOG(ARG)=ALOG(ARG) C ZERO=0.0 UNITY=1.0 TWO=2.0 C IR3=NB IRR=NA-NB IF(IRR)10,20,20 10 IR3=NA IRR=-IRR C TEST FOR SIGNP=1. 20 IF(IRR-1)30,30,40 30 PROB=UNITY GO TO 70 C 40 PROB=ZERO IF(IR3)70,70,50 50 N=NA+NB XN=N ANUM2=GAMLN(XN+UNITY) TWON=XN*ZLOG(TWO) C SET -PROB- EQUAL TO THE ZERO-TH TERM. PROB=TWO**(-N) DO 60 J=1,IR3 AJ=J DENOM=GAMLN(AJ+UNITY)+GAMLN(XN-AJ+UNITY)+TWON 60 PROB=PROB+ZEXP(ANUM2-DENOM) C DOUBLE PROB TO FORM A TWO-TAILED TEST. PROB=PROB*TWO 70 SIGNP=PROB RETURN C END SIGNP. END FUNCTION RUNPR(NA,NB,NRA,NRB) C C RUNPR 1.0 A.N.S.I. STANDARD BASIC FORTRAN JULY 1971 C C DOYLE HILL, DEPT. OF BIOCHEMISTRY, OKLAHOMA STATE UNIVERSITY C C RUNPR EVALUATES THE CONFIDENCE LEVEL FOR THE RUNS TEST OF RESIDUALS. C THIS IS THE PROBABILITY OF GETTING AS SMALL (NON-RANDOM) A NUMBER OF C RUNS AS WAS ACTUALLY OBTAINED, IF THE SIGNS OF THE RESIDUALS C WERE POSITIVE OR NEGATIVE WITH EQUAL PROBABILITY. C C REFERENCES..... C F. N. DAVID, -A CHI-SQUARE SMOOTH TEST FOR GOODNESS OF FIT-, C BIOMETRIKA 34 (1947) 299. C P. G. HOEL, -INTRODUCTION TO MATHEMATICAL STATISTICS- C (WILEY, 1954) PAGE 293. C C NA -- NUMBER OF POSITIVE RESIDUALS C NB -- NUMBER OF NEGATIVE RESIDUALS C NRA -- NUMBER OF RUNS OF POSITIVE RESIDUALS C NRB -- NUMBER OF RUNS OF NEGATIVE RESIDUALS C ZLOG(ARG)=ALOG(ARG) ZEXP(ARG)=EXP(ARG) C C INITIALIZE. ZERO=0.0 UNITY=1.0 TWO=2.0 N=NA+NB C IU ... TOTAL NUMBER OF RUNS IU=NRA+NRB XNA=NA XNB=NB XN=N XNRA=NRA XNRB=NRB FOFU=ZERO FUSUM=ZERO IF(IU-2)90,10,10 C C PERM IS THE LOG OF THE NUMBER OF PERMUTATIONS OF NA AND NB IN N. C 10 PERM=GAMLN(XNA+UNITY)+GAMLN(XNB+UNITY)-GAMLN(XN+UNITY) ANUM=GAMLN(XNA)+GAMLN(XNB)+PERM C C COMPUTE THE CONFIDENCE LEVEL. THIS IS THE SUM OF THE PROBABILITIES C OF OBTAINING 2 RUNS, 3 RUNS, ..., IU RUNS. C DO 80 I=2,IU U=I C C DETERMINE IF THE NUMBER OF RUNS IS ODD OR EVEN. C II=I/2 IF(2*II-I)20,70,20 C ODD NUMBER OF RUNS..... 20 UA=(U+UNITY)/TWO UB=(U-UNITY)/TWO CONST=ANUM-GAMLN(UA)-GAMLN(UB) C C IF (XNA-UA) OR (XNB-UA) IS NEGATIVE, THE DENOMINATOR CONTAINS C THE FACTORIAL OF A NEGATIVE INTEGER, AND THE TERM IS ZERO. C TERM1=ZERO IF(XNA-UA)40,30,30 30 TERM1=CONST-GAMLN(XNA-UA+UNITY)-GAMLN(XNB-UB+UNITY) TERM1=ZEXP(TERM1) C 40 TERM2=ZERO IF(XNB-UA)60,50,50 50 TERM2=CONST-GAMLN(XNA-UB+UNITY)-GAMLN(XNB-UA+UNITY) TERM2=ZEXP(TERM2) 60 FOFU=TERM1+TERM2 GO TO 80 C EVEN NUMBER OF RUNS..... C 70 DENOM=TWO*GAMLN(U/TWO)+GAMLN(XNA+UNITY-U/TWO)+ * GAMLN(XNB+UNITY-U/TWO) FOFU=ZEXP(ANUM+ZLOG(TWO)-DENOM) 80 FUSUM=FUSUM+FOFU C 90 RUNPR=FUSUM RETURN C END RUNPR. END FUNCTION GAMLN(X) C C GAMLN 2.0 A.N.S.I. STANDARD BASIC FORTRAN MARCH 1973 C C J. P. CHANDLER, COMPUTER SCIENCE DEPT., OKLAHOMA STATE UNIVERSITY C C NATURAL LOG OF THE GAMMA FUNCTION OF A POSITIVE REAL ARGUMENT. C C ASYMPOTIC EXPANSION FOR LN GAMMA .... C C GLNL(X,XSQ)=(X-0.5)*ALOG(X)-X+0.918938533205+ * (1.0/12.-(1.0/360.-(1.0/1260.-(1.0/1680.)/ * XSQ)/XSQ)/XSQ)/X C C VERSION FOR THE STUPID IBM 360.... C GLNL(X,XSQ)=(X-0.5)*ALOG(X)-X+0.9189385 + * (1.0/12.-(1.0/360.-(1.0/1260.-(1.0/1680.)/ * XSQ)/XSQ)/XSQ)/X C ASYMP=6.0 C Y=X IF(Y)10,20,20 10 Y=-Y 20 IF(Y-ASYMP)40,30,30 30 GAMLN=GLNL(Y,Y*Y) GO TO 60 C RECUR OUT TO THE ASYMPOTIC REGION. 40 N=1.0-(Y-ASYMP) TEMP=1.0 DO 50 J=1,N TEMP=TEMP*Y 50 Y=Y+1.0 GAMLN=GLNL(Y,Y*Y)-ALOG(TEMP) 60 RETURN C C END GAMLN C END SUBROUTINE DIFFUN (NEQ,T,Y,DYDT) C C DIFFUN EVALUATES THE DERIVATIVES FOR GEM, WHICH SOLVES THE C SYSTEM OF DIFFERENTIAL EQUATIONS (SEE P. 7, ADDENDUM OF 2-2-78). C C THIS IS THE FIVE-PARAMETER DIFFUN USED, E.G., TO GENERATE TABLES 1, C 2, AND 3 IN THE 1972 PAPER BY CHANDLER, HILL, AND SPIVEY. C C NOMENCLATURE...... C BD(J) -- K(J) (K SUB J) C Y(1) -- EP C Y(2) -- P C Y(3) -- EAB C YA=YDE(4) -- E C YB=YDE(5) -- A C CUSER(J) -- AT(J), J=1,(NCUSR-1) C CUSER(NCUSR) -- ET C FITJ -- PT C C YDE(4) AND YDE(5) ARE RETURNED IN CASE IT IS DESIRED TO PRINT C THE COMPLETE SET OF SOLUTIONS IN SOME OTHER ROUTINE. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION T,Y,DYDT,YA,YB,AK23,AK45 C DIMENSION Y(1),DYDT(1) C COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C YA=CUSER(NCUSR)-Y(1)-Y(3) YB=CUSER(JSET)-Y(3)-Y(1)-Y(2) AK23=BD(2)+BD(3) AK45=BD(4)+BD(5) DYDT(1)=BD(3)*Y(3)-AK45*Y(1) DYDT(2)=BD(5)*Y(1) DYDT(3)=BD(1)*YA*YB+BD(4)*Y(1)-AK23*Y(3) C C RETURN THE FIT FUNCTION. FITJ=Y(1)+Y(2) C RETURN THE VARIABLES DEFINED BY THE C CONSERVATION EQUATIONS. YDE(4)=YA YDE(5)=YB RETURN END $ENTRY GDH TRANSIENTS --- 25 JUNE 1971 --- 30 MIN CHARCOAL TREATMENT CRIC 5 10 .002 .000606 .00021 .004 .001223 .008 .002169 .020 .004142 .040 .005875 .068 .006958 .104 .008075 .132 .008840 .160 .009622 .200 .010625 13 .004 .000909 .008 .002001 .012 .002883 .016 .003622 .022 .004292 .028 .004975 .040 .005671 .060 .006562 .086 .007476 .112 .008225 .148 .009184 .168 .009772 .200 .010572 12 .003 .000654 .011 .001546 .019 .002230 .027 .002928 .043 .003881 .055 .004612 .071 .005360 .091 .006124 .119 .006904 .149 .007703 .175 .008246 .195 .008658 8 .012 .000605 .024 .001118 .044 .001954 .068 .002810 .096 .003688 .132 .004589 .176 .005513 .196 .005866 8 .011 .000607 .023 .001121 .039 .001748 .079 .003037 .115 .004148 .143 .004832 .175 .005530 .195 .005884 1 15 2000 6 .00001 .0001 10. .01 .0001 1.25 1.25 0.24 0.08 0.08 0.018 5 3 3 1 3 5 1. 5 0 0 0 0 0 0 0 0 0 0 150. 7.33 60.4 63.0 4.35 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2 SUBROUTINE DERIV (NACTV,BSOLV,NPTS,FITSL,P,LP) C C DERIV2 -- VERSION OF DERIV FOR OBTAINING THE DERIVATIVES OF THE C FITTED FUNCTION WITH RESPECT TO THE PARAMETERS BY INTEGRATING C ADDITIONAL DIFFERENTIAL EQUATIONS. C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION FITSL C DIMENSION BSOLV(1),FITSL(1),P(LP,NACTV) C COMMON /CRICK/ X(100),DFRAC,ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,KALC,NMAX,NPRNT,KW, * LOGB(30),MSK(30),NPAR,NFUPR COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C NSTPS=0 JPT=0 C LOOP OVER THE SETS OF DATA. DO 140 JS=1,NSETS JSET=JS MPTS=JPTS(JS) C LOOP OVER THE POINTS IN THE JS-TH SET. DO 140 J=1,MPTS KPSET=J JPT=JPT+1 C CALL FOFX TO COMPUTE DFDB AT THIS POINT. C CALL FOFX (X(JPT),BSOLV,NACTV,KPSET,NEQDE,NUSED) C C UNPACK THE DERIVATIVES INTO THE P ARRAY. K=0 DO 130 L=1,NPAR IF(MSK(L))130,100,130 100 K=K+1 C TRANSFORM THE DERIVATIVE IF LOGB.NE.0 . C IF(LOGB(L))110,120,110 110 DFDB(L)=DFDB(L)*BD(L)*2.302585 120 P(JPT,K)=DFDB(L) 130 CONTINUE 140 NSTPS=NSTPS+NUSED IF(NPRNT)170,170,150 150 WRITE(KW,160)NSTPS 160 FORMAT(///' DERIV2.... TOTAL NUMBER OF FUNCTION', * ' EVALUATIONS IN GEM =',I6) C END DERIV2. END SUBROUTINE DIFFUN (NEQ,T,Y,DYDT) C C DIFFUN2 -- VERSION OF DIFFUN FOR OBTAINING THE DERIVATIVES OF THE C FITTED FUNCTION WITH RESPECT TO THE PARAMETERS BY INTEGRATING C ADDITIONAL DIFFERENTIAL EQUATIONS. C THIS IS THE REGULAR FIVE-PARAMETER MODEL. C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION T,Y,DYDT,YA,YB,AK23,AK45 C DIMENSION Y(1),DYDT(1) C COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C YA=CUSER(NCUSR)-Y(1)-Y(3) YB=CUSER(JSET)-Y(3)-Y(1)-Y(2) AK23=BD(2)+BD(3) AK45=BD(4)+BD(5) DYDT(1)=BD(3)*Y(3)-AK45*Y(1) DYDT(2)=BD(5)*Y(1) DYDT(3)=BD(1)*YA*YB+BD(4)*Y(1)-AK23*Y(3) C C RETURN THE FIT FUNCTION. FITJ=Y(1)+Y(2) C C IF NEQ IS EQUAL TO 3, FOFX WAS CALLED FROM FUNC, NOT FROM DERIV. C IN THIS CASE THE EVALUATION OF THE REMAINING DERIVATIVES IS C UNECCESSARY, AND WE RETURN. C IF(NEQ-3)120,120,100 C C DYDT(4) IS THE DERIVATIVE OF DYDT(1) WITH RESPECT TO BD(1), ETC. C NOTE THAT THESE DERIVATIVES INCLUDE NOT ONLY THE EXPLICIT C DERIVATIVE BUT ALSO THE IMPLICIT CHAIN RULE TERMS THAT ENTER C VIA THE APPEARANCE OF Y(J) IN DYDT(K). C 100 DYDT(4)=BD(3)*Y(14)-AK45*Y(4) DYDT(5)=BD(3)*Y(15)-AK45*Y(5) DYDT(6)=Y(3)+BD(3)*Y(16)-AK45*Y(6) DYDT(7)=BD(3)*Y(17)-Y(1)-AK45*Y(7) DYDT(8)=BD(3)*Y(18)-Y(1)-AK45*Y(8) DYDT(9)=BD(5)*Y(4) DYDT(10)=BD(5)*Y(5) DYDT(11)=BD(5)*Y(6) DYDT(12)=BD(5)*Y(7) DYDT(13)=Y(1)+BD(5)*Y(8) DYDT(14)=YA*YB+BD(1)*((-Y(4)-Y(14))*YB+YA*(-Y(4)-Y(9)-Y(14)))+ * BD(4)*Y(4)-AK23*Y(14) DYDT(15)=BD(1)*((-Y(5)-Y(15))*YB+YA*(-Y(5)-Y(10)-Y(15)))+ * BD(4)*Y(5)-Y(3)-AK23*Y(15) DYDT(16)=BD(1)*((-Y(6)-Y(16))*YB+YA*(-Y(6)-Y(11)-Y(16)))+ * BD(4)*Y(6)-Y(3)-AK23*Y(16) DYDT(17)=BD(1)*((-Y(7)-Y(17))*YB+YA*(-Y(7)-Y(12)-Y(17)))+ * Y(1)+BD(4)*Y(7)-AK23*Y(17) DYDT(18)=BD(1)*((-Y(8)-Y(18))*YB+YA*(-Y(8)-Y(13)-Y(18)))+ * BD(4)*Y(8)-AK23*Y(18) C C RETURN THE PARTIAL DERIVATIVES. DFDB(L) C IS THE PARTIAL DERIVATIVE OF FITJ C WITH RESPECT TO B(L). DO 110 L=1,5 110 DFDB(L)=Y(L+3)+Y(L+8) C 120 RETURN C END DIFFUN2 5K. END $ENTRY GDH TRANSIENTS --- 25 JUNE 1971 --- 30 MIN CHARCOAL TREATMENT CRIC 5 10 .002 .000606 .00021 .004 .001223 .008 .002169 .020 .004142 .040 .005875 .068 .006958 .104 .008075 .132 .008840 .160 .009622 .200 .010625 13 .004 .000909 .008 .002001 .012 .002883 .016 .003622 .022 .004292 .028 .004975 .040 .005671 .060 .006562 .086 .007476 .112 .008225 .148 .009184 .168 .009772 .200 .010572 12 .003 .000654 .011 .001546 .019 .002230 .027 .002928 .043 .003881 .055 .004612 .071 .005360 .091 .006124 .119 .006904 .149 .007703 .175 .008246 .195 .008658 8 .012 .000605 .024 .001118 .044 .001954 .068 .002810 .096 .003688 .132 .004589 .176 .005513 .196 .005866 8 .011 .000607 .023 .001121 .039 .001748 .079 .003037 .115 .004148 .143 .004832 .175 .005530 .195 .005884 1 15 2000 6 .00001 .0001 10. .01 .0001 1.25 1.25 0.24 0.08 0.08 0.018 5 3 18 1 0 0 0. 5 0 0 0 0 0 0 0 0 0 0 150. 7.33 60.4 63.0 4.35 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2 SUBROUTINE GEM(NEQ,KN,INT,X,XF,H,HMAX,Y,ERR,RERR,FUNCT,NMAX,NUSED) C C -STIFF- GEM 2.1 A.N.S.I. STANDARD FORTRAN APRIL 1973 C C CONTROL ROUTINE FOR INTEGRATING A SYSTEM OF STIFF FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS. C FOR INSTRUCTIONS ON USAGE, SEE THE COMMENTS IN THE DECK OF THE C NON-STIFF VERSION OF GEM, AND THE WRITE-UP FOR IT. C EXTERNAL FUNCT C DOUBLE PRECISION X,XF,H,HMAX,Y,ERR,RERR,YT,YU, * UPFAC,DNFAC,SCALE,RATST,RUNIT,RCLOS,XX,FF,HT,HH,HHMAX, * XT,ABHH,XL,RATFC,XA,DY,TT,ABERR,ARERR C DIMENSION Y(1),ERR(1),RERR(1) DIMENSION YT(40),YU(40) C COMMON /MINK/ NSW C DATA JFRST/7/ C MOD FUNCTION (NOT IN BASIC FORTRAN).... MOD(J,K)=J-(J/K)*K C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C NDBUG ... =1 TO OBTAIN DEBUG PRINTOUT NDBUG=1 NDBUG=0 C NSW ... =0 FOR RUNGE-KUTTA, C =1 FOR TREANOR-S METHOD, C =2 FOR TREANOR-S MODIFIED METHOD NSW=1 C KW ... LOGICAL UNIT NUMBER OF THE PRINTER C (FOR DEBUG PRINT STATEMENTS) KW=6 C NEQMX ... MAXIMUM VALUE OF NEQ, AND THE C DIMENSION OF THE ARRAYS YT AND YU NEQMX=40 C C KEXTR=1 TO OBTAIN POLYNOMIAL EXTRAPOLATION, KEXTR=0 OTHERWISE. C THE USE OF EXTRAPOLATION MAY IMPROVE THE ACCURACY. IT DEGRADES C THE STABILITY. C KEXTR=1 KEXTR=0 C M IS THE ORDER OF THE INTEGRATION METHOD. M=4 C UPFAC ... FACTOR FOR INCREASING RATFC UPFAC=1.2 C DNFAC ... FACTOR FOR DECREASING RATFC DNFAC=.98 C SCALE ... FACTOR FOR CHANGING STEP SIZE SCALE=2. C NROT ... =1 TO -ROTATE- THE METHOD (I.E. C TO CHANGE THE INTEGRATION METHOD C WHENEVER THE STEP SIZE IS DECREASED) NROT=0 NROT=1 C RCLOS ... CRITERION FOR TAKING LAST STEP RCLOS=0.99 RUNIT=1.0 RATST=SCALE**M C 110 NDONE=0 IQUIT=0 XX=X FF=XF HT=H C MM=MAX0(1,MIN0(NEQMX,NEQ)) MM=NEQ IF(MM-NEQMX)130,130,120 120 MM=NEQMX 130 IF(MM)140,140,150 140 MM=1 C HHMAX=ABS(HMAX) 150 HHMAX=HMAX IF(HHMAX)160,170,170 160 HHMAX=-HHMAX C HH=SIGN(AMIN1(ABS(H),HHMAX),XF-X) 170 HH=HT IF(HH)180,190,190 180 HH=-HH 190 IF(HH-HHMAX)210,210,200 200 HH=HHMAX 210 IF(FF-XX)220,230,230 220 HH=-HH C CHECK THAT HH IS NOT NEGLIGIBLE. 230 XT=XX+HH IF(XT-XX)240,290,240 C SEE WHETHER WE CAN FINISH IN ONE STEP. C IF(ABS(HH)-RCLOS*ABS(XF-X)) 240 ABHH=HH IF(ABHH)250,260,260 250 ABHH=-ABHH 260 XL=RCLOS*(FF-XX) IF(XL)270,280,280 270 XL=-XL 280 IF(ABHH-XL)300,290,290 290 HH=XF-X IQUIT=1 300 ITCH=0 C H=AMIN1(ABS(HH),HHMAX) HT=ABHH IF(HT-HHMAX)320,320,310 310 HT=HHMAX C CHECK WHETHER THIS STEP IS A CONTINUATION. 320 IF(KN)330,340,330 330 IF(XX-XSAVE)340,350,340 340 RATFC=RATST C 350 XT=XX+HH IF(XT-XX)360,810,360 360 XA=XX CALL HORDE (MM,XA,XA+HH,HH,Y,YT,FUNCT,NUT) NDONE=NDONE+NUT XA=XX CALL HORDE (MM,XA,XA+HH,HH/SCALE,Y,YU,FUNCT,NUT) NDONE=NDONE+NUT INK=1 DO 580 K=1,MM DY=(YT(K)-YU(K))/(RATST-RUNIT) C C SHOULD WE EXTRAPOLATE.... IF(KEXTR)380,380,370 370 YU(K)=YU(K)-DY 380 TT=DY C TT=ABS(TT) IF(TT)390,400,400 390 TT=-TT C IF(TT-ABS(ERR(K))) 400 ABERR=ERR(K) IF(ABERR)410,420,420 410 ABERR=-ABERR C IF(TT-ABS(YT(K)*RERR(K))) 420 ARERR=YT(K)*RERR(K) IF(ARERR)430,440,440 430 ARERR=-ARERR 440 IF(TT-ABERR)550,550,450 450 IF(TT-ARERR)550,550,460 C C DECREASE THE STEP SIZE. 460 HH=HH/SCALE C H=ABS(HH) HT=HH IF(HH)470,480,480 470 HT=-HH 480 IQUIT=0 IF(ITCH)500,500,490 490 RATFC=UPFAC*RATFC 500 ITCH=-1 C SHOULD WE ROTATE THE METHOD.... IF(NROT-1)800,510,800 510 CONTINUE C USE RUNGE-KUTTA OR ONE OF TREANOR-S TWO C METHODS. C NSW=MOD(NSW+1,3) C USE ONE OF TREANOR-S TWO METHODS. NSW=1+MOD(NSW,2) C DEBUG PRINT STATEMENTS........ IF(NDBUG)800,800,520 520 WRITE(KW,530)NSW,HT,RATFC,NDONE,K,ABERR,ERR(K),ARERR,RERR(K) 530 FORMAT(' NSW =',I3,' IN STIFF GEM.',5X,'HT =', * 1PE11.3,5X,'RATFC =', * E11.3,5X,'NDONE =',I6/5X,'K =',I3,5X, * 'ABERR =',E11.3,5X,'ERR(K) =',E11.3/5X, * 'ARERR =',E11.3,5X,'RERR(K) =',E11.3) WRITE(KW,540)XX,Y(K),YT(K),YU(K),DY 540 FORMAT(' XX =',1PE15.7,5X,'Y(K) =',E15.7/5X, * 'YT(K) =',E15.7,5X,'YU(K) =',E15.7,5X,'DY =',E15.7) C GO TO 800 C CHECK WHETHER STEP SIZE CAN BE INCREASED. C 550 IF(RATFC*TT-ABERR)580,580,560 560 IF(RATFC*TT-ARERR)580,580,570 570 INK=0 580 CONTINUE DO 590 K=1,MM 590 Y(K)=YU(K) IF(IQUIT)600,600,820 600 XX=XX+HH IF(INK)630,630,610 610 IF(ITCH)620,660,660 620 RATFC=UPFAC*RATFC GO TO 650 C RATFC=AMAX1(DNFAC*RATFC,RATST) 630 RATFC=DNFAC*RATFC IF(RATFC-RATST)640,650,650 640 RATFC=RATST 650 ITCH=0 GO TO 740 660 ITCH=+1 C INCREASE THE STEP SIZE. C HH=SIGN(AMIN1(ABS(SCALE*HH),HHMAX),HH) ABHH=SCALE*HH IF(ABHH)670,680,680 670 ABHH=-ABHH 680 IF(ABHH-HHMAX)700,700,690 690 ABHH=HHMAX 700 IF(HH)710,720,720 710 ABHH=-ABHH 720 HH=ABHH C H=ABS(HH) HT=HH IF(HH)730,740,740 730 HT=-HH C SEE WHETHER WE CAN FINISH IN ONE STEP. C IF(ABS(HH)-RCLOS*ABS(XF-X)) 740 ABHH=HH IF(ABHH)750,760,760 750 ABHH=-ABHH 760 XL=RCLOS*(FF-XX) IF(XL)770,780,780 770 XL=-XL 780 IF(ABHH-XL)800,790,790 790 HH=FF-XX IQUIT=1 GO TO 350 800 IF(NDONE-NMAX)350,350,810 C C INTEGRATION UNSUCCESSFUL... XF NOT REACHED IN NMAX STEPS. C 810 NDONE=-NDONE GO TO 830 C C INTEGRATION SUCCESSFUL... XF REACHED IN NOT MORE THAN NMAX STEPS. C 820 XX=FF XSAVE=XX 830 X=XX H=HT NUSED=NDONE RETURN C END STGEM. END SUBROUTINE HORDE(NE,TA,TF,H,YA,Y,FUNCT,NUSED) C C HORDE 2.1 A.N.S.I. STANDARD FORTRAN APRIL 1973 C C HORDE INTEGRATES THE SYSTEM OF ORDINARY DIFFERENTIAL C EQUATIONS FROM TA TO TF USING THE FIXED STEP SIZE H. C HORDE CALLS SUBROUTINE -ONCE- ONCE FOR EACH INTEGRATION STEP. C C J. P. CHANDLER, COMPUTER SCIENCE DEPT., OKLAHOMA STATE UNIVERSITY C C NE IS THE NUMBER OF DIFFERENTIAL EQUATIONS IN THE SYSTEM. C TA IS THE INITIAL VALUE OF THE INDEPENDENT VARIABLE. C TF IS THE FINAL VALUE OF THE INDEPENDENT VARIABLE. C H IS THE STEP SIZE. IT IS FIXED FOR EACH ENTRY TO HORDE. C YA IS THE VECTOR OF INITIAL VALUES OF THE DEPENDENT VARIABLE. C Y RETURNS THE VECTOR OF FINAL VALUES OF THE DEPENDENT VARIABLE. C FUNCT IS THE NAME OF THE SUBROUTINE THAT COMPUTES THE DERIVATIVES. C NUSED RETURNS THE NUMBER OF CALLS MADE TO FUNCT. C EXTERNAL FUNCT C DOUBLE PRECISION TA,TF,H,YA,Y,X,XF,HH,XA,ABHH,ABXD C DIMENSION YA(1),Y(1) C C INITIALIZE. JGO=0 NEQ=NE C IF(NEQ.LT.1.OR.NEQ.GT.NEQMAX) CALL Q8QERROR(0,15H NEQ TOO LARGE.) X=TA XF=TF NUT=0 C HH=SIGN(H,XF-X) IF(H)100,110,110 100 IF(XF-X)120,130,130 110 IF(XF-X)130,120,120 120 HH=H GO TO 140 130 HH=-H 140 DO 150 K=1,NEQ 150 Y(K)=YA(K) C CHECK FOR THE END OF THE INTEGRATION. 160 XA=X+HH ABHH=HH IF(ABHH)170,180,180 C IF(ABS(HH)-.99*ABS(XF-X)) 170 ABHH=-ABHH 180 ABXD=XF-X IF(ABXD)190,200,200 190 ABXD=-ABXD 200 IF(ABHH-.99*ABXD)220,210,210 210 XA=XF JGO=1 C COMPUTE Y(J) AT THE END OF THE INTERVAL. C 220 CALL ONCE (NEQ,X,XA,Y,Y,FUNCT,NA) NUT=NUT+NA C IF FINISHED, BUG OUT. IF NOT, DO IT AGAIN. IF(JGO)230,230,240 230 X=XA GO TO 160 240 NUSED=NUT TA=XF RETURN C END HORDE. END SUBROUTINE ONCE (NE,XX,XF,Y,YF,FUNCT,NA) C C ONCE 2.1 A.N.S.I. STANDARD FORTRAN APRIL 1973 C C INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL C EQUATIONS FROM XX TO XF IN ONE STEP. USES TREANOR-S METHOD FOR C STIFF EQUATIONS. C C J. P. CHANDLER, COMPUTER SCIENCE DEPT., OKLAHOMA STATE UNIVERSITY C C REFERENCE.... C. E. TREANOR, MATH. OF COMPUTATION 20 (1966) 39 C DOUBLE PRECISION XX,XF,Y,YF,YT,F,P,FF,FT,YY,EXMIN,DEXP,ZEXP, * Q,RHALF,RUNIT,RTWO,RTHRE,RFOUR,THRSH,PHLIM,X,H,DENOM, * PH,FZ,DYZ,DYB C DIMENSION Y(1),YF(1) DIMENSION YT(40,4),F(40,4),P(40),FF(40,3) DIMENSION FT(40),YY(40) C C NSW -- IF .EQ.0, USE THE CLASSICAL FOURTH ORDER RUNGE-KUTTA METHOD C -- IF .EQ.1, USE TREANOR-S METHOD (FIRST METHOD) C -- IF .EQ.2, USE TREANOR-S SECOND METHOD C COMMON /MINK/ NSW C C SET LIBRARY FUNCTION FOR SINGLE PRECISION C OR DOUBLE PRECISION. ZEXP(Q)=DEXP(Q) C C KBIGP=1 IMPLEMENTS THE METHOD AS GIVEN BY LAPIDUS AND SEINFELD. C THIS USES THE MAXIMUM P VALUE IN ALL EQUATIONS. IN SOME CASES C AT LEAST, THIS IS MUCH SLOWER AND NO MORE ACCURATE. C IT IS NOT RECOMMENDED.......... C KBIGP=0 C C PHLIM ... VALUE OF P*H ABOVE WHICH C TREANOR-S MODIFIED METHOD IS USED C REGARDLESS OF THE VALUE OF NSW PHLIM=10.0 PHLIM=2. C THRSH ... VALUE OF P*H BELOW WHICH C THE RUNGE-KUTTA METHOD IS USED THRSH=.01 C C EXMIN IS THE MOST NEGATIVE NUMBER SUCH THAT EXP(EXMIN) DOES NOT C UNDERFLOW. ON THE IBM 360, THIS IS -179.4 . C EXMIN=-179.4 C RZERO=0.0 RUNIT=1.0 RTWO=2. RHALF=RUNIT/RTWO RTHRE=3. RFOUR=4. C 110 NEQ=NE C RETURN THE NUMBER OF DERIVATIVE C EVALUATIONS PER STEP. NA=4 X=XX H=XF-X IF(H)140,120,140 120 DO 130 J=1,NEQ 130 YF(J)=Y(J) C C COMPUTE THE Y-S. F IS TREANOR-S LOWER-CASE F. C FF IS TREANOR-S CAPITAL F. C GO TO 390 140 CALL FUNCT (NEQ,X,Y,FT) DO 150 J=1,NEQ F(J,1)=FT(J) YT(J,1)=Y(J) YT(J,2)=YT(J,1)+F(J,1)*H/RTWO 150 YY(J)=YT(J,2) CALL FUNCT (NEQ,X+H/RTWO,YY,FT) DO 160 J=1,NEQ F(J,2)=FT(J) YT(J,3)=YT(J,1)+F(J,2)*H/RTWO 160 YY(J)=YT(J,3) CALL FUNCT (NEQ,X+H/RTWO,YY,FT) C C MOVE THE F-S. DO 170 J=1,NEQ 170 F(J,3)=FT(J) C TEST NSW TO SEE WHICH METHOD TO USE. IF(NSW)260,260,180 C COMPUTE THE P-S. 180 DO 210 J=1,NEQ DENOM=YT(J,3)-YT(J,2) IF(DENOM)190,200,190 C 190 P(J)=-(F(J,3)-F(J,2))/DENOM C C P(J)=AMAX1(P(J),0.0) IF(P(J))200,210,210 200 P(J)=RZERO 210 CONTINUE C BRANCH ON KBIGP. IF(KBIGP)220,280,220 C FIND THE GREATEST P(J). 220 AMAXP=P(1) DO 240 J=1,NEQ IF(P(J)-AMAXP)240,240,230 230 AMAXP=P(J) 240 CONTINUE C KBIGP=1. SET ALL P(J) = MAXIMUM P(J). DO 250 J=1,NEQ 250 P(J)=AMAXP GO TO 280 C NSW=0. SET ALL P(J)=0. 260 DO 270 J=1,NEQ 270 P(J)=RZERO C TAKE THE STEP. 280 DO 370 J=1,NEQ PH=P(J)*H DENOM=-PH IF(PH-THRSH)290,290,300 C USE ORDINARY RUNGE-KUTTA. 290 P(J)=RZERO FF(J,1)=RUNIT FF(J,2)=RHALF FF(J,3)=RHALF/RTHRE PH=RZERO GO TO 330 C COMPUTE THE CAPITAL F-S. 300 FZ=RZERO IF(DENOM-EXMIN)320,320,310 310 FZ=ZEXP(DENOM) 320 FF(J,1)=(FZ-RUNIT)/DENOM FF(J,2)=(FF(J,1)-RUNIT)/DENOM FF(J,3)=(FF(J,2)-RHALF)/DENOM C C CHOOSE WHICH OF TREANOR-S METHODS TO USE. 330 IF(NSW-1)350,340,360 340 IF(PH-PHLIM)350,350,360 C C USE TREANOR-S FIRST (UNMODIFIED) METHOD. C 350 YT(J,4)=YT(J,1)+H*F(J,3) GO TO 370 C USE THE MODIFIED METHOD (P. 42). C 360 YT(J,4)=YT(J,1)+H*(RTWO*F(J,3)*FF(J,2)+ * F(J,1)*(FF(J,1)-RTWO*FF(J,2))+ * F(J,2)*P(J)*H*FF(J,2)) C C COMPUTE THE YF-S. 370 YY(J)=YT(J,4) CALL FUNCT (NEQ,X+H,YY,FT) DO 380 J=1,NEQ F(J,4)=FT(J) DYA=H*(F(J,1)*FF(J,1)+(-RTHRE*(F(J,1)+P(J)*YT(J,1))+ * RTWO*(F(J,2)+P(J)*YT(J,2))+RTWO*(F(J,3)+P(J)*YT(J,3))- * (F(J,4)+P(J)*YT(J,4)))*FF(J,2)) DYB=H*RFOUR*(F(J,1)+P(J)*YT(J,1)-(F(J,2)+P(J)*YT(J,2))- * (F(J,3)+P(J)*YT(J,3))+(F(J,4)+P(J)*YT(J,4)))*FF(J,3) 380 YF(J)=Y(J)+DYA+DYB 390 RETURN C END ONCE. END SUBROUTINE DIFFUN (NEQ,T,Y,DYDT) C C THIS IS THE SEVEN-PARAMETER DIFFUN USED, E.G., TO GENERATE TABLE 4 C IN THE 1972 PAPER BY CHANDLER, HILL, AND SPIVEY. C C YDE(5) -- A C YDE(6) -- E C YDE(7) -- B C CUSER(1),CUSER(3),CUSER(5),... -- AT(1),AT(2),AT(3),... C CUSER(2),CUSER(4),CUSER(6),... -- BT(1),BT(2),BT(3),... C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION T,Y,DYDT C DIMENSION Y(4),DYDT(4) C COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C JS=2*JSET-1 YDE(5)=CUSER(JS)-Y(1)-Y(2)-Y(3)-Y(4) YDE(6)=CUSER(NCUSR)-Y(1)-Y(3)-Y(4) YDE(7)=CUSER(JS+1)-Y(1)-Y(2)-Y(4) C DYDT(1)=BD(5)*Y(4)-(BD(6)+BD(7))*Y(1) DYDT(2)=BD(7)*Y(1) DYDT(3)=BD(1)*YDE(6)*YDE(5)+BD(4)*Y(4)-Y(3)*(BD(2)+BD(3)*YDE(7)) DYDT(4)=BD(3)*Y(3)*YDE(7)+BD(6)*Y(1)-(BD(4)+BD(5))*Y(4) FITJ=Y(1)+Y(2) RETURN C END DIFFUN 7K. END $ENTRY DATA FOR SEVEN-PARAMETER STIFF PROBLEM CRICF TEST PROBLEM 4 19 .002 .000606 .00021 .004 .001223 .008 .002169 .012 .002883 .016 .003622 .020 .004142 .022 .004292 .028 .004975 .040 .005875 .060 .006562 .068 .006958 .086 .007476 .104 .008075 .112 .008225 .132 .008840 .148 .009184 .160 .009622 .168 .009772 .200 .010625 19 .004 .001078 .012 .002523 .020 .003451 .028 .004089 .036 .005200 .044 .005100 .052 .005335 .060 .005773 .068 .006520 .076 .006369 .084 .006462 .092 .006663 .100 .007117 .108 .006596 .116 .007197 .124 .007446 .132 .007668 .140 .007501 .144 .008498 19 .004 .000426 .016 .002081 .024 .002226 .028 .002880 .036 .002521 .044 .003360 .052 .004073 .060 .004211 .068 .004027 .076 .004162 .084 .004422 .092 .004437 .100 .004564 .108 .004884 .116 .004953 .124 .004849 .132 .005392 .140 .005511 .148 .005392 19 .004 .000061 .012 .001083 .020 .001446 .028 .002211 .036 .002179 .044 .002404 .052 .002688 .060 .002731 .068 .002787 .076 .002729 .084 .003180 .092 .003053 .100 .003182 .108 .003466 .116 .003561 .124 .003468 .132 .003916 .140 .003785 .148 .003815 1 2 2000 9 .00001 .0001 10. 1. .0001 . 1.25 25. .08 25. 1.25 5. .08 5. .018 7 4 4 2 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 19813.51 1146.206 337.9353 2633.707 25.00938 36.19896 4.803050 100. 10. 3. 10. 1. 1. 1. 1000000. 100000. 30000. 100000. 15000. 10000. 100. 2 SUBROUTINE FOFX (XX,BSOLV,NACTV,KPTIN,NEQ,NUSED) C C FOFX COMPUTES AND RETURNS THE VALUE OF THE FITTED FUNCTION AT XX. C THIS VERSION IS USED WITH MODELS THAT ARE NOT THE SOLUTIONS OF C DIFFERENTIAL EQUATIONS. C DOUBLE PRECISION YDE,FITJ,BD,CUSER C DIMENSION BSOLV(1) C COMMON /CRICK/ X(100),DFRAC,ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,KALC,NMAX,NPRNT,KW, * LOGB(30),MSK(30),NPAR,NFUPR COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C C UNPACK THE PARAMETER VECTOR, EXPONENTIATING IT IF NECESSARY. C NUSED=0 J=0 DO 120 K=1,NPAR IF(MSK(K))120,100,120 100 J=J+1 BD(K)=BSOLV(J) IF(LOGB(K))110,120,110 110 BD(K)=10.0**BSOLV(J) 120 CONTINUE C C CALL CALCF TO COMPUTE THE FITTED VALUE. IT IS RETURNED THROUGH C THE VARIABLE FITJ IN COMMON/CDFUN/. IF ANALYTIC DERIVATIVES ARE C USED (I.E. IF DERIV2 IS USED), THEY ARE RETURNED THROUGH DFDB. C CALL CALCF (XX) RETURN C END FOFXC. END SUBROUTINE CALCF (XX) C C CALCF FOR FITTING B(1)+B(2)*EXP(B(3)*XX) TO DATA. C NOTE THAT -CUSER- IS NOT USED HERE. C DOUBLE PRECISION YDE,FITJ,BD,CUSER DOUBLE PRECISION DQ C COMMON /CDFUN/ YDE(40),FITJ,BD(30),CUSER(60),DFDB(30),JSET,NCUSR C C AVOID MIXED MODE ARITHMETIC.... Q=BD(3) DQ=EXP(Q*XX) FITJ=BD(1)+BD(2)*DQ RETURN C END CALCF. END CRICF TEST PROBLEM 4 .... LEAST SQUARE FIT OF B(1)+B(2)*EXP(B(3)*T) T $ENTRY 1 12 0. 167. 13. 1. 109. 11. 2. 91. 9. 3. 60. 8. 4. 56. 7. 5. 47. 7. 6. 32. 6. 7. 31. 6. 8. 28. 5. 9. 22. 5. 10. 25. 5. 11. 19. 4.5 1 .0001 .0 3 2 0 0 0 0 0 0 0. 100. -.5 0. 0. 0. 0. 0. 0. 2 //