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