C ******************************************************************* C **** RDFIT Was written by Michael Birnbaum in 1977 C THIS PROGRAM WILL FIT EITHER THE 1 OR 2 OPERATION THEORY C TO THE RATIO AND DIFFERENCE EXPERIMENTS C FOR FURTHER REFERENCE SEE BIRNBAUM JEP-G 1980, 109 (3), PP 304-319 C C The program and method is described in the following: C C Birnbaum, M. H. (1980). Comparison of two theories C of "ratio" and "difference" judgments. C Journal of Experimental Psychology: General, 109, 304-319. C **** RDFIT - MODIFIED FOR VAX/VMS 1-OCT-1992 **** C **** thanks to Steven McCormick for assistance **** C ******************************************************************* C C DATA FORMAT C CARD 1 I3 # OF PROBLEMS C 2 I3 OPERATION (1 OR 2 OPERATION MODEL) C 3 2I3 # OF ROWS, # OF COLUMNS C 4 (15F6.2) DIFFERENCE DATA - ONE LINE FOR EACH ROW C 5 (15F6.2) RATIO DATA ONE LINE FOR EACH ROW C 6 (15F6.0) WEIGHT MATRIX C WT(I,J)=0. IF THE VALUE IS MISSING C WT(I,J)<>0. TO WEIGHT THE MATRIX C 7 (7F6.2) SCALE VALUES SC(2),SC(3), ... , SC(NCOL), C (8) [SR(2),SR(3), ... , SR(NROW),] for nonsymmetric C AD,BD,CD,AR,BR,CR(=M) C 8-9 I3 MASK FOR CR(=M) 0=FREE TO VARY; <>0=FIXED C CARDS 2-8 REPEATED FOR # OF PROBLEMS C C ******************************************************************* C PROGRAM RDFIT C IMPLICIT REAL*8(A-H,O-Z) <--some compilers may prefer IMPLICIT DOUBLEPRECISION(A-H,O-Z) C COMMON/STPIT/CHISQ,X(50),XMAX(50),XMIN(50),DELTAX(50),ERR(50,50) +,DELMIN(50),MASK(50),NV,NTRACE,MATRIX COMMON/TOFUNK/D(15,15),R(15,15),PRED(15,15),PRER(15,15),RESD(15, +125),RESR(15,15),CHID,CHIR,SR(15),SC(15),AD,BD,CD,AR,BR,CR,SSD, +SSR,NROW,NCOL,NSW,WT(15,15) COMMON/FASTER/JVARY COMMON/ICALL/ICALLS,IMAX C C *** FILES TO BE USED *** C OPEN (5, FILE = 'RDFIT.DAT', STATUS = 'OLD') OPEN (6, FILE = 'RDFIT.OUT', STATUS = 'NEW') C C READ IN NUMBER OF PROBLEMS READ(5,99)NPROB DO 998NP=1,NPROB ICALLS=0 IMAX=100000 C READ IN OPERATION (NSW=1 FOR 1 OPERATION MODEL ANALYSIS C 2 2 ) READ(5,99)NSW IF (NSW .GT. 2 .OR. NSW .LT. 1) STOP 'OPERATION ERROR' C READ IN NO.OF ROWS (NROW),NO.OF COLUMN (NCOL) 99 FORMAT(2I3) READ(5,99)NROW,NCOL C READ IN DIFFERENCES DO 100I=1,NROW READ(5,101)(D(I,J),J=1,NCOL) 101 FORMAT(15F6.2) 100 CONTINUE C READ IN RATIOS DO 110I=1,NROW READ(5,101)(R(I,J),J=1,NCOL) 110 CONTINUE C READ IN MATRIX OF MISSING VALUES C IF MISSING, WT(I,J)=0. C FOR EQUAL WEIGHT, LET W(I,J)=1. WRITE (6,112) NP 112 FORMAT(////' PROBLEM #',I3,/,/,' WEIGHT MATRIX') DO 115 I=1,NROW READ(5,116)(WT(I,J),J=1,NCOL) WRITE(6,117) (WT(I,J),J=1,NCOL) 115 CONTINUE 116 FORMAT(15F6.0) 117 FORMAT(' ',15F5.0) WTOT=0. C FIND MEAN AND SUM OF SQUARED DEVIATIONS SUMD=0.0 SUMR=0.0 DO 120I=1,NROW DO 120J=1,NCOL SUMD=SUMD+D(I,J)*WT(I,J) WTOT=WT(I,J)+WTOT C THIS OPTION WAS ADDED TO MINIMIZE THE SQUARED DEVIATIONS IN THE LOG Z=R(I,J) SUMR=SUMR+DLOG10(Z)*WT(I,J) 120 CONTINUE T=WTOT C FIND MEANS,XBARD(DIFF) AND XBARR(RATIO) XBARD=SUMD/T XBARR=SUMR/T C FIND SUM SQUARED DEVIATIONS,SSD AND SSR SSD=0.0 SSR=0.0 DO 130I=1,NROW DO 130J=1,NCOL SSD=SSD+WT(I,J)*(D(I,J)-XBARD)**2 Z=R(I,J) Z=DLOG10(Z) SSR=SSR+WT(I,J)*(Z-XBARR)**2 130 CONTINUE C NV IS THE NO.OF ESTIMATED PARAMETERS NV=NCOL+5 C IF ASYMMETRIC DESIGN THEN EST SCALE VALUES FOR ROWS SEPARATELY. IF(NROW.NE.NCOL)NV=NROW+NCOL+5 C NTRACE=1 MEANS PRINT OUT OF MINIMIZATION C TO MAKE STEPIT HAPPY: NTRACE=0 MATRIX=105 C READ IN INITIAL ESTIMATES OF THE X VALUES C SCALE VALUE ESTIMATES ARE ENTERED IN ARRAY X IN THE FOLLOWING ORDER C (FOR SYMMETRIC RATIO AND DIFFERENCE EXPERIMENTS) C SC(2), SC(3), . . . , SC(NCOL), AD, BD, CD, AR, BR, CR(=M) C X(1) X(2) X(NCOL-1) X(NV) READ(5,151)(X(J),J=1,NV) 151 FORMAT(7F6.2) WRITE(6,155)(X(J),J=1,NV) 155 FORMAT(/,' INITIAL ESTIMATES',8(/,' ',7F9.2)) C ALL PARAMETERS ARE FREE AT FIRST DO 140I=1,NV XMAX(I)=0.0 XMIN(I)=0.0 DELMIN(I)=0.0 DELTAX(I)=X(I)*(.01) MASK(I)=0.0 140 CONTINUE C LIMIT SCALE VALUES TO POSITIVE VALUES DO 160I=1,NV XMIN(I)=.001 XMAX(I)=10000 160 CONTINUE C BUT SOME PARAMETERS ARE FIXED AT LAST C TO FIX A PARAMETER, SET MASK=VALUE C MASK=0 FOR FREE PARAMETERS C MASK<>0 FOR FIXED PARAMETERS C TRY THESE ESTIMATES FOR M (=CR) C X(NV)=1.47 C VALUE OF M IS FIXED C MASK(NV)=1 C READ IN THE MASK FOR M (=CR) READ (5,99) MASK(NV) C VALUE OF CD IS MASKED (NOT USED IN EITHER 1 OR 2 OP MODEL) L=NV-3 MASK(L)=1.0 C NOW CALL STEPIT(FUNK) CALL STEPIT C PRINT OUT RESULTS WRITE(6,801) 801 FORMAT('1',T36,'FIT OF RATIO-DIFFERENCE THEORY') IF (NSW .EQ. 1) WRITE(6,802) 802 FORMAT(T36,'ONE OPERATION MODEL') IF (NSW .EQ. 2) WRITE(6,803) 803 FORMAT(T36,'TWO OPERATION MODEL') C PRINT OUT DIFFERENCES,PREDICTED,RESIDUALS WRITE(6,804)NROW,NCOL 804 FORMAT(T35,//,' DIFFERENCES FOR ',I3,' X',I3,' DESIGN') DO 810I=1,NROW WRITE(6,811)(D(I,J),J=1,NCOL) 811 FORMAT(//,T5,'DATA=',5X,15F9.2) WRITE(6,812)(PRED(I,J),J=1,NCOL) 812 FORMAT(/,T5,'PRED=',5X,15F9.2) WRITE(6,813)(RESD(I,J),J=1,NCOL) 813 FORMAT(/,T5,'RESIDUAL=',1X,15F9.2) 810 CONTINUE WRITE(6,901)NROW,NCOL 901 FORMAT(//,T38,'RATIOS FOR',I3,' X ',I3,' DESIGN') DO 910I=1,NROW WRITE(6,811)(R(I,J),J=1,NCOL) WRITE(6,812)(PRER(I,J),J=1,NCOL) WRITE(6,813)(RESR(I,J),J=1,NCOL) 910 CONTINUE C PRINT SUMMARY STATISTIC WRITE(6,911)CHID,CHIR,CHISQ 911 FORMAT(//,' PROPORTION OF VARIANCE DEVIATIONS',//,10X, +'DIFF TASK= ',F8.4,/,10X,'RATIO TASK=',F8.4,/,10X,'SUM=',7X,F8. +4) WRITE(6,921)SSD,SSR 921 FORMAT(//' SUMS OF SQUARES FOR',/,10X,' DIFFS=',F10.2,/,10X, +' RATIOS=',F10.2) WRITE(6,912)(SC(I),I=1,NCOL) 912 FORMAT(//,' COLUMN SCALE VALUES',/,10X,11F8.3) WRITE(6,929)(SR(I),I=1,NROW) 929 FORMAT(///,' ROW SCALE VALUES',/,10X,11F8.3) WRITE(6,913)AD,BD,CD,AR,BR,CR 913 FORMAT(///,' J-FUNCTION PARAMETERS',/,10X,'AD=',F9.3,/,10X, +'BD=',F9.3,/,10X,'CD=',F9.3,/,10X,'AR=',F9.3,/,10X,'BR=',F9.3,/, +10X,'CR(=M)=',F9.3) IF (MASK(NV) .NE. 0.0) WRITE(6,997) 997 FORMAT(' (CR=M IS FIXED)') 998 CONTINUE WRITE(6,999) 999 FORMAT(///,T15,'NORMAL END') STOP END C C ******************************************************************** C **** **** C **** HERE IS SUBROUTINE ' FUNK ' **** C **** **** C ******************************************************************** C C HERE IS FUNK--WHICH FITS RATIO-DIFFERENCE THEORIES FOR C THE ONE OR TWO OPERATION MODEL SUBROUTINE FUNK C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT DOUBLEPRECISION(A-H,O-Z) COMMON/STPIT/CHISQ,X(50),XMAX(50),XMIN(50),DELTAX(50),ERR(50,50) +,DELMIN(50),MASK(50),NV,NTRACE,MATRIX COMMON/TOFUNK/D(15,15),R(15,15),PRED(15,15),PRER(15,15),RESD(15, +125),RESR(15,15),CHID,CHIR,SR(15),SC(15),AD,BD,CD,AR,BR,CR,SSD, +SSR,NROW,NCOL,NSW,WT(15,15) COMMON/FASTER/JVARY COMMON/ICALL/ICALLS,IMAX C DEFINE SCALE VALUES SC(1)=1.0 DO 501I=2,NCOL L=I-1 SC(I)=X(L) 501 CONTINUE C THIS TIME ASSUME SYMMETRIC DESIGN DO 502I=1,NROW SR(I)=SC(I) 502 CONTINUE IF(NROW.EQ.NCOL)GOTO511 DO 512I=1,NROW L=L+1 SR(I)=X(L) 512 CONTINUE 511 CONTINUE C DEFINE J-FUNCTION PARAMETERS L=L+1 AD=X(L) L=L+1 BD=X(L) L=L+1 CD=X(L) L=L+1 AR=X(L) L=L+1 BR=X(L) L=L+1 CR=X(L) C DEFINE ROW SCALE VALUES C FIND SUM OF SQUARED RESIDUALS FOR DIFFERENCES CHIDIF=0.0 DO 510I=1,NROW DO 510J=1,NCOL PRED(I,J)=AD*(SC(J)-SR(I))+BD RESD(I,J)=D(I,J)-PRED(I,J) C IF WT(I,J)=0., THEN RESIDUAL=0. CHIDIF=CHIDIF+WT(I,J)*(RESD(I,J))**2 510 CONTINUE CHID=CHIDIF/SSD C FIND SUM OF SQUARED RESIDUALS FOR RATIOS CHIRAT=0.0 DO 520I=1,NROW DO 520J=1,NCOL C ONE OR TWO OPERATION MODEL?? IF (NSW .EQ. 2) GOTO 815 C PREDICTION BASED ON ONE OPERATION MODEL TEMP=(SC(J)-SR(I)) PRER(I,J)=AR*2.0**TEMP+BR GOTO 817 C PREDICTIONS BASED ON 2 OPERATION MODEL 815 CONTINUE PRER(I,J)=AR*(SC(J)/SR(I))**CR+BR 817 CONTINUE RESR(I,J)=R(I,J)-PRER(I,J) Z1=PRER(I,J) Z2=R(I,J) CHIRAT=CHIRAT+WT(I,J)*(DLOG10(Z1)-DLOG10(Z2))**2 520 CONTINUE CHIR=CHIRAT/SSR CHISQ=CHID+CHIR ICALLS=ICALLS+1 RETURN END C********************************************************************* C C ****************************************************************** C **** **** C **** THIS IS SUBROUTINE STEPIT **** C **** **** C ****************************************************************** C SUBROUTINE STEPIT C COPYRIGHT 1965 -- J. P. CHANDLER, PHYSICS DEPT., INDIANA UNIVERSITY. C C STEPIT 5.2 AUGUST 1967 C AVAILABLE FROM.... QUANTUM CHEMISTRY PROGRAM EXCHANGE C I.U. CHEMISTRY DEPT., BLOOMINGTON, INDIANA. C C MODIFIED FOR FORTRAN G OS/360 MVT C C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT DOUBLEPRECISION(A-H,O-Z) DIMENSION VEC(50),TRIAL(50),XSAVE(50),CHI(50),DX(50),SECOND(2,2) DIMENSION OLDVEC(50),SALVO(50),XOSC(50,15),CHIOSC(15),NFLAT(50) C COMMON/STPIT/CHISQ,X(50),XMAX(50),XMIN(50),DELTAX(50),ERR(50,50) +,DELMIN(50),MASK(50),NV,NTRACE,MATRIX COMMON/FASTER/JVARY COMMON/ICALL/ICALLS,IMAX NVMAX=50 MOSQUE=15 KW=61 RATIO=10.0 COLIN=0.99 NCOMP=5 ACK=2.0 SIGNIF=2.E8 HUGE=1.E37 C JVARY=0 40 IF(NV)290,290,50 50 NACTIV=0 DO 150I=1,NV IF(MASK(I))150,60,150 60 IF(SIGNIF*DABS(DELTAX(I))-DABS(X(I)))70,70,100 70 IF(X(I))90,80,90 80 DELTAX(I)=0.01 GOTO 100 90 DELTAX(I)=0.01*X(I) 100 IF(DELMIN(I))120,110,120 110 DELMIN(I)=DELTAX(I)/SIGNIF 120 IF(XMAX(I)-XMIN(I))130,130,140 130 XMAX(I)=HUGE XMIN(I)=-HUGE 140 NACTIV=NACTIV+1 VACD=DMIN1(XMAX(I),X(I)) X(I)=DMAX1(XMIN(I),VACD) 150 CONTINUE COMPAR=0.0 IF(NACTIV-1)160,190,180 160 DO 170J=1,NV 170 MASK(J)=0 GOTO 50 180 A=NACTIV SUB=2.0/(A-1.0) P=2.0*(1.0/DSQRT(A)/(1.0-0.5**SUB)-1.0) COMPAR=DMIN1(9.99D-1,DABS((1.0-(1.0-COLIN)**SUB)*(1.0+P*(1.- +COLIN)))) 190 IF(NTRACE)280,200,200 200 WRITE(6,210) 210 FORMAT( +'1ENTER SUBROUTINE STEPIT. COPYRIGHT 1965 J. P. CHANDLER, PHYSICS + DEPT., INDIANA UNIVERSITY.'//' INITIAL VALUES....'/) WRITE(6,220)(MASK(J),J=1,NV) 220 FORMAT(/' MASK = '10(I6,6X)/(4X10I12)) WRITE(6,230)(X(J),J=1,NV) 230 FORMAT(/' X = '10E12.4/(10X10E12.4)) WRITE(6,240)(XMAX(J),J=1,NV) 240 FORMAT(/' XMAX = '10E12.4/(10X10E12.4)) WRITE(6,250)(XMIN(J),J=1,NV) 250 FORMAT(/' XMIN = '10E12.4/(10X10E12.4)) WRITE(6,260)(DELTAX(J),J=1,NV) 260 FORMAT(/' DELTAX = '10E12.4/(10X10E12.4)) WRITE(6,270)(DELMIN(J),J=1,NV) 270 FORMAT(/' DELMIN = '10E12.4/(10X10E12.4)) 280 CALL FUNK NF=1 JOCK=1 290 IF(NTRACE)320,300,300 300 WRITE(6,310)NV,NACTIV,MATRIX,NCOMP,RATIO,ACK,COLIN,COMPAR,CHISQ 310 FORMAT(//' 'I3,' VARIABLES,'I3,' ACTIVE.'10X8HMATRIX =I4, +10X7H NCOMP=I2//' RATIO ='F5.1,10X5HACK =F5.1,10X7HCOLIN =F6.3, +10X8HCOMPAR =F6.3///' CHISQ ='E15.8) 320 IF(NV)2150,2150,330 330 IF(NTRACE)360,360,340 340 WRITE(6,350) 350 FORMAT(//60(1X1H*)//10X37HTRACE MAP OF THE MINIMIZATION PROCESS//) C 360 DO 370I=1,NV DX(I)=DELTAX(I) VEC(I)=0. DO 370J=1,20 370 ERR(I,J)=0. CHIOLD=CHISQ NOSC=0 C 380 NCIRC=0 NZIP=0 C C MAIN DO LOOP FOR CYCLING THROUGH THE VARIABLES. C FIRST TRIAL STEP WITH EACH VARIABLE IS SEPARATE. 390 NACK=0 DO 1350I=1,NV OLDVEC(I)=VEC(I) VEC(I)=0.0 TRIAL(I)=0.0 IF(MASK(I))400,410,400 400 VEC(I)=-0.0 NFLAT(I)=1 GOTO 1350 410 NACK=NACK+1 SAVE=X(I) IF(SIGNIF*DABS(DX(I))-DABS(X(I)))580,580,420 420 X(I)=SAVE+DX(I) JVARY=0 IF(JOCK)440,440,430 430 JOCK=0 JVARY=I 440 NFLAG=1 IF(X(I)-XMIN(I))460,460,450 450 IF(X(I)-XMAX(I))470,460,460 460 NFLAG=NFLAG+3 GOTO 490 470 CALL FUNK IF(ICALLS.LE.IMAX)GOTO10 RETURN 10 CONTINUE NF=NF+1 JVARY=I CHIME=CHISQ IF(CHISQ-CHIOLD)620,480,490 480 NFLAG=NFLAG+1 490 X(I)=SAVE-DX(I) IF(X(I)-XMIN(I))590,590,500 500 IF(X(I)-XMAX(I))510,590,590 510 CALL FUNK NF=NF+1 JVARY=I IF(CHISQ-CHIOLD)610,520,530 520 NFLAG=NFLAG+1 530 IF(NFLAG-3)540,580,590 540 IF((CHISQ-CHIME)*(CHIME-2.*CHIOLD+CHISQ))550,590,550 550 TRIAL(I)=.5*DX(I)*(CHISQ-CHIME)/(CHIME-2.*CHIOLD+CHISQ) VEC(I)=TRIAL(I)/DABS(DX(I)) NFLAT(I)=0 X(I)=SAVE+TRIAL(I) CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)560,570,570 560 CHIOLD=CHISQ JOCK=1 GOTO 600 570 TRIAL(I)=0.0 VEC(I)=0.0 GOTO 590 580 VEC(I)=-0.0 NFLAT(I)=1 590 X(I)=SAVE 600 NCIRC=NCIRC+1 IF(NCIRC-NACTIV)690,1430,1430 610 DX(I)=-DX(I) C C A LOWER VALUE HAS BEEN FOUND. HENCE THIS VARIABLE WILL CHANGE. C 620 NCIRC=0 DEL=DX(I) 630 CHIME=CHIOLD CHIOLD=CHISQ VEC(I)=VEC(I)+DEL/DABS(DX(I)) NFLAT(I)=0 TRIAL(I)=TRIAL(I)+DEL DEL=ACK*DEL SAVE=X(I) X(I)=SAVE+DEL IF(X(I)-XMIN(I))680,680,640 640 IF(X(I)-XMAX(I))650,680,680 650 CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)630,660,660 660 ZZZ=ACK*CHIME-(ACK+1.0)*CHIOLD+CHISQ IF(ZZZ.EQ.0.)GOTO661 CINDER=(0.5/ACK)*(ACK**2*CHIME-(ACK**2-1.0)*CHIOLD-CHISQ)/ZZZ GOTO 662 661 CINDER=0.0 662 X(I)=SAVE+CINDER*DEL CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)670,680,680 670 CHIOLD=CHISQ TRIAL(I)=TRIAL(I)+CINDER*DEL VEC(I)=VEC(I)+CINDER*DEL/DABS(DX(I)) GOTO 690 680 X(I)=SAVE 690 IF(NZIP-1)1340,700,700 700 IF(DABS(VEC(I))-ACK)750,710,710 710 DX(I)=ACK*DABS(DX(I)) VEC(I)=VEC(I)/ACK OLDVEC(I)=OLDVEC(I)/ACK DO 720J=1,MOSQUE 720 ERR(I,J)=ERR(I,J)/ACK IF(NTRACE)750,750,730 730 WRITE(6,740)I,DX(I) 740 FORMAT(' STEP SIZE'I3,' INCREASED TO 'E12.5) 750 SUMO=0.0 SUMV=0.0 DO 760J=1,NV SUMO=SUMO+OLDVEC(J)**2 760 SUMV=SUMV+VEC(J)**2 IF(SUMO*SUMV)1340,1340,770 770 SUMO=DSQRT(SUMO) SUMV=DSQRT(SUMV) COSINE=0.0 DO 780J=1,NV 780 COSINE=COSINE+(OLDVEC(J)/SUMO)*(VEC(J)/SUMV) IF(NZIP-1)1340,790,800 790 IF(NACK-NACTIV)1340,820,820 800 IF(NACK-NACTIV)820,810,810 810 IF(NZIP-NCOMP)820,830,830 820 IF(COSINE-COMPAR)1340,830,830 C C SIMON SAYS, TAKE AS MANY GIANT STEPS AS POSSIBLE... C 830 IF(NTRACE)860,860,840 840 WRITE(6,850)CHIOLD,(VEC(J),J=1,I) 850 FORMAT(' CHISQ ='E15.8,5X14HNO. OF STEPS =10F9.2/(42X10F9.2)) 860 NGIANT=0 NTRY=0 NRETRY=0 KL=1 NOSC=NOSC+1 IF(NOSC-MOSQUE)890,890,870 870 NOSC=MOSQUE DO 880K=2,MOSQUE CHIOSC(K-1)=CHIOSC(K) DO 880J=1,NV XOSC(J,K-1)=XOSC(J,K) 880 ERR(J,K-1)=ERR(J,K) 890 DO 900J=1,NV XOSC(J,NOSC)=X(J) 900 ERR(J,NOSC)=VEC(J)/SUMV CHIOSC(NOSC)=CHIOLD IF(NOSC-3)960,910,910 C C SEARCH FOR A PREVIOUS SUCCESSFUL GIANT STEP IN A DIRECTION MORE C NEARLY PARALLEL TO THE DIRECTION OF THE PROPOSED STEP THAN WAS THE C IMMEDIATELY PREVIOUS ONE. C 910 COXCOM=0.0 DO 920J=1,NV 920 COXCOM=COXCOM+ERR(J,NOSC)*ERR(J,NOSC-1) NAH=NOSC-2 930 NTRY=0 DO 950K=KL,NAH NRETRY=NAH-K COSINE=0.0 DO 940J=1,NV 940 COSINE=COSINE+ERR(J,NOSC)*ERR(J,K) IF(COSINE-COXCOM)950,950,970 950 CONTINUE 960 CHIBAK=CHI(I) GOTO 1020 970 NTRY=1 KL=K+1 IF(NTRACE)1000,1000,980 980 NT=NOSC-K WRITE(6,990)NT 990 FORMAT(/1X8H********5X34HPOSSIBLEOSCILLATIONWITHPERIODI2, +' DETECTED.') 1000 DO 1010J=1,NV SALVO(J)=TRIAL(J) 1010 TRIAL(J)=(X(J)-XOSC(J,K))/ACK CHIBAK=CHIOLD+(CHIOSC(K)-CHIOLD)/ACK C 1020 DO 1040J=1,NV XSAVE(J)=X(J) TRIAL(J)=ACK*TRIAL(J) IF(MASK(J))1040,1030,1040 1030 TRAAL=X(J)+TRIAL(J) TRAIL=DMIN1(TRAAL,XMAX(J)) X(J)=DMAX1(TRAIL,XMIN(J)) 1040 CONTINUE JOCK=0 JVARY=0 CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)1050,1080,1080 1050 CHIBAK=CHIOLD CHIOLD=CHISQ NGIANT=NGIANT+1 IF(NTRACE)1020,1020,1060 1060 WRITE(6,1070)CHISQ,(X(J),J=1,NV) 1070 FORMAT(' CHISQ ='E15.8/' X(I)....'/(10(1XE12.5))) GOTO 1020 C 1080 IF(NRETRY)1100,1100,1090 1090 IF(NGIANT)1150,1150,1100 1100 CINDER=(0.5/ACK)*(ACK**2*CHIBAK-(ACK**2-1.0)*CHIOLD-CHISQ)/(ACK* +CHIBAK-(ACK+1.0)*CHIOLD+CHISQ) DO 1120J=1,NV IF(MASK(J))1120,1110,1120 1110 CANDER=XSAVE(J)+CINDER*TRIAL(J) XIJ=DMIN1(CANDER,XMAX(J)) X(J)=DMAX1(XIJ,XMIN(J)) 1120 CONTINUE JOCK=0 JVARY=0 CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)1280,1130,1130 1130 IF(NGIANT)1170,1140,1170 1140 IF(NTRY)1150,1170,1150 1150 DO 1160J=1,NV TRIAL(J)=SALVO(J) 1160 X(J)=XSAVE(J) GOTO 1190 1170 DO 1180J=1,NV TRIAL(J)=TRIAL(J)/ACK 1180 X(J)=XSAVE(J) 1190 IF(NTRACE)1240,1240,1200 1200 WRITE(6,1210)CHIOLD,NGIANT 1210 FORMAT(/' CHISQ ='E15.8,' AFTER'I3,' GIANT STEPS.') WRITE(6,1220)(X(J),J=1,NV) 1220 FORMAT(' X(I)....'/(10(1XE12.5))) WRITE(6,1230) 1230 FORMAT(/) 1240 IF(NGIANT)1250,1250,1310 1250 IF(NRETRY)1260,1260,930 1260 IF(NTRY)1270,1330,1270 1270 NTRY=0 GOTO 960 C 1280 CHIOLD=CHISQ JOCK=1 IF(NTRACE)1310,1310,1290 1290 STEPS=FLOAT(NGIANT)+CINDER WRITE(6,1300)CHIOLD,STEPS 1300 FORMAT(/' CHISQ ='E15.8,' AFTER'F6.1,' GIANT STEPS.') WRITE(6,1220)(X(J),J=1,NV) WRITE(6,1230) 1310 IF(NTRY)1320,380,1320 1320 NOSC=0 GOTO 380 1330 NOSC=MAX0(NOSC-1,0) 1340 CHI(I)=CHIOLD 1350 CONTINUE C C ANOTHER CYCLE THROUGH THE VARIABLES HAS BEEN COMPLETED. C PRINT ANOTHER LINE OF TRACES. C IF(NTRACE)1370,1370,1360 1360 WRITE(6,850)CHIOLD,(VEC(J),J=1,NV) 1370 CONTINUE 1380 IF(NZIP)1420,1390,1420 1390 IF(NTRACE)1420,1420,1400 1400 WRITE(6,1220)(X(J),J=1,NV) WRITE(6,1410) 1410 FORMAT(' ') 1420 NZIP=NZIP+1 GOTO 390 C C A MINIMUM HAS BEEN FOUND. PRINT THE REMAINING TRACES. C 1430 IF(NTRACE)1450,1450,1440 1440 WRITE(6,850)CHIOLD,(VEC(J),J=1,I) 1450 IF(NTRACE)1470,1470,1460 1460 WRITE(6,1220)(X(J),J=1,NV) WRITE(6,1230) C C DECREASE THE SIZE OF THE STEPS FOR ALL VARIABLES. C 1470 CONTINUE 1480 NOSC=0 NGATE=1 DO 1520J=1,NV IF(MASK(J))1520,1490,1520 1490 IF(NFLAT(J))1500,1500,1520 1500 IF(DABS(DX(J))-DABS(DELMIN(J)))1520,1520,1510 1510 NGATE=0 1520 DX(J)=DX(J)/RATIO IF(NGATE)1530,1530,1600 1530 IF(NTRACE)1570,1570,1540 1540 WRITE(6,1550)(DX(J),J=1,NV) 1550 FORMAT(60(1X1H*)//' STEP SIZES REDUCED TO....'//(10(1XE12.5))) WRITE(6,1560) 1560 FORMAT(//) 1570 GOTO 380 1580 WRITE(6,1590)(DX(J),J=1,NV) 1590 FORMAT(' OPERAT. TERM.'//10(1X,E12.5)) C 1600 CHISQ=CHIOLD IF(NTRACE)1630,1610,1610 1610 WRITE(6,1620)NF 1620 FORMAT(//1X,I5,' FUNCTION COMPUTATIONS ') 1630 CONTINUE 1640 IF(IABS(MATRIX-100)-50)1650,1650,2160 1650 IF(NACTIV-NV)2160,1660,2160 C C COMPUTE THE STANDARD ERRORS AND THE CORRELATIONS. C 1660 FAC=RATIO**(MATRIX-100) DO 1680I=1,NV DX(I)=DABS(FAC*DX(I)) XSAVE(I)=X(I) JVARY=0 DO 1670J=1,2 X(I)=XSAVE(I)+DX(I) CALL FUNK NF=NF+1 JVARY=I SECOND(1,J)=CHISQ 1670 DX(I)=-DX(I) X(I)=XSAVE(I) 1680 ERR(I,I)=(SECOND(1,1)-2.0*CHIOLD+SECOND(1,2))/DX(I)**2 DO 1710I=2,NV IM=I-1 DO 1710J=1,IM DO 1700K=1,2 X(I)=XSAVE(I)+DX(I) JVARY=0 DO 1690L=1,2 X(J)=XSAVE(J)+DX(J) CALL FUNK NF=NF+1 JVARY=J SECOND(K,L)=CHISQ X(J)=XSAVE(J) 1690 DX(J)=-DX(J) X(I)=XSAVE(I) 1700 DX(I)=-DX(I) ERR(I,J)=0.25*(SECOND(1,1)-SECOND(1,2)-SECOND(2,1)+SECOND(2,2))/ +DABS(DX(I)*DX(J)) 1710 ERR(J,I)=ERR(I,J) IF(NTRACE)1770,1720,1720 1720 WRITE(6,1730) 1730 FORMAT('1SIZES OF INCREMENTS TO BE USED BELOW....') WRITE(6,1740)(DX(J),J=1,NV) 1740 FORMAT(/(10(1XE12.5))) WRITE(6,1750) 1750 FORMAT(/////' MATRIX OF THE SECOND PARTIAL DERIVATIVES....'/) DO 1760I=1,NV 1760 WRITE(6,1740)(ERR(I,J),J=1,I) 1770 DO 1780I=1,NV DO 1780J=1,I IF(ERR(I,J))1780,1790,1780 1780 CONTINUE GOTO 1810 1790 WRITE(6,1800) 1800 FORMAT(//// +' THE ABOVE MATRIX CONTAINS ONE OR MORE ZEROES. A LARGER VALUE OF + -MATRIX- SHOULD BE TRIED, TO SEE IF THEY ARE LEGITIMATE. ') C C INVERT THE MATRIX USING SYMINV2 (COMM. OF THE A.C.M. 6, P. 67). C 1810 DET=1.0 DETLOG=0. DO 1820J=1,NV 1820 SALVO(J)=1.0 DO 1970I=1,NV BIGAJJ=0.0 DO 1850J=1,NV IF(SALVO(J))1830,1850,1830 1830 IF(DABS(ERR(J,J))-BIGAJJ)1850,1850,1840 1840 BIGAJJ=DABS(ERR(J,J)) K=J 1850 CONTINUE IF(BIGAJJ)1870,1860,1870 1860 DET=0.0 GOTO 1980 1870 SALVO(K)=0.0 IF(DET.GT.1.E35.OR.DET.LT.1.E-35)GOTO1871 GOTO 1875 1871 WRITE(6,1872) 1872 FORMAT(///,'DETERMINANT COMPUTATION WILL CAUSE OVERFLOW.') GOTO 2160 1875 DET=DET*ERR(K,K) DETLOG=DETLOG+DLOG(DABS(ERR(K,K)))/2.303 TRIAL(K)=1.0/ERR(K,K) ERR(K,K)=0.0 XSAVE(K)=1.0 M=K-1 IF(M)1910,1910,1880 1880 DO 1900J=1,M XSAVE(J)=ERR(K,J) TRIAL(J)=ERR(K,J)*TRIAL(K) IF(SALVO(J))1860,1900,1890 1890 TRIAL(J)=-TRIAL(J) 1900 ERR(K,J)=0.0 1910 M=K+1 IF(M-NV)1920,1920,1960 1920 DO 1950J=M,NV XSAVE(J)=ERR(J,K) IF(SALVO(J))1860,1930,1940 1930 XSAVE(J)=-XSAVE(J) 1940 TRIAL(J)=-ERR(J,K)*TRIAL(K) 1950 ERR(J,K)=0.0 1960 DO 1970J=1,NV DO 1970K=J,NV 1970 ERR(K,J)=ERR(K,J)+XSAVE(J)*TRIAL(K) IF(DET)2000,1980,2020 1980 WRITE(6,1990) 1990 FORMAT(//// +' ERROR MATRIX IS SINGULAR. -MATRIX- SHOULD PROBABLY BE INCREASED +. '/////) GOTO 2150 2000 WRITE(6,2010) 2010 FORMAT(//// +' ERROR MATRIX IS NEGATIVE DEFINITE. -MATRIX- SHOULD PROBABLY BE +DECREASED.') 2020 IF(NTRACE)2050,2030,2030 2030 WRITE(6,2040)DET,DETLOG 2040 FORMAT(////' DETERMINANT OF ABOVE MATRIX = 'E12.5, E12.5) C +10X14HLOG10F(DET)=E12.5) C ****** THE ABOVE LINE WAS COMMENTED OUT TO MAKE THE C ****** COMPILIER HAPPY. IT BELONGS WITH LINE 2040 C *************************************************** C 2050 DO 2090I=1,NV DO 2060J=1,I ERR(I,J)=ERR(I,J)*2.0 2060 ERR(J,I)=ERR(I,J) IF(ERR(I,I))2070,2070,2090 2070 WRITE(6,2080)ERR(I,I) 2080 FORMAT(///' NEGATIVE OR ZERO MEAN SQUARE ERROR ENCOUNTERED...', +3X,E15.8/' -MATRIX- SHOULD PROBABLY BE DECREASED.'///) 2090 XSAVE(I)=DSIGN(DSQRT(DABS(ERR(I,I))),ERR(I,I)) IF(NTRACE)2160,2100,2100 2100 WRITE(6,2110) 2110 FORMAT(/////' STANDARD ERRORS....') WRITE(6,1740)(XSAVE(J),J=1,NV) WRITE(6,2120) 2120 FORMAT(/////' LOWER TRIANGLE OF THE CORRELATION MATRIX....'/) DO 2140I=2,NV IM=I-1 DO 2130J=1,IM 2130 TRIAL(J)=ERR(I,J)/DABS(XSAVE(I)*XSAVE(J)) 2140 WRITE(6,1740)(TRIAL(J),J=1,IM) 2150 WRITE(6,1620)NF C 2160 CONTINUE 2190 JVARY=0 CALL FUNK IF(NTRACE)2230,2200,2200 2200 WRITE(6,2210)(X(J),J=1,NV) 2210 FORMAT(///10X24HFINAL VALUES OF X(I)....//(7(1XE16.9))) WRITE(6,2220)CHISQ 2220 FORMAT(//' FINAL VALUE OF CHISQ = 'E15.8//) 2230 RETURN END