C *******************************************************************
C **** ****
C **** RTFIT - MODIFIED FOR LS FORTRAN for MAC OCT-1997 ****
C **** ****
C *******************************************************************
C
C See BIRNBAUM, M. H., AND JOU,J. (1990).
C A Theory of Comparative Response Times and "Difference"
C Judgments. COGNITIVE PSYCHOLOGY, 22, 184-210.
C
C THIS PROGRAM WILL FIT CHOICE RT ("MORE") AND RT("LESS") TIMES
C AND "DIFFERENCE" JUDGMENTS
C This program was developed from RDFIT, which fits two data arrays,
C "Differences" and "Ratios"
C FOR FURTHER REFERENCE SEE BIRNBAUM JEP-G 1980, 109 (3), PP 304-319.
C
C DATA FORMAT INPUT SEQUENCE for FILE RDFIT.DAT
C SEQUENCE 1 I3 # OF PROBLEMS
C 2 2I3 # OF ROWS, # OF COLUMNS
C 3 (11F5.2) DIFFERENCE DATA - ONE LINE FOR EACH ROW
C 4 (11F5.2) REACTION TIME DATA - ONE LINE FOR EACH ROW
C 5 (11F5.0) WEIGHT MATRIX
C WT(I,J)=0. IF THE VALUE IS MISSING
C WT(I,J) not zero TO WEIGHT THE MATRIX; e.g., 1.
C 6 (13F5.2) SCALE VALUES SC(2),SC(3), ... , SC(NCOL),
C [SR(2),SR(3), ... , SR(NROW),] OPTIONAL
C AD,BD,BR,BRM,TRL(1),TRL(2),...TRM(1),TRM(2)...
C
C ITEMS 2-6 REPEATED FOR # OF PROBLEMS
C
C *******************************************************************
C
PROGRAM RTFIT
C
C IMPLICIT REAL*8(A-H,O-Z)
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
C
COMMON/TOFUNK/D(12,12),R(12,12),PRED(12,12),PRER(12,12),RESD(12,
+12),RESR(12,12),CHID,CHIR,SR(12),SC(12),AD,BD,CD,AR,BR,CR,SSD,
+SSR,NROW,NCOL,NSW,WT(12,12),RM(12,12),PRERM(12,12),RESRM(12,12),
+ARM,BRM,TRM(12),TRL(12),CHIRM,SSRM
C
COMMON/FASTER/JVARY
COMMON/ICALL/ICALLS,IMAX
C
C *** FILES TO BE USED ***
C
OPEN (5, FILE = 'RTFIT.DAT', STATUS = 'OLD')
OPEN (6, FILE = 'RTFIT.OUT', STATUS = 'NEW')
C
C READ IN NUMBER OF PROBLEMS
READ(5,99)NPROB
DO 998NP=1,NPROB
ICALLS=0
IMAX=25000
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(11F5.2)
100 CONTINUE
C READ IN REACTION TIMES (LESS)
DO 110I=1,NROW
READ(5,102)(R(I,J),J=1,NCOL)
102 FORMAT(11F5.0)
110 CONTINUE
C READ RT(MORE)
DO 111 I=1, NROW
READ (5,102)(RM(I,J),J=1,NCOL)
111 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('1 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(11F5.0)
117 FORMAT(' ',12F5.0)
WTOT=0.
C FIND MEAN AND SUM OF SQUARED DEVIATIONS
SUMD=0.0
SUMR=0.0
SUMRM=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
Z=R(I,J)
SUMR=SUMR+Z*WT(I,J)
SUMRM=SUMRM+WT(I,J)*RM(I,J)
120 CONTINUE
T=WTOT
C FIND MEANS,XBARD(DIFF) AND XBARR(RT)
XBARD=SUMD/T
XBARR=SUMR/T
XBARRM=SUMRM/T
C FIND SUM SQUARED DEVIATIONS,SSD AND SSR
SSD=0.0
SSR=0.0
SSRM=0.0
DO 130I=1,NROW
DO 130J=1,NCOL
SSD=SSD+WT(I,J)*(D(I,J)-XBARD)**2
Z=R(I,J)
SSR=SSR+WT(I,J)*(Z-XBARR)**2
SSRM=SSRM+WT(I,J)*(RM(I,J)-XBARRM)**2
130 CONTINUE
C NV IS THE NO.OF ESTIMATED PARAMETERS
NV=NCOL+2*NROW+3
C IF ASYMMETRIC DESIGN, ALSO ESTIMATE ROW SCALE VALUES
IF (NROW.EQ.NCOL) GOTO 133
NV=NV+NROW
133 CONTINUE
C INITIALIZE X ARRAY
DO 135 I=1, NV
X(I) = 1.0
135 CONTINUE
C IF ASYMMETRIC DESIGN THEN EST SCALE VALUES FOR ROWS SEPARATELY.
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 ADD ARM, BRM
C INITIAL ESTIMATES TO DATA
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, BR, BRM
C TRL(1)....,TRL(NROW),TRM(1),.... TRM(NROW)
C X(1) X(2) X(NCOL-1) X(NV)
READ(5,151)(X(J),J=1,NV)
151 FORMAT(9F8.2)
WRITE(6,155)(X(J),J=1,NV)
155 FORMAT(/,' INITIAL ESTIMATES',8(/,' ',11F7.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 BUT SOME PARAMETERS ARE FIXED AT LAST
C TO FIX A PARAMETER, SET MASK=VALUE
C MASK=0 FOR FREE PARAMETERS
C MASK not equal 0 FOR FIXED PARAMETERS
C READ MASK VALUES
READ(5,152)(MASK(J),J=1,NV)
152 FORMAT(9I8)
C MASK T VALUES POSITIVE
LIMV=NV-2*NROW+1
DO 153I=LIMV,NV
XMIN(I)=.000001
XMAX(I)=10000.0
153 CONTINUE
C NOW CALL STEPIT
CALL STEPIT
C PRINT OUT RESULTS
WRITE(6,801)
801 FORMAT(//,T3,'RTFIT-FITS MODEL OF BIRNBAUM & JOU (1990,
+COGNITIVE PSYCHOLOGY,22,184-210)')
802 FORMAT(//,T3,'FIT TO DIFFERENCE JUDGMENTS')
C PRINT OUT DIFFERENCES,PREDICTED,RESIDUALS
WRITE(6,804)NROW,NCOL
804 FORMAT(///,T3,'Difference JUDGMENTS for ',I3,' X',I3,' DESIGN',/)
DO 810I=1,NROW
WRITE(6,811)(D(I,J),J=1,NCOL)
811 FORMAT(//,T5,'DATA=',5X,11F9.2)
WRITE(6,812)(PRED(I,J),J=1,NCOL)
812 FORMAT(/,T5,'PRED=',5X,11F9.2)
WRITE(6,813)(RESD(I,J),J=1,NCOL)
813 FORMAT(/,T5,'RESIDUAL=',1X,11F9.2)
810 CONTINUE
WRITE(6,901)NROW,NCOL
901 FORMAT(///,T3,'Response TIMES-LESS 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
WRITE(6,981)
981 FORMAT(///,T3,'Response TIMES-MORE',/)
DO 989 I=1, NROW
WRITE(6,811)(RM(I,J),J=1,NCOL)
WRITE(6,812)(PRERM(I,J),J=1,NCOL)
WRITE(6,813)(RESRM(I,J),J=1,NCOL)
989 CONTINUE
C PRINT SUMMARY STATISTIC
WRITE(6,911)CHID,CHIR,CHIRM,CHISQ
911 FORMAT(//,T3,'PROPORTION OF VARIANCE DEVIATIONS',//,T22,
+'DIFF TASK = ',F8.4,/,T22,'LESS RT TASK= ',F8.4,/,T22,
+'MORE RT TASK= ',F8.4,/,T22,
+' SUM= ',F8.4)
WRITE(6,921)SSD,SSR,SSRM
921 FORMAT(///,T3,'SUMS OF SQUARES FOR',//,T22,' DIFFS =',F15.2,/,T22,
+' LESS RTS=',F15.2,/,T22,' MORE RTS=',F15.2)
WRITE(6,912)(SC(I),I=1,NCOL)
912 FORMAT(///,T10,'COLUMN SCALE VALUES',//,T15,11F9.3)
WRITE(6,929)(SR(I),I=1,NROW)
929 FORMAT(///,T10,'ROW SCALE VALUES',//,T15,11F9.3)
WRITE(6,913)AD,BD,BR,BRM
913 FORMAT(///,T3,'J-FUNCTION PARAMETERS',//,T15,'AD=',F9.3,/,T15,
+'BD=',F9.3,/,T15,'BR=',F9.3,/,
+T15,'BRM=',F9.3)
WRITE(6,955)
955 FORMAT(///,T3,'VALUES OF BIAS PARAMETERS',/)
WRITE(6,956)(TRL(J),J=1,NCOL)
WRITE(6,957)(TRM(J),J=1,NCOL)
956 FORMAT(//,T7,' LOWER BIAS',//,12(10X,F12.6,/))
957 FORMAT(//,T7,' UPPER BIAS',//,12(10X,F12.6,/))
WRITE(6,996)ICALLS,LIMV,NV
996 FORMAT(//,T15, ' ICALLS=', I8,//,T15,' LIMV=',I4,' NV=',I4)
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 Birnbaum-Jou model to
C Choice RTs and "difference" judgments
C Converted to Fit Birnbaum & Jou Model to "Difference" Judgments
C and Comparative RTs.
C
SUBROUTINE FUNK
C
C IMPLICIT REAL*8(A-H,O-Z)
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
C
COMMON/TOFUNK/D(12,12),R(12,12),PRED(12,12),PRER(12,12),RESD(12,
+12),RESR(12,12),CHID,CHIR,SR(12),SC(12),AD,BD,CD,AR,BR,CR,SSD,
+SSR,NROW,NCOL,NSW,WT(12,12),RM(12,12),PRERM(12,12),RESRM(12,12),
+ARM,BRM,TRM(12),TRL(12),CHIRM,SSRM
C
COMMON/FASTER/JVARY
COMMON/ICALL/ICALLS,IMAX
C
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
BR=X(L)
L=L+1
BRM=X(L)
DO 600 I=1, NROW
L=L+1
TRL(I)=X(L)
600 CONTINUE
DO 610 I=1,NROW
L=L+1
TRM(I)=X(L)
610 CONTINUE
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 RESPONSE TIMES
CHIRAT=0.0
DO 520I=1,NROW
DO 520J=1,NCOL
C PREDICTION BASED FOR "LESS" RTIMES
TEMP=(SC(J)-SR(I))
TEMP=DABS(TEMP)
IF (TEMP .EQ. 0) TEMP=TEMP+0.00001
PRER(I,J)=(TRL(I)*TRL(J))/TEMP+BR
RESR(I,J)=R(I,J)-PRER(I,J)
CHIRAT=CHIRAT+WT(I,J)*(RESR(I,J))**2
520 CONTINUE
CHIR=CHIRAT/SSR
CHIRATM=0.0
C PREDICTIONS FOR "MORE" TIMES
DO 530 I=1, NROW
DO 530 J=1, NCOL
TEMP=(SC(J)-SR(I))
TEMP=DABS(TEMP)
IF (TEMP .EQ. 0) TEMP=TEMP + .00001
PRERM(I,J)=(TRM(I)*TRM(J))/TEMP+BRM
RESRM(I,J)=RM(I,J)-PRERM(I,J)
CHIRATM=CHIRATM+WT(I,J)*(RESRM(I,J))**2
530 CONTINUE
CHIRM=CHIRATM/SSRM
CHISQ=CHIRM+CHID+CHIR
C
ICALLS=ICALLS+1
RETURN
END
C
C
C ******************************************************************
C **** ****
C **** THIS IS SUBROUTINE STEPIT ****
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