C        ANALYSIS OF COVARIANCE WITH MULTIPLE COVARIATES   APRIL 1, 1966
C     MAIN ROUTINE
C        THIS IS A SIFTED VERSION OF BMD04V ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
C
C     CODE=PROBLEM IDENTIFICATION
C     KVF=NO. OF VARIABLE FORMAT CARDS
C     LYNN=NO. OF TREATMENT DEFINITION CARDS
C     NADD=NO. OF VARIABLES TO BE ADDED AFTER TG
C     NCARD=NO. OF CARDS/CASE
C     NCASE=NO. OF CASES
C     NCOV=NO. OF COVARIATES
C     NDV=NO. OF THE DEPENDENT VARIABLE
C     NGRPS=NO. OF TREATMENT GROUPS
C     NTAPE=LOGICAL TAPE NO. FOR INPUT DATA
C     NTG=NO. OF TG CARDS
C     REP=REPROCESS SAME DATA MATRIX
C
      DOUBLE PRECISION TMP
      DIMENSION NAMES(36),KKK(100),TMP(8),INDX(8)
      DOUBLE PRECISION NAMES,KKK
      DIMENSION SL(26),NEW(99),JUMP(99),NA(99),BN(99),NSAM(99),FMT(180),
     1DATA(100),NSUB(99,6),CONT(99,6),REL(99,6),OPER(99,6),NG(99),X(36),
     2SSVAR(99,36),TOTXX(36,36),TOTXY(36),TREXX(36,36),TREXY(36),
     3ERRXX(36,36),ERRXY(36),LLL(35),B(35),C(35),D(35),
     4SEB(35),SEC(35),SED(35),ADY(99),SEADY(99)
     5,INDEX1(100),INDEX2(36)
      COMMON  DATA   , NSUB   , CONT   , REL    , OPER   , NG
      COMMON  LH     , NEW    , JUMP   , NA     , BN     , IERROR
      COMMON  NVAR   , NTG    , NCASE
      COMMON /WT/B,C,D,SEB,SEC,SSVAR,ADY,SEADY
 102  FORMAT('1BMD04V - ANALYSIS OF COVARIANCE - MULTIPLE COVARIATES -
     1REVISED OCTOBER 31, 1968'/
     241H HEALTH SCIENCES COMPUTING FACILITY, UCLA//)
      ITAPE=2
	CALL USAGEB('BMD04V')
      REWIND ITAPE
      IT3=3
      NTAPE=5
	DOUBLE PRECISION ZQQQ
	DOUBLE PRECISION A124,A125,A126,A321,A421,A521,P123
      DATA A124/6HFINISH/,A125/6HSAMSIZ/,A126/6HGRPDEF/,A321/6HTRNGEN/
     1,A421/6HCVRSEL/,A521/6HLAB   /,P123/6HPROBLM/
       DOUBLE PRECISION CODE,TODE,KODE
      DATA A123/3HYES/
  999 READ (5,100)TODE,CODE,NGRPS,NCASE,NCARD,NCOV,NADD,NDV,OUT,INAM,REP
     1,LYNN,NTG,MTAPE,KVF
      NITP=0
      IF(TODE .NE.A124) GO TO 1
    2 IF(NTAPE .LE. 5) GO TO 4
      REWIND NTAPE
    4 STOP
    1 ZQQQ=P123
      IF(TODE.NE.P123) GO TO 446
      IF(NGRPS.LE.0) GO TO 456
      IF(NCOV.LT.2)GO TO 496
      GO TO 400
  446 WRITE (6,4446) ZQQQ,TODE
 4446 FORMAT(1H0,A6,60H CARD EXPECTED. THE FIRST SIX COLUMNS OF THE CARD
     1 READ WERE    ,A6)
      GO TO 2
  441 WRITE (6,4441)
 4441 FORMAT(40H0ILLEGAL NUMBER OF TRANSGENERATION CARDS)
      GO TO 2
 1441 WRITE(6,2441)
 2441 FORMAT(37H0VIOLATING RULE(1) ON PAGE 525 OF BMD)
      GO TO 2
 1442 WRITE(6,2442)
 2442 FORMAT(44H0ILLEGAL NUMBER IN COL.24-26 ON PROBLM CARD.)
      GO TO 2
  448 WRITE(6,4488)
 4488 FORMAT(55H0VIOLATING EITHER RULE(2) OR RULE(5) ON PAGE 525 OF BMD)
      GO TO 2
  491 WRITE(6,4910)
 4910 FORMAT(41H0NUMBER OF GROUP DEFINITION CARDS ILLEGAL)
      GO TO 2
  496 WRITE(6,4969)
 4969 FORMAT(35H0ILLEGAL NUMBER OF VARIATES READ IN)
      GO TO 2
  456 WRITE(6,4569)
 4569 FORMAT(35H0ILLEGAL NUMBER OF TREATMENT GROUPS)
      GO TO 2
  443 WRITE (6,4443)
 4443 FORMAT(50H0THE NUMBER OF THE INPUT TAPE CANNOT BE 2 OR 3       )
      GO TO 2
  401 WRITE (6,403)
      GO TO 2
 400  WRITE (6,102)
      WRITE (6,103)CODE,NGRPS,NCASE,NCARD,NCOV,NADD,NDV,OUT, INAM,REP,LY
     1NN,NTG,MTAPE,KVF
      IF((NCOV+NADD)*(NCOV+NADD-100).GE.0)GO TO 1441
      IF((NADD + 98) * (NADD - 35) .GE. 0) GO TO 1442
      IF((MTAPE-ITAPE)*(MTAPE-IT3) .EQ. 0) GO TO 443
      CALL TPWD(MTAPE,NTAPE)
      NVAR=NCOV
      KIT3=0
      REWIND IT3
      IERROR=0
      DO 332 I=1,99
  332 DATA(I)=0.0
      DO 333 I=1,36
      NAMES(I)=0
      INDEX1(I)=0
 333  INDEX2(I)=0
      LH=NVAR+1+NADD
      NVAR1=LH-1
      LR=NVAR
      IF(LH .GE. 38 .OR. LH .LE. 2) GO TO 441
      IF(KVF.GT.0.AND.KVF.LE.10)GO TO 410
      WRITE (6,9999)
      KVF = 1
 410  ASSIGN 50 TO NSKIP
      IF(NTG) 441,10,11
 11   ASSIGN 49 TO NSKIP
      IF(NTG.GT.99) GO TO 441
      WRITE (6,104)
      WRITE (6,105)
      ZQQQ=A321
      DO 411 I=1,NTG
      READ (5,106)TODE,NEW(I),JUMP(I),NA(I),BN(I)
      IF(TODE.NE.A321) GO TO 446
      WRITE (6,107)I,NEW(I),JUMP(I),NA(I),BN(I)
      IF((JUMP(I)-18)*(JUMP(I)-19)*(JUMP(I)-40).EQ. 0) GO TO 4115
      IF(JUMP(I)-41)4111,411,4115
 4111 IF((JUMP(I)-25)*(39-JUMP(I)) .LT. 0) GO TO 411
 4115 WRITE (6,150)I
      LYNN=-99
  411 CONTINUE
 10   IF(LYNN)491,14,13
   13 WRITE (6,110)
      DO 21 I=1,LYNN
      READ (5,111)TODE,NG(I),(NSUB(I,J),REL(I,J),CONT(I,J),OPER(I,J),J=1
     1,6)
      ZQQQ=A126
      IF(TODE.NE.A126) GO TO 446
      WRITE (6,112)NG(I)
      DO 22 J=1,6
      WRITE (6,113)NSUB(I,J),REL(I,J),CONT(I,J),OPER(I,J)
      DATA STAR/4H*   /
      IF(OPER(I,J) .EQ. STAR) GO TO 21
   22 CONTINUE
   21 CONTINUE
      GO TO 15
 14   NCARDS=(NGRPS+21)/22
      N2=0
      DO 1450 J=1,NCARDS
      N1=N2+1
      N2=N2+22
      IF(N2 .LE. 99) GO TO 1430
      N2=99
 1430 READ (5,108)TODE,(NSAM(I),I=N1,N2)
      ZQQQ=A125
      IF(TODE.NE.A125) GO TO 446
 1450 CONTINUE
 15   IF(INAM)401,17,501
  501 N1=0
  503 READ(5,504)TODE,(INDX(I),TMP(I),I=1,7)
  504 FORMAT(A3,3X,7(I4,A6))
      IF(TODE.NE.A521)GO TO 446
      DO 505 I=1,7
      IF(INDX(I).EQ.0)GO TO 505
      N1=N1+1
      INDEX1(N1)=INDX(I)
      KKK(N1)=TMP(I)
  505 CONTINUE
      IF(N1.LT.INAM)GO TO 503
   17 KVF=KVF*18
      WRITE(6,1090)
      READ (5,109)(FMT(I),I=1,KVF)
      WRITE(6,209)(FMT(I),I=1,KVF)
 209  FORMAT(1X,18A4)
      IF( REP .EQ.A123) GO TO 37
      K1=NCARD*NCASE
      KIT3=1
C
      IF(NTAPE.NE.5) GO TO 3777
      DO 35 I=1,K1
      READ (NTAPE,114) (X(J),J=1,20)
   35 WRITE (IT3,114) (X(J),J=1,20)
      END FILE IT3
      NLAST=IT3
   37 NTAPE=NLAST
 3777 KNTP=KIT3
      IF(KNTP.EQ.0) GO TO 36
      REWIND NTAPE
   36 IF(LYNN .GT. 0) GO TO 26
      NITP=1
      KNTP=1
      DO 27 I=1,NGRPS
      NN=NSAM(I)
      DO 27 J=1,NN
      READ (NTAPE,FMT)(DATA(K),K=1,NVAR)
      DATA(100)=I
 416  WRITE(ITAPE)(DATA(K),K=1,NVAR),DATA(100)
   27 CONTINUE
      GO TO 28
   26 DO 29 I=1,99
   29 NSAM(I)=0
      KNTP=1
      DO 30 I=1,NCASE
      DATA(100)=0.0
      READ (NTAPE,FMT)(DATA(J),J=1,NVAR)
	J=1
311      CALL GRPSEL(J)
      IF(J .GT. 0) GO TO 299
      READ (5,151)KODE,NCO
      IF(NITP.EQ.0) GO TO 88
      REWIND ITAPE
   88 IF(NCO .GE. 33) GO TO 999
      READ (5,151)KODE
      GO TO 999
 299  IF(DATA(100) .LE. 0.) GO TO 31
      KK=NG(J)
      NSAM(KK)=NSAM(KK)+1
      NITP=1
      WRITE(ITAPE)(DATA(K),K=1,NVAR),DATA(100)
      GO TO 30
31	J=J+1
	IF(J.LE.LYNN)GO TO 311
   30 CONTINUE
 28   END FILE ITAPE
      REWIND ITAPE
      IF(LYNN .LE. 0) GO TO 43
      JJ=0
      DO 39 I=1,99
   39 NG(I)=0
      DO 40 I=1,99
      IF(NSAM(I) .LE. 0) GO TO 40
      JJ=JJ+1
      NG(JJ)=NSAM(I)
      NSAM(I)=JJ
   40 CONTINUE
      JX=MAX0(JJ,NGRPS)
      NGRPS=JJ
      WRITE (6,115)
      DO 45 I=1,JX
   45 WRITE (6,116)I,NSAM(I),NG(I)
      GO TO 252
   43 DO 251 I=1,NGRPS
  251 NG(I)=NSAM(I)
  252 NVAR=NVAR+NADD
      NCOV=NCOV+NADD-1
      DO 46 I=1,NVAR
      TOTXY(I)=0.0
      TREXY(I)=0.0
      ERRXY(I)=0.0
      DO 46 J=1,NVAR
      TOTXX(I,J)=0.0
      TREXX(I,J)=0.0
   46 ERRXX(I,J)=0.0
      DO 47 I=1,NGRPS
      DO 47 J=1,NVAR
   47 SSVAR(I,J)=0.0
      ASSIGN 466 TO LSKIP
      ASSIGN 490 TO MUD
      N1=1
      N3=32
 463  READ (5,151)TODE,NCO,(INDEX2(I),I=N1,N3)
      ZQQQ=A421
      IF(TODE.NE.A421) GO TO 446
      IF(NCO.GT.35.OR.(NGRPS+NCO+1).GT.NCASE) GO TO 448
      GO TO LSKIP,(466,468)
  466 IF(NCO .LT. 33) GO TO 470
 467  ASSIGN 468 TO LSKIP
      N1=33
      N3=NCO
      GO TO 463
 468  INDEX2(36)=INDEX2(35)
      INDEX2(35)=INDEX2(34)
      INDEX2(34)=INDEX2(33)
      INDEX2(33)=NCO
      NCO=N3
 470  NCOV=NCO
      NVAR=NCO+1
 471  INDEX2(NVAR)=NDV
      IF( OUT .NE.A123) GO TO 475
      ASSIGN 417 TO MUD
      MESS=0
      WRITE (6,415)
      NITP=1
 475  DO 48 II=1,NGRPS
      LG=NG(II)
      DO 48 M=1,LG
      READ(ITAPE)(DATA(L),L=1,LR  ),DATA(100)
      GO TO MUD,(490,417)
  417 MESS=MESS+1
      WRITE (6,418)MESS
      WRITE (6,250)(DATA(K),K=1,LR),DATA(100)
 490  GO TO NSKIP,(49,50)
   49 CALL TG41(II,M)
 50   K=DATA(100)
   51 DO 48 J=1,NVAR
      I=INDEX2(J)
      SSVAR(K,J)=SSVAR(K,J)+DATA(I)
      DO 48 N=1,NVAR
      L=INDEX2(N)
      TOTXX(J,N)=TOTXX(J,N)+(DATA(I)*DATA(L))
   48 CONTINUE
      REWIND ITAPE
      NT=0
      DO 54 I=1,NGRPS
      NT=NT+NG(I)
   54 CONTINUE
      DO 55 J=1,NVAR
      X(J)=0.0
      DO 55 I=1,NGRPS
   55 X(J)=X(J)+SSVAR(I,J)
      BIG=NT
      DO 58 I=1,NVAR
      DO 58 J=1,NVAR
      TOTXX(I,J)=TOTXX(I,J)-(X(I)*X(J))/BIG
      DO 56 K=1,NGRPS
      T=NG(K)
      TREXX(I,J)=TREXX(I,J)+(SSVAR(K,I)*SSVAR(K,J))/T
   56 CONTINUE
   57 TREXX(I,J)=TREXX(I,J)-(X(I)*X(J))/BIG
   58 ERRXX(I,J)=TOTXX(I,J)-TREXX(I,J)
      DO 59 I=1,NGRPS
      D(1)=NG(I)
      DO 59 J=1,NVAR
   59 SSVAR(I,J)=SSVAR(I,J)/D(1)
      WRITE (6,117)
      DO 63 LT=1,NVAR
      I=INDEX2(LT)
      DO 615 J=1,INAM
      IF(I .EQ. INDEX1(J)) GO TO 625
 615  CONTINUE
      NITP=1
      WRITE(ITAPE,888)I
      REWIND ITAPE
      READ(ITAPE,777)NAMES(LT)
      REWIND ITAPE
      GO TO 63
 625  NAMES(LT)=KKK(J)
 63   CONTINUE
   61 DO 60 I=1,NGRPS
      WRITE (6,118)I
   60 WRITE (6,119)(NAMES(J),SSVAR(I,J),J=1,NVAR)
      WRITE (6,120)
      WRITE (6,121)
      CALL PATTY2(TOTXX,NVAR,NAMES)
      WRITE (6,120)
      WRITE (6,122)
      CALL PATTY2(TREXX,NVAR,NAMES)
      WRITE (6,120)
      WRITE (6,123)
      CALL PATTY2(ERRXX,NVAR,NAMES)
      DO 64 I=1,NVAR
      TOTXY(I)=TOTXX(I,NVAR)
      TREXY(I)=TREXX(I,NVAR)
   64 ERRXY(I)=ERRXX(I,NVAR)
      TOTYY=TOTXY(NVAR)
      TREYY=TREXY(NVAR)
      ERRYY=ERRXY(NVAR)
      CALL MATINV(TOTXX,NCOV,LLL,KKK)
      CALL MATINV(TREXX,NCOV,LLL,KKK)
      CALL MATINV(ERRXX,NCOV,LLL,KKK)
      WRITE (6,124)
      WRITE (6,121)
      CALL PATTY2(TOTXX,NCOV,NAMES)
      WRITE (6,124)
      WRITE (6,122)
      CALL PATTY2(TREXX,NCOV,NAMES)
      WRITE (6,124)
      WRITE (6,123)
      CALL PATTY2(ERRXX,NCOV,NAMES)
      DO 65 I=1,NCOV
      B(I)=0.0
      C(I)=0.0
      D(I)=0.0
      DO 65 J=1,NCOV
      B(I)=B(I)+ERRXX(I,J)*ERRXY(J)
      C(I)=C(I)+TOTXX(I,J)*TOTXY(J)
   65 D(I)=D(I)+TREXX(I,J)*TREXY(J)
      SSDTRE=0.0
      SSDTOT=0.0
      SSDERR=0.0
      DO 66 I=1,NCOV
      SSDERR=SSDERR+B(I)*ERRXY(I)
      SSDTRE=SSDTRE+D(I)*TREXY(I)
   66 SSDTOT=SSDTOT+C(I)*TOTXY(I)
      SSDTRE=SSDTOT-SSDERR
      SSATRE=TREYY-SSDTRE
      SSATOT=TOTYY-SSDTOT
      SSAERR=ERRYY-SSDERR
      N1=NGRPS-1
      N2=BIG
      N2=N2-NGRPS
      N3=BIG-1.0
      N4=N1
      N5=N2-NCOV
      N6=N3-NCOV
      TEST=SSATOT-SSAERR
      DN4=N4
      DN5=N5
      DN6=N6
      DN1=N1
      SMATRE=SSATRE/DN4
      SMAERR=SSAERR/DN5
      ADJTES=TEST/DN1
      FTEST=ADJTES/SMAERR
      SMATOT=SSATOT/DN6
      CALL WRT(N1,N2,N3,N5,N6,TREYY,ERRYY,TOTYY,SSDERR,SSAERR,SMAERR,
     *SSDTOT,SSATOT,TEST,FTEST,ADJTES,NCOV,SMATOT,BIG,NGRPS,NAMES,ERRXX,
     XTOTXX)
      NLAST=NTAPE
      GO TO 999
  403 FORMAT(48H0NEGATIVE NUMBER IS ILLEGAL FOR NUMBER OF LEBALS)
 1090 FORMAT(26H0THE VARIABLE FORMAT CARDS)
  100 FORMAT(2A6,I2,I5,2I2,I3,I2, A3,I3,A3,I2,27X,3I2)
  101 FORMAT(A3,3X,7(I4,A6))
  103 FORMAT(33H PROBLEM CODE....................,A6/
     1  33H NUMBER OF TREATMENT GROUPS......,I6/
     2  33H TOTAL NUMBER OF CASES...........,I6/
     3  33H NUMBER OF CARDS PER CASE........,I6/
     4  33H NUMBER OF ENTERING VARIABLES....,I6/
     5  33H NUMBER OF VARIABLES ADDED.......,I6/
     6  33H INDEX OF THE DEPENDENT VARIABLE.,I6/
     7  33H PRINT DATA INPUT................,3X,A3/
     8  33H NUMBER OF VARIABLES NAMED.......,I6/
     9  33H USE DATA CARDS OF LAST PROBLEM..,3X,A3/
     *  33H NUMBER OF GROUP DEFINITION CARDS,I6/
     1  33H NUMBER OF TRANSGENERATION CARDS.,I6/
     2  33H SPECIFIED BCD TAPE NUMBER.......,I6/
     3  33H NUMBER OF VARIABLE FORMAT CARDS.,I6)
  104 FORMAT(23H TRANS-GENERATION CARDS//)
  105 FORMAT(30H CARD NEW  TRANS OLD    CONST./33H  NO. VAR.  CODE VAR(A
     1) OR VAR(B)//)
  106 FORMAT(A6,I3,I2,I3,F6.0)
  107 FORMAT(I4,I5,I6,I5,4X,F9.4)
  108 FORMAT(A6,22I3)
  109 FORMAT(18A4)
  110 FORMAT(23H0GROUP DEFINITION CARDS)
 111  FORMAT(A6,I2,6(I2,A2,F5.0,A1))
  112 FORMAT(/20H0A CASE IS IN GROUP I2,4H IF,)
  113 FORMAT(5H VAR(I2,5H) IS A2,1H F9.4,1H A1)
  114 FORMAT(20A4)
  115 FORMAT(/26H     GROUP NO.      SAMPLE/25H SPECIFIED ASSIGNED* SIZE
     1//)
  116 FORMAT(5X,I2,7X,I2,4X,I5)
  117 FORMAT(///30H VARIABLE MEANS FOR EACH GROUP//)
  118 FORMAT(/11H GROUP NO. I2//)
 119  FORMAT(6(' (',A6,1X,F10.4,')')/)
  120 FORMAT(41H-SUMS OF SQUARES AND CROSS PRODUCT MATRIX)
  121 FORMAT(11H FOR TOTAL.//)
  122 FORMAT(15H FOR TREATMENT.//)
  123 FORMAT(11H FOR ERROR.//)
  124 FORMAT(26H-INVERSE OF THE COVARIATES/21H CROSS PRODUCT MATRIX)
 150  FORMAT(1H029X20HTRANSGENERATION CARDI3,35H HAS INVALID CODE FOR TH
     1IS PROGRAM.)
 151  FORMAT(A6,33I2)
C
  250 FORMAT(1H 10F11.4)
  415 FORMAT(//64H DATA INPUT (LAST VALUE FOR EACH CASE IS TREATMENT GRO
     1UP NUMBER))
  418 FORMAT(13H0CASE NUMBER I5)
  777 FORMAT(A6)
  888 FORMAT(3H X(,I2,1H))
 9876 FORMAT(20A4)
 9999 FORMAT(1H0,23X,71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPEC
     1IFIED, ASSUMED TO BE 1.)
C
      END
C             SUBROUTINE MATINV FOR BMD04V           APRIL  1, 1966
C        MATINV INVERTS AN NXN SYMMETRIC MATRIX
      SUBROUTINE MATINV(A,N,IN,V)
      DIMENSION A(36,36),IN(36),V(36),U(36)
      DO 10 I = 1,N
      V(I)=A(I,I)
   10 IN(I)=0
      D=1.0
      K=1
      DO 50 L = 1,N
      DO 20 I = 1,K
      U(I)=A(K,I)
   20 A(K,I)=0.0
      P=U(K)
      DO 30 I = K,N
      U(I)=A(I,K)
   30 A(I,K)=0.0
      T=H
      H=-1.E20
      IN(K)=1
      U(K)=-1.0
      D=D*P
      DO 50 I = 1,N
      Y=U(I)/P
      DO 40 J = 1,I
   40 A(I,J)=A(I,J)-Y*U(J)
      IF(IN(I) .GT. 0) GO TO 50
      IF(A(I,I) / V(I) .LE. H) GO TO 50
      H=A(I,I)/V(I)
      K=I
   50 CONTINUE
      DO 60 I = 1,N
      DO 60 J = 1,I
      A(I,J)=-A(I,J)
   60 A(J,I)=A(I,J)
      RETURN
      END
C             SUBROUTINE GRPSEL FOR BMD04V           APRIL  1, 1966
      SUBROUTINE GRPSEL(M)
      DIMENSIONDATA(100),NSUB(99,6),CONT(99,6),REL(99,6),OPER(99,6),B(6)
     1,NG(99),E(7)
      COMMON  DATA   , NSUB   , CONT   , REL    , OPER   , NG
	DATA A,O,STAR/'A','O','*'/
      DATA E/2HGT,2HLT,2HLE,2HGE,
     1 2HEQ,2HNE,2H  /
      FB=1.0
      ASSIGN 105 TO ISKIP
 500  DO 90 I=1,6
      IS=NSUB(M,I)
      IF(REL(M,I) .EQ. E(1)) GO TO 2
      IF(REL(M,I) .EQ. E(2)) GO TO 4
      IF(REL(M,I) .EQ. E(3)) GO TO 6
      IF(REL(M,I) .EQ. E(4)) GO TO 8
      IF(REL(M,I) .EQ. E(5)) GO TO 10
      IF(REL(M,I) .EQ. E(6)) GO TO 11
      IF(REL(M,I) .EQ. E(7)) GO TO 100
 96   WRITE (6,4000)REL(M,I)
      M=-6
      GO TO 104
    2 IF(DATA(IS)-CONT(M,I))50,50,20
    4 IF(DATA(IS)-CONT(M,I))20,50,50
    6 IF(DATA(IS)-CONT(M,I))20,20,50
    8 IF(DATA(IS)-CONT(M,I))50,20,20
   10 IF(DATA(IS)-CONT(M,I))50,20,50
   11 IF(DATA(IS)-CONT(M,I))20,50,20
   20 B(I)=1.0
      GO TO 90
   50 B(I)=0.0
   90 CONTINUE
  100 GO TO ISKIP,(105,125)
  125 IF(OPER(M-1,6) .NE. A) GO TO 150
      FB=FB*B(1)
      GO TO 175
  105 IF(OPER(M,1) .EQ. A) GO TO 147
      IF(OPER(M,1) .NE. O) GO TO 170
 147  FB=B(1)
      GO TO 175
  150 IF(OPER(M-1,6) .EQ. O) GO TO 152
      I = 6
      REL(M,6)=OPER(M-1,6)
      GO TO 96
 170  FB=B(1)
      GO TO 350
  152 IF(FB .LE. 0.) GO TO 147
 175  DO 300 I=1,5
      IF(OPER(M,I) .NE. A) GO TO 222
      FB = FB*B(I+1)
      GO TO 300
  222 IF(OPER(M,I) .NE. O) GO TO 333
      IF(FB .GT. 0.) GO TO 300
      FB=B(I+1)
  300 CONTINUE
      IF(OPER(M,6) .EQ. STAR) GO TO 350
      M=M+1
      ASSIGN 125 TO ISKIP
      GO TO 500
  333 IF(OPER(M,I) .EQ. STAR) GO TO 350
      REL(M,I)=OPER(M,I)
      GO TO 96
  350 IF(FB .LE. 0.) GO TO 104
      DATA(100) = NG(M)
  104 RETURN
 4000 FORMAT(1H014X30HILLEGAL OPERATION OR RELATION A2,57H. PROGRAM WILL
     1 PROCEED TO NEXT JOB, IF ANY, OR TERMINATE.)
      END
C             SUBROUTINE PATTY2 FOR BMD04V           APRIL  1, 1966
      SUBROUTINE PATTY2 (A,N,NAMES)
      DIMENSION NAMES(36),NN(8),A(36,36)
      DOUBLE PRECISION NAMES,NN
      IT=1
      KK=0
      K1=IT
      K2=MIN0(8,N)
   10 KK=KK+8
      IF(N .LE. KK) GO TO 20
      IT=IT+1
      GO TO 10
   20 DO 70 JX=1,IT
      LLL=K2-K1+1
      LL=0
      DO 30 JJ=K1,K2
      LL=LL+1
   30 NN(LL)=NAMES(JJ)
      WRITE (6,50) (NN(II),II=1,LLL)
      DO 40 I=1,N
   40 WRITE (6,60)NAMES(I),(A(I,J),J=K1,K2)
      K1=K2+1
      K2=K1+7
      K2=MIN0(K2,N)
   50 FORMAT(1H0,13X,A6,7(8X,A6)/)
   60 FORMAT(1H ,A6,1X,8F14.4)
   70 CONTINUE
      RETURN
      END
      SUBROUTINE TG41(II,M)
      DIMENSIONDATA(100),NSUB(99,6),CONT(99,6),REL(99,6),OPER(99,6),NG(9
     19),NEW(99),JUMP(99),NA(99),BN(99)
      COMMON  DATA   , NSUB   , CONT   , REL    , OPER   , NG
      COMMON  LH     , NEW    , JUMP   , NA     , BN     , IERROR
      COMMON  NVAR   , NTG    , NCASE
      EXTERNAL SIGN
      ASN(X) = ATAN(X/SQRT(1.-X**2))
      SAMPX=NCASE
      PI2=1.57079633
      DO 1010 I=1,NTG
      JOKE=JUMP(I)
      NW=NEW(I)
      IF(NW .GE. 100) GO TO 97
      NAA=NA(I)
      IF(NAA .GE. 100) GO TO 97
  100 NB=BN(I)
      BBN=BN(I)
      D1=DATA(NAA)
      IF(JOKE .LE. 17) GO TO 33
      JOKE=JOKE-2
      IF(JOKE .LE. 22) GO TO 33
      JOKE=23
   33 GO TO (1,2,3,4,5,6,7,8,9,10,111,111,111,111,150,111,170,180,190,20
     1,21,22,230),JOKE
    1 IF(   D1    )99,11,12
   11 D2=0.0
      GO TO 1000
   12 D2=SQRT(D1)
      GO TO 1000
    2 IF(   D1    )99,13,14
   13 D2=1.0
      GO TO 1000
   14    D2   =SQRT(   D1    )+SQRT(   D1    +1.0)
      GO TO 1000
    3 IF(   D1     .LE. 0.) GO TO 99
      D2=ALOG10(D1)
      GO TO 1000
    4 D2=EXP(D1)
      GO TO 1000
    5 IF(   D1    )99,11,17
   17 IF(   D1    -1.0)18,19,99
   19 D2=PI2
      GO TO 1000
   18 A=SQRT(D1)
      D2=ASN(A)
      GO TO 1000
    6 A=   D1    /(SAMPX+1.0)
      B=A+1.0/(SAMPX+1.0)
      IF(A)99,27,24
   27    D2   =ASN(SQRT(B))
      GO TO 1000
   24 IF(B)99,28,29
 28   IF(A-1.0)285,287,99
 285     D2   =ASN(SQRT(A))*2.0
      GO TO 1000
 287     D2   =3.14159265
      GO TO 1000
   29 A=SQRT(A)
      B=SQRT(B)
      IF(A .GT. 1. .OR. B .GT. 1.) GO TO 99
         D2   =ASN(A)+ASN(B)
      GO TO 1000
    7 IF(   D1     .EQ. 0.) GO TO 99
         D2   =1.0/D1
      GO TO 1000
    8 D2=D1-BBN
      GO TO 1000
    9 D2=D1*BBN
      GO TO 1000
   10 IF(   D1    )99,11,32
   32 D2=D1**BBN
      GO TO 1000
  111 IF(NB .GE. 100) GO TO 97
      NSKIP=JOKE-10
      GO TO(110,120,130,140,150,160),NSKIP
  110    D2   =   D1    +DATA(NB)
      GO TO 1000
  120    D2   =   D1    -DATA(NB)
      GO TO 1000
  130    D2   =   D1    *DATA(NB)
      GO TO 1000
  140 IF(DATA(NB) .EQ. 0.) GO TO 98
         D2   =   D1    /DATA(NB)
      GO TO 1000
  150 IF(   D1    -BBN)11,13,13
  160 IF(   D1    -DATA(NB))11,13,13
  170 IF(   D1     .LE. 0.) GO TO 99
         D2   =ALOG(D1)
      GO TO 1000
  180    D2   =SIN(D1)
      GO TO 1000
  190    D2   =COS(D1)
      GO TO 1000
   20 IF((   D1    -PI2)*(-PI2-   D1    ) .GT. 0.) GO TO 99
      D2=ATAN(D1)
      GO TO 1000
   21 IF(   D1    )99,11,38
   38 IF(NB .GE. 100) GO TO 97
      D2=D1**DATA(NB)
      GO TO 1000
   22 IF(BBN)99,11,40
   40 D2=BBN**D1
      GO TO 1000
  230 IF(   D1     .NE. 0.) GO TO 1000
      IF(SIGN(10.0,   D1    ))43,1000,1000
   43 D2=BBN
      GO TO 1000
 97   WRITE (6,4001)I
      GO TO 991
 98   NAA=NB
 99   WRITE (6,4000)NAA,JUMP(I),II,M
 991  WRITE (6,4002)
      GO TO 1010
 1000 DATA(NW)=D2
 1010 CONTINUE
  999 RETURN
 4000 FORMAT(17H0DATA OF VARIABLEI4,30H VIOLATES RESTRICTION FOR CODEI3,
     130H. THIS WAS FOR TREATMENT GROUPI3,7H SAMPLEI4)
 4001 FORMAT(20H0TRANSGENERATOR CARDI3,44H HAS VARIABLE NUMBERS TOO LARG
     1E FOR PROGRAM.)
 4002 FORMAT(32H THE DATA WILL REMAIN UNCHANGED.)
      END
C             SUBROUTINE TPWD FOR BMD04V             APRIL  1, 1966
      SUBROUTINE TPWD(NT1,NT2)
      IF(NT1)40,10,12
 10   NT1=5
   12 IF(NT1 .EQ. NT2) GO TO 19
      IF(NT2.EQ.5)GO TO 18
   15 REWIND NT2
   19 IF(NT1 .EQ. 5) GO TO 24
   18 IF(NT1 .EQ. 6) GO TO 40
      REWIND NT1
 24   NT2=NT1
      RETURN
 40   WRITE (6,49)
      STOP
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      END
      SUBROUTINE WRT(N1,N2,N3,N5,N6,TREYY,ERRYY,TOTYY,SSDERR,SSAERR,
     1SMAERR,SSDTOT,SSATOT,TEST,FTEST,ADJTES,NCOV,SMATOT,BIG,NGRPS,NAMES
     2,ERRXX,TOTXX)
      DOUBLE PRECISION NAMES
      DIMENSION DATA(100),NSUB(99,6),CONT(99,6),REL(99,6),OPER(99,6),
     *NG(99),X(36),SSVAR(99,36),TOTXX(36,36),ERRXX(36,36),B(35),C(35),
     1D(35),SEB(35),SEC(35),ADY(99),SEADY(99),NAMES(36)
      COMMON DATA,NSUB,CONT,REL,OPER,NG
      COMMON LH,NEW(99),JUMP(99),NA(99),BN(99),IERROR
      COMMON NVAR,NTG,NCASE
      COMMON /WT/B,C,D,SEB,SEC,SSVAR,ADY,SEADY
      WRITE(6, 125)
      WRITE(6, 126)
      WRITE(6, 127)
      WRITE(6, 128)
      WRITE(6, 129)
      WRITE(6, 127)
      WRITE(6, 130)
      WRITE(6, 127)
      WRITE(6, 131)
      WRITE(6, 132)N1,TREYY
      WRITE(6, 127)
      WRITE(6, 130)
      WRITE(6, 127)
      WRITE(6, 133)
      WRITE(6, 134)N2,ERRYY,SSDERR,SSAERR,N5,SMAERR
      WRITE(6, 130)
      WRITE(6, 127)
      WRITE(6, 131)
      WRITE(6, 136)
      WRITE(6, 137)N3,TOTYY,SSDTOT,SSATOT,N6
      WRITE(6, 127)
      WRITE(6, 130)
      WRITE(6, 127)
      WRITE(6, 138)TEST,N1,ADJTES
      WRITE(6, 127)
      WRITE(6, 126)
      WRITE(6, 139)
      WRITE(6, 140)N1,N5,FTEST
      WRITE(6, 135)
      WRITE(6, 142)
      WRITE(6, 152)
  125 FORMAT(1H135X,28HANALYSIS OF COVARIANCE TABLE/)
  126 FORMAT(1H 93(1H*))
  127 FORMAT(2H *,11X,1HI,7X,1HI,3(15X,1HI),7X,1HI,15X,1H*)
  128 FORMAT(2H *,11X,1HI,7X,1HI,15X,33HI   SUM-SQUARES I   SUM-SQUARES  
     *I, 7X,1HI,15X,1H*)
  129 FORMAT(94H *   SOURCE  I   DF  I       YY      I     (DUE)     I  
     *  (ABOUT)    I   DF  I   MEAN-SQUARE *)
  130 FORMAT(2H *,91(1H-),1H*)
 131  FORMAT(14H * TREATMENT I,7X,1HI,3(15X,1HI),7X,1HI,15X,1H*)
  132 FORMAT(15H * (BETWEEN) I ,I5,1H ,2HI ,F13.4,1H ,2(2HI ,13X,1H ),2H
     *I ,5X,3H I ,13X,2H *)
  133 FORMAT(14H *     ERROR I,7X,1HI,3(15X,1HI),7X,1HI,15X,1H*)
  134 FORMAT(15H *  (WITHIN) I ,I5,1X,3(2HI F13.4,1H ),2HI I5,3H I F13.4
     X,2H *)
  135 FORMAT(/22H0TABLE OF COEFFICIENTS//)
  136 FORMAT(14H * + ERROR   I,7X,1HI,3(15X,1HI),7X,1HI,15X,1H*)
  137 FORMAT(15H *   (TOTAL) I ,I5,1X,3(2HI F13.4,1H ),1HII6,2H I15X1H*)
 138  FORMAT(55H *DIFFERENCE FOR TESTING ADJUSTED TREATMENT MEANS....I F
     *13.4,3H I ,I5,3H I ,F13.4,2H *)
  139 FORMAT(//54H0NULL HYPOTHESIS. NO DIFFERENCE AMONG TREATMENTS AFTER
     */18X,26HADJUSTING WITH COVARIATES.//)
  140 FORMAT(18X,2HF(,I5,1H,I5,2H)=,F8.3//)
  142 FORMAT(1H 16X,19HTREATMENT (BETWEEN))
  152 FORMAT(1H 8X,11HCOEFFICIENT//)
      DO 68 I=1,NCOV
   68 WRITE(6, 144)NAMES(I),D(I)
      WRITE(6, 141)
      WRITE(6, 145)
      WRITE(6, 143)
      DO 69 I=1,NCOV
      SEB(I)=SMAERR*ERRXX(I,I)
      SEB(I)=SQRT(SEB(I))
      T=B(I)/SEB(I)
   69 WRITE(6, 144)NAMES(I),B(I),SEB(I),T
      WRITE(6, 141)
      WRITE(6, 146)
      WRITE(6, 143)
      DO 70 I=1,NCOV
      SEC(I)=SMATOT*TOTXX(I,I)
      SEC(I)=SQRT(SEC(I))
      T=C(I)/SEC(I)
  70  WRITE(6, 144)NAMES(I),C(I),SEC(I),T
      DO 821 I=1,NCOV
      X(I)=0.0
      DO 821 J=1,NGRPS
      T1=NG(J)
  821 X(I)=X(I)+SSVAR(J,I)*T1/BIG
      DO 72 J=1,NGRPS
      T=0.0
      T1=NG(J)
      T1=1.0/T1
      DO 71 I=1,NCOV
 71   T=T+B(I)*(SSVAR(J,I)-X(I))
      ADY(J)=SSVAR(J,NVAR)-T
      DO 73 I=1,NCOV
      DO 73 K=1,NCOV
   73 T1=T1+(SSVAR(J,I)-X(I))*ERRXX(I,K)*(SSVAR(J,K)-X(K))
 72   SEADY(J)=SQRT(SMAERR*T1)
      WRITE(6, 147)
      WRITE(6, 148)
      DO 76 I=1,NGRPS
   76 WRITE(6, 149)I,SSVAR(I,NVAR),ADY(I),SEADY(I)
  141 FORMAT(/60H0TABLE OF COEFFICIENTS,STANDARD ERRORS AND COMPUTED T-V
     *ALUES//)
  143 FORMAT(1H 8X,39HCOEFFICIENT   STAND. ERROR      T-VALUE//)
 144  FORMAT(1H ,A6,3(1X,F13.4))
  145 FORMAT(1H 17X,14HERROR (WITHIN))
  146 FORMAT(1H 11X,23HTREATMENT+ERROR (TOTAL))
  147 FORMAT(44H1TABLE OF ADJUSTED MEANS AND STANDARD ERRORS//)
  148 FORMAT(48H TREATMENT  TREATMENT    ADJUSTED    SE ADJUSTED/6H   NO
     *.7X,4HMEAN9X,4HMEAN/)
 149  FORMAT(1H ,I5,F16.4,2F13.4)
      RETURN
      END







































































