//U11701AA JOB (*),TIME=(,59),CLASS=2,MSGCLASS=Y, // NOTIFY=* /*JOBPARM ROOM=H,FORMS=9021 // EXEC VSF2CLG, // REGION.FORT=2000K, // PARM.FORT='XREF,LANGLVL(77)', // REGION.GO=2000K, // PARM.GO='SIZE=2000K' //FORT.SYSIN DD * C C CRICFM3.CNTL CURRENT VERSION OF CRICFM -- JPC C SUBROUTINE DIFFUN (NEQ,T,Y,DYDT) IMPLICIT REAL*8(A-H,O-Z) 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(3),DYDT(3) C COMMON /CDFNC/ YDE(40),BD(20),CUSER(60),FITJ,NCUSR,JSET 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 C C END DIFFUN (5K MODEL) C END C PROGRAM CRICFM (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT) IMPLICIT REAL*8(A-H,O-Z) C C CRICFM 1.1 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 CRICFM IS A SHORT VERSION OF CRICF FOR USE ON SMALL COMPUTERS. 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 EXTERNAL CRFUN CHARACTER*1 KHINT,KHTEM,KHDF,KHTM2 DOUBLE PRECISION YDE,FIT,BD,CUSER,FITJ,X,DHSAV DOUBLE PRECISION Y,YSIG DIMENSION KHTEM(10),KHTM2(2) COMMON /CRICM/ T(300),DHSAV(10),ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,NMAX,NPRNT,LOGB(20),NFUPR,NACTV, * NVMAX COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20),DELMN(20), * ERR(20,21),FOBJ,NV,NTRAC,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW COMMON /NLLS4/ FLAMB,FNU,RELDF,RELMN,METHD,KALCP,KORDF,MAXIT, * LEQU,MXSUB,MXUPD COMMON /CHINT/ KHINT(10),KHDF(2) DATA KHTEM /'0','1','2','3','4','5','6','7','8','9'/ DATA KHTM2 /'D','F'/ C DO 10 J=1,10 10 KHINT(J)=KHTEM(J) KHDF(1)=KHTM2(1) KHDF(2)=KHTM2(2) C 20 CALL CRFMA CALL CRCON CALL STEPT(CRFUN) CALL CRFMB GO TO 20 C C END CRICFM C END SUBROUTINE CRFMA IMPLICIT REAL*8(A-H,O-Z) C C THIS SUBROUTINE READS DATA FOR THE CRICFM PACKAGE, AND INITIALIZES. C DOUBLE PRECISION YDE,FIT,BD,CUSER,FITJ,X,DHSAV,ZDABS,DARG DOUBLE PRECISION Y,YSIG C C THE DIMENSIONS OF THE MATRICES ARE C Y(NPTS),YSIG(NPTS),FIT(NPTS),T(NPTS), C LOGB(NV),MASK(NV),X(NV),XMAX(NV),XMIN(NV),DELMN(NV),BD(NV), C YDE(MAX(NEQFU,NEQDE)),CUSER(NCUSR),KTITL(80) C DIMENSION KTITL(80) C COMMON /CRICM/ T(300),DHSAV(10),ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,NMAX,NPRNT,LOGB(20),NFUPR,NACTV, X NVMAX COMMON /CDFNC/ YDE(40),BD(20),CUSER(60),FITJ,NCUSR,JSET COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20),DELMN(20), * ERR(20,21),FOBJ,NV,NTRAC,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW COMMON /NLLS4/ FLAMB,FNU,RELDF,RELMN,METHD,KALCP,KORDF,MAXIT, * LEQU,MXSUB,MXUPD COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS C ZDABS(DARG)=DABS(DARG) ZLOGT(ARG)=DLOG(ARG)/DLOG(10.0D0) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INITIALIZE. C C THE STEP SIZE FOR INTEGRATING OVER THE FIRST INTERVAL, FOR EACH C JSET (EACH SET OF DATA), IS SAVED IN DHSAV. C NHSAV=10 DO 10 J=1,NHSAV 10 DHSAV(J)=0.0 C CALL STSET NTRAC=1 C KR ... LOGICAL UNIT NUMBER OF THE READER KR=5 C NVMAX ... MAXIMUM NUMBER OF PARAMETERS, C AND THE DIMENSION OF B, XMIN, ETC. NVMAX=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 RELDF ... RELATIVE STEP SIZE FOR DERIV RELDF=.002 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,20)KTITL 20 FORMAT(80A1) WRITE(KW,30)KTITL 30 FORMAT('1 CRICF PROBLEM......'///10X,80A1) READ(KR,40)NSETS 40 FORMAT(I5) C CHECK FOR DEFAULT. IF(NSETS)90,90,60 90 STOP 60 WRITE(KW,70)NSETS 70 FORMAT(/////1X,I5,' SET(S) OF DATA') NPTS=0 YSDEF=1.0 DO 160 JSET=1,NSETS READ(KR,40)MPTS WRITE(KW,80)MPTS 80 FORMAT(/' SET OF',I4,' DATA POINTS....'/' ') IF(MPTS)90,90,100 100 JPTS(JSET)=MPTS DO 140 J=1,MPTS NPTS=NPTS+1 READ(KR,110)T(NPTS),Y(NPTS),YS 110 FORMAT(3E10.5) IF(YS)130,130,120 120 YSDEF=YS 130 YS=YSDEF YSIG(NPTS)=YS 140 WRITE(KW,150)NPTS,T(NPTS),Y(NPTS),YS 150 FORMAT(6X,'J =',I4,6X,'T =',1PE15.7,6X,'Y =',E15.7,6X, * 'YSIG =',E11.3) 160 CONTINUE C C READ IN THE OTHER INPUT QUANTITIES, SET DEFAULT VALUES, AND PRINT. C READ(KR,170)NPRNT,MAXIT,NMAX,NCUSR,ER,RER,FNU,FLAMB,RLCNV 170 FORMAT(4I5,4E10.5,10X,E10.5) IF(MAXIT)180,180,190 180 MAXIT=100 190 IF(NMAX)200,200,210 200 NMAX=5000 210 IF(ER)240,220,230 220 IF(RER)240,240,250 230 IF(RER)240,250,250 240 ER=0.0 RER=.0001 250 IF(FNU)260,260,270 260 FNU=10.0 270 IF(FLAMB)280,280,290 280 FLAMB=.1 290 WRITE(KW,300)NPRNT,MAXIT,NMAX,NCUSR,ER,RER,FNU,FLAMB,RLCNV 300 FORMAT(////'1'////' NPRNT =',I2,5X,'MAXIT =',I4,5X, * 'NMAX =',I6,5X,'NCUSR =',I3///5X,'ER =',1PG13.5,5X, * 'RER =',G13.5,5X,'FNU =',G13.5//5X,'FLAMB =',G13.5, * 5X,'RLCNV =',G13.5) IF(NCUSR)410,350,310 310 IF(NCUSR-MXNCU)320,320,410 320 READ(KR,330)(CUSER(J),J=1,NCUSR) 330 FORMAT(8D10.5) WRITE(KW,340)(J,CUSER(J),J=1,NCUSR) 340 FORMAT(/////' USER CONSTANTS....'// * (5X,'CUSER(',I2,') =',1PG15.7)) 350 READ(KR,360)NV,NEQFU,NEQDE,KORDF,MCSIM,KVSUP,SPFAC,NFUPR 360 FORMAT(6I5,E10.5,I5) IF(KORDF-2)370,380,370 370 KORDF=1 380 WRITE(KW,390)NV,NEQFU,NEQDE,KORDF,MCSIM,KVSUP,SPFAC,NFUPR 390 FORMAT(////' NV =',I3,5X,'NEQFU =',I3,5X,'NEQDE =',I3// * 5X,'KORDF =',I2,5X,'MCSIM =',I5,5X,'KVSUP =',I3,5X// * 5X,'SPFAC =',1PG13.5,5X,'NFUPR =',I4) IF(NV)410,410,400 400 IF(NV-NVMAX)430,430,410 410 WRITE(KW,420)NV,NVMAX,NCUSR,MXNCU 420 FORMAT(////' ERROR IN INPUT. NV =',I11,5X,'NVMAX =', * I11/5X,'NCUSR =',I11,5X,'MXNCU =',I11///' ') 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 430 CONTINUE READ(KR,440)(MASK(J),J=1,NV) 440 FORMAT(8I10) IF(MASK(1)-1)460,460,450 450 STOP 460 WRITE(KW,470)(MASK(J),J=1,NV) 470 FORMAT(////' MASK(J)....',I9,4I12/(9X,5I12)) READ(KR,440)(LOGB(J),J=1,NV) DO 490 J=1,NV IF(MASK(J))480,490,480 480 LOGB(J)=0 490 CONTINUE WRITE(KW,500)(LOGB(J),J=1,NV) 500 FORMAT(/' LOGB(J)....',I9,4I12/(9X,5I12)) READ(KR,510)(X(J),J=1,NV) 510 FORMAT(8E10.5) READ(KR,510)(XMIN(J),J=1,NV) READ(KR,510)(XMAX(J),J=1,NV) C C DEFAULT XMIN AND XMAX. DO 530 J=1,NV IF(XMAX(J)-XMIN(J))520,520,530 520 XMAX(J)=HUGE XMIN(J)=-HUGE 530 CONTINUE WRITE(KW,540)(X(J),J=1,NV) 540 FORMAT(/' X(J)....... ',1PE12.4,4E12.4/(10X,5E12.4)) WRITE(KW,550)(XMIN(J),J=1,NV) 550 FORMAT(/' XMIN(J).... ',1PE12.4,4E12.4/(10X,5E12.4)) WRITE(KW,560)(XMAX(J),J=1,NV) 560 FORMAT(/' XMAX(J).... ',1PE12.4,4E12.4/(10X,5E12.4)) C C MAP X, XMAX, AND XMIN INTO LOGS (BASE 10) C IF REQUESTED. DO 680 J=1,NV IF(LOGB(J))570,640,570 570 IF(XMIN(J))580,580,590 580 XMIN(J)=1.0/HUGE 590 XMIN(J)=ZLOGT(XMIN(J)) IF(XMAX(J))600,600,610 600 XMAX(J)=HUGE 610 XMAX(J)=ZLOGT(XMAX(J)) IF(X(J))620,620,630 620 X(J)=1.0 630 X(J)=ZLOGT(X(J)) 640 IF(X(J)-XMAX(J))660,660,650 650 X(J)=XMAX(J) 660 IF(X(J)-XMIN(J))670,680,680 670 X(J)=XMIN(J) 680 BD(J)=X(J) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C WRITE(KW,690) 690 FORMAT(////'1'///' FIT THE DATA..........'//' ') C WRITE(KW,700)(J,X(J),XMIN(J),XMAX(J),MASK(J),LOGB(J),J=1,NV) 700 FORMAT(////' PARAMETER',7X,'STARTING',9X,'LOWER',11X, * 'UPPER',13X,'MASK',7X,'LOGB'/' NUMBER',10X,'VALUE', * 2(11X,'LIMIT')//(1X,I6,7X,1PE16.7,2E16.7,2I10)) C NACTV=0 DO 720 J=1,NV IF(MASK(J))705,704,705 704 NACTV=NACTV+1 705 DELMN(J)=ZDABS(RLCNV*X(J)) IF(DELMN(J))710,710,720 710 DELMN(J)=RLCNV 720 CONTINUE RETURN END SUBROUTINE CRCON IMPLICIT REAL*8(A-H,O-Z) C C CRCON 1.2 A.N.S.I. STANDARD FORTRAN MARCH 1980 C C THIS SUBROUTINE CONTROLS THE OPERATION OF THE CRICFM PACKAGE. C THIS VERSION OF CRCON IS A DUMMY, FOR BATCH RUNS ONLY. C ANOTHER VERSION OF CRCON IS USED FOR INTERACTIVE OPERATION. C RETURN END SUBROUTINE CRFMB IMPLICIT REAL*8(A-H,O-Z) C C PRINTS ERROR BARS FOR CRICFM FITS. C DOUBLE PRECISION X,RTEN COMMON /CRICM/ T(300),DHSAV(10),ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,NMAX,NPRNT,LOGB(20),NFUPR,NACTV, * NVMAX COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20),DELMN(20), * ERR(20,21),FOBJ,NV,NTRAC,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C ZSQRT(ARG)=DSQRT(ARG) ZABS(ARG)=DABS(ARG) C RTEN=10.0 WRITE(KW,10) 10 FORMAT(/////' EXPONENTIATED RESULTS FROM THE FIT .....'/// * 83X,'ONE-STD.-ERROR'/8X,'PARAMETER',6X,'MASK',3X, * 'LOGB',10X,'XMIN',11X,'XMAX',22X,'ENDPOINTS'/' ') DO 70 K=1,NV IF(MASK(K))20,40,20 20 WRITE(KW,30)X(K),MASK(K),LOGB(K) 30 FORMAT(5X,1PE14.7,2I7,5X,2E15.5,5X,2E15.5) GO TO 70 40 BERR=ZSQRT(ZABS(ERR(K,K))) BB=X(K) BBPL=BB+BERR BBMI=BB-BERR BXMIN=XMIN(K) BXMAX=XMAX(K) IF(LOGB(K))50,60,50 50 BB=RTEN**BB BXMAX=RTEN**BXMAX BXMIN=RTEN**BXMIN BBPL=RTEN**BBPL BBMI=RTEN**BBMI 60 WRITE(KW,30)BB,MASK(K),LOGB(K),BXMIN,BXMAX,BBMI,BBPL 70 CONTINUE RETURN END SUBROUTINE FRPIN (KH,NKH,JX) IMPLICIT REAL*8(A-H,O-Z) C C FRPIN -- FREE FORM INPUT OF AN INTEGER C C ON ENTRY, THE ARRAY (KH(J),J=1,NKH) CONTAINS A STRING OF CHARACTERS C IN A1 FORM. THIS SUBROUTINE CONSTRUCTS THE INTEGER JX FROM THE C NUMERIC CHARACTERS, IF ANY, IN THIS STRING. C C BEFORE CALLING THIS SUBROUTINE, THE NUMERIC CHARACTERS C 1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9 C MUST HAVE BEEN STORED IN THE ARRAY KHINT(*). (THIS CANNOT BE DONE C IN THIS SUBROUTINE USING 1966 A.N.S.I. STANDARD FORTRAN.) C CHARACTER*1 KH,KHINT,KHDF DIMENSION KH(2) COMMON /CHINT/ KHINT(10),KHDF(2) C JX=0 DO 30 J=1,NKH DO 20 K=1,10 IF(KH(J).NE.KHINT(K)) GO TO 20 10 JX=10*JX+(K-1) GO TO 30 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE CRFUN (FITM) IMPLICIT REAL*8(A-H,O-Z) C C CRFUN 1.0 JULY 1979 C C CRFUN COMPUTES AND RETURNS THE VALUES OF THE FITTED FUNCTION, FITM(J) C THIS VERSION IS USED WITH MODELS THAT ARE THE SOLUTIONS OF C SYSTEMS OF DIFFERENTIAL EQUATIONS. THE EQUATIONS ARE INTEGRATED C BY SUBROUTINE GEM. C GEM CALLS SUBROUTINE DIFFUN TO COMPUTE THE VALUES OF THE DERIVATIVES. C CRFUN CALLS SUBROUTINE DINIT TO INITIALIZE TZERO AND YDE(J) FOR C THE INTEGRATION TO THE FIRST TT OF EACH DATA SET. C EXTERNAL DIFFUN C DOUBLE PRECISION YDE,FITJ,BD,CUSER,FITM,X,RTEN DOUBLE PRECISION ABERR,RERR,DTT,DTOLD,DH,DHMAX,DHSAV,DYDT C DIMENSION FITM(1) DIMENSION ABERR(40),RERR(40),DYDT(40) C COMMON /CRICM/ T(300),DHSAV(10),ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,NMAX,NPRNT,LOGB(20),NFUPR,NACTV, X NVMAX COMMON /CDFNC/ YDE(40),BD(20),CUSER(60),FITJ,NCUSR,JSET COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20),DELMN(20), * ERR(20,21),FOBJ,NV,NTRAC,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C RTEN=10.D0 NEQ=NEQFU IF(NEQDE-NEQ)20,20,10 10 NEQ=NEQDE C C MOVE THE PARAMETERS, EXPONENTIATING THEM IF NECESSARY. C 20 DO 40 K=1,NV BD(K)=X(K) IF(LOGB(K))30,40,30 30 BD(K)=RTEN**BD(K) 40 CONTINUE C DHMAX=1.0E30 DO 50 J=1,NEQ ABERR(J)=ER 50 RERR(J)=RER C IF(NPRNT-2)80,60,60 60 WRITE(KW,70) 70 FORMAT(//' INTEGRATED SOLUTION OF THE SYSTEM OF', * ' DIFFERENTIAL EQUATIONS....'/' ') C C LOOP OVER THE SETS OF DATA. 80 NSTPS=0 JPT=0 DO 270 JS=1,NSETS JSET=JS MPTS=JPTS(JS) C CALL DINIT (NV,BD,JSET,CUSER,NCUSR,NEQ,TZERO,YDE) DTOLD=TZERO C USE THE FIRST STEPSIZE FROM THE PREVIOUS C INTEGRATION FOR THIS SET OF DATA. DH=DHSAV(JSET) JSTEP=0 C LOOP OVER THE POINTS IN THE JS-TH SET. DO 260 KPTIN=1,MPTS JPT=JPT+1 TT=T(JPT) DTT=TT IF(TZERO-TT)90,110,90 90 IF(DH)140,100,140 100 DH=DTT-DTOLD GO TO 140 C 110 DH=0.0 NUSED=1 IF(NPRNT-2)170,120,120 120 WRITE(KW,130) 130 FORMAT(/' ') GO TO 170 C C CALL GEM TO INTEGRATE THE SYSTEM OF DIFFERENTIAL EQUATIONS TO XX, C THE ABSCISSA OF THE NEXT DATA POINT. C 140 CALL GEM (NEQ,DTOLD,DTT,DH,DHMAX,YDE,ABERR,RERR, * DIFFUN,NMAX,NUSED) IF(DHSAV(JSET))160,150,160 150 DHSAV(JSET)=DH 160 JSTEP=JSTEP+1 C CALL DIFFUN TO SET FITJ FOR THE FINAL C VALUES OF YDE. C 170 CALL DIFFUN (NEQ,DTOLD,YDE,DYDT) IF(NPRNT-2)180,190,190 180 IF(NUSED)190,250,250 190 NYPR=NEQ IF(NFUPR-NYPR)210,210,200 200 NYPR=NFUPR 210 WRITE(KW,220)JSTEP,TT,DH,FITJ,NUSED,(YDE(J),J=1,NYPR) 220 FORMAT(' STEP',I5,5X,'TT =',1PE15.7,5X,'HH =',E13.5/5X, * 'FITJ =',E15.7,5X,'NUSED =',I6,5X,'YDE(J)....'/ * (5X,6E15.7)) IF(NUSED)230,230,250 230 WRITE(KW,240)DTOLD,DTT,DH,NMAX,NUSED 240 FORMAT(//' GEM FAILURE.'/5X,'DTOLD =',1PG15.7,5X,'DTT =', * G15.7,5X,'DH =',G15.7,5X,'NMAX =',I8,5X,'NUSED =',I8) STOP 250 FITM(JPT)=FITJ 260 NSTPS=NSTPS+NUSED 270 CONTINUE C IF(NPRNT)300,300,280 280 WRITE(KW,290)NSTPS 290 FORMAT(' CRFUN. NSTPS =',I6) 300 RETURN END SUBROUTINE STEPT (FUNK) C C INTERFACE TO MAKE MARQ LOOK LIKE STEPT. C C J. P. CHANDLER, DEPARTMENT OF COMPUTER SCIENCE, C OKLAHOMA STATE UNIVERSITY C C TO USE THIS ROUTINE, SET THE VALUES OF LPDMA AND LPDMB, AND SET C THE DIMENSIONS OF THE ARRAYS P, FITSV, FIT, Y, AND YSIG. C NORMALLY THE DIMENSIONS ARE.... C P(LPDMA,LPDMB), FITSV(LPDMA), FIT(LPDMA), Y(LPDMA), C YSIG(LPDMA), WHERE LPDMA IS .GE. NPTS AND LPDMB IS .GE. NV. C C COMMON/CDAT/ DOES NOT APPEAR IN ANY ROUTINE OF THE MARQ PACKAGE C OTHER THAN THIS ONE, SO THIS IS THE ONLY ROUTINE THAT MUST BE C RECOMPILED WHEN THE DIMENSIONS OF THE ARRAYS ARE CHANGED. C C THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME COMPILERS C (WATFIV, FOR EXAMPLE) AND FORBIDDEN BY OTHERS (MODCOMP II). EXTERNAL FUNK C DOUBLE PRECISION FIT,FITSV ,Y,YSIG,P C INTEGER LPDMA,LPDMB,NPTS C DIMENSION P(300,20),FITSV(300) C COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS C LPDMA=300 LPDMB=20 C CALL MARQ (FUNK,Y,YSIG,NPTS,FIT,FITSV,P,LPDMA,LPDMB) C RETURN END SUBROUTINE MARQ (FUNK,Y,YSIG,NPTS,FIT,FITSV,P, * LPDMA,LPDMB) C C MARQ 3.2 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C COPYRIGHT (C) 1981, 1991 J. P. CHANDLER C C J. P. CHANDLER AND LEON W. JACKSON, C COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY, C STILLWATER, OKLAHOMA 74078 C C MARQ PERFORMS A NONLINEAR LEAST SQUARES FIT OF A C USER-SUPPLIED FUNCTION TO A GIVEN SET OF DATA, USING C MARQUARDT'S METHOD, THE GAUSS-NEWTON METHOD, OR THE C MODIFIED GAUSS-NEWTON METHOD. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INPUT QUANTITIES..... FUNK,X(*),XMAX(*),XMIN(*),DELMIN(*), C NV,NTRACE,MASK(*),MATRX,NFMAX, C NFLAT,KW, C Y(*),YSIG(*),NPTS,LPCOL, C FLAMBD,FNU,RELDIF,METHD,KALCP, C KORDIF,MAXIT,LEQU,MAXSUB C C OUTPUT QUANTITIES.... X(*),FOBJ,ERR(*,*),KFLAG,FIT(*) C C SCRATCH VECTOR....... DELTX(*) C C UNUSED QUANTITIES (INCLUDED FOR COMPATIBILITY WITH C STEPIT COMMON).... JVARY,NXTRA,NOREP,KERFL C C C FUNK -- THE NAME OF THE SUBROUTINE CALLED TO C OBTAIN THE FITTED VALUES IF KALCP=0 C C X(JX) -- THE JX-TH PARAMETER C C XMAX(JX) -- THE UPPER LIMIT FOR X(JX) C C XMIN(JX) -- THE LOWER LIMIT FOR X(JX) C C DELMIN(JX) -- THE CONVERGENCE TOLERANCE FOR X(JX) C C FOBJ -- RETURNS THE FINAL VALUE OF PHI, THE C WEIGHTED SUM OF SQUARES (PHI IS BEING C MINIMIZED AS A FUNCTION OF THE X(JX)) C C NV -- THE NUMBER OF PARAMETERS C C NTRACE -- USER PRINT CONTROL ... C C =-4 FOR NO OUTPUT EXCEPT FATAL ERROR C MESSAGE, C C =-3 FOR NO OUTPUT EXCEPT ERROR C MESSAGES, C C =-2 FOR NO OUTPUT EXCEPT DIAGNOSTIC C MESSAGES, C C =-1 FOR STANDARD OUTPUT EXCEPT NO FINAL C FIT VALUES, C C = 0 FOR STANDARD OUTPUT, C C = 1 TO PRINT ALSO THE RESULTS OF EACH C ITERATION, C C = 2 TO PRINT ALSO THE COEFFICIENT C MATRIX, QSAV(*,*), C C = 3 TO PRINT ALSO THE JACOBIAN C MATRIX, P(*,*) C C MATRX -- NONZERO IF STATISTICAL ERRORS ARE TO BE C COMPUTED C C MASK(JX) -- NONZERO IF X(JX) IS TO BE HELD FIXED C C NFMAX -- THE MAXIMUM NUMBER OF FUNCTION C COMPUTATIONS C C NFLAT -- NONZERO IF THE SEARCH IS TO TERMINATE WHEN C TWO ITERATIONS GIVE IDENTICAL VALUES C OF PHI. C THE RECOMMENDED VALUE OF NFLAT IS C USUALLY NFLAT=1 . C C KW -- THE LOGICAL UNIT NUMBER OF THE PRINTER C C ERR(JX,KX) -- RETURNS THE ERROR MATRIX C C KFLAG -- RETURNED .GT. ZERO FOR A NORMAL EXIT, C RETURNED .LT. ZERO FOR AN ABNORMAL EXIT C C Y(JPT) -- THE JPT-TH DATA ORDINATE C C YSIG(JPT) -- THE EXPECTED ERROR IN Y(JPT) C C NPTS -- THE NUMBER OF DATA OBSERVATIONS C C FIT(JPT) -- THE JPT-TH FITTED VALUE C C FITSV -- A SCRATCH VECTOR OF NPTS VALUES C C P(JPT,JX) -- THE FIRST PARTIAL DERIVATIVE C OF FIT(JPT) WITH RESPECT TO X(JX) C C LPDMA -- THE FIRST DIMENSION OF THE ARRAY C CONTAINING P C (LPCOL MUST BE .GE. NPTS IF KALCP IS C .GE. ZERO) C C LPDMB -- THE SECOND DIMENSION OF THE ARRAY C CONTAINING P (LPDMB MUST BE .GE. NV) C C FLAMBD -- MARQUARDT'S LAMBDA, THE RELATIVE AMOUNT BY C WHICH THE DIAGONAL COEFFICIENTS OF THE C NORMAL EQUATIONS ARE AUGMENTED. C MARQUARDT RECOMMENDED FLAMBD=0.01 BUT C SOME USERS PREFER FLAMBD=1.0, SLOWER C BUT SAFER. C C FNU -- MARQUARDT'S NU, THE FACTOR BY WHICH FLAMBD C IS CHANGED. C MARQUARDT RECOMMENDED FNU=10.0 . C C RELDIF -- DETERMINES THE MAGNITUDE OF THE C DIFFERENCING STEP. C IF KORDIF.EQ.1 , RELDIF SHOULD USUALLY C BE OF THE ORDER OF C SQRT(MACHINE EPSILON). C IF KORDIF.EQ.2 , RELDIF SHOULD USUALLY C BE OF THE ORDER OF C (MACHINE EPSILON)**(1.0/3.0) . C FOR DOUBLE PRECISION, YOU CAN USE C MACHINE EPSILON = 1.0D-16 . C C RELMIN -- RELATIVE CONVERGENCE CRITERION USED FOR C EACH X(J) FOR WHICH DELMIN(J) IS ZERO C C METHD -- DETERMINES THE METHOD USED... C C =-1, MODIFIED GAUSS-NEWTON METHOD, C C = 0, GAUSS-NEWTON METHOD, C C = 1, MODIFIED FORM OF MARQUARDT'S C METHOD, C C = 2, MARQUARDT'S METHOD C C THE RECOMMENDED VALUE OF METHD IS C USUALLY METHD=1 . C C KALCP -- DETERMINES WHICH ROUTINE IS CALLED C TO COMPUTE THE JACOBIAN MATRIX C OF FIRST PARTIAL DERIVATIVES... C C =-1, ONE ROW OF P RETURNED BY DERIV C OR CALCD. DERIV CALLS FOFX. C C = 0, ALL OF P RETURNED BY DERIV OR C CALCD. DERIV CALLS FUNK. C C = 1, ALL OF P RETURNED BY DERIV OR C CALCD. DERIV CALLS FOFX. C C THE RECOMMENDED VALUE OF KALCP IS C USUALLY KALCP=0 . C C KORDIF -- DETERMINES HOW P IS TO BE CALCULATED... C C = 1, FIRST ORDER APPROXIMATION USED BY C DERIV, C C = 2, SECOND ORDER APPROXIMATION USED BY C DERIV, C C = 3, ANALYTICAL DERIVATIVE SUPPLIED BY C CALCD C C THE RECOMMENDED VALUE OF KORDIF IS C KORDIF=1 FOR SPEED OR KORDIF=2 FOR C ACCURACY. KORDIF=3 IS ACCURATE AND C FAIRLY FAST, BUT A LOT MORE WORK. C C MAXIT -- THE MAXIMUM NUMBER OF ITERATIONS TO BE C PERFORMED C C LEQU -- SAVES STORAGE IF ALL YSIG(JPT) ARE C THE SAME... C C =1 IF ALL YSIG(JPT) ARE EQUAL C (IN THIS CASE ONLY YSIG(1) IS C REFERENCED), C C =0 OTHERWISE C C MAXSUB -- THE MAXIMUM NUMBER OF SUBITERATIONS C PERMITTED C C MAXUPD -- NOT USED AT PRESENT C C C MODIFICATIONS TO MARQUARDT'S METHOD... C C 1. THE QUANTITY (1+LAMBDA) IS NOT ALLOWED TO INCREASE BY C MORE THAN A FACTOR OF FNULM, IN ANY SINGLE CHANGE. C C 2. WHEN CUTSTEPS ARE USED IN A LINE SEARCH (WHEN THE C ANGLE GAMMA .LT. GAMMA-SUB-ZERO), THE VALUE OF LAMBDA C IS INCREASED. C C IN BOTH THE MODIFIED MARQUARDT'S METHOD AND THE MODIFIED C GAUSS-NEWTON METHOD, A HALF STEP IS ATTEMPTED FOLLOWING EACH C SUCCESSFUL STEP. THIS AVOIDS ONE FORM OF VERY SLOW C CONVERGENCE. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME C COMPILERS (WATFIV, FOR EXAMPLE) AND IT IS FORBIDDEN BY C SOME OTHERS (MODCOMP II, FOR EXAMPLE). C EXTERNAL FUNK C DOUBLE PRECISION P,Y,YSIG,H,XMAX,XMIN,DELTX,DELMIN, * ERR,FOBJ,FLAMBD,FNU,RELDIF,ZSQRT,ARG,CRIT,FLDEF, * RELMIN,RELTOL,FACCL,HUGE,RZERO,UNITR,RTWO,DELN,FMGN, * FCUT,STFAC,PTERM,SCALJ,SA,PIVOT,EM,SUM,COSIN,SB,SC,HH, * FRMIN,XMX,XMN,DENOM,FRAC,UPFAC,RLFAC,RSFAC,RMSDV, * SDVMX,ZMAX1,ZMIN1,ZABS,ARGB,DOWNU,FNULM,UPNU,YY, * RLAMBD DOUBLE PRECISION X,XSAVE,XTEMP,GRAD,FIT,FITSV, * XSAV,SIG,RTERM,PHI,PHNEW,PHALF,XLIM C INTEGER ITER,J,JBND,JINV,JPT,JPU,JQ,JSMAX,JSUB, * JT,JVARY,JX,JXLIM,K,KALCP,KBND,KERFL,KFLAG, * KORDIF,KPT,KQ,KRANK,KT,KW,KX,L,LEQU,LPDMA, * LPDMB,MASK,MASKT,MATRX,MAXIT,METHD,MRANK, * MAXSUB,MAXUPD,NACSV,NACT,NACTIV,NF,NFLAT,NFMAX, * NLOOP,NMU,NOREP,NPTS,NSMAL,NTRACE,NV,NVMAX, * NXTRA C C THE DIMENSIONS OF THE VECTORS AND MATRICES ARE.... C P(NPTS,NACTIV) (OR P(1,NACTIV) IF KALCP.EQ.-1), C FITSV(NPTS) (OR FITSV(1) IF KALCP.EQ.-1), C X(NV),XMAX(NV),XMIN(NV),DELMIN(NV),ERR(NV,NV+1), C XSAVE(NV),H(NV),MASKT(NV), C GRAD(NACTIV),SCALE(NACTIV), C Y(NPTS),FIT(NPTS),YSIG(NPTS) (OR YSIG(1) IF LEQU.NE.0), C WHERE NACTIV IS THE NUMBER OF ACTIVE (UNMASKED) X(J). C DIMENSION P(LPDMA,LPDMB) DIMENSION Y(1),YSIG(1),FIT(1),FITSV(1) DIMENSION XSAVE(20),H(20),GRAD(20),MASKT(20),XTEMP(20) C C USER COMMON..... COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C C MARQ COMMON..... COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C C SET THE LIBRARY FUNCTIONS FOR DOUBLE PRECISION C (DSQRT, DABS, DMAX1, DMIN1) C OR FOR SINGLE PRECISION C (SQRT, ABS, AMAX1, AMIN1). C IT IS RECOMMENDED THAT AN ARITHMETIC WITH AT LEAST C TWELVE SIGNIFICANT DECIMAL DIGITS BE USED. C ZSQRT(ARG)=DSQRT(ARG) ZABS(ARG)=DABS(ARG) ZMAX1(ARG,ARGB)=DMAX1(ARG,ARGB) ZMIN1(ARG,ARGB)=DMIN1(ARG,ARGB) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C SET FIXED QUANTITIES .... C C NVMAX IS THE MAXIMUM PERMISSIBLE VALUE OF NV. IT IS ALSO C THE DIMENSION OF THE ARRAYS X, XMAX, XMIN, MASK, DELMIN, C XSAVE, H, MASKT, GRAD, AND SCALE, AND THE FIRST DIMENSION C OF ERR. THE SECOND DIMENSION OF ERR IS NVMAX+1. C NVMAX=20 C C CRIT = COSINE OF MARQUARDT'S CRITICAL ANGLE, GAMMA-SUB-ZERO C CRIT=0.70711D0 C C FLDEF = DEFAULT VALUE FOR LAMBDA C FLDEF=1.0D0 C C RELTOL = TOLERANCE FOR A MESSAGE WARNING OF ILL-CONDITIONING C RELTOL=1.0D-4 C C FACCL = ACCELERATION FACTOR FOR FCUT C FACCL=2.0D0 C C FNULM = LIMIT ON THE FACTOR BY WHICH (1+FLAMBD) MAY CHANGE C FNULM=2.0D0 C C HUGE = A LARGE REAL NUMBER, THE DEFAULT VALUE FOR XMAX C AND (-XMIN) C HUGE=1.0D30 C RZERO=0.0D0 UNITR=1.0D0 RTWO=2.0D0 C C NO FLOATING POINT CONSTANTS ARE USED BEYOND THIS POINT. C KFLAG=0 NOREP=0 KERFL=0 C RLAMBD=FLAMBD ITER=0 PHI=HUGE SC=HUGE C IF(NTRACE.GE.-1) WRITE(KW,10) 10 FORMAT('1MARQ.... BEGIN NONLINEAR LEAST SQUARES', * ' SOLUTION') C IF(NV.GT.0 .AND. NV.LE.NVMAX .AND. NV.LE.LPDMB .AND. * NPTS.GE.1 .AND. (KALCP.LT.0 .OR. NPTS.LE.LPDMA)) * GO TO 30 C WRITE(KW,20)NV,NVMAX,NPTS,LPDMA,LPDMB,KALCP 20 FORMAT(//' ****** ILLEGAL INPUT VALUE IN SUBROUTINE', * ' MARQ.'/6X,'NV =',I14,6X,'NVMAX =',I14,6X,'NPTS =', * I14/6X,'LPDMA =',I14,6X,'LPDMA =',I14,6X,'KALCP =', * I14) STOP C C CHECK SOME INPUT QUANTITIES, AND SET THEM TO DEFAULT VALUES C IF DESIRED. C C NACTIV = NUMBER OF ACTIVE X(J) C 30 NACTIV=0 C DO 60 JX=1,NV DELN=ZABS(DELMIN(JX)) IF(MASK(JX).NE.0) GO TO 50 NACTIV=NACTIV+1 IF(DELN.EQ.RZERO) DELN=ZABS(RELMIN*X(JX)) IF(DELN.EQ.RZERO) DELN=RELMIN IF(XMAX(JX).GT.XMIN(JX)) GO TO 40 XMAX(JX)=HUGE XMIN(JX)=-HUGE C 40 X(JX)=ZMAX1(XMIN(JX),ZMIN1(XMAX(JX),X(JX))) C 50 DELMIN(JX)=DELN 60 CONTINUE C IF(NTRACE.LT.-1) GO TO 130 WRITE(KW,70)(MASK(J),J=1,NV) 70 FORMAT(/' MASK =',I7,4I13/(2X,5I13)) WRITE(KW,80)(X(J),J=1,NV) 80 FORMAT(/' X =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,90)(XMAX(J),J=1,NV) 90 FORMAT(/' XMAX =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,100)(XMIN(J),J=1,NV) 100 FORMAT(/' XMIN =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,110)(DELMIN(J),J=1,NV) 110 FORMAT(/' DELMIN=',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,120)NV,NPTS,LPDMA,LPDMB,NTRACE,METHD,KALCP, * KORDIF,NFLAT,NFMAX,MAXIT,MAXSUB,CRIT,RELDIF,RELMIN 120 FORMAT(//' NV =',I11,8X,'NPTS =',I11,8X,'LPDMA =',I11// * ' LPDMB =',I11,8X,'NTRACE =',I4,8X,'METHD =',I3// * ' KALCP =',I3,8X,'KORDIF =',I2,8X,'NFLAT =',I2,8X, * 'NFMAX =',I11//' MAXIT =',I8,8X,'MAXSUB =',I6// * ' CRIT =',1PG12.5,8X,'RELDIF =',G12.5,8X,'RELMIN =', * G12.5) C 130 JVARY=0 C C SET FMGN. IF NECESSARY, RESET RLAMBD AND/OR FACCL. C FMGN=UNITR UPNU=FNU DOWNU=FNU IF(RLAMBD.LE.RZERO) RLAMBD=FLDEF IF(METHD.LE.0) RLAMBD=RZERO IF(METHD.NE.1) FACCL=UNITR C C COMPUTE THE INITIAL GOODNESS OF FIT OF THE MODEL TO THE C DATA. CALL FUNC TO CALCULATE THE VECTOR OF FITTED VALUES. C CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) C C NF = EQUIVALENT NUMBER OF CALLS TO FUNC C NF=1 IF(NTRACE.GE.-1) WRITE(KW,140)PHI,RLAMBD 140 FORMAT(//' PHI (THE SUM OF SQUARES) =',1PG15.8,10X, * 'LAMBDA =',G12.5//' ') C IF(NACTIV.GT.0) GO TO 160 KFLAG=-1 IF(NTRACE.GE.-3) WRITE(KW,150) 150 FORMAT(//' ****** WARNING... MASK(J).NE.0 FOR ALL J,', * ' IN SUBROUTINE MARQ.'//6X, * 'PHI WILL BE EVALUATED BUT NOT MINIMIZED.') GO TO 1210 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C BEGIN THE NEXT ITERATION. C C THIS IS THE ENTRY POINT AFTER A SUCCESSFUL STEP, IF THE C CONVERGENCE CRITERION IS NOT MET. C 160 JSUB=0 FCUT=RTWO ITER=ITER+1 IF(NTRACE.GE.1) WRITE(KW,170)ITER,FMGN,RLAMBD 170 FORMAT(//' BEGIN ITERATION',I5,8X,'FMGN =',1PG9.2, * 8X,'LAMBDA =',G9.2) IF(NTRACE.GE.3) WRITE(KW,180) 180 FORMAT(/' P(*,*) (THE JACOBIAN MATRIX) ='/' ') C C INITIALIZE FOR THIS ITERATION. C STFAC=UNITR C DO 200 JX=1,NACTIV GRAD(JX)=RZERO C DO 190 KX=1,JX ERR(JX,KX)=RZERO 190 CONTINUE C 200 CONTINUE C C CALL DERIV (OR CALCD) TO COMPUTE THE JACOBIAN MATRIX, C P(*,*). DERIV IS CALLED NPTS TIMES IF KALCP=-1 . C SIG=YSIG(1) C DO 270 JPT=1,NPTS KPT=JPT IF(KALCP.LT.0) GO TO 210 IF(JPT.NE.1) GO TO 230 C 210 KPT=1 IF(KORDIF.LE.2) GO TO 220 CALL CALCD (JPT,P,LPDMA,LPDMB) GO TO 230 C 220 CALL DERIV (JPT,FUNK,NPTS,FIT,FITSV,P,LPDMA,LPDMB) C 230 IF(NTRACE.GE.3) WRITE(KW,240)JPT, * (P(KPT,JX),JX=1,NACTIV) 240 FORMAT(1X,I3,2X,1PE13.5,4E13.5/(6X,5E13.5)) C C COMPUTE QSAV(*,*) AND GRAD(*). C QSAV(*,*), WHICH IS STORED IN THE UPPER TRIANGLE OF C THE ARRAY ERR(*,*), IS PT*P . C GRAD(*) IS HALF THE GRADIENT OF PHI. C IF(LEQU.EQ.0) SIG=YSIG(JPT) RTERM=(FIT(JPT)-Y(JPT))/SIG**2 C DO 260 JX=1,NACTIV GRAD(JX)=GRAD(JX)+P(KPT,JX)*RTERM PTERM=P(KPT,JX)/SIG**2 C DO 250 KX=1,JX ERR(JX,KX)=ERR(JX,KX)+P(KPT,KX)*PTERM 250 CONTINUE C 260 CONTINUE C 270 CONTINUE C C RESTORE FIT(*) IF IT WAS DESTROYED IN DERIV. C IF(KORDIF.NE.2) GO TO 290 IF(KALCP.NE.0) GO TO 280 CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) NF=NF+1 C 280 NF=NF+NACTIV C 290 NF=NF+NACTIV C C COMPUTE THE SCALE FACTORS AND STORE THEM IN DELTX(*). C SCALE GRAD(*). C DO 300 JX=1,NACTIV SCALJ=ZSQRT(ERR(JX,JX)) IF(SCALJ.EQ.RZERO) SCALJ=UNITR DELTX(JX)=SCALJ GRAD(JX)=GRAD(JX)/SCALJ 300 CONTINUE C IF(NTRACE.GE.1) WRITE(KW,310)(GRAD(JX),JX=1,NACTIV) 310 FORMAT(/' SCALED GRADIENT VECTOR GRAD(*) ='/ * (1X,1PG15.5,4G15.5)) C C SCALE QSAV(*,*). C THE DIAGONAL ELEMENTS OF QSAV(*,*) ARE SCALED TO UNITY. C DO 370 JX=1,NACTIV C DO 360 KX=1,JX SA=ERR(JX,KX)/(DELTX(JX)*DELTX(KX)) IF(KX.NE.JX) GO TO 320 IF(SA.EQ.RZERO) GO TO 330 SA=UNITR GO TO 350 C 320 IF(ZABS(SA).LT.UNITR-RELTOL) GO TO 350 C 330 IF(NTRACE.GE.-2) WRITE(KW,340)JX,KX,SA,ITER 340 FORMAT(' ******** POSSIBLY DANGEROUS VALUE OF', * ' COEFFICIENT....'/6X,'QSAV(',I3,',',I3,') =', * 1PG15.7,10X,'ITERATION',I5) C 350 ERR(JX,KX)=SA 360 CONTINUE C 370 CONTINUE C IF(NTRACE.LT.2) GO TO 400 WRITE(KW,380) 380 FORMAT(/' QSAV(*,*) (PT*P, SCALED, WHERE P IS THE', * ' JACOBIAN MATRIX) ='/' ') C DO 390 JX=1,NACTIV WRITE(KW,240)JX,(ERR(JX,KX),KX=1,JX) 390 CONTINUE C 400 DO 410 JX=1,NV XSAVE(JX)=X(JX) 410 CONTINUE C C INITIALIZE MASKT AND NACT. C 420 NACT=NACTIV C DO 430 JX=1,NV MASKT(JX)=MASK(JX) 430 CONTINUE C C COPY QSAV(*,*) INTO Q(*,*) AND COPY GRAD(*) INTO H(*), C AND SET THE DIAGONAL ELEMENTS OF Q(*,*). C Q(*,*) IS STORED IN THE UPPER TRIANGLE OF ERR(*,*) . C C THIS IS THE ENTRY POINT FOR SUBITERATIONS IN WHICH RLAMBD IS C INCREASED OR CONSTRAINTS ARE IMPOSED. C 440 KRANK=0 JQ=0 JT=0 C DO 470 JX=1,NV IF(MASK(JX).NE.0) GO TO 470 JQ=JQ+1 IF(MASKT(JX).NE.0) GO TO 470 JT=JT+1 H(JT)=-GRAD(JQ) KQ=0 KT=0 C DO 460 KX=1,JX IF(MASK(KX).NE.0) GO TO 460 KQ=KQ+1 IF(MASKT(KX).NE.0) GO TO 460 KT=KT+1 SA=ERR(JQ,KQ) IF(KX.NE.JX .OR. SA.EQ.RZERO) GO TO 450 SA=UNITR+RLAMBD KRANK=KRANK+1 C 450 ERR(KT,JT+1)=SA 460 CONTINUE C 470 CONTINUE C C SOLVE THE NORMAL EQUATIONS FOR H(*), THE CORRECTION VECTOR. C NSMAL=0 NMU=NACT-1 IF(NMU.EQ.0) GO TO 520 C C REDUCE THE SYSTEM TO TRIANGULAR FORM, UTILIZING THE SYMMETRY C OF THE MATRIX. C DO 510 J=1,NMU PIVOT=ERR(J,J+1) IF(PIVOT.EQ.RZERO) GO TO 510 JPU=J+1 C DO 500 K=JPU,NACT EM=ERR(J,K+1)/PIVOT IF(EM.EQ.RZERO) GO TO 490 C DO 480 L=K,NACT ERR(K,L+1)=ERR(K,L+1)-ERR(J,L+1)*EM 480 CONTINUE C 490 H(K)=H(K)-H(J)*EM 500 CONTINUE C 510 CONTINUE C C DO THE BACK SOLUTION. C 520 DO 560 JINV=1,NACT J=(NACT+1)-JINV PIVOT=ERR(J,J+1) IF(PIVOT.LE.RZERO) NSMAL=NSMAL+1 IF(PIVOT.NE.RZERO) GO TO 530 H(J)=RZERO GO TO 560 C 530 SUM=RZERO IF(J.EQ.NACT) GO TO 550 JPU=J+1 C DO 540 K=JPU,NACT SUM=SUM+ERR(J,K+1)*H(K) 540 CONTINUE C 550 H(J)=(H(J)-SUM)/PIVOT 560 CONTINUE C C IF THE COEFFICIENT MATRIX WAS RANK DEFICIENT, PRINT A C MESSAGE. C MRANK=NACT-NSMAL IF(MRANK.EQ.NACT) GO TO 590 COSIN=HUGE IF(NTRACE.GE.-2) WRITE(KW,570)MRANK,NACT,ITER 570 FORMAT(/' RANK-DEFICIENT NORMAL EQUATIONS IN MARQ.'/8X, * 'RANK =',I4,8X,'ORDER OF MATRIX =',I4,8X, * 'ITERATION',I5) IF(MRANK.GT.0) GO TO 580 KFLAG=-4 GO TO 1210 C 580 IF(METHD.GT.0 .AND. MRANK.LT.KRANK) GO TO 870 C C UNPACK AND DE-SCALE THE CORRECTION VECTOR H(*). C COMPUTE THE INNER PRODUCTS SA, SB, AND SC. C 590 SA=RZERO SB=RZERO SC=RZERO KX=NV KQ=NACTIV KT=NACT C DO 620 JX=1,NV HH=RZERO IF(MASK(KX).NE.0) GO TO 610 IF(MASKT(KX).NE.0) GO TO 600 HH=H(KT) SA=SA-HH*GRAD(KQ) SB=SB+HH*HH SC=SC+GRAD(KQ)**2 HH=HH*FMGN/DELTX(KQ) KT=KT-1 C 600 KQ=KQ-1 C 610 H(KX)=HH KX=KX-1 620 CONTINUE C C ADD THE CORRECTION VECTOR TO THE PARAMETER VECTOR AND C CHECK FOR CONSTRAINT VIOLATIONS. C C THIS IS THE ENTRY POINT FOLLOWING A CUTSTEP. C 630 IF(NTRACE.GE.1) WRITE(KW,640)(H(JX),JX=1,NV) 640 FORMAT(/' CORRECTION VECTOR H(*) ='/(1X,1PG15.5,4G15.5)) NACSV=NACT FRMIN=UNITR NLOOP=0 JXLIM=0 C 650 DO 740 JX=1,NV IF(MASKT(JX).NE.0) GO TO 740 XSAV=XSAVE(JX) XMX=XMAX(JX) XMN=XMIN(JX) HH=H(JX) IF((XSAV.LT.XMX .OR. HH.LE.RZERO) .AND. * (XSAV.GT.XMN .OR. HH.GE.RZERO)) GO TO 670 MASKT(JX)=1 NACT=NACT-1 X(JX)=XSAV IF(NTRACE.GE.-1) WRITE(KW,660)JX,XSAV,ITER 660 FORMAT(/' FIX X(',I3,') = ',1PG15.8, * ' TEMPORARILY'/5X,'TO AVOID VIOLATING A', * ' CONSTRAINT.',8X,'ITERATION',I5) GO TO 740 C 670 XLIM=XSAV+HH*FRMIN IF(JX.NE.JXLIM) GO TO 680 IF(KBND.LT.0) GO TO 700 C 680 X(JX)=XLIM IF(XLIM.LE.XMX) GO TO 690 C X(JX)=XMX JBND=1 GO TO 710 C 690 IF(XLIM.GE.XMN) GO TO 740 C 700 X(JX)=XMN JBND=-1 C 710 IF(JX.NE.JXLIM) GO TO 730 IF(NTRACE.LT.-1) GO TO 740 WRITE(KW,720)JX,X(JX),FRMIN,ITER 720 FORMAT(/' CONSTRAINT VIOLATED BY X(',I3, * ').'/' VALUE RESET TO',1PG15.7/6X, * ' USING CUTSTEP FACTOR =',G12.5,8X,'ITERATION',I5) GO TO 740 C 730 IF(NLOOP.NE.0) GO TO 740 DENOM=XLIM-XSAV IF(DENOM.EQ.RZERO) GO TO 740 FRAC=(X(JX)-XSAV)/DENOM IF(FRAC.GE.FRMIN) GO TO 740 FRMIN=FRAC JXLIM=JX KBND=JBND 740 CONTINUE C C IF THE PROPOSED STEP WOULD VIOLATE ANY ALREADY ACTIVE C CONTRAINTS, FIX THOSE COMPONENTS OF H(*) EQUAL TO ZERO AND C RECOMPUTE THE OTHER COMPONENTS. C IF(NACT.GT.0) GO TO 760 KFLAG=3 IF(NTRACE.LT.-2) GO TO 1210 WRITE(KW,750) 750 FORMAT(//' APPARENT CONSTRAINED OPTIMUM LIES IN A', * ' CORNER.') GO TO 1210 C 760 IF(NACT.LT.NACSV) GO TO 440 IF(NLOOP.NE.0) GO TO 770 NLOOP=1 IF(JXLIM.NE.0) GO TO 650 C 770 IF(NTRACE.GE.1) WRITE(KW,780)(X(JX),JX=1,NV) 780 FORMAT(/' X(*) ='/(1X,1PG15.7,4G15.7)) C C CALCULATE THE NEW FITTED VALUES. C CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHNEW) NF=NF+1 IF(PHNEW.LT.PHI) GO TO 940 IF(PHNEW.GT.PHI) GO TO 800 C C THE NEW VALUE OF PHI IS EXACTLY EQUAL TO THE OLD VALUE. C CHECK FOR CONVERGENCE UNDER THE NFLAT OPTION. C IF(NFLAT.EQ.0) GO TO 940 KFLAG=2 IF(NTRACE.LT.-1) GO TO 840 WRITE(KW,790) 790 FORMAT(//' CONVERGENCE ACHIEVED UNDER THE NFLAT OPTION', * ' IN MARQ.') GO TO 840 C C THE NEW VALUE OF PHI IS GREATER THAN THE OLD VALUE. C 800 IF(NTRACE.GE.1) WRITE(KW,810)PHI,PHNEW 810 FORMAT(/' OLD PHI =',1PG15.8,8X,'NEW PHI =',G15.8) C C CHECK WHETHER JSUB HAS EXCEEDED MAXSUB. C JSUB=JSUB+1 IF(JSUB.GT.MAXSUB) GO TO 820 IF(METHD.EQ.0) GO TO 1100 IF(METHD.LT.0) GO TO 890 GO TO 860 C 820 KFLAG=-1 IF(NTRACE.GE.-1) WRITE(KW,830)MAXSUB 830 FORMAT(//' EXCEEDED MAXIMUM NUMBER OF SUBITERATIONS =', * I3,' IN MARQ.') C C RESTORE X(*) TO THE BASE POINT. C 840 DO 850 JX=1,NV X(JX)=XSAVE(JX) 850 CONTINUE C CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) GO TO 1210 C C THE NEW FIT IS WORSE THAN THE OLD FIT. C COMPUTE COSIN, THE COSINE OF THE ANGLE BETWEEN THE SCALED C GRADIENT AND THE SCALED CORRECTION VECTOR. C 860 DENOM=SB*SC IF(DENOM.LE.RZERO) GO TO 870 COSIN=SA/ZSQRT(DENOM) IF(COSIN.GT.CRIT) GO TO 890 C C COSIN IS NOT GREATER THAN CRIT. C INCREASE THE VALUE OF LAMBDA. C 870 UPFAC=UPNU UPNU=ZMIN1(UPNU*RTWO,FNU) IF(METHD.EQ.1) UPFAC=ZMIN1(UPFAC,FNULM+UNITR/RLAMBD) RLAMBD=RLAMBD*UPFAC IF(NTRACE.GE.1) WRITE(KW,880)JSUB,COSIN,RLAMBD 880 FORMAT(/' **** SUBITERATION',I3,8X,'INCREASE LAMBDA.'/8X, * 'COS(GAMMA) =',1PG12.5,8X,'LAMBDA =',G12.5) C C GO BACK AND FORM THE NORMAL EQUATIONS USING A LARGER VALUE C OF LAMBDA. THIS WILL CAUSE A SHORTER STEP VECTOR H(*) TO BE C COMPUTED. C GO TO 420 C C COSIN IS GREATER THAN CRIT. C CUT THE MAGNITUDE OF THE STEP VECTOR H(*). C 890 STFAC=STFAC/FCUT IF(METHD.GE.0) GO TO 900 FMGN=FMGN/FCUT GO TO 910 C 900 IF(METHD.EQ.1) RLAMBD=RLAMBD*RTWO C 910 DO 920 JX=1,NV H(JX)=(X(JX)-XSAVE(JX))/FCUT 920 CONTINUE C FCUT=FCUT*FACCL IF(NTRACE.GE.1) WRITE(KW,930)JSUB,COSIN,STFAC 930 FORMAT(/' **** SUBITERATION',I3,8X,'TAKE CUT STEPS.'/5X, * 'COS(GAMMA) =',1PG12.5,8X,'CUTSTEP FACTOR =',G12.5) C C GO BACK AND TRY A SMALLER CUTSTEP. C GO TO 630 C C THE VALUE OF PHI HAS DECREASED. TRY A HALF STEP. C 940 IF(METHD.EQ.0 .OR. METHD.EQ.2) GO TO 1100 C DO 950 JX=1,NV XTEMP(JX)=X(JX) IF(MASK(JX).NE.0) GO TO 950 X(JX)=XSAVE(JX)+(X(JX)-XSAVE(JX))/RTWO X(JX)=ZMAX1(XMIN(JX),ZMIN1(XMAX(JX),X(JX))) 950 CONTINUE C DO 960 JPT=1,NPTS FITSV(JPT)=FIT(JPT) 960 CONTINUE C CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHALF) NF=NF+1 C C USE QUADRATIC INTERPOLATION, IN ORDER TO TRY TO REFINE THE C POSITION OF THE MINIMUM OF PHI. C RLFAC=UNITR DENOM=RTWO*((PHNEW-PHALF)-(PHALF-PHI)) STFAC=RZERO IF(DENOM.LE.RZERO) GO TO 970 STFAC=(PHI-PHNEW)/DENOM RSFAC=(UNITR+STFAC)/RTWO C C DO NOT EXTRAPOLATE. C IF(STFAC.GE.UNITR) STFAC=RZERO C 970 DO 980 JX=1,NV H(JX)=X(JX) X(JX)=X(JX)+(XTEMP(JX)-X(JX))*STFAC 980 CONTINUE C IF(PHALF.GE.PHNEW) GO TO 1020 RLFAC=UNITR/RTWO JSUB=JSUB+1 C DO 990 JX=1,NV XTEMP(JX)=H(JX) 990 CONTINUE C DO 1000 JPT=1,NPTS FITSV(JPT)=FIT(JPT) 1000 CONTINUE C IF(NTRACE.GE.1) WRITE(KW,1010)PHNEW,PHALF 1010 FORMAT(/' HALF STEP SUCCEEDED.'/8X,'PHNEW =',1PG15.8,8X, * 'PHALF =',G15.8) PHNEW=PHALF C 1020 IF(STFAC.EQ.RZERO) GO TO 1030 CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) NF=NF+1 IF(PHI.LT.PHNEW) GO TO 1070 C 1030 DO 1040 JX=1,NV X(JX)=XTEMP(JX) 1040 CONTINUE C DO 1050 JPT=1,NPTS FIT(JPT)=FITSV(JPT) 1050 CONTINUE C IF(STFAC.NE.RZERO .AND. NTRACE.GE.1) WRITE(KW,1060)RSFAC, * PHI 1060 FORMAT(/' QUADRATIC INTERPOLATION.'/8X,'RSFAC =',1PG12.5, * 8X,'PHI =',G15.8) GO TO 1090 C 1070 RLFAC=RSFAC PHNEW=PHI IF(NTRACE.GE.1) WRITE(KW,1080)RLFAC,PHI 1080 FORMAT(/' QUADRATIC INTERPOLATION SUCCEEDED.'/8X, * 'RLFAC =',1PG12.5,8X,'PHI =',G15.8) C 1090 IF(RLFAC.LE.RZERO) GO TO 1100 RLAMBD=RLAMBD/RLFAC IF(METHD.LT.0) FMGN=FMGN*RLFAC C C THE STEP IS NOW ACCEPTED. C C TEST FOR CONVERGENCE, PROVIDED THAT NO CONSTRAINT BECAME C ACTIVE DURING THIS ITERATION. C 1100 IF(NTRACE.GE.1) WRITE(KW,1110)ITER,PHNEW 1110 FORMAT(/' END OF ITERATION',I5,8X,'PHI =',1PG15.8) PHI=PHNEW IF(JXLIM.GT.0) GO TO 1140 C DO 1120 JX=1,NV IF(MASK(JX).NE.0) GO TO 1120 IF(ZABS(X(JX)-XSAVE(JX)).GT.DELMIN(JX)) GO TO 1140 1120 CONTINUE C KFLAG=1 IF(NTRACE.LT.-1) GO TO 1210 WRITE(KW,1130) 1130 FORMAT(//' MARQ CONVERGED WHEN THE STEP BECAME SMALL.') GO TO 1210 C C THE ITERATION HAS NOT YET CONVERGED. C 1140 IF(ITER.LT.MAXIT) GO TO 1160 KFLAG=-6 IF(NTRACE.GE.-3) WRITE(KW,1150)MAXIT 1150 FORMAT(//' MAXIMUM NUMBER OF ITERATIONS REACHED IN', * ' MARQ.',8X,'MAXIT =',I5) GO TO 1210 C C IF SUBITERATIONS WERE NOT PERFORMED THIS ITERATION, DECREASE C THE VALUE OF LAMBDA. C 1160 IF(NF.GE.NFMAX) GO TO 1190 IF(JSUB.GT.0) GO TO 1170 FMGN=ZMIN1(FMGN*RTWO,UNITR) SCALJ=UNITR+RLAMBD IF(SCALJ.GT.UNITR) RLAMBD=RLAMBD/DOWNU UPNU=DOWNU DOWNU=ZMIN1(DOWNU*RTWO,FNU) GO TO 1180 C 1170 IF(METHD.NE.1) GO TO 1180 UPNU=FNULM DOWNU=FNULM C 1180 CONTINUE C C GO BACK AND DO ANOTHER ITERATION. C GO TO 160 C 1190 KFLAG=-7 IF(NTRACE.GE.-3) WRITE(KW,1200)NFMAX 1200 FORMAT(//' NF HAS REACHED NFMAX =',I8,'IN MARQ.') C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE ITERATION HAS TERMINATED. C C PRINT OUT THE DATA, FITTED VALUES, AND RESIDUALS. C COMPUTE AND PRINT THE STANDARD DEVIATION OF THE DATA FROM C THE FIT. C 1210 SC=ZSQRT(SC) IF(NTRACE.LT.-1) GO TO 1300 WRITE(KW,1220)ITER,NF,PHI,FMGN,RLAMBD,SC 1220 FORMAT(/1X,I4,' ITERATIONS',8X,'NF =',I5,8X,'PHI =', * 1PG15.8//8X,'FMGN =',G12.5,8X,'LAMBDA =',G12.5// * ' NORM OF SCALED GRADIENT VECTOR =',G12.5) WRITE(KW,780)(X(JX),JX=1,NV) IF(NTRACE.GE.0) WRITE(KW,1230) 1230 FORMAT(///4X,'J',5X,'Y(J)',10X,'FIT(J)',6X, * 'Y(J)-FIT(J)',2X,'YSIG(J)',3X,'(Y-FIT)/YSIG'/' ') SIG=YSIG(1) RMSDV=RZERO SDVMX=RZERO JSMAX=1 C DO 1250 JPT=1,NPTS IF(LEQU.EQ.0) SIG=YSIG(JPT) YY=Y(JPT) PTERM=YY-FIT(JPT) RTERM=PTERM/SIG IF(NTRACE.GE.0) WRITE(KW,1240)JPT,YY,FIT(JPT),PTERM, * SIG,RTERM 1240 FORMAT(1X,I4,1PE15.7,E15.7,3E12.4) RMSDV=RMSDV+RTERM**2 IF(ZABS(RTERM).LE.SDVMX) GO TO 1250 SDVMX=ZABS(RTERM) JSMAX=JPT 1250 CONTINUE C DENOM=NPTS-NACTIV WRITE(KW,1260)DENOM 1260 FORMAT(//' NUMBER OF DEGREES OF FREEDOM = ',1PG12.5) IF(DENOM.LE.RZERO) GO TO 1280 RMSDV=ZSQRT(RMSDV/DENOM) WRITE(KW,1270)RMSDV 1270 FORMAT(/' R.M.S. SCALED DEVIATION OF DATA FROM FIT =', * 1PG12.5) C 1280 WRITE(KW,1290)SDVMX,JSMAX 1290 FORMAT(/' MAXIMUM SCALED DEVIATION 0F',1PG12.5/ * 5X,'OCCURRED AT DATA POINT NUMBER',I6) C IF(NTRACE.GE.0) CALL RPLOT (KW,NPTS,Y,YSIG,LEQU,FIT,0) C C CALL FUNC TO SET THE FINAL VALUES OF FIT(*). C 1300 CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) FOBJ=PHI C C CALL MQERR TO PRINT THE PARAMETER ERRORS AND CORRELATIONS. C A DUMMY ROUTINE MAY BE SUBSTITUTED FOR MQERR IF THESE ARE C NOT NEEDED. C IF(MATRX.GT.0) CALL MQERR (NACTIV,NPTS) C RETURN C C END MARQ C END SUBROUTINE MQERR (NACTIV,NPTS) C C MQERR 3.2 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C MQERR IS CALLED BY MARQ TO COMPUTE AND PRINT C APPROXIMATE VALUES OF THE PARAMETER ERRORS AND CORRELATIONS. C C FOR THE MEANING OF THE "MAXIMUM VARIANCE INFLATION FACTOR", C SEE... D. W. MARQUARDT AND R. D. SNEE, C "RIDGE REGRESSION IN PRACTICE", C THE AMERICAN STATISTICIAN 29 (1975) 3-20 C C INPUT QUANTITIES..... KW,ERR(*,*),NACTIV,DELTX(*),NPTS,NV, C NTRACE,MASK(*),FOBJ C C OUTPUT QUANTITIES.... ERR(*,*) C DOUBLE PRECISION SCALJ,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * ZSQRT,ARG,RZERO,UNITR,HUGE,PIVOT,Q,VIFMX,ER,TEMP, * DENOM,SCFAC,RESCL ,X C INTEGER JV,JVARY,JX,K,KERFL,KFLAG,KV,KW,KX,L,LINV, * M,MASK,MATRX,NACTIV,NDF,NFLAT,NFMAX,NOREP,NPTS, * NRANK,NSMAL,NTRACE,NV,NVPLU,NXTRA C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C ZSQRT(ARG)=DSQRT(ARG) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C RZERO=0.0D0 UNITR=1.0D0 HUGE=1.0D30 C C PRINT QSAV(*,*). C IF(NTRACE.LT.-1) GO TO 40 WRITE(KW,10) 10 FORMAT(///' SUBROUTINE MQERR.'//' QSAV(*,*) (PT*P,', * ' SCALED, WHERE P IS THE JACOBIAN) =') DO 30 JX=1,NACTIV WRITE(KW,20)JX,(ERR(JX,KX),KX=1,JX) 20 FORMAT(/1X,I3,2X,1PG13.5,4G13.5/(6X,5G13.5)) 30 CONTINUE C C COMPUTE THE SCALED ERROR MATRIX, WHICH IS THE INVERSE OF C QSAV. C C INVERT QSAV USING THE GAUSS-JORDAN METHOD WITHOUT PIVOTING. C F. L. BAUER AND C. REINSCH, P. 45 IN "LINEAR ALGEBRA" C BY J. H. WILKINSON AND C. REINSCH (SPRINGER-VERLAG, 1971). C ERR(*,NV+1) IS USED AS A SCRATCH VECTOR. C 40 NVPLU=NV+1 NSMAL=0 C DO 140 LINV=1,NACTIV L=(NACTIV+1)-LINV PIVOT=ERR(1,1) IF(PIVOT.LE.RZERO) NSMAL=NSMAL+1 IF(NACTIV.LT.2) GO TO 100 C DO 90 K=2,NACTIV Q=ERR(K,1) IF(PIVOT.NE.RZERO) GO TO 50 ERR(K,NVPLU)=RZERO GO TO 70 C 50 IF(K.GT.L) GO TO 60 ERR(K,NVPLU)=-Q/PIVOT GO TO 70 C 60 ERR(K,NVPLU)=Q/PIVOT C 70 DO 80 M=2,K C ERR(K-1,M-1)=ERR(K,M)+Q*ERR(M,NVPLU) 80 CONTINUE C 90 CONTINUE C 100 IF(PIVOT.NE.RZERO) GO TO 110 ERR(NACTIV,NACTIV)=RZERO GO TO 120 C 110 ERR(NACTIV,NACTIV)=UNITR/PIVOT C 120 IF(NACTIV.LT.2) GO TO 140 C DO 130 K=2,NACTIV ERR(NACTIV,K-1)=ERR(K,NVPLU) 130 CONTINUE C 140 CONTINUE C NRANK=NACTIV-NSMAL IF(NRANK.LT.NACTIV .AND. NTRACE.GE.-3) * WRITE(KW,150)NRANK,NACTIV 150 FORMAT(//' THE SECOND DERIVATIVE MATRIX IS SINGULAR', * ' IN MQERR.'/8X,'RANK =',I3,8X,'ORDER =',I3/ * ' THEREFORE ALL PARAMETER ERRORS ARE INFINITE.') C C UNPACK THE ERROR MATRIX INTO THE UPPER C TRIANGLE OF ERR(*,*), DE-SCALING IT. C JV=0 VIFMX=RZERO C DO 210 JX=1,NV IF(MASK(JX).EQ.0) JV=JV+1 KV=0 C DO 200 KX=1,JX ER=RZERO IF(MASK(JX).NE.0 .OR. MASK(KX).NE.0) GO TO 190 KV=KV+1 TEMP=RZERO DENOM=DELTX(JV)*DELTX(KV) IF(DENOM.EQ.RZERO) GO TO 160 TEMP=ERR(JV,KV) ER=TEMP/DENOM C 160 IF(JX.NE.KX) GO TO 190 IF(TEMP.GT.RZERO) GO TO 180 IF(NTRACE.GE.-2) WRITE(KW,170)JX,KX,TEMP 170 FORMAT(/' THE (',I3,',',I3, * ') ELEMENT OF QSAV**-1 =',1PG13.5/ * ' . THEREFORE ALL PARAMETER ERRORS ARE', * ' INFINITE.') TEMP=-TEMP C 180 IF(TEMP.GT.VIFMX) VIFMX=TEMP C 190 ERR(KX,JX+1)=ER C 200 CONTINUE C 210 CONTINUE C C COMPUTE AND PRINT THE STANDARD ERRORS. C NDF=NPTS-NACTIV SCFAC=HUGE IF(NDF.LE.0) GO TO 220 SCFAC=NDF SCFAC=ZSQRT(FOBJ/SCFAC) C 220 RESCL=NDF+NDF IF(RESCL.GT.RZERO) RESCL=ZSQRT(RESCL) IF(NTRACE.LT.-1) GO TO 350 WRITE(KW,230)NDF,NDF,RESCL,FOBJ,SCFAC 230 FORMAT(///' NUMBER OF DEGREES OF FREEDOM (N.D.F.) =', * ' (NPTS-NACTIV) =',I5//' EXPECTED VALUE OF PHI =', * ' N.D.F. PLUS OR MINUS SQRT(2*N.D.F.)'/23X,'=',I5, * ' PLUS OR MINUS',1PG12.5// * ' ACTUAL VALUE OF PHI =',G12.5// * ' RESCALING FACTOR = SQRT(PHI/N.D.F.) =',G12.5) WRITE(KW,240)VIFMX 240 FORMAT(///' MAXIMUM VARIANCE INFLATION FACTOR =',1PG12.5, * ' = CONDITION NUMBER'/// * ' APPROXIMATE STANDARD ERRORS....'/ * 57X,'RESCALED'/3X,'J',6X,'MASK(J)',7X,'X(J)', * 13X,'ERROR',12X,'ERROR') RESCL=HUGE C DO 280 JX=1,NV SCALJ=UNITR ER=ERR(JX,JX+1) IF(ER.GE.RZERO) GO TO 250 ER=-ZSQRT(-ER) SCALJ=-ER GO TO 260 C 250 ER=ZSQRT(ER) SCALJ=ER C 260 IF(NDF.GT.0 .AND. NRANK.EQ.NACTIV) RESCL=SCFAC*ER DELTX(JX)=SCALJ WRITE(KW,270)JX,MASK(JX),X(JX),ER,RESCL 270 FORMAT(/1X,I3,I10,4X,1PG16.8,4X,G13.5,4X,G13.5) 280 CONTINUE C C COMPUTE AND PRINT THE CORRELATIONS. C IF(NV.LT.2) GO TO 350 WRITE(KW,290)(K,K=1,NV) 290 FORMAT(///' LOWER TRIANGLE OF THE CORRELATION MATRIX....' * //11X,'K =',I3,4I12/(6X,5I12)) WRITE(KW,300)(MASK(K),K=1,NV) 300 FORMAT(/6X,'MASK(K) =',I3,4I12/(6X,5I12)) WRITE(KW,310) 310 FORMAT(/' J MASK(J)') C DO 340 JX=1,NV C DO 320 KX=1,JX DENOM=DELTX(JX)*DELTX(KX) ERR(KX,1)=RZERO IF(DENOM.NE.RZERO) ERR(KX,1)=ERR(KX,JX+1)/DENOM 320 CONTINUE C WRITE(KW,330)JX,MASK(JX),(ERR(KX,1),KX=1,JX) 330 FORMAT(/1X,I3,I6,2X,1PG12.4,4G12.4/(12X,5G12.4)) 340 CONTINUE C C RESCALE ERR AND SYMMETRIZE IT. C 350 SCFAC=SCFAC**2 IF(SCFAC.LT.UNITR) SCFAC=UNITR C DO 370 JX=1,NV C DO 360 KX=1,JX ERR(JX,KX)=ERR(KX,JX+1)*SCFAC ERR(KX,JX)=ERR(JX,KX) 360 CONTINUE C 370 CONTINUE C RETURN C C END MQERR C END SUBROUTINE FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) C C FUNC CALLS FUNK OR FOFX TO COMPUTE THE ARRAY OF FITTED C VALUES, FITM(*), FOR THE MARQ PACKAGE. . C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C DOUBLE PRECISION Y,YSIG,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * FLAMBD,FNU,RELDIF,RELMIN,RZERO DOUBLE PRECISION X,FIT,F,PHI,SIG C INTEGER JPT,JVARY,KALCP,KERFL,KFLAG,KORDIF,KW, * LEQU,MASK,MATRX,MAXIT,METHD,MAXSUB,MAXUPD,NFLAT, * NFMAX,NOREP,NPTS,NTRACE,NV,NXTRA C DIMENSION Y(1),YSIG(1),FIT(1) C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C RZERO=0.0D0 C IF(KALCP.NE.0) GO TO 10 C CALL FUNK (FIT) GO TO 30 C 10 DO 20 JPT=1,NPTS CALL FOFX(JPT,NV,X,F) FIT(JPT)=F 20 CONTINUE C 30 PHI=RZERO SIG=YSIG(1) C DO 60 JPT=1,NPTS IF(LEQU.EQ.0) SIG=YSIG(JPT) C C CHECK FOR AN ILLEGAL VALUE OF SIG. C IF(SIG.GT.RZERO) GO TO 50 WRITE(KW,40)LEQU,JPT,SIG 40 FORMAT(/' ERROR IN MARQ.... LEQU =',I2,6X,'JPT =', * I5,6X,'YSIG =',1PG13.5,' IS .LE. ZERO.') STOP C 50 PHI=PHI+((FIT(JPT)-Y(JPT))/SIG)**2 60 CONTINUE C RETURN C C END FUNC C END SUBROUTINE DERIV (JPT,FUNK,NPTS,FIT,FITSV,P,LPDMA,LPDMB) C C DERIV 4.3 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C DERIV APPROXIMATES THE JACOBIAN MATRIX P(*,*) FOR THE C MARQ PACKAGE, USING DIFFERENCE QUOTIENTS. C C P(J,K) IS THE PARTIAL DERIVATIVE OF FIT(J) WITH RESPECT TO C X(K) . C IF KORDIF.EQ.1, DERIV USES A NONCENTRAL DIFFERENCE FORMULA. C IF KORDIF.EQ.2, DERIV USES A CENTRAL DIFFERENCE FORMULA. C KORDIF.EQ.1 IS ABOUT TWICE AS FAST AS KORDIF.EQ.2, BUT LESS C ACCURATE. C DOUBLE PRECISION P,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * FLAMBD,FNU,RELDIF,RELMIN,RZERO DOUBLE PRECISION X,FIT,FITSV,DEL,ABSX,XSAVE,FX0,FX1, * XPLUS,XMINUS C INTEGER J,JPT,JVARY,JX,KALCP,KERFL,KFLAG,KORDIF, * KW,KX,LEQU,LPDMA,LPDMB,MASK,MATRX,MAXIT,METHD, * MAXSUB,MAXUPD,NFLAT,NFMAX,NOREP,NPTS,NTRACE,NV, * NXTRA C DIMENSION FIT(1),FITSV(1),P(LPDMA,LPDMB) C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C RZERO=0.0D0 C JVARY=0 C C SAVE FIT IF KALCP.GE.0 . C IF(KALCP.LT.0) GO TO 20 C DO 10 J=1,NPTS FITSV(J)=FIT(J) 10 CONTINUE C C LOOP OVER THE ACTIVE PARAMETERS X(JX). C 20 KX=0 C DO 160 JX=1,NV IF(MASK(JX).NE.0) GO TO 160 KX=KX+1 XSAVE=X(JX) ABSX=XSAVE IF(ABSX.LT.RZERO) ABSX=-ABSX DEL=RELDIF*ABSX IF(DEL.EQ.RZERO) DEL=RELDIF XPLUS=XSAVE+DEL X(JX)=XPLUS XMINUS=XSAVE-DEL IF(XPLUS.NE.XSAVE .AND. XPLUS.NE.XMINUS) GO TO 40 WRITE(KW,30)RELDIF,JX,XSAVE,DEL,XPLUS,XMINUS 30 FORMAT(/' ****** UNABLE TO COMPUTE A SUITABLE VALUE', * ' FOR DEL IN SUBROUTINE DERIV.'/6X,'RELDIF =', * 1PG15.7,6X,'JX =',I11,6X,'XSAVE =',G15.7,6X,'DEL =' * ,G15.7/6X,'XPLUS =',G15.7,6X,'XMINUS =',G15.7) STOP C 40 IF(KALCP.LT.0) GO TO 120 IF(KALCP.GT.0) GO TO 90 C C KALCP.EQ.0 . COMPUTE P(*,*), ONE COLUMN AT A TIME. C CALL FUNK (FIT) IF(KORDIF.EQ.2) GO TO 60 C DO 50 J=1,NPTS P(J,KX)=(FIT(J)-FITSV(J))/(XPLUS-XSAVE) 50 CONTINUE C GO TO 150 C C KALCP.EQ.0 AND KORDIF.EQ.2 . C IN THIS CASE, THE INPUT VALUES OF FIT(J) WILL BE DESTROYED. C 60 X(JX)=XMINUS C DO 70 J=1,NPTS FITSV(J)=FIT(J) 70 CONTINUE C JVARY=JX CALL FUNK (FIT) JVARY=0 C DO 80 J=1,NPTS P(J,KX)=(FITSV(J)-FIT(J))/(XPLUS-XMINUS) 80 CONTINUE C GO TO 150 C C KALCP.GT.0 . COMPUTE P(*,*), ONE ELEMENT AT A TIME. C 90 DO 110 J=1,NPTS CALL FOFX (J,NV,X,FX1) IF(KORDIF.EQ.2) GO TO 100 P(J,KX)=(FX1-FITSV(J))/(XPLUS-XSAVE) GO TO 110 C 100 X(JX)=XSAVE-DEL CALL FOFX (J,NV,X,FX0) P(J,KX)=(FX1-FX0)/(XPLUS-XMINUS) X(JX)=XPLUS 110 CONTINUE C GO TO 150 C C KALCP.LT.0 . C COMPUTE ONE ROW OF P(*,*), ONE ELEMENT AT A TIME. C 120 FITSV(1)=FIT(JPT) CALL FOFX (JPT,NV,X,FX1) IF(KORDIF.EQ.2) GO TO 130 P(1,KX)=(FX1-FITSV(1))/(XPLUS-XSAVE) GO TO 140 C 130 X(JX)=XMINUS CALL FOFX (JPT,NV,X,FX0) P(1,KX)=(FX1-FX0)/(XPLUS-XMINUS) C 140 FIT(JPT)=FITSV(1) C C RESTORE X(JX). C 150 X(JX)=XSAVE 160 CONTINUE C IF(KALCP.LT.0) RETURN C DO 170 J=1,NPTS FIT(J)=FITSV(J) 170 CONTINUE C RETURN C C END DERIV C END SUBROUTINE FOFX (JPT,NV,X,F) C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C THIS IS A DUMMY VERSION OF SUBROUTINE FOFX, FOR THE MARQ C PACKAGE. C A NON-DUMMY VERSION OF FOFX MAY BE USED (OPTIONALLY) TO C SUPPLY TO SUBROUTINE MARQ THE VALUES OF THE FUNCTION BEING C FITTED, INSTEAD OF USING A 'FUNK' SUBROUTINE TO DO THIS. C THE USE OF FOFX REQUIRES SUBSTANTIALLY MORE OVERHEAD TIME C DURING EXECUTION, BUT SAVES CONSIDERABLE STORAGE BY NOT C REQUIRING THAT THE JACOBIAN MATRIX P(*,*) BE STORED. C DOUBLE PRECISION X,F C INTEGER JPT,NV C DIMENSION X(20) C RETURN END SUBROUTINE CALCD (JPT,P,LPDMA,LPDMB) C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C THIS IS A DUMMY VERSION OF SUBROUTINE CALCD, FOR THE MARQ C FOR THE MARQ PACKAGE. C A NON-DUMMY VERSION OF CALCD MAY BE USED (OPTIONALLY) TO C SUPPLY TO SUBROUTINE MARQ ANALYTIC VALUES OF THE ELEMENTS C OF THE JACOBIAN MATRIX, INSTEAD OF APPROXIMATION THE C ELEMENTS USING FINITE DIFFERENCES. C HOWEVER, MOST USERS PREFER TO USE FINITE DIFFERENCES. C DOUBLE PRECISION P C INTEGER JPT,LPDMA,LPDMB C DIMENSION P(LPDMA,LPDMB) C RETURN END SUBROUTINE STSET C C STSET 3.2 DECEMBER 1991 C C STSET SETS SOME INPUT QUANTITIES TO DEFAULT VALUES, FOR C SUBROUTINES STEPIT, SIMPLEX, MARQ, STP, MINF, OR KAUPE. C C J. P. CHANDLER, DEPARTMENT OF COMPUTER SCIENCE, C OKLAHOMA STATE UNIVERSITY C C C USAGE..... C C CALL STSET. C THEN SET SOME INPUT QUANTITIES (NV, AT LEAST) AND RESET ANY C OF THOSE SET IN STSET (BETTER VALUES OF X(J), ETC.) BEFORE C CALLING STEPIT OR SIMPLEX OR MARQ, ETC. C C DOUBLE PRECISION XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ,FLAMBD, * FNU,RELDIF,RELMIN,RZERO,HUGE, X C INTEGER JVARY,JX,KALCP,KERFL,KFLAG,KORDIF,KW,LEQU,MASK, * MATRX,MAXIT,METHD,MAXSUB,MAXUPD,NFLAT,NFMAX,NOREP, * NTRACE,NV,NVMAX,NXTRA C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C C HUGE=1.0D30 C RZERO=0.0D0 C C KW = LOGICAL UNIT NUMBER OF THE PRINTER C KW=6 C C NVMAX IS THE MAXIMUM PERMISSIBLE VALUE OF NV, GIVEN THE C PRESENT DIMENSIONS OF ARRAYS. C NVMAX IS THE DIMENSION OF THE ARRAYS X(*), XMAX(*), XMIN(*), C DELTX(*), DELMIN(*), AND MASK(*). NVMAX IS ALSO THE FIRST C DIMENSION OF ERR(*,*). THE SECOND DIMENSION OF ERR(*,*) C IS NVMAX+1. C NVMAX=20 C C THE USER MUST SET NV AFTER CALLING STSET. C NV=-1 C NTRACE=0 NFMAX=1000000 MAXIT=50 MAXSUB=30 METHD=1 KALCP=0 LEQU=0 NFLAT=1 MATRX=105 NXTRA=0 FLAMBD=1.0D0 FNU=10.0D0 C KORDIF=1 RELDIF=1.0D-8 RELMIN=1.0D-6 C C FOR GREATER ACCURACY BUT LESS SPEED IN MARQ, SET... C C KORDIF=2 C RELDIF=O.4D-5 C RELMIN=1.0D-9 C DO 10 JX=1,NVMAX X(JX)=RZERO XMAX(JX)=HUGE XMIN(JX)=-HUGE DELTX(JX)=RZERO DELMIN(JX)=RZERO MASK(JX)=0 10 CONTINUE C RETURN C C END STSET C END SUBROUTINE RPLOT (LP,NPTS,Y,YSIG,LEQU,FIT,MLARGE) C C RPLOT 2.2 JANUARY 1992 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C RPLOT PRINTS A PLOT OF DATA AND FITTED VALUES, SIDE BY SIDE C WITH A PLOT OF NORMALIZED RESIDUALS. C C THE X AXIS RUNS FROM TOP TO BOTTOM ON THE PAPER, AND THE Y C AXIS RUNS FROM LEFT TO RIGHT. C THE POINTS ARE PLOTTED ONE PER LINE, REGARDLESS OF THE C VALUES OF THE ABSCISSAS. C C C INPUT QUANTITIES.... C C LP -- LOGICAL UNIT NUMBER OF THE PRINTER C C NPTS -- NUMBER OF POINTS C C Y(J) -- ORDINATE OF THE J-TH POINT C C YSIG(J) -- EXPECTED ERROR IN Y(J) C C LEQU -- =1 IF ALL YSIG(J) ARE EQUAL (IN THIS CASE C ONLY YSIG(1) IS REFERENCED), C C =0 OTHERWISE C C FIT(J) -- FITTED VALUE AT THE J-TH POINT C C MLARGE -- =0 TO PRINT IN 72 COLUMNS, C =1 TO PRINT IN 120 COLUMNS C C CHARACTER*1 KBL,KY,KF,KBOTH,KR,KI,KPLUS,KMINUS,KYLINE, * KRLINE C DOUBLE PRECISION Y,YSIG,FIT, * RZERO,UNITR,R15, * YMAX,YMIN,RMAX,RMIN,YY,F,SIG,R,DY,DR C INTEGER LP,NPTS,LEQU,MLARGE, * NR,NY,JPT,J,JYZERO,JRZERO,K,L,JY,JF,JR C DIMENSION Y(1),YSIG(1),FIT(1) DIMENSION KYLINE(81),KRLINE(21) C KBL=' ' KY='Y' KF='F' KBOTH='2' KR='R' KI='I' KPLUS='+' KMINUS='-' C RZERO=0.0D0 UNITR=1.0D0 R15=1.5D0 C MAXMSG=3 NR=21 C IF(MLARGE.EQ.0) THEN NY=41 ELSE NY=81 ENDIF C C COMPUTE THE EXTREMES OF Y(*) AND FIT(*), AND OF ABS(R(*)) . C YMAX=Y(1) YMIN=YMAX RMAX=RZERO JPT=1 NBADSG=0 C DO 10 J=1,NPTS YY=Y(J) IF(YY.GT.YMAX) YMAX=YY IF(YY.LT.YMIN) YMIN=YY F=FIT(J) IF(F.GT.YMAX) YMAX=F IF(F.LT.YMIN) YMIN=F IF(LEQU.EQ.0) JPT=J SIG=YSIG(JPT) C IF(SIG.EQ.RZERO) THEN SIG=UNITR NBADSG=NBADSG+1 IF(NBADSG.LE.MAXMSG) WRITE(LP,110)J,JPT 110 FORMAT(/' ERROR IN SUBROUTINE RPLOT ...'/ * 5X,'SIG = ZERO FOR J =',I6,' AND JPT =',I6/ * 5X,'SIG = 1.0 WILL BE USED INSTEAD.') IF(NBADSG.EQ.MAXMSG+1) WRITE(LP,115) 115 FORMAT(/' FURTHER ERROR MESSAGES SUPPRESSED.') ENDIF C R=(YY-FIT(J))/SIG IF(R.GT.RMAX) RMAX=R IF(-R.GT.RMAX) RMAX=-R 10 CONTINUE C IF(YMAX.LT.RZERO) YMAX=RZERO IF(YMIN.GT.RZERO) YMIN=RZERO C IF(YMAX.EQ.RZERO .AND. YMIN.EQ.RZERO) THEN YMAX=UNITR YMIN=-UNITR ENDIF C IF(RMAX.EQ.RZERO) RMAX=UNITR RMIN=-RMAX C C PRINT THE HEADING. C IF(MLARGE.EQ.0) THEN WRITE(LP,20) 20 FORMAT(/13X,'Y = DATA',8X,'F = FIT',14X, * 'NORMALIZED RESIDUALS') WRITE(LP,30)YMIN,YMAX,RMIN,RMAX 30 FORMAT(1X,1PE9.2,28X,E9.2,1X,E8.1,9X,E8.1) ELSE WRITE(LP,40) 40 FORMAT(31X,'Y = DATA',12X,'F = FIT',32X, * 'NORMALIZED RESIDUALS') WRITE(LP,50)YMIN,YMAX,RMIN,RMAX 50 FORMAT(1X,1PE9.2,68X,E9.2,1X,E8.1,9X,E8.1) ENDIF C DY=NY-1 DY=(YMAX-YMIN)/DY DR=NR-1 DR=(RMAX-RMIN)/DR JYZERO=-YMIN/DY+R15 C DO 150 JJ=1,2 C DO 60 K=2,NY KYLINE(K)=KMINUS 60 CONTINUE C KYLINE(1)=KPLUS KYLINE(NY)=KPLUS KYLINE(JYZERO)=KPLUS C IF(MLARGE.EQ.0) THEN WRITE(LP,70)(KYLINE(K),K=1,NY) 70 FORMAT(4X,41A1,5X,'+---------+---------+') ELSE WRITE(LP,80)(KYLINE(K),K=1,NY) 80 FORMAT(4X,81A1,5X,'+---------+---------+') ENDIF C IF(JJ.EQ.2) GO TO 150 C C PLOT THE TWO GRAPHS SIDE BY SIDE. C JRZERO=(NR+1)/2 C DO 90 K=1,NY KYLINE(K)=KBL 90 CONTINUE C DO 100 K=1,NR KRLINE(K)=KBL 100 CONTINUE C JPT=1 L=0 C DO 140 J=1,NPTS L=L+1 IF(L.GT.999) L=0 YY=Y(J) JY=(YY-YMIN)/DY+R15 IF(JY.LT.1) JY=1 IF(JY.GT.NY) JY=NY KYLINE(JYZERO)=KI F=FIT(J) JF=(F-YMIN)/DY+R15 IF(JF.LT.1) JF=1 IF(JF.GT.NY) JF=NY IF(LEQU.EQ.0) JPT=J SIG=YSIG(JPT) IF(SIG.EQ.RZERO) SIG=UNITR R=(YY-F)/SIG JR=JRZERO IF(DR.NE.RZERO) JR=(R-RMIN)/DR+R15 KYLINE(JY)=KY KYLINE(JF)=KF IF(JY.EQ.JF) KYLINE(JY)=KBOTH KRLINE(JRZERO)=KI KRLINE(JR)=KR C IF(MLARGE.EQ.0) THEN WRITE(LP,120)(KYLINE(K),K=1,NY),L, * (KRLINE(K),K=1,NR) 120 FORMAT(4X,41A1,I4,1X,21A1) ELSE WRITE(LP,130)(KYLINE(K),K=1,NY),L, * (KRLINE(K),K=1,NR) 130 FORMAT(4X,81A1,I4,1X,21A1) ENDIF C C RESET THE PRINT LINE TO BLANKS. C KYLINE(JF)=KBL KYLINE(JY)=KBL KRLINE(JR)=KBL C 140 CONTINUE C 150 CONTINUE C IF(MLARGE.EQ.0) THEN WRITE(LP,30)YMIN,YMAX,RMIN,RMAX ELSE WRITE(LP,50)YMIN,YMAX,RMIN,RMAX ENDIF C WRITE(LP,160) 160 FORMAT(' ') C RETURN C C END RPLOT C END SUBROUTINE GEM(NEQ,T,TF,H,HMAX,Y,ERR,RERR,FUNCT,NMAX,NUSED) IMPLICIT REAL*8(A-H,O-Z) C C GEM 3.2 A.N.S.I. STANDARD FORTRAN MAY 1980 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 T, 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,T,TF,H,HMAX,Y,ERR,RERR,FUNCT,NMAX C OUTPUT QUANTITIES.... T,H,Y,NUSED C C NEQ -- NUMBER OF EQUATIONS IN THE SYSTEM C T -- THE INDEPENDENT VARIABLE. SET T TO THE INITIAL C VALUE BEFORE THE FIRST CALL TO GEM. IT WILL C BE RETURNED EQUAL TO THE FINAL VALUE BY GEM. C TF -- THE DESIRED FINAL VALUE OF T 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 ACCEPTED IF FOR ALL J EITHER THE C ABSOLUTE ERROR PER STEP IN Y(J) IS .LE. C ERR(J) OR THE RELATIVE ERROR PER STEP IN C Y(J) IS .LE. RERR(J). C FUNCT -- THE NAME OF THE SUBROUTINE THAT EVALUATES THE C DERIVATIVES. CALL FUNCT(T,Y,DYDT) MUST RETURN C THE VECTOR OF DERIVATIVES DYDT(J) OF Y(J) WITH C RESPECT TO T. 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. C C EXAMPLE OF USAGE.... D2Y/DT2=-Y , Y(0)=0 , DY/DT(0)=1 C WRITE AS TWO FIRST ORDER EQUATIONS... DY(1)/DT=Y(2) , DY(2)/DT=-Y(1) C DIMENSION Y(2),ERR(2),RERR(2) C EXTERNAL TD C T=0.0 $ Y(1)=0.0 $ Y(2)=1.0 C ERR(1)=ERR(2)=RERR(1)=RERR(2)=1.0E-5 $ H=0.5 C DO 1 J=1,50 C CALL GEM(2,KN,0,T,0.5*J,H,1.0,Y,ERR,RERR,TD,1000,NUSED) C IF(NUSED.LT.0) CALL EXIT C D=SIN(T)-Y(1) C 1 PRINT 2,T,H,Y(1),D,NUSED C 2 FORMAT(' T =',1PE13.5,' H =',E13.5,' Y =',E15.7, C * ' ERROR =',E13.5,' NUSED =',I5) C END C SUBROUTINE TD (NEQ,T,Y,DYDT) C DIMENSION Y(2),DYDT(2) C DYDT(1)=Y(2) $ DYDT(2)=-Y(1) $ RETURN C END C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE FOLLOWING STATEMENTS CONVERT GEM TO DOUBLE PRECISION. C DOUBLE PRECISION T,TF,H,HMAX,Y,ERR,RERR,YT,FT,F DOUBLE PRECISION TT,FF,HT,HHMAX,HH,TEMP,ABHH,RATFC,ERT,ABERR, * ABER2,A,B,C,D,ALPHA,BETA,GAMMA,DELTA,RUNIT,RTWO,RTHRE, * FRAC,RATIO,UPFAC,DWNFC,SCALE,U,V,EPS,R,VMR,S,QT,CST DOUBLE PRECISION REIGH,RNINE,ZABS,ARG C DIMENSION Y(1),ERR(1),RERR(1) DIMENSION YT(40),FT(40),F(40,4) C ZABS(ARG)=DABS(ARG) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C 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 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=1.5 RUNIT=1.0 RTWO=2. RTHRE=3. C RATIO ... INITIAL VALUE FOR RATFC C (RTWO IS A SAFETY FACTOR) RATIO=RTWO*SCALE**4 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 QT=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,FAC,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 QT=(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 IQUIT=0 C MOVE THE NON-SUBSCRIPTED INPUT VARIABLES. TT=T FF=TF HT=H IF(TT-FF)10,410,10 10 MM=NEQ C CHECK VALIDITY OF INPUT VALUES. IF(MM)40,40,20 20 IF(MM-NEQMX)30,30,40 30 N=NMAX IF(N)40,40,50 40 NDONE=-2 GO TO 420 50 HHMAX=ZABS(HMAX) C HH=SIGN(AMIN1(ABS(H),HHMAX),TF-T) HH=ZABS(HT) IF(HH-HHMAX)70,70,60 60 HH=HHMAX 70 HT=HH IF(FF-TT)80,90,90 80 HH=-HH 90 ITCH=0 NDONE=0 RATFC=RATIO CALL FUNCT (MM,TT,Y,FT) NDONE=1 DO 100 K=1,MM 100 F(K,1)=FT(K) C START OF EACH STEP.... C CHECK THAT HH IS NOT NEGLIGIBLE. 110 TEMP=TT+HH IF(TEMP-TT)120,400,120 C SEE WHETHER WE CAN FINISH IN ONE STEP. C IF(ABS(HH)-0.99*ABS(TF-T)) 120 IF(ZABS(HH)-FRAC*ZABS(FF-TT))140,130,130 130 HH=FF-TT IQUIT=1 140 DO 150 K=1,MM 150 YT(K)=Y(K)+U*HH*F(K,1) CALL FUNCT (MM,TT+U*HH,YT,FT) DO 160 K=1,MM F(K,2)=FT(K) 160 YT(K)=Y(K)+HH*(VMR*F(K,1)+R*F(K,2)) CALL FUNCT (MM,TT+V*HH,YT,FT) DO 170 K=1,MM F(K,3)=FT(K) 170 YT(K)=Y(K)+HH*(CST*F(K,1)+S*F(K,2)+QT*F(K,3)) CALL FUNCT (MM,TT+HH,YT,FT) DO 180 K=1,MM F(K,4)=FT(K) 180 YT(K)=Y(K)+HH*(A*F(K,1)+B*F(K,2)+C*F(K,3)+D*F(K,4)) CALL FUNCT (MM,TT+HH,YT,FT) NDONE=NDONE+4 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. INK=1 DO 260 K=1,MM ERT=ALPHA*F(K,1)+BETA*F(K,2)+GAMMA*F(K,3)+DELTA*F(K,4)+EPS*FT(K) ERT=ZABS(HH*ERT) ABERR=ZABS(ERR(K)) ABER2=ZABS(RERR(K)*YT(K)) IF(ERT-ABERR)230,230,190 190 IF(ERT-ABER2)230,230,200 C C DECREASE THE STEP SIZE. 200 HH=HH/SCALE HT=ZABS(HH) IQUIT=0 C WAS THE STEP SIZE INCREASED LAST STEP.... IF(ITCH)220,220,210 C YES. RATFC IS TOO SMALL. INCREASE IT. 210 RATFC=UPFAC*RATFC C THE STEP SIZE IS BEING DECREASED (ITCH=-1). 220 ITCH=-1 GO TO 390 C CHECK WHETHER STEP SIZE CAN BE INCREASED. 230 IF(RATFC*ERT-ABERR)260,260,240 240 IF(RATFC*ERT-ABER2)260,260,250 C C NO IT CANNOT. COMPONENT K FAILS. 250 INK=0 260 CONTINUE C END OF TRUNCATION ERROR TESTS. C THE STEP WAS SUCCESSFUL. MOVE YT AND FT. DO 270 K=1,MM Y(K)=YT(K) 270 F(K,1)=FT(K) C WAS THIS THE FINAL STEP.... IF(IQUIT)280,280,410 C NO. INCREMENT T. 280 TT=TT+HH C ARE THE ERRORS SMALL ENOUGH TO INCREASE C THE STEP SIZE.... IF(INK)310,310,290 C YES. WAS IT DECREASED ON THE PREVIOUS C STEP.... 290 IF(ITCH)300,340,340 C YES. RATFC IS TOO SMALL. INCREASE IT. 300 RATFC=UPFAC*RATFC GO TO 330 C ALL IS WELL. DECREASE RATFC SLIGHTLY. C RATFC=AMAX1(DWNFC*RATFC,RATIO) 310 RATFC=DWNFC*RATFC IF(RATFC-RATIO)320,330,330 320 RATFC=RATIO C THE STEP SIZE IS NOT BEING CHANGED. 330 ITCH=0 GO TO 390 C THE STEP SIZE IS BEING INCREASED (ITCH=1). 340 ITCH=+1 C INCREASE THE STEP SIZE. C HH=SIGN(AMIN1(ABS(SCALE*HH),HHMAX),HH) ABHH=ZABS(SCALE*HH) IF(ABHH-HHMAX)360,360,350 350 ABHH=HHMAX 360 IF(HH)370,380,380 370 ABHH=-ABHH 380 HH=ABHH HT=ZABS(HH) C CHECK NDONE, AND START THE NEXT STEP. 390 IF(NDONE-N)110,110,400 C TOO MANY STEPS -- INTEGRATION UNSUCCESSFUL. 400 NDONE=-NDONE GO TO 420 C THE INTEGRATION WAS SUCCESSFUL. 410 TT=FF C RETURN MODIFIED VALUES OF T, H, AND NUSED. 420 T=TT H=HT NUSED=NDONE RETURN C C END GEM C END SUBROUTINE DINIT (NV,BD,JSET,CUSER,NCUSR,NEQ,TZERO,YDE) IMPLICIT REAL*8(A-H,O-Z) C C DINIT 1.2 FEBRUARY 1974 C C INITIALIZES TZERO 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 TZERO=0.0 DO 10 J=1,NEQ 10 YDE(J)=0.0 RETURN C C END DINIT C END SUBROUTINE FRAM2 (NSETS,X1,Y1,X2,Y2,NPTS,KHAR,KW) C C FRAM2 2.1 JULY 1992 C IMPLICIT REAL*8(A-H,O-Z) C C FRAMES ONE OR TWO SETS OF POINTS, AND CALLS PLOTUM TO PLOT C THE POINTS ON A LINE PRINTER OR A MONITOR. C CHARACTER*1 KAREA(23,53),KHAR(2) C DIMENSION X1(1),Y1(1),X2(1),Y2(1) C NCOLS=51 NROWS=21 MRIGHT=0 MMULTI=1 LDIM=23 C XMAX=X1(1) XMIN=XMAX YMAX=Y1(1) YMIN=YMAX DO 10 J=1,NPTS XMAX=DMAX1(XMAX,X1(J)) XMIN=DMIN1(XMIN,X1(J)) YMAX=DMAX1(YMAX,Y1(J)) YMIN=DMIN1(YMIN,Y1(J)) IF(NSETS.LE.1) GO TO 10 XMAX=DMAX1(XMAX,X2(J)) XMIN=DMIN1(XMIN,X2(J)) YMAX=DMAX1(YMAX,Y2(J)) YMIN=DMIN1(YMIN,Y2(J)) 10 CONTINUE C CALL PLOT2 (XMAX,XMIN,YMAX,YMIN,NCOLS,NROWS,MRIGHT, * MMULTI,KW,KAREA,LDIM) C CALL PLOT3 (KHAR(1),X1,Y1,NPTS,KAREA,LDIM) C CALL PLOT3 (KHAR(2),X2,Y2,NPTS,KAREA,LDIM) C CALL PLOT4 (KAREA,LDIM) C RETURN C C END FRAM2 (NEW) C END SUBROUTINE PLOT2 (XXRITE,XXLEFT,YYTOP,YYBOTM,NNX,NNY, * MRIGHT,MMULTI,LLP,KAREA,LDIM) C C PLOTUM 1.0 JANUARY 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C SUBROUTINES PLOT2, PLOT3, AND PLOT4 COMPRISE "PLOTUM", A C PORTABLE FORTRAN PACKAGE FOR PLOTTING GRAPHS ON A LINE C PRINTER OR ON THE SCREEN OF A TERMINAL OR PERSONAL COMPUTER. C C PLOTUM USES THE "STORED GRID" METHOD, ALLOWING SEVERAL C PLOTTING SYMBOLS TO BE USED IN ONE GRAPH. C C AS FAR AS POSSIBLE, WHILE MAINTAINING PORTABILITY, PLOTUM IS C PATTERNED AFTER THE "UM PLOT" PACKAGE WRITTEN AT THE C UNIVERSITY OF MICHIGAN AND FORMERLY DISTRIBUTED BY THE C "SHARE" ORGANIZATION OF IBM USERS. C C SUBROUTINE PLOT2 SETS UP THE PLOTTING AREA AND PUTS IN THE C BORDERS AND THE COORDINATE AXES. C C SUBROUTINE PLOT3 PUTS ONE SET OF DATA POINTS INTO THE PLOT C ARRAY. C C SUBROUTINE PLOT4 PRINTS OUT THE GRAPH. C C TO PRINT ONE GRAPH, CALL PLOT2 ONCE, CALL PLOT3 AS OFTEN AS C DESIRED (ONCE FOR EACH SET OF POINTS), AND CALL PLOT4 TO C PRINT THE GRAPH. C C XXRITE -- THE VALUE OF X AT THE RIGHTHAND SIDE OF THE C GRAPH C C XXLEFT -- THE VALUE OF X AT THE LEFTHAND SIDE OF THE C GRAPH C C YYTOP -- THE VALUE OF Y AT THE TOP OF THE GRAPH C C YYBOTM -- THE VALUE OF Y AT THE BOTTOM OF THE GRAPH C C XXRITE IS USUALLY GREATER THAN XXLEFT, AND YYTOP IS USUALLY C GREATER THAN YYBOTM, BUT THIS IS NOT AT ALL NECESSARY. C C NNX -- NUMBER OF INCREMENTS (PRINT POSITIONS) IN THE C X DIRECTION. C NNX SHOULD BE ONE GREATER THAN A MULTIPLE C OF 20 OR 25, FOR EXAMPLE C NNX=101 FOR A PRINTER PLOT, OR C NNX=51 FOR A PLOT ON A MONITOR. C C NNY -- NUMBER OF INCREMENTS (PRINT POSITIONS) IN THE C Y DIRECTION. C NNY SHOULD BE ONE GREATER THAN A MULTIPLE C OF 10, FOR EXAMPLE C NNY=51 FOR A PRINTER PLOT, OR C NNY=21 FOR A PLOT ON A MONITOR. C C MRIGHT -- =1 TO PUT A BORDER ON ALL FOUR SIDES OF THE C GRAPH, C =0 TO OMIT THE BORDER ON THE RIGHTHAND SIDE. C MRIGHT=0 PLOTS FASTER WHEN RUNNING VIA A C MODEM. C C MMULTI -- =1 TO PLOT A '2' WHEN TWO POINTS FALL IN THE C SAME PRINT POSITION, ETC. C ('9' SIGNIFIES "MORE THAN 8 POINTS".) C =0 TO PLOT ONLY THE LAST POINT INSERTED INTO C EACH PRINT POSITION C C LLP -- LOGICAL UNIT NUMBER FOR OUTPUT OF THE GRAPH C C KAREA(*,*) -- CHARACTER ARRAY FOR THE GRAPH. C THE FIRST DIMENSION OF KAREA(*,*), LDIM, C MUST BE AT LEAST AS LARGE AS NNY+2, AND C THE SECOND DIMENSION MUST BE AT LEAST AS C LARGE AS NNX+2 . C C LDIM -- THE FIRST DIMENSION OF THE ARRAY KAREA(*,*) C IN THE CALLING MODULE. C LDIM MUST BE AT LEAST AS LARGE AS NNY+2 . C C CHARACTER*1 KAREA,KBLANK,KPLUS,KMINUS,KI,NUM C DOUBLE PRECISION XXRITE,XXLEFT,YYTOP,YYBOTM, * XRITE,XLEFT,YTOP,YBOTM,DX,DY,RZERO,UNITR,R2P5 C INTEGER NNX,NNY,MRIGHT,MMULTI,LLP,LDIM, * NX,NY,NDIVX,NDIVY,MULTIP,LP, * NXPLU,NYPLU,MX,MY,J,K,JYZERO,KXZERO C DIMENSION KAREA(LDIM,1) C COMMON /CPLOTU/ * XRITE,XLEFT,YTOP,YBOTM,DX,DY, * NX,NY,NDIVX,NDIVY,MULTIP,LP, * NUM(9),KBLANK,KPLUS,KMINUS,KI C C COPY PARAMETER VALUES INTO COMMON. C XRITE=XXRITE XLEFT=XXLEFT YTOP=YYTOP YBOTM=YYBOTM NX=NNX NY=NNY MULTIP=MMULTI LP=LLP C C CHECK FOR BAD INPUT VALUES. C IF(NX.LT.2 .OR. NY.LT.2 .OR. LDIM.LT.NY+2) THEN WRITE(LP,10)NX,NY,LDIM 10 FORMAT(/' ERROR IN PLOT2 INPUT VALUES:'/5X, * 'NX =',I11,5X,'NY =',I11,5X,'LDIM =',I11) STOP ENDIF C C SET THE VALUES OF SOME CONSTANTS. C NDIVY=10 NDIVX=20 IF((NX-1)/NDIVX*NDIVX.NE.NX-1 .AND. * (NX-1)/25*25.EQ.NX-1) NDIVX=25 C KBLANK=' ' KPLUS='+' KMINUS='-' KI='I' NUM(2)='2' NUM(3)='3' NUM(4)='4' NUM(5)='5' NUM(6)='6' NUM(7)='7' NUM(8)='8' NUM(9)='9' C RZERO=0.0D0 UNITR=1.0D0 R2P5=2.5D0 C NXPLU=NX+1 NYPLU=NY+1 MX=NX+2 MY=NY+2 C C BLANK OUT THE PLOTTING AREA. C DO 30 K=1,MX C DO 20 J=1,MY KAREA(J,K)=KBLANK 20 CONTINUE C 30 CONTINUE C C ADJUST XRITE AND XLEFT IF THEY ARE EQUAL, C AND SIMILARLY FOR YTOP AND YBOTM. C IF(XRITE.EQ.XLEFT) THEN IF(XRITE.LT.RZERO) XRITE=RZERO IF(XLEFT.GT.RZERO) XLEFT=RZERO C IF(XRITE.EQ.XLEFT) THEN XRITE=UNITR XLEFT=-UNITR ENDIF C ENDIF C IF(YTOP.EQ.YBOTM) THEN IF(YTOP.LT.RZERO) YTOP=RZERO IF(YBOTM.GT.RZERO) YBOTM=RZERO C IF(YTOP.EQ.YBOTM) THEN YTOP=UNITR YBOTM=-UNITR ENDIF C ENDIF C DX=(XRITE-XLEFT)/(NX-1) DY=(YTOP-YBOTM)/(NY-1) C C PUT IN THE TOP BORDER OF THE GRAPH, C THE BOTTOM BORDER, AND THE X AXIS. C JYZERO=-YBOTM/DY+R2P5 KXZERO=-XLEFT/DX+R2P5 C DO 40 K=1,NX KAREA(1,K+1)=KMINUS KAREA(MY,K+1)=KMINUS IF(JYZERO.GE.2 .AND. JYZERO.LE.NYPLU) * KAREA(JYZERO,K+1)=KMINUS 40 CONTINUE C DO 50 K=1,NX,NDIVX KAREA(1,K+1)=KPLUS KAREA(MY,K+1)=KPLUS IF(JYZERO.GE.2 .AND. JYZERO.LE.NYPLU) * KAREA(JYZERO,K+1)=KPLUS 50 CONTINUE C C PUT IN THE LEFT BORDER, THE RIGHT BORDER, AND THE Y AXIS. C DO 60 J=1,NY KAREA(J+1,1)=KI IF(MRIGHT.EQ.1) KAREA(J+1,MX)=KI IF(KXZERO.GE.2 .AND. KXZERO.LE.NXPLU) * KAREA(J+1,KXZERO)=KI 60 CONTINUE C DO 70 J=1,NY,NDIVY KAREA(J+1,1)=KPLUS IF(MRIGHT.EQ.1) KAREA(J+1,MX)=KPLUS IF(KXZERO.GE.2 .AND. KXZERO.LE.NXPLU) * KAREA(J+1,KXZERO)=KPLUS 70 CONTINUE C RETURN C C END PLOT2 C END SUBROUTINE PLOT3 (KHAR,X,Y,NPTS,KAREA,LDIM) C C SUBROUTINE PLOT3 PUTS ONE SET OF POINTS INTO THE PLOTTING C ARRAY FOR THE PLOTUM PACKAGE. C C KHAR -- PLOTTING SYMBOL. C SUGGESTION: DO NOT USE '+', '-', OR 'I' FOR C THE PLOTTING SYMBOL. C IF MULTIP=1, DO NOT USE '2', '3', ETC. FOR C THE PLOTTING SYMBOL. C C X(J) -- ABSCISSA OF THE J-TH POINT IN THE SET C C Y(J) -- ORDINATE OF THE J-TH POINT IN THE SET C C NPTS -- NUMBER OF POINTS IN THIS SET C C CHARACTER*1 KHAR,KAREA,NUM,KBLANK,KPLUS,KMINUS,KI,KK C DOUBLE PRECISION X,Y,XRITE,XLEFT,YTOP,YBOTM, * DX,DY,R2P5,XMAX,XMIN,YMAX,YMIN C INTEGER NPTS,LDIM,NX,NY,NDIVX,NDIVY, * MULTIP,LP,J,K,NXPLU,NYPLU,JX,JY C DIMENSION X(1),Y(1),KAREA(LDIM,1) C COMMON /CPLOTU/ * XRITE,XLEFT,YTOP,YBOTM,DX,DY, * NX,NY,NDIVX,NDIVY,MULTIP,LP, * NUM(9),KBLANK,KPLUS,KMINUS,KI C R2P5=2.5D0 C XMAX=XRITE IF(XLEFT.GT.XMAX) XMAX=XLEFT XMIN=XLEFT IF(XRITE.LT.XMIN) XMIN=XRITE YMAX=YTOP IF(YBOTM.GT.YMAX) YMAX=YBOTM YMIN=YBOTM IF(YTOP.LT.YMIN) YMIN=YTOP C NXPLU=NX+1 NYPLU=NY+1 C DO 20 J=1,NPTS C IF(X(J).GE.XMIN .AND. X(J).LE.XMAX .AND. * Y(J).GE.YMIN .AND. Y(J).LE.YMAX) THEN C JX=(X(J)-XLEFT)/DX+R2P5 JY=(Y(J)-YBOTM)/DY+R2P5 C IF(JX.GE.2 .AND. JX.LE.NXPLU .AND. * JY.GE.2 .AND. JY.LE.NYPLU) THEN C C KK=KAREA(JY,JX) IF(KK.EQ.KBLANK .OR. KK.EQ.KPLUS .OR. * KK.EQ.KMINUS .OR. KK.EQ.KI .OR. * MULTIP.EQ.0) THEN KAREA(JY,JX)=KHAR C ELSE C DO 10 K=2,9 C IF(KK.EQ.NUM(K)) THEN IF(K.LE.8) KAREA(JY,JX)=NUM(K+1) GO TO 20 ENDIF C 10 CONTINUE C KAREA(JY,JX)=NUM(2) ENDIF C ENDIF C ENDIF 20 CONTINUE C RETURN C C END PLOT3 C END SUBROUTINE PLOT4 (KAREA,LDIM) C C SUBROUTINE PLOT4 PLOTS (PRINTS) A GRAPH FOR THE PLOTUM C PLOTTING PACKAGE. C C CHARACTER*1 KAREA,NUM,KBLANK,KPLUS,KMINUS,KI C DOUBLE PRECISION XX,XRITE,XLEFT,YTOP,YBOTM,DX,DY,YY C INTEGER LDIM,NX,NY,NDIVX,NDIVY,MULTIP,LP, * LXXDIM,MX,MY,JJ,J,JMOD,NXPLU,L,K C DIMENSION KAREA(LDIM,1) DIMENSION XX(11) C COMMON /CPLOTU/ * XRITE,XLEFT,YTOP,YBOTM,DX,DY, * NX,NY,NDIVX,NDIVY,MULTIP,LP, * NUM(9),KBLANK,KPLUS,KMINUS,KI C C LXXDIM IS THE DIMENSION OF THE ARRAY XX(*) ABOVE. C LXXDIM=11 C C PRINT THE GRAPH WITH THE Y VALUES. C MX=NX+2 MY=NY+2 C DO 30 JJ=1,MY J=MY+1-JJ JMOD=(J-2)-(J-2)/NDIVY*NDIVY C IF(JMOD.EQ.0) THEN YY=YBOTM+(J-2)*DY IF(JJ.EQ.2) YY=YTOP WRITE(LP,10)YY,(KAREA(J,K),K=1,MX) 10 FORMAT(1X,1PG12.5,1X,103A1) ELSE WRITE(LP,20)(KAREA(J,K),K=1,MX) 20 FORMAT(14X,103A1) ENDIF C 30 CONTINUE C C COMPUTE AND PRINT THE X VALUES. C NXPLU=NX+1 L=0 C DO 40 K=2,NXPLU,NDIVX L=L+1 XX(L)=XLEFT+(K-2)*DX IF(K.EQ.NXPLU) XX(L)=XRITE IF(L.GE.LXXDIM) GO TO 50 40 CONTINUE C 50 IF(NDIVX.EQ.25) THEN WRITE(LP,60)(XX(K),K=1,L) 60 FORMAT(11X,1PG12.5,4G25.5) ELSE WRITE(LP,70)(XX(K),K=1,L) 70 FORMAT(11X,1PG12.5,5G20.5) ENDIF C RETURN C C END PLOT4 C END //GO.SYSIN DD * GDH TRANSIENTS --- 25 JUNE 1971 --- 30 MIN CHARCOAL TREATMENT 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 .0 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. TERMINATE THIS CRICFM RUN. 0 // 000000000111111111122222222223333333333444444444455555555556666666666777 123456789012345678901234567890123456789012345678901234567890123456789012 SUBROUTINE CRCON C C CRCON 1.3 JULY 1992 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C THIS SUBROUTINE CONTROLS THE OPERATION OF THE CRICFM PACKAGE. C THIS VERSION OF CRCON IS FOR INTERACTIVE USE ONLY. C A DUMMY VERSION OF CRCON IS USED FOR BATCH RUNS. C IMPLICIT REAL*8 (A-H,O-Z) C EXTERNAL CRFUN C CHARACTER*1 KHINT,KHDF,KH C DOUBLE PRECISION YDE,FIT,BD,CUSER,FITJ,X,DHSAV C DIMENSION KH(2) DIMENSION TTEMA(300),YTEMP(300),TTEMB(300),FTEMP(300) C COMMON /CRICM/ T(300),DHSAV(10),ER,RER, * NSETS,JPTS(10),NEQFU,NEQDE,NMAX,NPRNT,LOGB(20),NFUPR,NACTV, * NVMAX COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS COMMON /CDFNC/ YDE(40),BD(20),CUSER(60),FITJ,NCUSR,JSET COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20),DELMN(20), * ERR(20,21),FOBJ,NV,NTRAC,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW COMMON /NLLS3/ FLAMB,FNU,RELDF,METHD,KALCP,KORDF,MAXIT,LEQU, * MXSUB,MXUPD COMMON /CHINT/ KHINT(10),KHDF(2) C C KEYBD ... LOGICAL UNIT NUMBER OF THE C TERMINAL KEYBOARD KEYBD=5 C KDISP ... LOGICAL UNIT NUMBER OF THE C TERMINAL OUTPUT DISPLAY KDISP=6 C NCOMD ... NUMBER OF LEGAL COMMANDS NCOMD=6 C KFITF ... INTERNAL FLAG FOR REFITTING KFITF=0 C RETYPE THE MENU. 10 WRITE(KDISP,20)(MASK(J),J=1,NV) 20 FORMAT(/' MASK(J)...'/(1X,I8,4I13)) WRITE(KDISP,30)(X(J),J=1,NV) 30 FORMAT(/' X(J)...'/(1X,5E13.5)) WRITE(KDISP,40) 40 FORMAT(//' MENU.... PLEASE TYPE'/ * ' 1 TO PLOT ALL DATA AND FITTED VALUES'/ * ' 2 TO PLOT ONE SET OF DATA AND FITTED VALUES'/ * ' 3 TO FIT THE DATA, AND THEN PLOT'/ * ' 4 TO CHANGE ONE MASK(J) VALUE'/ * ' 5 TO CHANGE ONE X(J) VALUE'/ * ' 6 TO TERMINATE THE JOB') READ(KEYBD,50)KH 50 FORMAT(2A1) CALL FRPIN (KH,2,KCHUZ) WRITE(KDISP,55)KCHUZ 55 FORMAT(/' YOU TYPED....',I11) IF(KCHUZ.GE.1 .AND. KCHUZ.LE.NCOMD) GO TO 70 WRITE(KDISP,60)KCHUZ 60 FORMAT(/1X,I3,5X,' ILLEGAL CHOICE') GO TO 10 70 GO TO (220,220,270,90,150,80),KCHUZ C C MENU CHOICE 6 ... TERMINATE THE JOB 80 STOP C MENU CHOICE 4 ... ALTER MASK(J) 90 WRITE(KDISP,100) 100 FORMAT(/' TYPE THE VALUE OF J FOR WHICH'/ * 5X,' YOU WISH TO CHANGE MASK(J).') READ(KEYBD,50)KH CALL FRPIN (KH,2,JX) WRITE(KDISP,110)JX 110 FORMAT(' J =',I3) IF(JX.LT.1 .OR. JX.GT.NV) GO TO 200 IF(MASK(JX).NE.0) GO TO 120 MASK(JX)=1 GO TO 130 120 MASK(JX)=0 130 WRITE(KDISP,140)JX,MASK(JX) 140 FORMAT(/' THE VALUE OF MASK(',I2,' HAS BEEN RESET TO ',I1) GO TO 10 C MENU CHOICE 5 ... ALTER X(J) 150 WRITE(KDISP,160) 160 FORMAT(' TYPE THE VALUE OF J FOR WHICH'/ * 5X,' YOU WISH TO CHANGE X(J).') READ(KEYBD,50)KH CALL FRPIN (KH,2,JX) WRITE(KDISP,110)JX IF(JX.LT.1 .OR. JX.GT.NV) GO TO 200 WRITE(KDISP,170)JX 170 FORMAT(/' TYPE (IN F20 FORMAT WITH A DECIMAL POINT)'/ * 5X,' THE NEW VALUE OF X(',I2,').') READ(KEYBD,180)X(JX) 180 FORMAT(F20.10) WRITE(KDISP,190)JX,X(JX) 190 FORMAT(/' THE VALUE OF X(',I2,') HAS BEEN RESET TO'/8X,E15.7) KFITF=0 GO TO 10 C 200 WRITE(KDISP,210)NV 210 FORMAT(//' ILLEGAL VALUE OF J. NV = ',I3) GO TO 10 C MENU CHOICE 1 OR 2 ... RECOMPUTE FIT(J) C AND THEN PLOT 220 IF(KFITF.NE.0) GO TO 230 CALL CRFUN (FIT) KFITF=1 230 IF(KCHUZ.EQ.1) GO TO 280 C MENU CHOICE 2 ... PLOT ONE SET OF DATA WRITE(KDISP,240) 240 FORMAT(/' WHICH SET OF DATA WOULD YOU LIKE TO PLOT...') READ(KEYBD,50)KH CALL FRPIN (KH,2,JSET) WRITE(KDISP,250)JSET 250 FORMAT(' JSET = ',I3) IF(JSET.GE.1 .AND. JSET.LE.NSETS) GO TO 280 WRITE(KDISP,260)NSETS 260 FORMAT(/' ILLEGAL VALUE OF JSET. NSETS = ',I3) GO TO 10 C MENU CHOICE 3 ... FIT, AND THEN PLOT 270 CALL STEPT (CRFUN) CALL CRFMB KFITF=1 C MENU CHOICES 1, 2, OR 3 ... PLOT 280 LPTS=0 JPT=0 DO 300 KSET=1,NSETS JMAX=JPTS(KSET) DO 290 J=1,JMAX JPT=JPT+1 IF(KCHUZ.EQ.2 .AND. KSET.NE.JSET) GO TO 290 LPTS=LPTS+1 TTEMA(LPTS)=T(JPT) YTEMP(LPTS)=Y(JPT) TTEMB(LPTS)=T(JPT) FTEMP(LPTS)=FIT(JPT) 290 CONTINUE 300 CONTINUE C CALL FRAM2 (2,TTEMA,YTEMP,TTEMB,FTEMP,LPTS,KHDF,KDISP) C C USE A "PAUSE" STATEMENT HERE ON SOME COMPUTER SYSTEMS, C IF DESIRED. C C PAUSE C GO TO 10 C C END CRCON (INTERACTIVE VERSION) C END 3 DATA FOR SIMULATED INTERACTIVE EXECUTION IN BATCH 6