//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....