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