      PROGRAM QAQISH
C
C  Main program module 
C  
      CHARACTER FMT(2)*72,TITLE(2)*72,AM*4,LBL[ALLOCATABLE](:)*16
C
      REAL*8 W1A[ALLOCATABLE,HUGE](:,:),
     * W3A[ALLOCATABLE,HUGE](:,:),W2[ALLOCATABLE,HUGE](:),
     * UIJKA[ALLOCATABLE,HUGE](:),UIJKLA[ALLOCATABLE,HUGE](:)
      REAL*4 YX[ALLOCATABLE,HUGE](:,:)
C
      REAL*8 S1[ALLOCATABLE](:,:),S1A[ALLOCATABLE](:,:),
     * S3[ALLOCATABLE](:,:),RVAR[ALLOCATABLE](:,:),U[ALLOCATABLE](:),
     * V1[ALLOCATABLE](:),V2[ALLOCATABLE](:),RZ[ALLOCATABLE](:),
     * B[ALLOCATABLE](:),S2[ALLOCATABLE](:),TOL,TEMP1,TEMP2,TEMP3,PS,
     * PS1,W3I,W3J1,W3K1,C0,C1,C2,C3,C4,C5,C6,C8,C9,C10,C12,DX,Z,P1,P2,
     * P3,P4,U1,U2,U12,U3,U13,U23,U123,U4,U14,U24,U34,U124,U134,U234,
     * UK1J2,UK1K2
      REAL*4 Y1[ALLOCATABLE](:),XX[ALLOCATABLE](:),XX1[ALLOCATABLE](:),
     * X[ALLOCATABLE](:)
      INTEGER*4 OFF3[ALLOCATABLE](:),IJSYM,MRAM,IJK,NDIAG,L1,L2
      INTEGER*2 TYPTBL[ALLOCATABLE](:,:),REG[ALLOCATABLE](:,:),
     * SUM2[ALLOCATABLE](:),SUM3[ALLOCATABLE](:),OFF2[ALLOCATABLE](:),
     * IT1[ALLOCATABLE](:),SIZ[ALLOCATABLE](:),IBXN[ALLOCATABLE](:),
     * METH[ALLOCATABLE](:)
      LOGICAL USED[ALLOCATABLE](:),CONVERGE,LASTPASS,MONBFLAG
C
C The INDUA function is used to locate appropriate 3rd/4th product moment
C calculation values in the UIJKA and UIJKLA vectors
C
      INDUA(I1,I2,I3,N)=
     * (I2-I1-1)*(N-I1-1) - (I2-I1-1)*(I2-I1)/2 + I3-I1-1
C
C  Maxtype= maximum number of classes
C  Maxdimb= maximum number of parameters
C  Maxsize= maximum cluster size
C  Maxnclu= maximum number of clusters
C 
      MAXTYPE= 10
      MAXDIMB= 10
      MAXSIZE= 50
      MAXNCLU= 1000
C 
C  Call subroutines to read/store/output all specifications 
C     
      CALL GETSPEC1(TOL,TITLE,FMT,AM,MAXTYPE,MAXDIMB,NTYP,NDIMB,NDIMB1,
     * IGEE,MAXITER,MONBFLAG)
      ALLOCATE (LBL(NDIMB),B(NDIMB),USED(NDIMB),TYPTBL(0:NTYP,0:NTYP),
     * REG(NTYP+NTYP*(NTYP+1)/2,2*NDIMB+4),IBXN(2*NDIMB),XX(NDIMB1),
     * XX1(NDIMB1),X(NDIMB),Y1(MAXSIZE),IT1(MAXSIZE),SIZ(MAXNCLU),
     * STAT=IK)
      IF (IK.NE.0) THEN
        WRITE(*,'(/A)') ' PROGRAM TERMINATED: Array allocation 1 failed'
        WRITE(9,'(/A)') ' PROGRAM TERMINATED: Array allocation 1 failed'
        STOP
      ENDIF
      CALL GETSPEC2(LBL,B,NDIMB)
      CALL GETSPEC3(USED,TYPTBL,REG,IBXN,NTYP,NDIMB,NDIMB1)
C
C  Call subroutine to read data and create expanded rows in scratch data file;
C  cluster sizes also calculated and checked and checked in subroutine
C  
      OPEN(10,FORM='UNFORMATTED') 
      OPEN(11,FORM='UNFORMATTED')
      CALL GETDATA(FMT,TYPTBL,REG,XX,XX1,X,Y1,IT1,SIZ,NTYP,NDIMB,NDIMB1,
     * MAXSIZE,MAXNCLU,IMAX,NCLUS)
C
      DEALLOCATE (TYPTBL,USED,REG,XX,XX1,Y1,X,IT1,IBXN,STAT=IK)
      ALLOCATE (METH(IMAX),S1(NDIMB,NDIMB),S1A(NDIMB,NDIMB),
     * S3(NDIMB,NDIMB),S2(NDIMB),V2(NDIMB),IT1(NDIMB),STAT=IK)
      IF (IK.NE.0) THEN
        WRITE(*,'(/A)') ' PROGRAM TERMINATED: Array allocation 2 failed'
        WRITE(9,'(/A)') ' PROGRAM TERMINATED: Array allocation 2 failed'
        STOP
      ENDIF
C
C Initialize some values as needed
C
      DO 12 I=1,IMAX
   12 METH(I)= 0

      IF (AM.EQ.'VMBQ'.AND.IMAX.GT.2) THEN
        ALLOCATE(SUM2(IMAX-2),SUM3(IMAX-2),OFF2(IMAX-2),OFF3(IMAX-2),
     *   STAT=IK)
        IF (IK.NE.0) THEN
          WRITE(*,'(/A)')
     *     ' PROGRAM TERMINATED: Array allocation 3 failed'
          WRITE(9,'(/A)')
     *     ' PROGRAM TERMINATED: Array allocation 3 failed'
          STOP
        ENDIF
        SUM2(1)= 1
        SUM3(1)= 1
        DO 14 I=2,IMAX-2
        SUM2(I)= SUM2(I-1) + I
   14   SUM3(I)= SUM3(I-1) + SUM2(I)
      ENDIF
         
      ITER=0
      CONVERGE= (ITER.GT.0)
      LASTPASS= CONVERGE
      WRITE(9,'(A)')
C
C  Top of loop for convergence of estimates; ends with statement 114
C
      DO 114 WHILE (.NOT.LASTPASS)
      REWIND 10
      ITER= ITER+1
      IF (MONBFLAG) WRITE(9,'(A,7F9.4,(/7X,8F9.4))') ' Beta estimates:',
     * B
      DO 16 I=1,NDIMB
      S2(I)= 0.0
      DO 15 J=1,NDIMB
   15 S1A(I,J)= 0.0
      DO 16 J=I,NDIMB
      S1(I,J)= 0.0
   16 S3(I,J)= 0.0
C
C  Top of loop for individual cluster calculations; ends with statement 100
C
      DO 100 IC=1,NCLUS

      N = SIZ(IC)
      N1= N*(N+1)/2
      ALLOCATE (W1A(NDIMB,N1),YX(N1,NDIMB+1),U(N1),V1(N1),STAT=IK)
      IF (IK.NE.0) THEN
        WRITE(*,'(/2A,I4)') ' PROGRAM TERMINATED: Array allocation 4 ',
     *   'failed for cluster of size',N
        WRITE(9,'(/2A,I4)') ' PROGRAM TERMINATED: Array allocation 4 ',
     *   'failed for cluster of size',N
        STOP
      ENDIF

C 
C  This block of statements calculates the array W1.  It ends with the
C  DO 22 structure.
C
      DO 18 I=1,N1
   18 READ(10) (YX(I,J),J=1,NDIMB+1)

      J1=1
      K1=1
      DO 22 I=1,N1

C  calclp,calcpsi

      PS= 0.
      DO 20 J=1,NDIMB
   20 PS= PS + YX(I,J+1)*B(J)
      IF (ABS(PS).LT..000001) THEN
        PS= 1.0
      ELSE
        PS= EXP(PS)
      ENDIF
 
C calcu, calcr, calcg, calcg11, calcg21, calcg22

      IF (I.LE.N) THEN
        U(I)= PS/(1.0+PS)
        W3I = U(I)*(1.0-U(I)) 
        W3J1= 0.
        W3K1= 0.
      ELSE
        IF (K1.EQ.N) THEN
          J1= J1+1
          K1= J1+1
        ELSE
          K1=K1+1
        ENDIF
        IF (PS.EQ.1.0) THEN
          U(I)= U(J1)*U(K1)
          W3I = 0.
          W3J1= U(J1) * (1.0-U(J1)) * U(K1) 
          W3K1= U(K1) * (1.0-U(K1)) * U(J1) 
        ELSE
          PS1= PS - 1.
          TEMP1= 1.0 + (U(J1)+U(K1)) * PS1
          TEMP2= SQRT(TEMP1*TEMP1-4.*PS*PS1*U(J1)*U(K1))
          U(I)= .5 * (TEMP1-TEMP2) / PS1
          TEMP3= (TEMP1*(U(J1)+U(K1))-2.*U(J1)*U(K1)*(PS+PS1))/TEMP2
          TEMP3= 2.*(U(J1)+U(K1)-TEMP3)*PS1
          W3I = PS*(TEMP3-2.*(TEMP1-TEMP2)) / (4.*PS1*PS1) 
          W3J1= U(J1)*(1.-U(J1)) * 0.5*(TEMP2-TEMP1+2.*PS*U(K1))/TEMP2
          W3K1= U(K1)*(1.-U(K1)) * 0.5*(TEMP2-TEMP1+2.*PS*U(J1))/TEMP2 
        ENDIF
      ENDIF
      V1(I)= YX(I,1) - U(I)

C calcxtgt

      DO 22 J=1,NDIMB
   22 W1A(J,I)= W3I*YX(I,J+1) + W3J1*YX(J1,J+1) + W3K1*YX(K1,J+1)  

      DEALLOCATE (YX,STAT=IK)

C
C  If method VMBQ is used, determine sizes of arrays needed for 3rd/4th
C  product moment calculations
C
      MRAM= 0
      IF (AM.EQ.'VMBQ'.AND.N.GT.2) THEN
        OFF2(1)= 0
        OFF3(1)= 0
        DO 24 I=2,N-2
        OFF2(I)= OFF2(I-1) + SUM2(N-I)
   24   OFF3(I)= OFF3(I-1) + SUM3(N-I-1)
        IJK= MAX(OFF3(N-2),1)
        ALLOCATE (UIJKA(SUM3(N-2)),UIJKLA(IJK),STAT=IK)
        IF (IK.NE.0) THEN
          WRITE(*,'(/2A,I4)') ' PROGRAM TERMINATED: Array ',
     *     'allocation 5 failed for cluster of size',N
          WRITE(9,'(/2A,I4)') ' PROGRAM TERMINATED: Array ',
     *     'allocation 5 failed for cluster of size',N
          STOP
        ENDIF
        MRAM= SUM3(N-2) + IJK
C  
C  This block of statements calculates the values for the 3rd/4th product
C  moment calculations.  It ends with the DO 30 structure. 
C
        IJK= 0
        I1= 0
        DO 28 I=1,N-2
        U1= U(I)
        I1= I1 + N-I 
        J1= I1
        DO 28 J=I+1,N-1
        U2 = U(J)
        U12= U(I1+J)
        J1= J1 + (N-J)
        DO 28 K=J+1,N
        U3 = U(K)
        U13= U(I1+K)
        U23= U(J1+K)

        C0= 1. + U13 - U1 + U12 - U2 + U23 - U3 
        C1= U1 - U12 - U13
        C2= U2 - U12 - U23     
        C4= U3 - U13 - U23
        Z= (MIN(U12,U13,U23) + MIN(0.,C1,C2,C4)) * 0.5
        DX   = 9999.
        ITER1= 0
        DO WHILE (ITER1.LT.100.AND.ABS(DX).GT.1E-6)
        ITER1= ITER1+1
        DX= ((C0-Z)*(U12-Z)*(U13-Z)*(U23-Z) - (C1+Z)*(C2+Z)*(C4+Z)*Z) /
     *   ((U12+C0-2.*Z)*(U13-Z)*(U23-Z) + (C0-Z)*(U12-Z)*(U13+U23-2.*Z)  
     *   + (C1+C2+2.*Z)*(C4+Z)*Z + (C1+Z)*(C2+Z)*(C4+2.*Z))
        Z= Z + DX
        END DO

        IJK= IJK+1
   28   UIJKA(IJK)= Z   
C
        IJK= 0
        I1= 0
        DO 30 I=1,N-3
        U1= U(I)
        I1= I1 + N-I
        J1= I1
        DO 30 J=I+1,N-2
        U2 = U(J)
        U12= U(I1+J)
        J1= J1 + N-J
        K1= J1
        DO 30 K=J+1,N-1
        U3  = U(K)
        U13 = U(I1+K)
        U23 = U(J1+K)
        U123= UIJKA(OFF2(I)+INDUA(I,J,K,N))
        K1= K1 + N-K
        DO 30 L=K+1,N
        U4  = U(L)
        U14 = U(I1+L)
        U24 = U(J1+L)
        U34 = U(K1+L)
        U124= UIJKA(OFF2(I)+INDUA(I,J,L,N))
        U134= UIJKA(OFF2(I)+INDUA(I,K,L,N))
        U234= UIJKA(OFF2(J)+INDUA(J,K,L,N))

        C0 = 1. + U34 - U1 + U13 + U12 - U124 + U14 - U123 - U134 - U2 - 
     *   U3 + U24 + U23 - U234 - U4
        C1 = U1 - U13 - U12 + U124 - U14 + U123 + U134
        C2 = U2 + U124 - U12 + U123 - U24 + U234 - U23
        C3 = U12 - U124 - U123
        C4 = U3 - U34 - U13 + U123 - U23 + U134 + U234
        C5 = U13 - U123 - U134
        C6 = U23 - U123 - U234
        C8 = U4 - U34 - U14 - U24 + U124 + U134 + U234
        C9 = U14 - U124 - U134
        C10= U24 - U124 - U234
        C12= U34 - U134 - U234 
        Z= (MIN(U123,U124,U134,U234) + MIN(0.,C12,C10,C6,C9,C5,C3))*0.5
        DX   = 9999.
        ITER1 = 0  
        DO WHILE (ITER1.LT.100.AND.ABS(DX).GT.1E-6)
        ITER1= ITER1+1
        P1= (C0+Z)*(C3+Z)*(C5+Z)*(C6+Z)   
        P2= (C9+Z)*(C10+Z)*(C12+Z)*Z
        P3= (C1-Z)*(C2-Z)*(C4-Z)*(U123-Z)
        P4= (C8-Z)*(U124-Z)*(U134-Z)*(U234-Z)
        DX= (P3*P4 - P1*P2) /
     *   (P1*((C10+C9+2.*Z)*(C12+Z)*Z + (C9+Z)*(C10+Z)*(C12+2.*Z)) +
     *    P2*((C3+C0+2.*Z)*(C5+Z)*(C6+Z) + (C0+Z)*(C3+Z)*(C6+C5+2.*Z)) +
     *    P3*((U124+C8-2.*Z)*(U134-Z)*(U234-Z) +
     *    (C8-Z)*(U124-Z)*(U134+U234-2.*Z)) +
     *    P4*((C2+C1-2.*Z)*(C4-Z)*(U123-Z) +
     *    (C1-Z)*(C2-Z)*(U123+C4-2.*Z)))
        Z= Z + DX
        END DO

        IJK= IJK+1
   30   UIJKLA(IJK)=Z

      ENDIF
C
C  Computational structures for rest of loop depends on RAM available
C  and size of cluster.  One of two methods will be used, if possible.
C  The difference is that the first method keeps everything in RAM
C  while the second method makes use of a disk scratch file.  The code
C  for the method used for a cluster size is stored in METH.  If the METH
C  element is zero, the method must be calculated for that cluster size.
C
C  Modules called based on METH values cover procedures calcb, calcb11,
C  calcb12, and calcb22.
C
      IF (METH(N).EQ.0) THEN
        ALLOCATE (W2(N1*(N1+1)/2 + NDIMB*N1 - MRAM),STAT=IK)
        IF (IK.EQ.0) THEN
          METH(N)=1
        ELSE
          ALLOCATE (W2(N1*(N1+1)/2),STAT=IK)
          IF (IK.EQ.0) THEN
            METH(N)=2
          ELSE
            WRITE(*,'(/2A,I4)') ' PROGRAM TERMINATED: Array ',
     *       'allocation 6 failed for cluster of size',N
            WRITE(9,'(/2A,I4)') ' PROGRAM TERMINATED: Array ',
     *       'allocation 6 failed for cluster of size',N
            STOP
          ENDIF
        ENDIF
        DEALLOCATE (W2,STAT=IK)
      ENDIF

C

      ALLOCATE (W2(N1*(N1+1)/2),STAT=IK)
      IJK= 0

      IF (N.EQ.1) THEN
        W2(1)= 1./(U(1)*(1.-U(1)))
        GOTO 76
      ENDIF

      K=N
      II=0
      DO 44 I=1,N

      IJK=IJK+1
      W2(IJK)= U(I)*(1.-U(I))
      DO 35 J=I+1,N
      K= K+1
      IJK=IJK+1
   35 W2(IJK)= U(K)-U(I)*U(J)
     
      IF (IGEE.EQ.1) THEN
        DO 37 J=N+1,N1
        IJK=IJK+1
   37   W2(IJK)= 0.
        GOTO 44
      ENDIF

      J1=1
      K1=2
      II=II + N-I
      JJ= N-1
      KK= 2*N-3
      DO 43 J=N+1,N1
      IF (I.EQ.J1.OR.I.EQ.K1) THEN
        TEMP1= U(JJ+ K1) * (1.-U(I))
      ELSEIF (AM.EQ.'VMZP') THEN
        IF (I.LT.J1) THEN       
          TEMP1= (U(II+J1) - U(I)*U(J1)) * U(K1)
     *     + (U(II+K1) - U(I)*U(K1)) * U(J1)
        ELSEIF (I.LT.K1) THEN
          TEMP1= (U(JJ+I) - U(I)*U(J1)) * U(K1)
     *     + (U(II+K1) - U(I)*U(K1)) * U(J1)
        ELSE
          TEMP1= (U(JJ+I) - U(I)*U(J1)) * U(K1) 
     *     + (U(KK+I) - U(I)*U(K1)) * U(J1)
        ENDIF
      ELSEIF (I.LT.J1) THEN
        TEMP1= UIJKA(OFF2(I)+INDUA(I,J1,K1,N)) - U(I)*U(JJ+K1)
      ELSEIF (I.LT.K1) THEN
        TEMP1= UIJKA(OFF2(J1)+INDUA(J1,I,K1,N)) - U(I)*U(JJ+K1)
      ELSE
        TEMP1= UIJKA(OFF2(J1)+INDUA(J1,K1,I,N)) - U(I)*U(JJ+K1)
      ENDIF
      IF (K1.EQ.N) THEN
        J1=J1+1
        K1=J1+1
        JJ=JJ + N-J1
        KK=JJ + N-K1
      ELSE
        K1=K1+1
        KK=KK + N-K1
      ENDIF
      IJK=IJK+1
   43 W2(IJK)= TEMP1
   44 CONTINUE
C
      J1=1
      K1=2
      JJ= N-1
      KK= 2*N-3
      DO 48 I=N+1,N1-1

      IJK=IJK+1
      W2(IJK)= U(I)*(1.-U(I))
      
      J2= J1
      K2= K1
      LL= JJ
      DO 46 J=I+1,N1
      IF (K2.EQ.N) THEN
        J2=J2+1
        K2=J2+1
        LL= LL + N-J2
      ELSE
        K2=K2+1
      ENDIF
      IF (K1.EQ.J2) THEN
        IF (AM.EQ.'VMZP') THEN
          TEMP1= (U(KK+K2) - U(K1)*U(K2)) * U(J1) +
     *           (U(JJ+K1) - U(J1)*U(K1)) * U(K2) +
     *           U(JJ+K2)*U(K1) - U(JJ+K1)*U(KK+K2)
        ELSE
          TEMP1= UIJKA(OFF2(J1)+INDUA(J1,K1,K2,N)) - U(JJ+K1)*U(KK+K2)
        ENDIF   
      ELSEIF (K1.EQ.K2) THEN
        IF (AM.EQ.'VMZP') THEN
          TEMP1= (U(LL+K2) - U(J2)*U(K2)) * U(J1) +
     *           (U(JJ+K2) - U(J1)*U(K2)) * U(J2) +
     *           U(JJ+J2)*U(K2) - U(JJ+K2)*U(LL+K2)
        ELSE
          TEMP1= UIJKA(OFF2(J1)+INDUA(J1,J2,K2,N)) - U(JJ+K2)*U(LL+K2)
        ENDIF 
      ELSEIF (J1.EQ.J2) THEN
        IF (AM.EQ.'VMZP') THEN
          IF (K1.LT.K2) THEN
            TEMP1= (U(JJ+K2) - U(J1)*U(K2)) * U(K1) + (U(JJ+K1) - 
     *       U(J1)*U(K1)) * U(K2) + U(KK+K2)*U(J1) - U(JJ+K1)*U(JJ+K2)
          ELSE
            TEMP1= (U(JJ+K2) - U(J1)*U(K2)) * U(K1) + (U(JJ+K1) - 
     *       U(J1)*U(K1)) * U(K2) +
     *       U(N-K2+(K2-1)*(2*N-K2)/2 + K1)*U(J1) - U(JJ+K1)*U(JJ+K2)
          ENDIF
        ELSE
          IF (K1.LT.K2) THEN
            TEMP1= UIJKA(OFF2(J1)+INDUA(J1,K1,K2,N)) - U(JJ+K1)*U(JJ+K2)
          ELSE
            TEMP1= UIJKA(OFF2(J1)+INDUA(J1,K2,K1,N)) - U(JJ+K1)*U(JJ+K2)
          ENDIF
        ENDIF
      ELSEIF (AM.EQ.'VMZP') THEN
        IF (K1.LT.J2) THEN
          UK1J2= U(KK+J2)
          UK1K2= U(KK+K2)
        ELSE
          UK1J2= U(LL+K1)
          IF (K1.LT.K2) THEN
            UK1K2= U(KK+K2)
          ELSE
            UK1K2= U(N-K2+(K2-1)*(2*N-K2)/2 + K1)
          ENDIF
        ENDIF
        TEMP1= (U(JJ+J2) - U(J1)*U(J2)) * U(K1)*U(K2) +
     *           (U(JJ+K2) - U(J1)*U(K2)) * U(K1)*U(J2) +
     *           (UK1J2-U(K1)*U(J2))*U(J1)*U(K2) +
     *           (UK1K2-U(K1)*U(K2))*U(J1)*U(J2)
      ELSEIF (K1.LT.J2) THEN
        TEMP1= UIJKLA(OFF3(J1)+OFF2(K1)-OFF2(J1+1)+INDUA(K1,J2,K2,N)) -
     *   U(JJ+K1)*U(LL+K2)
      ELSEIF (K1.GT.K2) THEN
        TEMP1= UIJKLA(OFF3(J1)+OFF2(J2)-OFF2(J1+1)+INDUA(J2,K2,K1,N)) -
     *   U(JJ+K1)*U(LL+K2)
      ELSE
        TEMP1= UIJKLA(OFF3(J1)+OFF2(J2)-OFF2(J1+1)+INDUA(J2,K1,K2,N)) -
     *   U(JJ+K1)*U(LL+K2)
      ENDIF
      IJK=IJK+1
   46 W2(IJK)= TEMP1
      IF (K1.EQ.N) THEN
        J1=J1+1
        K1=J1+1
        JJ=JJ + N-J1
        KK=JJ + N-K1
      ELSE
        K1=K1+1
        KK=KK + N-K1 
      ENDIF
   48 CONTINUE

      IJK=IJK+1
      W2(IJK)= U(N1)*(1.-U(N1))

C  
C  This block of statements inverts the array whose upper triangular elements
C  are stored in the vector W2.  It ends with the DO 7205 structure.
C
      NDIAG= -N1
      DO 7204 K=1,N1 
      NDIAG= NDIAG + N1 + 2 - K
      TEMP1= W2(NDIAG)

      IF (TEMP1.EQ.0.0) THEN
        WRITE(*,'(2A,I5,A,I5)') ' PROGRAM TERMINATED: W2 inversion',
     *   ' failed on cluster/column',IC,'/',K
        WRITE(9,'(2A,I5,A,I5)') ' PROGRAM TERMINATED: W2 inversion',
     *   ' failed on cluster/column',IC,'/',K
        STOP
      ENDIF

      IJK=0
      DO 7202 I=1,N1
      IF (I.EQ.K) THEN
        L1 = IJK+2
        L2 = L1+N1-I-1
        IJK= L2
      ELSE     
        TEMP2= W2(IJSYM(N1,I,K))
        DO 7201 J=I,N1
        IJK=IJK+1
        IF (J.EQ.K) THEN
          W2(IJK)= W2(IJK)/TEMP1
        ELSE
          W2(IJK)= W2(IJK) - TEMP2*W2(IJSYM(N1,J,K))/TEMP1
        ENDIF
 7201   CONTINUE
      ENDIF          
 7202 CONTINUE 
      DO 7203 IJK=L1,L2
 7203 W2(IJK)= W2(IJK)/TEMP1
      W2(NDIAG)= -1.0/TEMP1
 7204 CONTINUE
      DO 7205 IJK=1,NDIAG
 7205 W2(IJK)= -W2(IJK)

      IF (METH(N).EQ.2) THEN
        DO 75 I=1,N1
        DO 74 J=1,N1
   74   U(J)= W2(IJSYM(N1,I,J))
   75   WRITE(11) (U(J),J=1,N1)
        ENDFILE 11
        DEALLOCATE (W2,STAT=IK)
        REWIND 11
      ENDIF

      IF (MRAM.GT.0) DEALLOCATE (UIJKA,UIJKLA,STAT=IK)
                
   76 ALLOCATE (W3A(NDIMB,N1),STAT=IK)

      IF (METH(N).EQ.1) THEN
        DO 85 I=1,NDIMB
        DO 85 J=1,N1
        W3A(I,J)= 0.
        DO 83 K=1,N
   83   W3A(I,J)= W3A(I,J) + W1A(I,K)*W2(IJSYM(N1,K,J))
        IF (IGEE.NE.1.OR.I.GT.NDIMB1) THEN
          DO 84 K=N+1,N1
   84     W3A(I,J)= W3A(I,J) + W1A(I,K)*W2(IJSYM(N1,K,J))
        ENDIF
   85   CONTINUE
        DEALLOCATE (U,W2,STAT=IK)
      ELSE
        DO 88 J=1,N1
        READ(11) (U(I),I=1,N1)
        DO 88 I=1,NDIMB
        W3A(I,J)= 0.
        DO 86 K=1,N
   86   W3A(I,J)= W3A(I,J) + W1A(I,K)*U(K)
        IF (IGEE.NE.1.OR.I.GT.NDIMB1) THEN
          DO 87 K=N+1,N1
   87     W3A(I,J)= W3A(I,J) + W1A(I,K)*U(K)
        ENDIF
   88   CONTINUE
        REWIND 11
        DEALLOCATE (U,STAT=IK)
      ENDIF

      DO 91 I=1,NDIMB
      V2(I)= 0.
      DO 90 J=1,N1
   90 V2(I)= V2(I) + W3A(I,J)*V1(J)
   91 S2(I)= S2(I) + V2(I)

      IF (IGEE.EQ.2) THEN
        DO 92 I=1,NDIMB
        DO 92 J=I,NDIMB
        DO 92 K=1,N1         
   92   S1(I,J)= S1(I,J) + W3A(I,K)*W1A(J,K)
      ELSE
        DO 96 I=1,NDIMB
        DO 95 J=I,NDIMB
        DO 93 K=1,N
   93   S1(I,J)= S1(I,J) + W3A(I,K)*W1A(J,K)
        IF (J.GT.NDIMB1) THEN
          DO 94 K=N+1,N1
   94     S1(I,J)= S1(I,J) + W3A(I,K)*W1A(J,K)
        ENDIF 
   95   CONTINUE
        DO 96 J=1,NDIMB
        DO 96 K=1,N1
   96   S1A(I,J)= S1A(I,J) + W3A(I,K)*W1A(J,K)
      ENDIF
   
      IF (CONVERGE) THEN
        DO 97 I=1,NDIMB
        DO 97 J=I,NDIMB
   97   S3(I,J)= S3(I,J) + V2(I)*V2(J)
      ENDIF

      DEALLOCATE(W1A,V1,W3A,STAT=IK)
  100 CONTINUE

C  End of cluster loop; final convergence loop section below 

      DO 104 I=2,NDIMB
      DO 104 J=1,I-1
  104 S1(I,J)= S1(J,I)
      
      IF (CONVERGE) THEN
        DO 106 I=2,NDIMB
        DO 106 J=1,I-1
  106   S3(I,J)= S3(J,I)
      ENDIF

      CALL SINV(S1,IT1,NDIMB,IK)
      IF (IK.NE.0) THEN
        WRITE(*,'(2A,I3)') ' PROGRAM TERMINATED: S1 inversion failed ',
     *   'on column',IK
        WRITE(9,'(2A,I3)') ' PROGRAM TERMINATED: S1 inversion failed ',
     *   'on column',IK
        STOP
      ENDIF
     
      TEMP1= 0.
      DO 112 I=1,NDIMB
      TEMP2= 0.
      DO 110 J=1,NDIMB
  110 TEMP2= TEMP2 + S1(I,J)*S2(J)
      TEMP1= TEMP1 + ABS(TEMP2)
  112 B(I)= B(I) + TEMP2

      IF (CONVERGE) THEN
        LASTPASS= CONVERGE
      ELSE
        CONVERGE= (TEMP1.LT.TOL)
        LASTPASS= (.NOT.CONVERGE.AND.ITER.GE.MAXITER)
      ENDIF
 
  114 CONTINUE

C  End of convergence loop; if converged, CALCVSE procedure calculations are
C  followed by procedure CALCCORR calculations.  Print results. (V2 used as
C  temporary storage vector of NDIMB elements)

      IF (CONVERGE) THEN
        ALLOCATE (RVAR(NDIMB,NDIMB),RZ(NDIMB),STAT=IK)

        IF (IGEE.EQ.2) THEN
          DO 120 I=1,NDIMB
          DO 118 J=1,NDIMB
          V2(J)= 0.
          DO 118 K=1,NDIMB
  118     V2(J)= V2(J) + S1(I,K)*S3(K,J)
          DO 120 J=1,NDIMB
          RVAR(I,J)= 0.
          DO 120 K=1,NDIMB
  120     RVAR(I,J)= RVAR(I,J) + V2(K)*S1(K,J)
        ELSE       
          CALL SINV(S1A,IT1,NDIMB,IK)
          IF (IK.NE.0) WRITE(9,'(A,I3/9X,A)')
     *     ' WARNING: S1A INVERSION FAILED ON COLUMN',IK,
     *     ' OUTPUT FOR ROBUST DATA LISTED BELOW IS INVALID'
          DO 124 I=1,NDIMB
          DO 122 J=1,NDIMB
          V2(J)= 0.
          DO 122 K=1,NDIMB
  122     V2(J)= V2(J) + S1A(I,K)*S3(K,J)
          DO 124 J=1,NDIMB
          RVAR(I,J)= 0.
          DO 124 K=1,NDIMB
  124     RVAR(I,J)= RVAR(I,J) + V2(K)*S1A(J,K)
        ENDIF

        DO 128 I=1,NDIMB  
        IF (RVAR(I,I).GT.0.) THEN
          RZ(I)= B(I)/SQRT(RVAR(I,I))
        ELSE
          RZ(I)= 0.
        ENDIF
        DO 128 J=1,I-1
        IF (RVAR(I,I)*RVAR(J,J).LE.0.) THEN
          RVAR(I,J)= 999.0
        ELSE
          RVAR(I,J)= RVAR(I,J)/SQRT(RVAR(I,I)*RVAR(J,J))
        ENDIF
        IF (S1(I,I)*S1(J,J).LE.0.) THEN
          S1(I,J)= 999.0
        ELSE
          S1(I,J)= S1(I,J)/SQRT(S1(I,I)*S1(J,J))
        ENDIF
  128   CONTINUE

        WRITE(9,'(/A,I4,A//2A)') ' Converged after',ITER-1,
     *   ' iterations',' Final beta estimates:             Beta    ',
     *   'Robust SE     Naive SE     Robust Z'
        WRITE(9,'(5X,A,5X,4F13.5)') (LBL(I),B(I),SQRT(RVAR(I,I)),
     *   SQRT(S1(I,I)),RZ(I),I=1,NDIMB) 
        WRITE(9,'(/2A)') ' Robust estimates matrix:  ',
     *   'lower = correlations \ upper = variances'
        DO 132 I=1,NDIMB
  132   WRITE(9,'(A,I2,A,6F11.5,(/7X,6F11.5))') ' (',I,')  ',
     *   (RVAR(I,J),J=1,NDIMB)
        WRITE(9,'(/2A)') ' Naive estimates matrix:  ',
     *   'lower = correlations \ upper = variances' 
        DO 134 I=1,NDIMB
  134   WRITE(9,'(A,I2,A,6F11.5,(/7X,6F11.5))') ' (',I,')  ',
     *   (S1(I,J),J=1,NDIMB)

C  No convergence
   
      ELSE
        WRITE(9,'(/A,I4,A/A,6F8.5,2(/9F8.5))') ' No convergence after',
     *   ITER,' iterations',' Last estimate of betas: ',B
      ENDIF
C
      STOP
      END
