Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
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