C*************************************************************                                                       
C ****       TAXFIT - FOR LS FORTRAN   15-NOV-1996        ****
C ****   COPYRIGHT Michael Birnbaum, all rights reserved
C ****        MODIFIED 3-MAY-1998                         ****
C *** No warranties or guarantees.
C *** This program may be used for noncommercial, scholarly
C *** purposes without charge.  Publications should cite following:
C
C     Birnbaum, M. H., & Chavez, A. (1997). Tests of Theories of 
C       Decision Making: Violations of Branch Independence 
C       and Distribution Independence. Organizational Behavior 
C       and Human Decision Processes, 71(2), 161-194.
C     Birnbaum, M. H., & Stegner, S. E. (1979). Source credibility
C       in social judgment:  Bias, expertise, and the judge's 
C       point of view. Journal of Personality and Social Psychology,
C       37, 48-74.
C
C***  Note:  This version includes Chandler's STEPIT Subroutine.
C***  If you use STEPIT, also cite Chandler (1969).
C*************************************************************
C
C THIS PROGRAM WILL FIT SIMPLIFIED BIRNBAUM-STEGNER MODEL TO EXPT              
C WITH Judged Strength of pref TO GET ONE GAMBLE OR ANOTHER
C FOR  CHOICES BETWEEN TWO, THREE, OR FOUR OUTCOME GAMBLES  
C  MODEL ALLOWS ONE CONFIGURAL WEIGHT PARAMETER, HERE CALLED C
C  THIS VERSION FITS SINGLE S DATA AS A FUNCTION OF UTILITY DIFF  
C     THERE ARE TWO INDICES OF FIT: 
C              1. -LOG LIKELIHOOD OF CHOICES, GIVEN MODEL
C              2. SUM OF SQUARED DEVIATIONS between Data and Model
C  DATA FORMAT    (Each CARD is a line of 80 columns )                                               
C  CARD 1  2I3, F5.0      # OF TRIALS TO BE FIT, MODEL ID, 1-h 
C                    1 - h is the weight of -Log(likelihood| Model)
C                     h is the weight of sum of squared deviations
C             MODEL ID = 1 JUDGMENTS, MODEL ID = 2 is not active here
C       2  I3  No. of Subjects (I3)
C       3  5F5.2  NEXT CARDS (5F5.2)  INITIAL ESTIMATES, 
C
C  Parameters are GAMMA, BETA, A, B, and DELTA (here called C)
C  Warning: Your design may not allow separation of parameters
C    In that case, you should fix one or more parameters
C    A good parameter to fix is beta = 1, which has been
C       found as a reasonable approximation for small outcomes.
C    You should also constrain delta to avoid negative weights
C UNLIKE Birnbaum & Chavez, B is the scaling parameter for choice
C                           A is the scaling parameter for judgment
C
C       4-6   5F5.2 MASKS, MIN, MAX-These constrain parameter values
C             EACH SET STARTS ON A NEW CARD-Format is 5F5.2
C             MASK=0 is free, MASK=1 is fixed to initial estimate
C             XMIN =0 XMAX =0 for free parameters,
C       7    FORMAT for mean judgments in this file
C       8    FORMAT for raw data of subjects in file = SUBJECTS.DAT
C       9    FORMAT for output of predictions and residuals
C      10    2I3       # OF OUTCOMES LEFT, # OF OUTCOMES RIGHT       
C      11    (8F5.0)   PROBABILITIES OF LEFT AND RIGHT OUTCOMES
C      12    (8F5.0)  OUTCOMES OF LEFT AND RIGHT GAMBLES
C             NOTE: lowest outcomes are on the left, highest on right
C      13    (1F7.0)  DATA VALUE (E.G., MEAN JUDGMENT OF STR preference)
C
C  CARDS 10-13 REPEATED FOR # OF PROBLEMS                              
C      Note: This file MUST have the problems in the same order
C       as the order of judgments in the subjects data array.
C       Data coded: negative numbers = chose gamble on the left
C                   positive numbers = chose gamble on the right
C
C **************************************************************
C
      PROGRAM TAXFIT
      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(120,300),PD(120,300),RESD(120,300),
     1GAMMA(120),BETA(120),A(120),B(120),NPAIRS,WWLIK,
     2PX(300,5),PY(300,5),XO(300,5),YO(300,5),NOL(300),NOR(300),
     3W(5),C(110),CWT(5,5),UTL(300),UTR(300),MID,ISUB,CHISS,CHILK
      COMMON/FASTER/JVARY                                               
      COMMON/ICALL/ICALLS,IMAX
      DIMENSION FMTDAT(20),FMT2(20),FMT3(20)
	  DIMENSION XS(50),D1(300)
C
C *** FILES TO BE USED **********************************  FILES
C  
      OPEN (5, FILE = 'TAXFIT.DAT', STATUS = 'OLD')
      OPEN (6, FILE = 'TAXFIT.OUT', STATUS = 'NEW')
	  OPEN (7, FILE = 'SUBJECTS.DAT', STATUS = 'OLD')
	  OPEN (8, FILE = 'TAXFIT.PAR', STATUS = 'NEW')
	  OPEN (9, FILE = 'TAXFIT.RES', STATUS = 'NEW')
	  OPEN (10, FILE = 'TAXFIT.PRE', STATUS = 'NEW')
C                                                      
C  READ IN THE NUMBER OF TRIALS, GAMBLE PAIRS = NPAIRS
      READ(5,99)NPAIRS,MID,WWLIK  
	  READ(5,99)NSUBS
	  NSP1 = NSUBS +1
   99 FORMAT(2I3,F5.0)            
C NV IS THE NO.OF ESTIMATED PARAMETERS   
      NV=5                                
C NTRACE=1 MEANS PRINT OUT OF MINIMIZATION                             
C TO MAKE STEPIT HAPPY:                                                
      NTRACE=0                                                         
      MATRIX=110                                                       
C  READ IN INITIAL ESTIMATES OF THE X VALUES                           
C                           TRANSFER TABLE   
C                           GAMMA, BETA,  A,    B,    C = DELTA  
C                           X(1)   X(2)  X(3)  X(4)   X(5)    
      READ(5,151)(XS(J),J=1,NV)                                        
  151 FORMAT(10F5.2)                                                   
      WRITE(6,155)(XS(J),J=1,NV)                                       
  155 FORMAT(/,'         INITIAL ESTIMATES',8(/,'     ',7F9.2))        
C SOME PARAMETERS ARE FIXED OR CONSTRAINED                     
C TO FIX A PARAMETER, SET MASKS TO 1 --TO CONSTRAIN USE XMIN AND XMAX         
C       MASK = 0 FOR FREE PARAMETERS   XMIN = XMAX = 0 = UNCONSTRAINED  
C       MASK = 1 FOR FIXED PARAMETERS                                     
C  IF YOU WANT EVERYTHING FREE, LEAVE ALL ZEROS
      READ(5,161) (MASK(J),J=1,NV)
  161 FORMAT(5I5)
      READ(5,151) (XMIN(J),J=1,NV)
      READ(5,151) (XMAX(J),J=1,NV)
      READ(5,167) FMTDAT
	  READ(5,167) FMT2
	  READ(5,167) FMT3
  167 FORMAT(20A4)
C *************************** START DATA LOOP HERE *********

      DO 110 IP=1,NPAIRS                                              
      READ(5,99)NOL(IP),NOR(IP)                                           
C   NOL = NO OF OUTCOMES LEFT, NOR = NO OF OUTCOMES ON RIGHT
C READ IN PROBABILITIES
      NL = NOL(IP)
      NR = NOR(IP)
      READ(5,101)(PX(IP,J),J=1,NL),(PY(IP,K),K=1,NR)   
  101 FORMAT(8F5.0)                                                    
      READ(5,101)(XO(IP,J),J=1,NL),(YO(IP,K),K=1,NR)
      READ(5,FMTDAT) D1(IP)
C READ IN THE DATA AVERAGED DATA, D1(IP), USING SPECIFIED FORMAT
  102 FORMAT(F7.0)
  110 CONTINUE  
      TOTSS=0.0
	  TOTLK=0.0
C ******************************+++++++++++++++************ SUBJECT LOOP
C START SUBJECT LOOP HERE
C *********************************************************
      DO 990 ISUB=1,NSUBS
C READ SUBJECTS DATA
	  READ(7,FMT2) (D(ISUB,IP),IP=1,NPAIRS),ISID
C SET INITIAL VALUES FOR EACH SUBJECT
      DO 120 I =1,NV
	  X(I) = XS(I)
  120 CONTINUE
C ICALLS IS NO. OF CALLS TO FUNK, IMAX IS THE UPPER LIMIT FOR ICALLS
      ICALLS=0                                                          
      IMAX=2500   
C SET VALUES OF DELMIN AND DELTAX 
      DO 140I=1,NV                                                     
      DELMIN(I)=0.0                                                     
      DELTAX(I)=X(I)*(.01)                                              
  140 CONTINUE                                                                                                                    
C **********************************************************
      WRITE(6,141) ISUB
  141 FORMAT(' SUBJECT=',I5)
C NOW CALL STEPIT                                                     
      CALL STEPIT      
C **********************************************************
C PRINT OUT RESULTS 
C PRINTOUTS PER SUBJECT
C                        FILE OF PARAMETERS AND INDICES OF FIT
      WRITE(8,703) GAMMA(ISUB),BETA(ISUB),A(ISUB),B(ISUB),C(ISUB),
     1 CHISS,CHILK,ISID
  703 FORMAT(5F8.3,2F12.3,2X,I5)
C FIND WEIGHTS FOR EACH SUBJECT OF 4 EQUALLY LIKELY OUTCOMES
C CALCULATE W(P)
      TTT = 4*(.25)**GAMMA(ISUB)
	  V1=(.25)**GAMMA(ISUB)
	  V2=C(ISUB)
      W1 = V1*(1.0-(.6)*V2)/TTT
	  W2 = (V1-(.2)*C(ISUB)*V1)/TTT
	  W3 = (V1+(.2)*C(ISUB)*V1)/TTT
	  W4 = (V1+(.6)*C(ISUB)*V1)/TTT	
C PRINT WEIGHTS FOR EACH SUBJECT
      WRITE(8,723) W1,W2,W3,W4,ISID
  723 FORMAT(4F8.3,34X,I5)
C CUMULATE TOTAL INDICES OF FIT
	  TOTSS=TOTSS+CHISS
	  TOTLK=TOTLK+CHILK
C PRINT FILE OF RESIDUALS
	  WRITE(9,FMT3) (RESD(ISUB,IP),IP=1,NPAIRS),ISID
	  WRITE(10,FMT3) (PD(ISUB,IP),IP=1,NPAIRS),ISID	
C ****************************************************** 
  990 CONTINUE
C
C  ********************************* SUBJECT LOOP *******
C END PRINTOUT PER SUBJECT HERE ************************* SUBJECT LOOP ENDS
C
      WRITE(6,803) NSUBS,TOTSS,TOTLK
  803 FORMAT(' FIT OF BIRNBAUM-STEGNER MODEL FOR',I5,2X,' SUBJECTS',
     1 /, '  TOTAL SUM OF SQUARED DEVIATIONS =',F20.2,
     2 /, '  TOTAL NEGATIVE LOG LIKELIHOOD   =',F20.2)
C FIND AVERAGE DEVIATION
	  AVSS = TOTSS/(NSUBS*NPAIRS)
	  AVLK = TOTLK/(NSUBS*NPAIRS)
	  RMS = AVSS**.5
	  WRITE (6,825) AVSS,RMS,AVLK
  825 FORMAT(/,10X,'   AVERAGE SQUARED DEVIATION=',F10.3,
     1     /,10X,  ' ROOT MEAN SQUARED DEVIATION=',F10.3,
     2     /,10X,  ' AVERAGE NEG. LOG LIKELIHOOD=',F10.3)
C 
C *************************************************************
C NOW FIT THE AVERAGED DATA  --These results will appear as the last Subject
      NSP1 = NSUBS + 1
      ISUB = NSP1
	  DO 901 IP = 1,NPAIRS
	  D(ISUB,IP)=D1(IP)
  901 CONTINUE
C**************************-----------
C SET INITIAL VALUES 
      DO 920 I =1,NV
	  X(I) = XS(I)
  920 CONTINUE
C ICALLS IS NO. OF CALLS TO FUNK, IMAX IS THE UPPER LIMIT FOR ICALLS
      ICALLS=0                                                          
      IMAX=2500   
C SET VALUES OF DELMIN AND DELTAX 
      DO 940I=1,NV                                                     
      DELMIN(I)=0.0                                                     
      DELTAX(I)=X(I)*(.01)                                              
  940 CONTINUE                                                                                                                    
C **********************************************************
      WRITE(6,941) 
  941 FORMAT(' FIT FOR AVERAGED DATA')
C NOW CALL STEPIT                                                     
      CALL STEPIT      

C*************************----------------
      WRITE(6,801)                                                     
  801 FORMAT('1',T36,'FIT OF BIRNBAUM-STEGNER TAX MODEL-ONE OMEGA = C')               
C PRINT OUT GAMBLES FIT
      WRITE(6,804)
  804 FORMAT(T15,/,T15, ' GAMBLE ON THE LEFT    GAMBLE ON THE RIGHT') 
C----------BEGIN LOOP OF PAIRS HERE--------------
      DO 810 IP=1,NPAIRS                                                 
      NL = NOL(IP)
      NR = NOR(IP)
      WRITE(6,821) IP,NOL(IP),NOR(IP)
  821 FORMAT(/,T5,' TRIAL=',I4,'  NO OUTS LEFT=',I4,8X,
     + ' NO OUTS RIGHT=',I4)
      WRITE(6,805) (PX(IP,J),J=1,NL),(PY(IP,K),K=1,NR)
  805 FORMAT(,T10, ' PROBS=',2X,8F8.2)
      WRITE(6,806) (XO(IP,J),J=1,NL),(YO(IP,K),K=1,NR)
  806 FORMAT(,T10,' OUTS =',2X,8F8.2)
      WRITE(6,822) UTL(IP),UTR(IP)
  822 FORMAT(,T15, ' CE EQUIV LEFT=', F7.2,5X,' CE EQUIV RIGHT =',F7.2)
      WRITE(6,807) D(NSP1,IP),PD(NSP1,IP),RESD(NSP1,IP)
  807 FORMAT(T5,' DATA =',3X,F10.2,3X,' PRED =',3X,F10.2,3X,
     + ' RESID=',3X,F10.2)
  808 FORMAT(/,T20,' PRED =',10X,F10.2)
  809 FORMAT(/T20,' RESID=',10X,F10.2)
  810 CONTINUE
      WRITE(8,703) GAMMA(ISUB),BETA(ISUB),A(ISUB),B(ISUB),C(ISUB),
     1 CHISS,CHILK,NSP1
C FIND WEIGHTS FOR GROUP DATA OF 4 EQUALLY LIKELY OUTCOMES
C CALCULATE W(P)
      TTT = 4*(.25)**GAMMA(ISUB)
	  V1=(.25)**GAMMA(ISUB)
	  V2=C(ISUB)
      W1 = V1*(1.0-(.6)*V2)/TTT
	  W2 = (V1-(.2)*C(ISUB)*V1)/TTT
	  W3 = (V1+(.2)*C(ISUB)*V1)/TTT
	  W4 = (V1+(.6)*C(ISUB)*V1)/TTT	  
      WRITE(8,723) W1,W2,W3,W4,ISUB	
	  WRITE(9,FMT3) (RESD(ISUB,IP),IP=1,NPAIRS),ISUB
	  WRITE(10,FMT3) (PD(ISUB,IP),IP=1,NPAIRS),ISUB
C PRINT WEIGHTS FOR GROUP DATA INTO SUBJECT FILE SUBJECT
C 
C 
C IF WE GET THIS MESSAGE, ALL IS WELL
      WRITE(6,999)                                                      
  999 FORMAT(///,T15,'NORMAL END')                                     
      STOP                                                              
      END                                                               
C
C ********************************************************************
C ****                                                            ****
C ****           HERE IS SUBROUTINE   '  FUNK   '                 ****
C ****              Copyright Michael Birnbaum                    ****
C ********************************************************************
C
C HERE IS FUNK--WHICH FITS BIRNBAUM-STEGNER MODEL WITH ONE PARAMETER
C   TO JUDGMENTS of Strength of Preferences between pairs of Gambles
C    Negative Numbers mean preference for the Gamble on the Left
C    Positive Numbers mean preference for the Gamble on the Right
C    0 = No preference.
C              
      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(120,300),PD(120,300),RESD(120,300),
     1GAMMA(120),BETA(120),A(120),B(120),NPAIRS,WWLIK,
     2PX(300,5),PY(300,5),XO(300,5),YO(300,5),NOL(300),NOR(300),
     3W(5),C(110),CWT(5,5),UTL(300),UTR(300),MID,ISUB,CHISS,CHILK
      COMMON/FASTER/JVARY                                               
      COMMON/ICALL/ICALLS,IMAX                                         
C DECODE THE TRANSFER TABLE                                            
      GAMMA(ISUB) = X(1)
      BETA(ISUB) = X(2)
      A(ISUB) = X(3)
      B(ISUB) = X(4)
      C(ISUB) = X(5)
C FIND SUM OF SQUARED RESIDUALS, CHISS, AND NEGATIVE LOG LIKELIHOOD, CHILK   
      CHISQ = 0.0
	  CHISS = 0.0
	  CHILK = 0.0
C ************************************************************
C--------------IN THIS SECTION, WE CALCULATE
C       THE PREDICTION FOR EACH DATA POINT, FIND THE RESIDUAL
C       SQUARE THE RESIDUALS, AND ADD THEM UP INTO CHISS 
C--------------WE ALSO CALCULATE HERE NEGATIVE LOG LIKELIHOOD
C                  OF CHOICES GIVEN THE PARAMETERS AND MODEL
C *************************************************************
      DO 550 IP =1,NPAIRS
C  COMPUTE UTILITIES OF THE GAMBLES ON THE LEFT AND THE RIGHT
      NL = NOL(IP)
      NR = NOR(IP)
      SUML = 0.0
	  TTL = 0.0
C GAMBLE ON THE LEFT *****************************       LEFT
C FIRST LOOP== FIND SUMS OF WEIGHTS AND WEIGHT-UTILITY PRODS
      DO 520, I = 1,NL
	  W(I) = PX(IP,I)**GAMMA(ISUB)
      SUML = SUML + W(I)*XO(IP,I)**BETA(ISUB)
	  TTL = TTL + W(I)
  520 CONTINUE
C SECOND SET OF LOOPS==FIND SUM OF CONFIGURAL TERMS
	  SUMC = 0.0
      DO 525 I = 2,NL
	  JL = I - 1
	  DO 525 J =1,JL
	  TN = NL+1
	  XDIF = (XO(IP,I)**BETA(ISUB)) - (XO(IP,J)**BETA(ISUB))
	  IF (C(ISUB) .LE. 0.0) THEN
	  SUMC = SUMC + C(ISUB)*(XDIF/TN)*(PX(IP,I)**GAMMA(ISUB))
	  ELSE
	  SUMC = SUMC + C(ISUB)*(XDIF/TN)*(PX(IP,J)**GAMMA(ISUB))
	  END IF
  525 CONTINUE
      UTL(IP) = (SUML + SUMC)/TTL
C  
C GAMBLE ON THE RIGHT *****************************         RIGHT
C FIRST LOOP== FIND SUMS OF WEIGHTS AND WEIGHT-UTILITY PRODS
      SUMR = 0.0
	  TTR = 0.0
      DO 530, I = 1,NR
	  W(I) = PY(IP,I)**GAMMA(ISUB)
      SUMR = SUMR + W(I)*YO(IP,I)**BETA(ISUB)
	  TTR = TTR + W(I)
  530 CONTINUE
C SECOND SET OF LOOPS==FIND SUM OF CONFIGURAL TERMS
	  SUMC = 0.0
      DO 535 I = 2,NR
	  JL = I - 1
	  DO 535 J =1,JL
	  TN = NR+1
	  YDIF = YO(IP,I)**BETA(ISUB) - YO(IP,J)**BETA(ISUB)
	  IF (C(ISUB) .LE. 0.0) THEN
	  SUMC = SUMC + C(ISUB)*(YDIF/TN)*(PY(IP,I)**GAMMA(ISUB))
	  ELSE
	  SUMC = SUMC + C(ISUB)*(YDIF/TN)*(PY(IP,J)**GAMMA(ISUB))
	  END IF
  535 CONTINUE
      UTR(IP) = (SUMR + SUMC)/TTR
C ******************************************* CALCULATE PREDICTED RESPONSE
      PD(ISUB,IP) = A(ISUB)*(UTR(IP) - UTL(IP))
C
C Here WE CALCULATE SUM OF SQUARED DEVIATIONS AND SUM OF -LOG(PROBS)
      RESD(ISUB,IP) = D(ISUB,IP) - PD(ISUB,IP)
      CHISS = CHISS + RESD(ISUB,IP)**2
C 
      PPP =  B(ISUB)*(UTR(IP) - UTL(IP))
      PP = (1.0/(1.0+EXP(PPP)))
	  IF (D(ISUB,IP) .LE. 0.0) THEN
	     CHILK = CHILK - LOG(PP)
	  ELSE
	  	PP = 1.0 - PP
		CHILK = CHILK - LOG(PP)
	  END IF
C
  550 CONTINUE 
      CHISQ = (WWLIK)*CHILK +(1.0-WWLIK)*CHISS
      ICALLS=ICALLS+1                                                   
      RETURN                                                            
      END                                                               
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