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