C	WESTERN MICHIGAN UNIVERSITY
C	TWOAOV.F4 (FILENAME ON LIBRARY DECTAPE)
C	TWOAOV, 1.9.2 (CALLING NAME, SUBLST #)
C	TWO WAY ANALYSIS OF VARIANCE
C	PROGRAMMED BY EVA GAINES AT WMU, LATER MODIFIED BY SAM ANEMA,
C	 RUSS BARR
C	STATISTICAL CONSULTANT - DR. M. STOLINE
C	FORWMU PROGRAMS USED:  TTYPTY, EXISTS, ALLCOR, PRINTS, DEVCHG
C	APLIB PROGRAMS USED:  GETFOR, IO, FISHER
C	LIBRARY DECTAPE PROGRAMS USED:  USAGE
C	INTERNAL SUBR. USED:  MNL, GAINE
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C---------------
C---------------MAX NO. OF OBS. IN ANY ONE CELL.  I.E., 400 IN METHOD
C--------------- 1 (SEE ST. 42 AND 42-1 IN SUBR. MNL)
C---------------ID(16) LIMITS OUTPUT IDENTIFICATION TO 80 CHARACTERS
C---------------JEFFMT LIMITS OBJECTIVE FORMAT TO 80 CHARACTERS
C---------------S(2) IS USED IN CALL MNL
	DIMENSION JEFFMT(16),Y(400),ID(16),S(1)
C---------------ASSIGNMENT OF I/O CHANNELS (B. GRANET)
      IUT=3
      INP=2
      NOUTD=-1
      IND=-4
C---------------PRINTS 5 LEVELS PER LINE ON FACTOR 2 FOR TERMINAL
C--------------- OR 12 FOR BATCH
      KS=5
C---------------TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB
      CALL TTYPTY(ICOD)
      IF(ICOD.EQ.-1)KS=12
      WRITE(NOUTD,9872)
9872  FORMAT(20X,'WMU TWO-WAY ANALYSIS OF VARIANCE'/)
C	CALL USAGE('TWOAOV')
C---------------IO DEALS WITH INPUT? AND OUTPUT? DIALOGUE
      CALL IO(IUT,NDV1,NOUTD,IND,1,ICOD)
1000  CALL IO(INP,NDV2,NOUTD,IND,0,ICOD)
C---------------DEALS WITH FORMAT DIALOGUE
      CALL GETFOR(NOUTD,IND,JEFFMT,ISTD,16,4)
      WRITE(NOUTD,8291)
8291  FORMAT(' ENTER OUTPUT IDENTIFICATION IF DESIRED, OTHERWISE RETURN.
     1'/)
      READ(IND,8292)ID
8292  FORMAT(16A5)
      WRITE(NOUTD,450)
450   FORMAT(' ENTER INPUT METHOD DESIRED.(1,2,OR 3)'/)
      READ(IND,1002)METHOD
      GO TO (1004,456,456),METHOD
1004  WRITE(NOUTD,1001)
1001  FORMAT(' NUMBER OF LEVELS ON FACTOR 1 = ',$)
      READ(IND,1002)N1
      IF(N1.LT.1)GO TO 7882
1002  FORMAT(12I)
      WRITE(NOUTD,1003)
1003  FORMAT(' NUMBER OF LEVELS ON FACTOR 2 = ',$)
      READ(IND,1002)N2
      IF(N2.LT.1)GO TO 7882
C	DYNAMIC ALLOCATION OF MEMORY
460	NM=MAX0(N1,N2)
	NMAX=4*N1+5*N2+2*N1*N2+N2*N2+NM*N2
	CALL ALLCOR(NMAX,IER,K1,S)
	IF(IER.LT.0)GO TO 7782
	IYS=K1
	IY1=IYS+NM*N2
	IY2=IY1+N1
	IS1=IY2+N2
	IS2=IS1+N1
	IB1=IS2+N2
	IYM1=IB1+N2
	IT=IYM1+N1
	IC=IT+N1*N2
	IAM=IC+N2
	IYT=IAM+N1
	IYM2=IYT+N1*N2
	IYR=IYM2+N2
	IXXX=IYR+N2*N2
	CALL MNL(N1,N2,NM,METHOD,IUT,INP,NOUTD,IND,KS,ICOD,NDV1,NDV2,
     1S(IYS),S(IY1),S(IY2),S(IS1),S(IS2),S(IB1),S(IYM1),S(IT),S(IC),
     2S(IAM),S(IYT),S(IYM2),S(IYR),JEFFMT,Y,ID,ISTD)
	GO TO 1000
7882	WRITE(NOUTD,7883)
7883	FORMAT(' IMPROPER VALUE',/)
	GO TO 1004
7782	WRITE(NOUTD,7783)
7783	FORMAT(' PROBLEM TOO LARGE',/)
	GO TO 1000
456	N1=12
	N2=12
	GO TO 460
	END
C---------------N1, N2, NM, METHOD, IUT, INP, NOUTD, INDP, KS,
C--------------- ICOD, NDV1, NDV2 ARE INPUT. OTHER ARGS. ARE SPACES
C--------------- RESERVED BY DYN. ALLOC.
	SUBROUTINE MNL(N1,N2,NM,METHOD,IUT,INP,NOUTD,IND,KS,ICOD,
     1NDV1,NDV2,YS,Y1,Y2,S1,S2,B1,YM1,T,C,AM,YT,YM2,YR,
     2JEFFMT,Y,ID,ISTD)
      DIMENSION YS(NM,N2),Y1(1),Y2(1),S1(1),S2(1),B1(1),YM1(1)
      DIMENSION T(N1,N2),C(1),AM(1),YT(N1,N2),JEFFMT(16),Y(400),ID(16)
      DIMENSION YM2(1),YR(N2,N2)
	N1S=N1
	N2S=N2
	NMS=NM
      DO 66 I=1,N1
      AM(I)=0.
      Y1(I)=0.
      S1(I)=0.
      DO 66 J=1,N2
      T(I,J)=0.
      YT(I,J)=0.
66    CONTINUE
	DO 68 J=1,N2
      Y2(J)=0.
      C(J)=0.
      B1(J)=0.
      S2(J)=0.
	DO 68 I=1,NM
68	YS(I,J)=0.
      RAB=0.
      SST=0.
      SS=0.
      G=0.
      SSC=0.
      SSAB=0.
      AN=0.
      SSBA=0.
C---------------INITIALIZE TO PRINT PRELIMINARY ANOVA (SEE LINE
C--------------- FOLLOWING "PRINT MEANS ETC." LOOP)
      ITYP=1
      TW=0.
      TC=0.
	GO TO(1004,456,456),METHOD
1004     WRITE(NOUTD,1020)
1020  FORMAT(' ENTER CELL SIZES'/)
      DO 5 I=1,N1
5     READ(IND,1008) (T(I,J),J=1,N2)
1008  FORMAT(12F)
456   IF(ISTD.EQ.1)JEFFMT(1)='(10F)'
457   GO TO (500,703,471),METHOD
7782  WRITE(NOUTD,7783)
7783   FORMAT(' PROBLEM IS TOO LARGE'/)
      GO TO 1004
C---------------METHOD 2 OF INPUT
703   IF(NDV2.EQ.'TTY')WRITE(NOUTD,700)
700   FORMAT(' ENTER DATA'/)
      IF(NDV2.NE.'TTY')WRITE(NOUTD,830)
830   FORMAT(' YOUR DATA IS BEING READ.'/)
      N1=1
      N2=1
      I=1
      J=1
701   READ(INP,702,END=1500)IAST
702   FORMAT(A1)
      IF(IAST.EQ.'*')GO TO 750
      REREAD JEFFMT,YIN
      YT(I,J)=YT(I,J)+YIN
      YS(I,J)=YS(I,J)+YIN*YIN
      T(I,J)=T(I,J)+1.0
      GO TO 701
750   REREAD 751,I,J
751   FORMAT(1X,2I)
      IF(I.LT.0.OR.I.GT.12.OR.J.LT.0.OR.J.GT.12)GO TO 752
      IF(I.EQ.0.OR.J.EQ.0)GO TO 1500
      IF(I.GT.N1)N1=I
      IF(J.GT.N2)N2=J
      GO TO 701
752   WRITE(NOUTD,753)
753   FORMAT(' INDEX OUT OF RANGE ENTRY IGNORED'//)
      GO TO 701
C---------------METHOD 3 OF INPUT
471   IF(NDV2.EQ.'TTY')WRITE(NOUTD,700) 
      IF(NDV2.NE.'TTY')WRITE(NOUTD,830)
      IF(ISTD.EQ.1)JEFFMT(1)='(2I,F'
      IF(ISTD.EQ.1)JEFFMT(2)=')    '
      N1=0
      N2=0
474   READ(INP,JEFFMT,END=1500)I,J,YIN
      IF(I.LT.0.OR.I.GT.12.OR.J.LT.0.OR.J.GT.12)GO TO 472
      IF(I.EQ.0.OR.J.EQ.0)GO TO 1500
      IF(I.GT.N1)N1=I
      IF(J.GT.N2)N2=J
      YT(I,J)=YT(I,J)+YIN
      YS(I,J)=YS(I,J)+YIN*YIN
      T(I,J)=T(I,J)+1.0
      GO TO 474
472   WRITE(NOUTD,753)
      GO TO 474
C---------------METHOD 1 OF INPUT
500   IF(NDV2.NE.'TTY')WRITE(NOUTD,830)
      DO 30 I=1,N1
      DO 25 J=1,N2
      IF(T(I,J))40,25,40
40    JK=T(I,J)
      IF(NDV2.EQ.'TTY')WRITE(NOUTD,1025)I,J
1025  FORMAT(' ENTER DATA FOR CELL (',I2,',',I2,')'/)
      READ (INP,JEFFMT)(Y(K),K=1,JK)
      DO 42 K=1,JK
      YT(I,J)=YT(I,J)+Y(K)
42    YS(I,J)=YS(I,J)+Y(K)*Y(K)
25    CONTINUE
30    CONTINUE
1500  IF(METHOD.NE.1.AND.(N1.GT.12.OR.N2.GT.12))GO TO 1005
C---------------PROCEED IF METHOD 2 OR 3 AND N1 AND N2 ARE NOT
C--------------- TOO LARGE
      DO 4250 I=1,N1
      DO 4251 J=1,N2
      IF(T(I,J).NE.0.0)GO TO 4250
4251  CONTINUE
      DO 4260 II=I,N1-1
      DO 4261 J=1,N2
      T(II,J)=T(II+1,J)
      YT(II,J)=YT(II+1,J)
4261  YS(II,J)=YS(II+1,J)
4260  CONTINUE
      N1=N1-1
4250  CONTINUE
      DO 4252 J=1,N2
      DO 4253 I=1,N1
      IF(T(I,J).NE.0.0)GO TO 4252
4253  CONTINUE
      DO 4262JJ=J,N2-1
      DO 4263 I=1,N1
      T(I,JJ)=T(I,JJ+1)
      YT(I,JJ)=YT(I,JJ+1)
4263  YS(I,JJ)=YS(I,JJ+1)
4262  CONTINUE
      N2=N2-1
4252  CONTINUE
      DO 31 I=1,N1
      DO 32 J=1,N2
      IF(T(I,J))41,32,41
41    TC=TC+1.0
      TW=TW+T(I,J)-1.0
      Y1(I)=Y1(I)+YT(I,J)
      S1(I)=S1(I)+T(I,J)
      SSC=SSC+(YT(I,J)*YT(I,J))/T(I,J)
      SST=SST+YS(I,J)
32    CONTINUE
      SSAB=SSAB+(Y1(I)*Y1(I))/S1(I)
      SS=SS+S1(I)
31    G=G+Y1(I)
      IWM=0
      ITYP=0
C---------------CHECK FOR EMPTY CELLS AND CALC. MEANS AND ST. DEV.
C--------------- FOR CELLS WITH MORE THAN ONE OBS.  SEE LINE #250 
C--------------- AND 825.
      DO 400 J=1,N2
      DO 300 I=1,N1
      IF(T(I,J))250,825,250
250   IF(T(I,J).NE.1.0)ITYP=1
      Y2(J)=Y2(J)+YT(I,J)
      S2(J)=S2(J)+T(I,J)
      YS(I,J)=YS(I,J)-(YT(I,J)*YT(I,J))/T(I,J)
      IF(T(I,J)-1.)15,300,15
15    YS(I,J)=YS(I,J)/(T(I,J)-1.)
      YS(I,J)=SQRT(YS(I,J))
      YT(I,J)=YT(I,J)/T(I,J)
      GO TO 300
C---------------WE HAVE AN EMPTY CELL
825   IWM=1
300   CONTINUE
400   SSBA=SSBA+(Y2(J)*Y2(J))/S2(J)
      DO 312 I=1,N1
312   YM1(I)=Y1(I)/S1(I)
      DO 313 I=1,N2
313   YM2(I)=Y2(I)/S2(I)
      G1=(G*G)/SS
      SST=SST-G1
      SSC=SSC-G1
      SSAB=SSAB-G1
      SSBA=SSBA-G1
      SSW=SST-SSC
      NP=0
      NP=NP+1
      IF(TC-1.)83,997,83
83    SMC=SSC/(TC-1.)
      NP=NP+1
      IF(TW)93,6780,93
93    SMW=SSW/TW
      NP=NP+1
6780  IF(SMW)53,6791,53
53    F=SMC/SMW
6791  KLL=FLOAT(N2)/FLOAT(KS)+.9
      IF(ID(1).NE.' ')WRITE(IUT,8293)ID
8293  FORMAT('1',16A5)
      K2=0
C---------------PRINT MEANS, SIZES, AND ST. DEV. OF CELLS
C--------------- AND MEANS OF LEVELS
      DO 760 KL=1,KLL
      K1=K2+1
      K2=MIN0(N2,K1+KS-1)
      WRITE(IUT,900)
900   FORMAT(//' ',35X,'FACTOR 2')
      WRITE(IUT,90)(I,I=K1,K2)
90    FORMAT(/6X,'LEVEL - - - - - - -',12(I4,5X))
      WRITE(IUT,91)
91    FORMAT(8X,'.')
      WRITE(IUT,92)(YM2(I),I=K1,K2)
92    FORMAT(8X,'.',4X,'MEAN - - -',12F9.2)
      WRITE(IUT,937)
937   FORMAT(8X,'.',5X,'.')
      WRITE(IUT,94)(T(1,J),J=K1,K2)
94    FORMAT(19X,'SIZE',12F9.2)
      WRITE(IUT,95)YM1(1),(YT(1,J),J=K1,K2)
95    FORMAT(8X,'1',1X,F8.2,' MEAN',12F9.2)
      WRITE(IUT,96)(YS(1,J),J=K1,K2)
96    FORMAT(18X,'ST DEV',12(F8.2,1X))
      IF(N1.EQ.1)GO TO 760
      WRITE(IUT,977)(T(2,J),J=K1,K2)
977    FORMAT(/1X,'FACTOR',12X,'SIZE',12F9.2)
      WRITE(IUT,98)YM1(2),(YT(2,J),J=K1,K2)
98    FORMAT(4X,'1',3X,'2',F9.2,' MEAN',12F9.2)
      WRITE(IUT,96)(YS(2,J),J=K1,K2)
      IF(N1.LT.3)GO TO 760
      DO 897 I=3,N1
      WRITE(IUT,898)(T(I,J),J=K1,K2)
898   FORMAT(/19X,'SIZE',12F9.2)
      WRITE(IUT,899)I,YM1(I),(YT(I,J),J=K1,K2)
899   FORMAT(7X,I2,F9.2,' MEAN',12F9.2)
897   WRITE(IUT,96)(YS(I,J),J=K1,K2)
760   CONTINUE
      IF(ITYP.EQ.0)GO TO 6781
C---------------PREL. ANOVA ONLY WHEN AT LEAST ONE CELL HAS
C--------------- MORE THAN ONE OBS.
80    WRITE (IUT,39)
39    FORMAT(//25X17HPRELIMINARY ANOVA)
      WRITE (IUT,49)
49    FORMAT(6X6HSOURCE,17X2HDF,8X2HSS,8X2HMS,8X1HF,6X,4HPROB)
      TP=TC-1.
      PROB=FISHER(IFIX(TP),IFIX(TW),F)
      WRITE (IUT,59)TP,SSC,SMC,F,PROB
59    FORMAT(6X5HCELLS,10X4F10.2,F9.3)
      TP=N1-1
      I=1
      J=2
      WRITE (IUT,69)I,J,TP,SSAB
69    FORMAT(4XI1,9H IGNORING,I2,5X2F10.2)
      TP=N2-1
      WRITE (IUT,69)J,I,TP,SSBA
      WRITE (IUT,79)TW,SSW,SMW
79    FORMAT(6X6HWITHIN,9X3F10.2)
      TP=SS-1.
      WRITE (IUT,89)TP,SST
89    FORMAT(6X5HTOTAL,10X2F10.2)
      IF(N1.EQ.1.OR.N2.EQ.1)GO TO 1000
C---------------START WORK ON LEAST SQ. ANOVA
C---------------PUT IDENTITY MATRIX IN YS FOR CALL GAINE BELOW,
C--------------- INIT. YR TO 0
6781  DO 33 I=1,N2
      DO 33 J=1,N2
      YR(I,J)=0.
      IF(I-J)133,37,133
37    YS(I,J)=1.
      GO TO 33
133   YS(I,J)=0.
33    CONTINUE
C---------------CALC. NORMAL (LEAST SQ.) EQS. COEFFS. IN YR AND C
      DO 180 I=1,N2
      DO 180 J=1,N2
      DO 180 K=1,N1
      IF(S1(K))180,999,180
999   WRITE(NOUTD,409)K,S1(K)
409   FORMAT(1X,3HS1(,I3,2H)=,F6.0)
180   YR(I,J)=YR(I,J)+(T(K,I)*T(K,J))/S1(K)
      DO 102 I=1,N2
102   YR(I,I)=YR(I,I)-S2(I)
      DO 120 J=1,N2
171   DO 110 K=1,N1
      IF(S1(K))110,992,110
992   WRITE(NOUTD,409)K,S1(K)
110   C(J)=C(J)+(T(K,J)*Y1(K))/S1(K)
      YR(N2,J)=1.
120   C(J)=C(J)-Y2(J)
      C(N2)=0.
C---------------CALC. INVERSE OF YR, RETURN RESULT IN YS
      CALL GAINE(YR,YS,N2,N2S,NMS,IUT)
C---------------SOLVE NORMAL (LEAST SQ.) EQS.
      DO 280 J=1,N2
      DO 280 K=1,N2
280   B1(J)=B1(J)+YS(J,K)*C(K)
      DO 210 I=1,N1
      DO 220 J=1,N2
220   AM(I)=AM(I)+T(I,J)*B1(J)
      AM(I)=Y1(I)-AM(I)
      IF(S1(I))13,993,13
993   WRITE(NOUTD,409)I,S1(I)
13    AM(I)=AM(I)/S1(I)
210   AN=AN+AM(I)
      P1=N1
      AN=AN/P1
      DO 230 I=1,N1
      A=AM(I)-AN
230   RAB=RAB+A*Y1(I)
      DO 260 J=1,N2
260   RAB=RAB+B1(J)*Y2(J)
      RAB=RAB+AN*G-G1
      SI=SSC-RAB
      P2=N2
      TI=TC-P1-P2+1.
      MP=0
      MP=MP+1
      IF(TI)77,998,77
77    SMI=SI/TI
      TAB=RAB-SSBA
      MP=MP+1
      IF(P1-1.)87,998,87
87    TMAB=TAB/(P1-1.)
      TBA=RAB-SSAB
      MP=MP+1
      IF(P2-1.)97,998,97
97    TMBA=TBA/(P2-1.)
      MP=MP+1
      IDF=IFIX(TW)
      IF(ITYP.EQ.0)SMW=SMI
      IF(ITYP.EQ.0)IDF=IFIX(TI)
      IF(SMW)107,998,107
107   FI=SMI/SMW
      FA=TMAB/SMW
      FB=TMBA/SMW
      WRITE (IUT,109)
109   FORMAT(//25X,'LEAST SQUARES ANOVA')
      WRITE (IUT,49)
      TP=TC-1.
      WRITE (IUT,59)TP,SSC
      TP=P1-1.
      I=1
      J=2
      PROB=FISHER(IFIX(TP),IDF,FA)
      WRITE (IUT,204)I,J,TP,TAB,TMAB,FA,PROB
204   FORMAT(3XI1,12H ELIMINATING,I2,3X4F10.2,F9.3)
      TP=P2-1.
      PROB=FISHER(IFIX(TP),IDF,FB)
      WRITE (IUT,304)J,I,TP,TBA,TMBA,FB,PROB
304   FORMAT(3XI1,12H ELIMINATING,I2,3X4F10.2,F9.3)
      PROB=FISHER(IFIX(TI),IDF,FI)
C---------------ITYP=0 MEANS EVERY CELL HAS ONE OBS., ITYP=1
C--------------- MEANS AT LEAST ONE CELL HAS MORE THAN ONE OBS.
      IF(ITYP.EQ.0)WRITE(IUT,404)I,J,TI,SI,SMI
      IF(ITYP.NE.0)WRITE (IUT,404)I,J,TI,SI,SMI,FI,PROB
404   FORMAT(6XI1,3H BY,I2,9X4F10.2,F9.3)
      IF(ITYP.EQ.0)GO TO 6783
      WRITE (IUT,79)TW,SSW,SMW
6783  TP=SS-1.
      WRITE (IUT,89)TP,SST
C
C	CHECK FOR BALANCED CASE
C
	T11=T(1,1)
	DO 821 I=1,N1
	DO 821 J=1,N2
	IF(T11.NE.T(I,J))GO TO 823
821	CONTINUE
	WRITE(IUT,822)
822	FORMAT(/,' THE DESIGN IS BALANCED. IN THIS CASE THE WEIGHTED',/,
     #' MEANS AND LEAST SQUARES ANOVA''S ARE IDENTICAL.',//)
	GO TO 1000
C
C	CHECK FOR PROPORTIONALITY
C
823	DO 824 I=2,N1
	TI1=T(I,1)
	DO 824 J=2,N2
	IF(T11*T(I,J).NE.T(1,J)*TI1)GO TO 826
824	CONTINUE
	WRITE(IUT,831)
831	FORMAT(/,' THE DESIGN IS PROPORTIONAL, BUT NOT BALANCED.',/,
     #' IN THIS CASE THE LEAST SQUARES ANOVA GIVES THE',/,
     #' CORRECT SUM OF SQUARES AND MEAN SQUARES, BUT NOT',/,
     #' CORRECT F-TESTS. SEE T.A. BANCROFT "TOPICS IN',/,
     #' INTERMEDIATE STATISTICAL METHODS(P. 14-15)" FOR',/,
     #' THE PROPER F-TEST PROCEDURES FOR THE PROPOR-',/,
     #' TIONAL CASE.',//)
	GO TO 1000
C
C        WEIGHTED MEANS ANALYSIS
C
C---------------WE HAVE NON-PROPORTIONALITY.  IWM=1 MEANS AT LEAST
C--------------- ONE EMPTY CELL, SEE LINE ABOVE #250
826      IF(IWM.EQ.1)GO TO 870
      DO 811 I=1,N1
      YM1(I)=0.0
811   S1(I)=0.0
      DO 812 J=1,N2
      YM2(J)=0.0
812   S2(J)=0.0
      DO 813 I=1,N1
      DO 814 J=1,N2
      YM1(I)=YM1(I)+YT(I,J)
      S1(I)=S1(I)+1.0/T(I,J)
      YM2(J)=YM2(J)+YT(I,J)
814   S2(J)=S2(J)+1.0/T(I,J)
      YM1(I)=YM1(I)/FLOAT(N2)
813   S1(I)=FLOAT(N2**2)/S1(I)
      DO 815 J=1,N2
      YM2(J)=YM2(J)/FLOAT(N1)
815   S2(J)=FLOAT(N1**2)/S2(J)
      A=0.0
      B=0.0
      D=0.0
      DO 816 I=1,N1
      A=A+S1(I)*YM1(I)**2
      B=B+S1(I)*YM1(I)
816   D=D+S1(I)
      SSR=A-B**2/D
      A=0.0
      B=0.0
      D=0.0
      DO 817 J=1,N2
      A=A+S2(J)*YM2(J)**2
      B=B+S2(J)*YM2(J)
817   D=D+S2(J)
      SSCO=A-B**2/D
      SMR=SSR/FLOAT(N1-1)
      SMCO=SSCO/FLOAT(N2-1)
C
C         OUTPUT ANOVA TABLE FOR WEIGHTED MEANS
C
      WRITE(IUT,818)
818   FORMAT(//25X,'WEIGHTED MEANS ANOVA')
      WRITE(IUT,49)
      TP=TC-1
      WRITE(IUT,59)TP,SSC,SMC
      TP=P1-1.0
      F=SMR/SMW
      PROB=FISHER(IFIX(TP),IDF,F)
      I=1
      J=2
      WRITE(IUT,819)I,TP,SSR,SMR,F,PROB
819   FORMAT(8X,I1,12X,4F10.2,F9.3)
      TP=P2-1.0
      F=SMCO/SMW
      PROB=FISHER(IFIX(TP),IDF,F)
      WRITE(IUT,819)J,TP,SSCO,SMCO,F,PROB
      PROB=FISHER(IFIX(TI),IDF,FI)
      IF(ITYP.EQ.0)WRITE(IUT,404)I,J,TI,SI,SMI
      IF(ITYP.NE.0)WRITE(IUT,404)I,J,TI,SI,SMI,FI,PROB
      IF(ITYP.EQ.0)GO TO 820
      WRITE(IUT,79)TW,SSW,SMW
 820  TP=SS-1.0
      WRITE(IUT,89)TP,SST
      GO TO 1000
870   WRITE(IUT,871)
871   FORMAT(//' NO WEIGHTED MEANS ANALYSIS CAN BE PERFORMED WITH'
     1,' EMPTY CELLS'//)
      GO TO 1000
1005  WRITE(NOUTD,1007)
1007  FORMAT(' PROBLEM TOO LARGE'/)
      IF(ICODE.EQ.-1)CALL EXIT
      GO TO 1000
C     VALUES OF NP MAY BE 1,2,3
C     1 IMPLIES TC-1=0,2 IMPLIES TW=0,3 IMPLIES SMW=0
997   WRITE(NOUTD,410)NP
410   FORMAT(1X,3HNP=,I2)
      GO TO(83,93,53),NP
C     VALUES OF MP MAY BE 1,2,3,4
C     1 IMPLIES TI=0,2 IMPLIES P1-1=0,3 IMPLIES P2-1=0,4 IMPLIES SMW=0
998   WRITE(NOUTD,1410)MP
1410  FORMAT(1X,3HMP=,I2)
      GO TO(77,87,97,107),MP
1000	RETURN
      END
C---------------A, N, IUT, N2 (USED IN DIM. ST.), NM (USED IN
C--------------- DIM. ST.) ARE INPUT.  B RETURNED
C--------------- CALCULATE INVERSE OF A, STORE INVERSE IN B
      SUBROUTINE  GAINE(A,B,N,N2,NM,IUT)
      DIMENSION A(N2,N2),B(NM,N2)
      I=0
100   I=I+1
      IF(A(I,I)-100.+100.)19,8,19
19    DIA=A(I,I)
C     MAKE DIAGONAL ELEMENT ONE
      DO 39 J=1,N
      A(I,J)=A(I,J)/DIA
39    B(I,J)=B(I,J)/DIA
      IF(I-N)29,18,29
C     MAKE ALL ELEMENTS BELOW THE DIAGONAL ELEMENT ZERO
29    M=I+1
      DO 11 L=M,N
      XA=A(L,I)
      DO 11 K=1,N
      A(L,K)=A(L,K)-XA*A(I,K)
11    B(L,K)=B(L,K)-XA*B(I,K)
      IF(I-1)18,100,18
C     MAKE ALL ELEMENTS ABOVE THE DIAGONAL ELEMENT ZERO
18    KP=0
      MN=I-1
      DO 40 J=1,MN
      KP=KP+1
      YA=A(KP,I)
      DO 40 L=1,N
      A(J,L)=A(J,L)-YA*A(I,L)
40    B(J,L)=B(J,L)-YA*B(I,L)
      IF(I-N)100,99,100
8     MN=I
48    MN=MN+1
      IF(MN-N)9,9,90
C     MAKE DIAGONAL ELEMENT NON-ZERO
9     DO 6 J=1,N
      A(I,J)=A(MN,J)+A(I,J)
6     B(I,J)=B(MN,J)+B(I,J)
      IF(A(I,I)-100.+100.)19,48,19
90    WRITE (IUT,80)
80    FORMAT(1X,26HTHE INVERSE DOES NOT EXIST)
99    RETURN
      END
