PDP-10 Archive: decus/20-0137/bmd/bmdx69.for from decuslib20-05

Trailing-Edge - PDP-10 Archives - decuslib20-05 - decus/20-0137/bmd/bmdx69.for
There is 1 other file named bmdx69.for in the archive. Click here to see a list.
C     MULTIVARIATE ANALYSIS OF VARIANCE AND COVARIANCE - OCTOBER 7,1966
C THIS IS A CONVERTED VERSION OF BMDX69 FOR THE 360/75 ORIGINALLY WRITTE
C IN FORTRAN TWO.  SUBROUTINES LWH, PMEANS AND SS ARE IDENTICAL TO THOSE
C IN BMD08V, AND ANOVA IS IDENTICAL EXCEPT FOR ABOUT A DOZEN STATEMENTS.
      DIMENSION IND(100)
      COMMON/ANARG/FF(180),NI,NCM,KM,MT,NV,ILL
C
      DIMENSION FIN(11),IN(11),AA(11),P(100),PL(11),A(99
     1),B(50),C(100),D(100),E(100),G(100),P1(100),P2(100),MN(100),IL(100
     2),S(5000),ID(100),J1(10)
      DOUBLE PRECISION PROB,FINI,PF,PC,Q003HL,Q005HL,Q006HL
       DOUBLE PRECISION P1,P2,Q007HL,Q009HL,Q010HL,AP,PA,Q008HL
      INTEGER PCOV
      DOUBLE PRECISION SPB
      INTEGER AA,A,B,C,D,E,U,UNION,H,Q012CT,Q013CT,Q014CT
      LOGICAL INCL
      DIMENSION EE(10)
	CALL USAGEB('BMDX69')
      L=10
      I=1
      DO 500 J=1,10
      AA(L)=I
      I=2*I
 500  L=L-1
      Q012CT=I
      Q013CT=2*I
      Q014CT=4*I
      DATA PROB,FINI/6HPROBLM,6HFINISH/
      DATA Q003HL/6HINDEX /
C
      DATA Q005HL/6HDESIGN/
      DATA Q006HL/6H,     /
      DATA Q007HL/6H      /
      DATA Q008HL/6H(     /
      DATA Q009HL/6H)     /
      DATA Q010HL/6HMEAN  /
      DATA QMINUS/4H-   /
      DATA SPB/6HSUBPRO/
  134  FORMAT ( '1BMDX69 - MULTIVARIATE ANALYSIS OF VARIANCE AND COVARIA
     *NCE - REVISED ',
     2'APRIL 14, 1969'
     * / ' HEALTH SCIENCES COMPUTING FACILITY, UCLA' /
     *  '0PROBLEM CODE 'A6)
      IIIII=1
      KLV=-1
      MTP=5
 101  READ (5,100)PF,PC,NV,NI,KM,NF,MT,PCOV
      IF(PF.EQ.PROB) GO TO 9143
      GO TO 20
  100 FORMAT(2A6,6I2)
 7777 IF(PF.EQ.PROB)GO TO 111
  110 IF(PF.EQ.FINI)GO TO 113
 20   IF(KLV)2111,2111,7766
2111  KLV=1
      GO TO (2101,2102,2103,2104,2105,2106,2107),IIIII
 7766 READ (5,100)PF
  118 FORMAT(13A6,A2)
      GO TO 7777
  113 IF(MTP.EQ.5)GO TO 117
 116  REWIND MTP
 117  STOP
 2000 WRITE (6,2001)
2001  FORMAT(55H0 THE NUMBER OF ANALYSIS OF VARIANCE INDICES IS TOO BIG)
      KLV=1
      GO TO 101
111   KLV=-1
 9143 IF(MT)102,102,103
 102  MT=5
 103  IF(MT-MTP) 9247,9248,9247
 9247 REWIND MT
 9248 IF(MT-MTP)104,107,104
  104 IF(MTP.EQ.5)GO TO 107
 106  REWIND MTP
 107  MTP=MT
 5013 IF(MT.EQ.5)GO TO 5114
 5116 REWIND MT
 5114 WRITE (6,134)PC
      READ (5,191)AP, (J1(I),I=1,10)
      IIIII=2
 191  FORMAT(A6,10I3)
      IF(AP.NE.Q003HL)GO TO 20
 193  N=NI
      IF(NI*(11-NI))2000,2000,2011
 989   FORMAT(36X,10A3)
 2011 DO 135 I=1,NI
       IN(N)=J1(I)
 135  N=N-1
      IF(NF*(11-NF))7000,7000,200
3000  WRITE(6,4000)
4000  FORMAT(82H0THE NUMBER OF COMPONENTS IN THE ANALYSIS OF VARIANCE PA
     1RT OF THE MODEL IS TOO BIG)
      KLV=1
      GO TO 101
5000  WRITE(6,6000)
6000  FORMAT(26H0VIOLATING RESTRICTION (B))
      KLV=1
      GO TO 101
7000  WRITE(6,9000)
9000  FORMAT(38H0THE NUMBER OF FORMAT CARDS IS TOO BIG)
      KLV=1
      GO TO 101
2101  WRITE(6,2121)
2121  FORMAT(32H0PROBLEM CARD IS OUT OF SEQUENCE)
      GO TO 7766
2102  WRITE(6,2122)
2122  FORMAT(30H0INDEX CARD IS OUT OF SEQUENCE)
      GO TO 7766
2103  WRITE(6,2123)
2123  FORMAT(31H0DESIGN CARD IS OUT OF SEQUENCE)
      GO TO 7766
2104  WRITE(6,2124)
2124  FORMAT(74H0SUBPROBLEM CARD OR THE NEXT PROBLEM CARD OR FINISH CARD
     1IS OUT OF SEQUENCE)
      GO TO 7766
2105  WRITE(6,2125)
2125  FORMAT(33H0ILLEGAL CHARACTER ON DESIGN CARD)
      GO TO 7766
2106  WRITE(6,2126)
2126  FORMAT(35H0WRONG SPECIFICATION ON DESIGN CARD)
      GO TO 7766
2107  WRITE(6,2127)
2127  FORMAT(43H0DISAGREEMENT OF DESIGN CARD AND INDEX CARD)
      GO TO 7766
200   NF=MAX0(1,NF)*20
      READ (5,179)PA,(P(I),I=1,66)
      IIIII=3
 179  FORMAT(A6,66A1)
      IF(PA.NE.Q005HL)GO TO 20
 152  L=NI
      M=0
      IIIII=5
      DO 31 I=1,10
      FIN(I)=IN(I)
      K=LWH(P(I))
      GO TO (31,33,20,20,20,20,20,20),K
 33   PL(L)=P(I)
      M=M+1
      EE(M)=P(I)
      L=L-1
 31   CONTINUE
      IIIII=7
      IF(L)20,32,20
 32   WRITE (6,722) (EE(I),I=1,NI)
      WRITE (6,723)(J1(I),I=1,NI)
 722  FORMAT(17H0INDEX           10(4X,A1))
 723  FORMAT(17H NUMBER OF LEVELS10I5)
      WRITE (6,93)(P(I),I=1,66)
 93   FORMAT(12H0DESIGN CARD 5X,72A1)
      P(10)=','
      CO=1.0
      LB=-1
      RB=-1.0
      N=0
      MO=10
      MI=66
 178  DO 4 I=MO,MI
      IIIII=6
      K=LWH(P(I))
      GO TO (4,36,37,8,9,10,11,20),K
 37   IF(CO)20,20,12
 12   DO=1.0
      CO=-1.0
      PE=-1.0
      LL=1
      NL=0
      N=N+1
      A(N)=0
      B(N)=0
      GO TO 4
 11   IF(DO)20,20,13
 13   A(N)=Q012CT
      DO=-1.0
      GO TO 4
 9    IF(RB)20,20,14
 14   NL=1
      CO=1.0
      PE=1.0
      RB=-1.0
      GO TO 4
 8    IF(LB)20,20,15
 15   NL=-1
      RB=1.0
      LB=-1
      CO=-1.0
      PE=-1.0
      GO TO 4
 36   IF(LL)20,20,16
 16   DO 17 M=1,NI
      IF(P(I).EQ.PL(M))  GO TO 19
 17   CONTINUE
      GO TO 20
 19   IF(NL) 21,22,23
 22   A(N)=A(N)+Q013CT+AA(M)
      NL=1
      DO=-1.0
      LB=1
      CO=1.0
      PE=1.0
      GO TO 4
 23   IF(B(N)/Q013CT.EQ.0) B(N)=B(N)+Q013CT
 21   B(N)=B(N)+AA(M)
 4    CONTINUE
      READ (5,130)(P(I),I=1,72)
      MO=1
      MI=72
      WRITE (6,471)(P(I),I=1,72)
 471  FORMAT(18X,72A1)
      GO TO 178
 10   IF(PE)20,20,24
 24   C(1)=Q013CT
      D(1)=0
      MN(1)=1
      J=1
      DO 51 K=1,N
      J0=J
      DO 51 I=1,J0
      IF(.NOT.INCL(B(K),Q014CT-1-C(I)) .OR.
     X .NOT.INCL(A(K),Q014CT-1-D(I))) GO TO 51
 52   J=J+1
      IF(J-100)2221,3000,3000
 2221 MN(J)=0
      D(J)=UNION(B(K),D(I))
      C(J)=UNION(A(K),C(I))
      E(J)=MOD(C(J)+D(J),Q012CT)
 51   CONTINUE
      NCM=J
      DO 25 I=1,N
      IF(MOD(A(I)/Q012CT,2).EQ.0) GO TO 25
 26   H=MOD(A(I)+B(I),Q012CT)
      DO 27 J=1,NCM
       IF(E(J).EQ.H) GO TO 28
 27   CONTINUE
      GO TO 20
 28   MN(J)=1
 25   CONTINUE
      E(1)=0
      DO 53 I=1,NCM
      C(I)=MOD(C(I),Q012CT)
      D(I)=MOD(D(I),Q012CT)
 53   G(I)=1024*(10*NBITS(E(I))+NBITS(C(I)))+C(I)
      DO 86 I=1,NCM
      X=1.E20
      DO 89 K=I,NCM
      IF(X-G(K))89,89,88
 88   J=K
      X=G(K)
 89   CONTINUE
      G(J)=G(I)
      U=C(J)
      C(J)=C(I)
      C(I)=U
      U=D(J)
      D(J)=D(I)
      D(I)=U
      NUU=MN(I)
      MN(I)=MN(J)
      MN(J)=NUU
 86   E(I)=C(I)+D(I)
      DO 122 I=2,NCM
      P2(I)=Q007HL
      P1(I)=Q007HL
      N=0
      L=NI
      DO 123 J=1,NI
      IF(MOD(C(I)/AA(L),2).EQ.0) GO TO 123
 124  N=N+1
      IF(N.LE.6) CALL PUTCHR(P1(I),N,PL(L))
      IF(N.GT.6) CALL PUTCHR(P2(I),N-6,PL(L))
 123  L=L-1
C
      IF(D(I)) 125,122,125
 125  N=N+1
      IF(N.LE.6) CALL PUTCHR(P1(I),N,Q008HL)
      IF(N.GT.6) CALL PUTCHR(P2(I),N-6,Q008HL)
      L=NI
      DO 126 J=1,NI
      IF(MOD(D(I)/AA(L),2) .EQ.0) GO TO 126
 127  N=N+1
      IF(N.LE.6) CALL PUTCHR(P1(I),N,PL(L))
      IF(N.GT.6) CALL PUTCHR(P2(I),N-6,PL(L))
 126  L=L-1
      N=N+1
      IF(N.LE.6) CALL PUTCHR(P1(I),N,Q009HL)
      IF(N.GT.6) CALL PUTCHR(P2(I),N-6,Q009HL)
 130  FORMAT(72A1)
 122  CONTINUE
      P1(1)=(+Q010HL)
      P2(1)=(+Q007HL)
 121  FORMAT(12A6)
       DO 602 I=2,NCM
       IF(MN(I))602,602,603
C
C
  603 DO  601 J=1,I
 601  IF(INCL(E(J),E(I))) MN(J)=1
 602   CONTINUE
C
  921 FORMAT(20A4)
       READ(5,921)(FF (I),I=1,NF)
      WRITE (6,5432)(FF(I),I=1,NF)
 5432 FORMAT(21H0VARIABLE FORMAT          20A4/(21X,20A4))
      M1P1=1
      DO 90 I=M1P1,5000
 90   S(I)=0.0
      ILL=5000-M1P1
      CALL ANOVA(AA,MN,E,IN,IL,S(M1P1),ID)
      IF(ILL)5000,5000,2888
 2888 NVV= (NV*(NV+1))/2
      NN=NVV*NCM
      DO 131 J=1,NV
      DO 131 I=1,NCM
      IF(MN(I))132,131,132
 132   L1=IL(I)+J+M1P1-1
       M=L1+MN(I)
 131  CONTINUE
        K=NVV+M1P1-1
      R=QMINUS
       IF(PCOV)8000,8001,8000
 8000 DO 8002 I=2,NCM
 8003 WRITE(6,8004)P1(I),P2(I)
 8004 FORMAT( ////' CROSS PRODUCT MATRIX FOR COMPONENT',2X,2A6)
      J2=0
 8045 JI=J2+1
      J2=MIN0(J2+10,NV)
      WRITE(6,8046)(J,J=JI,J2)
 8046 FORMAT(//' VARIABLE',10(I3, 9X))
      M1=K+JI*(JI+1)/2
      DO 8005 J=JI,NV
      M2=M1+MIN0(J-JI,9)
      WRITE(6,8006)J,(S(M),M=M1,M2)
 8005 M1=M1+J
      IF(J2.LT.NV)GO TO 8045
 8002 K=K+NVV
 8006 FORMAT(I3 ,10F12.5)
 8001 MM=IL(1)+M1P1
      DO 401 I=1,NV
      L=IL(1)+I+M1P1-1
      WRITE(6,400)I,S(L),I
  400 FORMAT('-MEAN FOR VARIABLE',I3,F14.5/'0CELL MEANS FOR VARIABLE',
     * I3)
      DO 401 J=2,NCM
      IF(MN(J))402,401,402
 402  L1=IL(J)+I+M1P1-1
      CALL PMEANS(E(J),AA,PL,IN,S(L1),NI,NV)
 401  CONTINUE
      KLV=0
 1    READ (5,1028) PF,PC,(IND(I),I=7,66)
 1028 FORMAT(2A6,60I1)
      IF(PF.EQ.SPB) GO TO 1029
      NV=10*IND(7)+IND(8)
      NI=10*IND(9)+IND(10)
      IIIII=4
      KM=10*IND(11)+IND(12)
      NF=10*IND(13)+IND(14)
      MT=10*IND(15)+IND(16)
      PCOV=10*IND(17)+IND(18)
      GO TO 7777
 1029 CALL CONV(PC,IND,6)
      CALL COVAR(S(M1P1),ID,MN,IL,E,AA,PL,IN,IND,P1,P2)
      GO TO 1
       END
      SUBROUTINE MULTIV(S,JND,JC,ID,T,MV,NC,P1,P2)
      COMMON/ANARG/FF(180),NI,NCM,KM,MT,NV,ILL
      DIMENSION S(2),JND(2),JC(2),ID(2),P1(2),P2(2),T(2),DT(100),IND(100
     1),NP(100),KKK(100),X(100),KND(100)
      DOUBLE PRECISIONP1,P2
      MV1=MV+1
      M6=NC+1
      IF(NC)15,15,16
   16 WRITE(6,17)(JC(I),I=M6,MV)
      WRITE(6,18)(JC(I),I=1,NC)
      GO TO 19
   15 WRITE(6,25)(JC(I),I=M6,MV)
 17   FORMAT(36H1MULTIVARIATE ANALYSIS OF COVARIANCE/20H0DEPENDENT VARIA
     1BLES   20I5/(20X,20I5))
 18   FORMAT(20H0COVARIATES           20I5/(20X,20I5))
 25   FORMAT(34H1MULTIVARIATE ANALYSIS OF VARIANCE/20H0DEPENDENT VARIABL
     1ES   20I5/(20X,20I5))
 19   NVV=(NV*(NV+1))/2
      WRITE(6,50)
   50 FORMAT(1H-,'SOURCE     LOG(GENERALIZED   U-STATISTIC     DEGREES
     1OF   APPROXIMATE   DEGREES OF'/ 18X,
     2    'VARIANCE)',22X,'FREEDOM    F- STATISTIC     FREEDOM'//)
C
      MM=NVV
      LL=(NCM-1)*NVV
      DO 1 K=2,NCM
      N=0
      DO 2 I=1,MV
      DO 2 J=1,MV
      I1=MAX0(JC(I),JC(J))
      I1=(I1*(I1-1))/2+JC(I)+JC(J)-I1
      M1=MM+I1
      L1=LL+I1
      N=N+1
 2    T(N)=S(M1)+S(L1)
      MM=MM+NVV
      NP(K)=0
      CALL SOLVIT(T,MV,JND,1.E-6,NP(K),D)
      L=1
      DO 3 I=1,MV
      X(I)=T(L)
      L=L+MV1
      IF(K-NCM)3,54,3
 54   KND(I)=-I
      CALL SOLVIT(T,MV,IND,1.E-6,KND(I),D)
 3    IND(I)=1-JND(I)
      KKK(K)=0
 1    CALL SOLVIT(T,MV,IND,1.E-6,KKK(K),DT(K))
      IP=MV-NC
      DE=DT(NCM)+IP*ALOG(.5)
      P=IP
      PP=P*P
      N=ID(NCM)-NP(NCM)
      TN=N
      NCM1=NCM-1
      DO 4 I=2,NCM1
      IQ=ID(I)-NP(I)+NP(NCM)
      IF(KKK(I)-IP)30,31,30
   30 WRITE(6,32)P1(I),P2(I)
 32   FORMAT(A7,A6,6X,9HUNDEFINED)
      GO TO 4
 31   Q=IQ
      QQ=Q*Q
      IF(PP+QQ-5.)51,60,51
 60   Z=1.
      GO TO 52
   51 Z=SQRT((PP*QQ-4.)/(PP+QQ-5.))
   52 U=EXP(DE-DT(I))
      Y=U**(1./Z)
      D2=(TN-(P-Q+1.)*.5)*Z-P*Q*.5+1.
      I1=P*Q
      F=(1.-Y)/Y*D2/(P*Q)
      WRITE(6,5)P1(I),P2(I),DT(I),U,IP,IQ,N,F,I1,D2
 4    CONTINUE
    5 FORMAT(1X,2A6,F10.5,6X,F10.6,4X,3I5,F12.4,I5,F8.2)
      L=1
      IF(NC)40,40,41
 41   DO 42 I=1,NC
      IF(KND(I))43,43,44
 43   I1=0
      IQ=0
      U=1.
      GO TO 45
 44   I1=IP
      IQ=1
      U=X(I)/T(L)
 45   D2=N-IP+1
      F=(1.-U)/U*D2/P
      GH=DE/U
      L=L+MV1
   42 WRITE(6,37)JC(I),GH,U,IP,IQ,N,F,I1,D2
 37   FORMAT(' COVARIATE',I3,F10.5,6X,F10.6,4X,3I5,F12.4,I5,F8.2)
   40 WRITE(6,5)P1(NCM),P2(NCM),DE
      RETURN
      END
      SUBROUTINE COVAR(S,ID,MN,IL,Z,AA,PL,IN,IND,P1,P2)
      COMMON/ANARG/FF(180),NI,NCM,KM,MT,NV,ILL
       DIMENSION JC(100),NP(100),S(2),ID(2),P1(2),P2(2),IND(100),MN(2),I
     1L(2),Z(2),JND(100)
      DIMENSION AA(2),PL(2),IN(2)
      DOUBLE PRECISION UQ,P1,P2,C,D,Q01,Q02,Q1,Q2
      DATA Q01,Q02/6HCOVARI,6HATES  /
      DOUBLE PRECISION Q03,Q04
      DATA Q03,Q04/6HFULL M,6HODEL  /
       Q1=P1(NCM)
       Q2=P2(NCM)
       NVV=(NV*(NV+1))/2
C
      IF(NV.GT.66) READ (5,1) (IND(I),I=67,NV)
 1    FORMAT(6X,66I1)
      L=0
      DO 250 I=1,NV
      IF(IND(I)-2)250,251,250
 251  L=L+1
      JND(L)=0
      JC(L)=I
 250  CONTINUE
      NC=L
      DO 252 I=1,NV
      IF(IND(I)-1)252,253,252
 253  L=L+1
      JND(L)=1
      JC(L)=I
 252  CONTINUE
      MV=L
      LLL=KM*NV+1
      JD=JC(L)
      IF(MV-NC-1)255,256,255
 255  CALL MULTIV(S,JND,JC,ID,S(LLL+1),MV,NC,P1,P2)
      RETURN
 256  MM=0
      IF(ILL-MV*(MV+NCM)-3) 2001,2001,2000
 2001 ILL=0
      WRITE(6,1111)UQ,JD,NC,(JC(I),I=1,NC)
 1111 FORMAT(1H0,A6,22I3/(7X,22I3))
      RETURN
 2000 P1(NCM)=Q03
      P2(NCM)=Q04
      LL=(NCM-1)*NVV
       NN=MV*MV+LLL
       L0=NN
      DO 2 NO=1,NCM
       K=LLL
       DO 3 I=1,MV
       DO 3 J=1,MV
       I1= MAX0 (JC(I),JC(J))
       N=(I1*(I1-1))/2+JC(I)+JC(J)-I1
       M1=MM+N
       L1=LL+N
       K=K+1
 3     S(K)=S(L1)+S(M1)
       R=S(K)/2.
       MM=MM+NVV
       NP(NO)=0
      CALL SOLVIT(S(LLL+1),MV,JND,1.E-6,NP(NO),D)
       K=MV+LLL
       N=NN+1
       NN=NN+MV
       DO 2 I=N,NN
       S(I)=S(K)
 2     K=K+MV
      WRITE(6,241)JD
  241 FORMAT(1H1,'ANALYSIS OF COVARIANCE FOR DEPENDENT VARIABLE',I3//)
      IF(NC.EQ.0)GO TO 12
      WRITE(6,9008)
 9008 FORMAT(11X,'REGRESSION COEFICIENTS UNDER EACH HYPOTHESIS')
C
      I1=0
      II=NCM
 7     I0=I1+1
       JJ= MIN0 (II,10)
       I1=I1+JJ
       II=II-JJ
      WRITE(6,8)(P1(I),P2(I),I=I0,I1)
 8     FORMAT(1H010X,20A6)
      WRITE(6,9)
 9     FORMAT(10H COVARIATE)
      K0=L0+(I0-1)*MV+1
      K1=L0+(I1-1)*MV+1
       DO 10 I=1,NC
      WRITE(6,11)JC(I),(S(K),K=K0,K1,MV)
 11    FORMAT(I5,10F12.5)
       K0=K0+1
 10    K1=K1+1
       IF(II)12,12,7
 12    L=L0+MV
       E=S(NN)/2.
       DO 5 I=L,NN,MV
 5     S(I)=S(I)-E
      WRITE(6,14)
 14    FORMAT(11H-    SOURCE13X,48HSUM OF SQUARES  DEGREES OF  MEAN SQUA
     1RE        F/41X,7HFREEDOM//)
       IDE=ID(NCM)-NP(NCM)
       EM=E/FLOAT (IDE)
       NM=NCM-1
      DO 15 I=1,NM
       IDI=ID(I)-NP(I)+NP(NCM)
       T=S(L)/FLOAT (IDI)
       F=T/EM
      WRITE(6,16)P1(I),P2(I),S(L),IDI,T,F
 16    FORMAT(6X,2A6,F18.4,I10,F16.4,F12.4)
 15    L=L+MV
      IF(NC.EQ.0)GO TO 9016
       R=R-E
      C=Q01
      D=Q02
       T=R/FLOAT (NP(NCM))
       F=T/EM
      WRITE(6,16)C,D,R,NP(NCM),T,F
       J=LLL+1
       K=LLL
       DO 17 I=1,NC
       K=K+MV
      R=-S(K)*S(K)/S(J)/2.
       J=J+MV+1
       IDF=-I
      CALL SOLVIT(S(LLL+1),MV,JND,1.E-6,IDF,D)
       IF(IDF)994,994,995
 994   R=0.
 995   F=R/EM
   17 WRITE(6,18)JC(I),R,IDF,R,F
 18    FORMAT(6X,9HCOVARIATEI3,F18.4,I10,F16.4,F12.4)
 9016 WRITE(6,16)Q1,Q2,E,IDE,EM
      IF(NC.EQ.0)RETURN
       P1(NCM)=Q1
       P2(NCM)=Q2
       K0=(NCM-1)*MV+L0
       L=IL(1)
       KKK=K0
       A=0.
       DO 50 I=1,NC
       KKK=KKK+1
       KK=L+JC(I)
 50    A=A+S(KK)*S(KKK)
      LZZ=0
       DO 139 I=2,NCM
       IF(MN(I))139,139,133
 133   K=IL(I)
      IF(LZZ)400,400,401
 400  LZZ=1
      WRITE(6,402)
 402  FORMAT(20H-ADJUSTED CELL MEANS)
 401  N=K0+MV
       KK=IL(I)+MN(I)-NV
       DO 135 L=K,KK,NV
       M=L+JD
       C=0.
       KKK=K0
       DO 136 II=1,NC
       J=L+JC(II)
       KKK=KKK+1
 136   C=C+S(J)*S(KKK)
      N=N+1
      IF(ILL-N)2007,135,135
 2007 WRITE(6,2008)P1(I),P2(I)
 2008 FORMAT(80H0THIS PROBLEM IS TOO LARGE TO ALLOW COMPUTATION OF ADJUS
     1TED MEANS FOR COMPONENT  2A6)
      GO TO 139
 135  S(N)=S(M)+A-C
      L1=K0+MV+1
        CALL PMEANS (Z(I),AA,PL,IN,S(L1),NI,1)
 139  CONTINUE
      RETURN
       END
C             SUBROUTINE ANOVA FOR BMD08V            MARCH  1, 1966
      SUBROUTINE ANOVA(AJ,MN,CW,IN,IL,S,ID)
      DIMENSIONAJ(2),MN(2),CW(2),IN(2),IL(2),S(2),FF(180),ID(2),ST(100),
     1IC(11,100),II(11),IJ(11,100),X(255),FP(100),SF(100),SG(100)
      COMMON/ANARG/FF,NI,NCW,MK,MT,NV,ILL
      LOGICAL INCL
      INTEGER CW,AJ,QCT
      QCT=2**10
       IF(MT)500,500,501
 500   DO 403 I=2,NCW
      AJ(NI1)=QCT
      CW(I)=CW(I)+QCT
       J1=I-1
 410   DO 404 J=1,J1
      IF(.NOT.INCL(CW(J),CW(I))) GO TO 404
 405  M=IL(I)-MT
      L=IL(J)-MT
      S(M)=S(M)-S(L)
 404   CONTINUE
       II(NI1)=-2
       DO 408 K=1,NI1
      IF(MOD(CW(I)/AJ(K),2).EQ.0) GO TO 408
 402   II(K)=II(K)+1
       IF(II(K))401,400,400
 400   II(K)=-IN(K)
 408   CONTINUE
 401   DO 407 J=1,I
      IF(.NOT.INCL(CW(J),CW(I))) GO TO 407
 418   IL(J)=IL(J)+IJ(K,J)
 407   CONTINUE
       IF(NI1-K)410,403,410
 403   CONTINUE
       RETURN
 501   KM=MK
      DO 228 I=1,100
 228  ST(I)=0.0
      NVV=(NV*(NV+1))/2
      NN=1
      DO 51 I=1,NI
 51   NN=NN*IN(I)
      NI1=NI+1
      DO 10 J=1,NCW
      ID(J)=1
      IL(J)=1
      IO=1
      DO 12 I=1,NI
      IF(MOD(CW(J)/AJ(I),2).EQ.0) GO TO 13
 52   IC(I,J)=1
      ID(J)=ID(J)*IN(I)
      GO TO 12
 13   IO=I
      IC(I,J)=0
 12   CONTINUE
      IF(MN(J))11,11,107
 11   DO 14 K=IO,NI
 14   IC(K,J)=-IC(K,J)
 107  FP(J)=NN/ID(J)
 10   IC(NI1,J)=-1
      IN(NI1)=2
      DO 7 I=1,NI1
      DO 15 J=1,NCW
      IF(IC(I,J)-1)16,17,16
 17   IJ(I,J)=NV
      IL(J)=IL(J)*IN(I)
      GO TO 15
 16   IJ(I,J)=(1-IL(J))*NV
 15   CONTINUE
 7    II(I)=-IN(I)
      N=(NVV*NCW)/NV+1
       DO 118 J=1,NCW
       IF(MN(J))118,118,110
 110   N=N+IL(J)
       IL(J)=(N-IL(J))*NV
 118   CONTINUE
       MK=N+1
       DO 119 J=1,NCW
       IF(MN(J))111,111,119
 111   N=N+IL(J)
       IL(J)=(N-IL(J))*NV
 119   CONTINUE
      KU=KM
      FQ=0.0
      IF(N*NV-ILL)19,2000,2000
 2000  ILL=0
       RETURN
 19   DO 23 M=1,NV
      KU=KU+1
      IF(KU-KM)22,22,21
 21   READ (MT,FF)(X(K),K=1,KM)
      KU=1
 22   XK=X(KU)
      ST(M)=ST(M)+XK
      DO 23 J=1,NCW
      N=IL(J)+M
 23   S(N)=S(N)+XK
      FQ=FQ+1.0
      DO 20 I=1,NI1
      II(I)=II(I)+1
      IF(II(I))24,20,20
 20   II(I)=-IN(I)
 24   DO 25 J=2,NCW
      IL(J)=IL(J)+IJ(I,J)
      IF(IC(I,J))26,25,25
 26   N=IL(J)
      DO 222 K=1,NV
 222  SG(K)=0.0
      LL=-IJ(I,J)/NV+1
      K0=(J-1)*NVV
      GR=FP(J)/FQ
      DO 27 L=1,LL
      M=K0
      DO 27 K=1,NV
      N=N+1
      SF(K)=S(N)-GR*ST(K)
      SG(K)=SG(K)+SF(K)
       DO 277 KK=1,K
      M=M+1
  277 S(M)=S(M)+SF(K)*SF(KK)
      IF(I-NI1)60,27,60
 60   S(N)=0.0
 27   CONTINUE
      GR=FQ/FP(J)-FLOAT(LL)
      IF(GR)25,25,372
 372  M=K0
      DO 327 K=1,NV
        DO 327 KK=1,K
      M=M+1
 327  S(M)=S(M)+SG(K)*SG(KK)/GR
 25   CONTINUE
      IF(I-NI1)19,30,19
 30   K1=0
      K=0
      L0=IL(1)+1
      L1=IL(1)+NV
      DO 373 L=L0,L1
        DO 373 LL=L0,L
      K=K+1
 373  S(K)=S(L)*S(LL)
      DO 31 J=1,NCW
      NM=NV-IJ(NI1,J)
      MN(J)=MN(J)*NM
      F=(NN*NV)/NM
      K2=IL(J)+1
      K3=IL(J)+NM
      DO 70 K=K2,K3
 70   S(K)=S(K)/F
      F=NN/ID(J)
      K0=K1+1
      K1=K1+NVV
      DO 31 K=K0,K1
 31   S(K)=S(K)/F
      L1=NVV
      DO 32 J=2,NCW
      J1=J-1
      L0=L1+1
      L1=L1+NVV
      DO 32 K=1,J1
      IF(.NOT.INCL(CW(K),CW(J))) GO TO 32
C   FUNCTION LWH FOR BMDX69  MAY 14, 1968
 33   IF(K-1)62,62,61
 61   M=(K-1)*NVV
      DO 34 L=L0,L1
      M=M+1
 34   S(L)=S(L)-S(M)
 62   ID(J)=ID(J)-ID(K)
 32   CONTINUE
      RETURN
      END
      FUNCTION LWH(P)
      DIMENSION A(17)
      DATA A/' ',' ',',','(',')','.','$','=','+','-','*','/',1H',
     *'=','+','(',')'/
      DO 1 I=1,17
      IF(P.EQ.A(I)) GO TO 2
 1    CONTINUE
      LWH=2
      RETURN
 2    IF(I.GE.16) I=I-12
      LWH=MIN0(8,I)
      RETURN
      END
C             SUBROUTINE PMEANS FOR BMD08V           MARCH  1, 1966
      SUBROUTINE PMEANS(E,A,PL,IN,S,NI,NV)
      DIMENSION A(2),PL(2),S(2),LL(11),P(11),JO(11),IN(2)
      INTEGER E,A
 40   M=NI+1
      DO 1 I=1,NI
      IF(MOD(E/A(I),2).EQ.0) GO TO 1
 2    M=M-1
      P(M)=PL(I)
      JO(M)=IN(I)
 1    CONTINUE
      NJ=NI-2
      N=NI-M+1
      JO1=JO(NI)
      JO2=JO(NI-1)
      L1=1-NV
      IF(N-2)3,88,5
 88   WRITE (6,127)P(NI),(I,I=1,JO1)
      GO TO 89
 5    DO 6 I=M,NI
 6    LL(I)=1
 11   WRITE (6,7)(P(K),LL(K),K=M,NJ)
 7    FORMAT(1H08(4X,A1,2H =I3))
      GO TO 8
 9    I=NJ
      DO 10 J=3,N
      LL(I)=LL(I)+1
      IF(LL(I)-JO(I))11,11,12
 12   LL(I)=1
 10   I=I-1
      RETURN
 8    WRITE (6,27)P(NI),(I,I=1,JO1)
 27   FORMAT(   5X,A1,2H =I7,9I12, /(3X,10I12))
 89   I=1
      L0=L1+NV
      L1=L1+JO1*NV
      WRITE (6,24)P(NI-1),I,(S(L),L=L0,L1,NV)
 24   FORMAT(1X,A1,2H =I3,10F12.5, /(7X,10F12.5))
      DO 29 I=2,JO2
      L0=L1+NV
      L1=L1+JO1*NV
 29   WRITE (6,25)I,(S(L),L=L0,L1,NV)
 25   FORMAT(4X,I3,10F12.5, /(7X,10F12.5))
      IF(N-2)37,37,9
 3    L0=L1+NV
      L1=L1+NV*JO1
      WRITE (6,127)P(NI),(I,I=1,JO1)
 127  FORMAT(1H04X,A1,2H =I7,9I12, /(3X,10I12))
      WRITE (6,35)(S(L),L=L0,L1,NV)
 35   FORMAT(7X,10F12.5)
 37   RETURN
      END
       SUBROUTINE SOLVIT(A,N,IND,T,IDF,DT)
       DIMENSION A(2),U(100),V(100),IN(100),IND(2)
      N1=N
      IF(IDF)20,30,30
 20   I=-IDF
       K0=(N+1)*I-N
      IF(IN(I))21,21,22
 22   DO 23 J=1,N1
      IF(IN(J))24,24,23
 24    K=(I-1)*N+J
       K1=(N+1)*J-N
       IF((A(K1)-A(K)*A(K)/A(K0))/V(I)-T)23,23,21
 23   CONTINUE
      IDF=1
      RETURN
 21   IDF=0
      RETURN
 30    IDF=0
      DT=0.
      K=0
       J=1
       DO 1 I=1,N1
       V(I)=A(J)
       J=J+N+1
      IF(IND(I))1,41,1
 41    K=I
 1    IN(I)=IND(I)
       G=1.
      IF(K)18,18,19
 19   T1=G
       IDF=IDF+1
      IN(K)=1
       L=K
 10   DO 11 I=1,K
       U(I)=A(L)
       A(L)=0.
 11    L=L+N
      X=U(K)
       L=L-N
      DO 12 I=K,N
       U(I)=A(L)
       A(L)=0.
 12    L=L+1
      U(K)=-1.
      DT=DT+ALOG(X)
       L=1
      DO 8 I=1,N
       P=U(I)/X
       DO 48 J=I,N
       A(L)=A(L)-U(J)*P
 48    L=L+1
 8     L=L+I
 14   G=T
       J=1
      DO 15 I=1,N1
      IF(IN(I))15,16,15
 16    H=A(J)/V(I)
      IF(H-G)15,15,17
 17   G=H
      K=I
 15    J=J+N+1
      IF(G-T)18,18,19
 18    L=N
       DO 31 I=1,N1
       IF(IN(I))32,32,31
 32    A(L)=0.
 31    L=L+N
       L=1
       DO 55 I=1,N
       K=L
       DO 50 J=I,N
       A(K)=A(L)
       L=L+1
 50    K=K+N
 55    L=L+I
      RETURN
      END
      FUNCTION NBITS(II)
      K=II
      NBITS=0
      DO 1 I=1,10
      IF(MOD(K,2).EQ.1) NBITS=NBITS+1
 1    K=K/2
      RETURN
      END
      INTEGER FUNCTION UNION(II,JJ)
      UNION=0
      L=1
      DO 1 I=1,12
      IF(MOD(II/L,2).NE.0 .OR. MOD(JJ/L,2).NE.0) UNION=UNION+L
 1    L=L*2
      RETURN
      END
      LOGICAL FUNCTION INCL(II,JJ)
      INCL=.FALSE.
      L=1
      DO 1 I=1,12
      IF(MOD(II/L,2).EQ.0) GO TO 1
      IF(MOD(JJ/L,2).EQ.0) RETURN
 1    L=L*2
      INCL=.TRUE.
      RETURN
      END
      SUBROUTINE CONV(A,IN,N)
      DOUBLE PRECISION A
      DIMENSION C(10),IN(10)
      DATA  C/'0','1','2','3','4','5','6','7','8','9'/
      DO 1 I=1,N
      CALL GETCHR(A,I,B)
      IN(I)=0
      DO 1 J=1,10
 1    IF(B.EQ.C(J)) IN(I)=J-1
      RETURN
      END