Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/bmd/bmd02s.for
There is 1 other file named bmd02s.for in the archive. Click here to see a list.
C CONTINGENCY TABLE ANALYSIS JUNE 15, 1966
C THIS IS A SIFTED VERSION OF BMD02S ORIGINALLY WRITTEN IN
C FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
DOUBLE PRECISION PPPP,PN,POP,PUP,YUV,XUV,PIP,PAP,Q004HL,VN,VV
DIMENSION X(351),P(2000),LPI(351),MP(351),JV(351),KV(351),VN(351),
1LIA(200),IA(200),IB(200),IC(15,200),LB(200),BT(200),IX(1024),F(162
2),NR(351),Q(351)
COMMON ICOL , RP , CP , TP
C
75 FORMAT(50H1BMD02S - CONTINGENCY TABLE ANALYSIS - VERSION OF ,
118HJUNE 15, 1966 ,/
2 41H HEALTH SCIENCES COMPUTING FACILITY, UCLA //
3 23H PROBLEM IDENTIFICATION 3X,A6)
DATA POP,PUP,PIP,PAP,ONO/6HPROBLM,6HFINISH,6HINTRVL,6HTABLES,2HNO/
NP=2000
MV1=351
MD=15
MNP = NP
MLL=200
CALL USAGEB('BMD02S')
C NPR=NO.OF WORDS PER PHYSICAL BINARY RECORD
NPR=255
MT=5
2 MTP=MT
READ (5,1)PPPP,PN,NV,NO,NS,NTC,ICOL,TP,RP,CP,MT,QV,NF
1 FORMAT(2A6,I3,I5,2I3,I2,3A3,I2,A2,29X,I2)
IF(MT)110,110,109
110 MT=5
109 IF(PPPP.EQ.POP)GO TO 101
IF(PPPP.EQ.PUP)GO TO 103
80 XUV=PPPP
YUV=POP
GO TO 91
90 MTP=MT
91 WRITE (30,72)YUV,XUV
72 FORMAT(21H0CONTROL CARD ERROR, A6,40H CARD EXPECTED, FIRST 6 CHARA
1CTERS WERE A6)
103 IF(MTP-5)106,115,106
85 WRITE (30,86)
86 FORMAT(28H0INTERVAL CARDS OUT OF ORDER )
GO TO 103
106 REWIND MTP
115 STOP
101 WRITE (30,75)PN
IF((MT-3)*(MT-6))202,201,202
201 WRITE (30,203)MT
203 FORMAT(5H0TAPEI3,65HMAY NOT BE USED AS AN ALTERNATE INPUT TAPE FOR
1 THIS PROGRAM )
GO TO 103
104 WRITE (30,4000) NS
4000 FORMAT (I4,' INTERVAL CARDS IS ILLEGAL')
GO TO 103
202 IF(MT-5)107,108,107
107 REWIND MT
108 IF(MT-MTP)111,114,111
111 IF(MTP.EQ.5)GO TO 114
112 REWIND MTP
114 II=0
IF((NF-1)*(9-NF))2345,3456,3456
2345 WRITE (30,4567)
4567 FORMAT(56H0NUMBER OF FORMAT CARDS INCORRECTLY SPECIFIED, 1 ASSUMED
1)
NF=1
3456 NF=18*NF
CALL OFILE(20,'DATA')
LPL=-7
LPN=-55
L=0
YUV=POP
IF(NV*(NV-MV1))7,7777,7777
7 DO 8 I=1,NV
8 KV(I)=0
YUV=PIP
11 IF(NS) 104,65,66
66 DO 13 J=1,NS
READ (5,14)XUV,V,VV,N,(X(I),I=1,13)
14 FORMAT(A6,F4.0,A6,I2,13F4.0)
IF(XUV.NE.PIP)GO TO 90
71 IF(X(3))15,16,15
15 IF(N-13)17,17,18
16 D=X(2)-X(1)
D = D-D*1.E-6
DO 19 I=3,N
19 X(I)=X(I-1)+D
GO TO 17
18 IF(N.GT.23) WRITE(30,9973)
9973 FORMAT(58H NO OF INTERVALS SPECIFIED BY INTERVAL CARD CANNOT EXCEE
1D '23')
IF(N.GT.23 ) GO TO 103
185 READ (5,20)(X(I),I=14,N)
17 II=II+1
C VN(II)=NAME OF INTERNAL VARIABLE II
VN(II)=VV
C I1=EXTERNAL VARIABLE NO.
I1=V
IF((I1-1)*(NV-I1))1357,2468,2468
1357 MTP=MT
GO TO 7779
2468 KV(I1)=1
C KV(I1)=1 INDICATES EXTERNAL VARIABLE I1 HAS AT LEAST ONE INTERVAL CAR
20 FORMAT(15F4.0)
C LPN AND LPM ARE FOR SEQUENCE CHECKING INTERVAL CARDS FOR THE SAME
C EXTERNAL VARIABLE
LPM=(V-FLOAT(I1))*10.0+.5
IF(I1-LPL)87,84,87
84 IF(LPN+1-LPM)85,22,85
87 IF(LPM)85,21,85
21 MP(I1)=II
C LPI(II)=LOCATION IN P PRECEEDING FIRST BOUNDARY POINT FOR INTERNAL
C VARIABLE II
22 LPI(II)=L
LPL=I1
LPN=LPM
C JV(II)=EXTERNAL VARIABLE CORRESPONDING TO INTERNAL VARIABLE II
JV(II)=I1
DO 13 I=1,N
L=L+1
IF(MNP-L) 12,12,13
12 WRITE (30,4001)
1MNP
4001 FORMAT(11H0MORE THAN I4,23H INTERVALS NOT ALLOWED )
GO TO 1549
C P=LINEAR LIST OF ALL BOUNDARY POINTS
13 P(L)=X(I)
65 LLG=L
A=1.0
DO 25 J=1,10
L=L+1
P(L)=A
25 A=A+1.0
DO 23 I=1,NV
IF(KV(I))24,24,23
C INSERT AUTOMATIC BOUNDARY POINTS FOR UNSPECIFIED VARIABLES
24 II=II+1
JV(II)=I
DATA Q004HL/6H /
VN(II)=(+Q004HL)
LPI(II)=LLG
MP(I)=II
23 CONTINUE
IF(LPI(II)-LLG)235,236,236
235 LPI(II+1)=LLG
GO TO 237
236 LPI(II+1)=L
237 IF(II.LE.( MV1-1))GO TO 1246
MV11 = MV1-1
1348 WRITE (30,1548)MV11
1548 FORMAT(10H0MORE THAN I4, 10H VARIABLES )
1549 MTP=MT
WRITE (30,4002)
4002 FORMAT(19H0PROGRAM TERMINATED)
GO TO 103
1246 NIV=II
1247 NOPR=((1024/NPR)*NPR)/NIV
C NOPR=NO.OF OBSERVATIONS TO WRITE PER LOGICAL BINARY RECORD
DO 57 I=1,NIV
57 NR(I)=0
READ (5,30)(F(I),I=1,NF)
30 FORMAT(18A4)
WRITE (30,7776) (F(I),I=1,NF)
7776 FORMAT (' VARIABLE FORMAT IS'/(1X,18A4))
NOB=NO
28 NN=MIN0(NOB,NOPR)
N=0
DO 34 K=1,NN
READ (MT,F)(X(I),I=1,NV)
C FORM INTERNAL VARIABLES (ONE FOR EACH CATEGORIZATION)
DO 29 I=1,NIV
L1=LPI(I)+1
L2=LPI(I+1)
L2=MAX0(L1-L2,0)*10+L2
J=JV(I)
DO 31 L=L1,L2
IF(X(J)-P(L))32,31,31
31 CONTINUE
L=L1
32 N=N+1
IX(N)=L-L1
IF(IX(N))29,33,29
C NR=REJECT COUNT
33 NR(I)=NR(I)+1
29 CONTINUE
34 CONTINUE
C WRITE INTERNAL VARIABLES ON SCRATCH TAPE
WRITE(20,8000) (IX(I),I=1,N)
NOB=NOB-NN
IF(NOB)35,35,28
35 END FILE 20
CALL IFILE(20,'DATA')
37 I=0
NCL=NTC
DO 150 I=1,NIV
I1=JV(I)
I2=I-MP(I1)
150 Q(I)=FLOAT(I1)+FLOAT(I2)/10.0
WRITE (30,120)
120 FORMAT(//3(32H VARIABLE NO. OF REJECTS 8X)//)
IF(NIV.LT.3) WRITE(30,121) (Q(I),VN(I),NR(I),I=1,NIV)
IF(NIV.LT.3) GO TO 122
NIV3=NIV/3
NOV=NIV-MOD(NIV,3)
WRITE (30,121)((Q(I),VN(I),NR(I),I=J,NOV,NIV3),J=1,NIV3)
121 FORMAT(3(F6.1,1X,A6,I8,19X))
IF(NIV-NOV)123,122,123
123 NOV=NOV+1
WRITE (30,124)(Q(I),VN(I),NR(I),I=NOV,NIV)
124 FORMAT(80X,F6.1,1X,A6,I8,11X)
122 LL=0
L=0
302 FORMAT(' COMPRESSED OUTPUT'/'1',16X,'VARIABLES',11X,'NUMBER',
12X,'OF',4X,'CHI SQUARE',3X,'DEGREES OF',3X,'CHI SQUARE/D.F.',
23X,'CONTINGENCY',5X,'-2 LOG (MAXIMUM'/36X,'EMPTY CELLS',18X,
3'FREEDOM',23X,'COEFFICIENT',3X,'LIKELIHOOD RATIO)')
IF(QV.NE.ONO)GO TO 301
300 TP=0.0
WRITE (30,302)
301 I=0
38 I=I+1
C READ TABLES CARDS TILL MLL(400) COLUMNS ARE SPECIFIED
READ (5,39)XUV,A,NB,(BT(J),J=1,NB)
39 FORMAT(A6,F4.0,I2,15F4.0)
IF(XUV.EQ.PAP)GO TO 73
82 YUV=PAP
GO TO 90
73 LIA(I)=L
I1=A
I2=(A-FLOAT(I1))*10.0+.5
C IA(I)=INTERNAL ROW VARIABLE FOR TABLE CARD
IA(I)=MP(I1)+I2
IF((IA(I)-1)*(NIV-IA(I)))1289,1290,1290
1289 WRITE (30,1291) A,I
1291 FORMAT (' ROW VARIABLE ',F5.0,' ON TABLES CARD ',I3,' IS ILLEGAL'
2)
MTP=MT
GO TO 103
1290 DO 40 J=1,NB
I1=BT(J)
I2=(BT(J)-FLOAT(I1))*10.0+.5
L=L+1
C IB(L)=INTERNAL COLUMN VARIABLE
IB(L)=MP(I1)+I2
IF((IB(L)-1)*(NIV-IB(L)))1292,40,40
1292 WRITE (30,1293) BT(J),I
1293 FORMAT (' COLUMN VARIABLE ',F5.0,' ON TABLES CARD ',I3,' IS ILLEGA
2L')
MTP=MT
GO TO 103
40 CONTINUE
L0=LIA(I)+1
41 DO 42 J=L0,L
K=IB(J)
LB(J)=LL
LLG=LPI(K+1)-LPI(K)
LL=LL+10*MAX0(1-LLG,0)+LLG
IF(LL-MLL)42,42,43
42 CONTINUE
NCL=NCL-1
IF(NCL)44,44,38
44 LSW=0
LIA(I+1)=L
GO TO 45
43 LS1=J
LS2=L
LIA(I+1)=J-1
LSW=1
45 NIA=I
DO 60 I=1,MD
DO 60 J=1,MLL
60 IC(I,J)=0
NOB=NO
46 NL=MIN0(NOB,NOPR)
NN=NL*NIV
READ(20,8000) (IX(I),I=1,NN)
8000 FORMAT(20G)
C NIA=NO.OF ROW VARIABLES NOW IN LIST
DO 47 I=1,NIA
JA=IA(I)
L1=LIA(I)+1
L2=LIA(I+1)
IF(L2-L1)47,48,48
48 DO 49 L=JA,NN,NIV
LL=IX(L)
C IX(L)=0 IMPLIES REJECT
IF(LL)49,49,50
50 LLL=L-JA
DO 51 N=L1,L2
JB=IB(N)
LLLL=LLL+JB
MM=IX(LLLL)
IF(MM)51,51,52
52 M=LB(N)+MM
C FORM FREQUENCY TABLES IN IC
IC(LL,M)=IC(LL,M)+1
51 CONTINUE
49 CONTINUE
47 CONTINUE
NOB=NOB-NL
IF(NOB)53,53,46
C SET TAPE FOR NEXT PASS IF NEEDED
53 CALL OFILE(20,'DATA')
DO 54 I=1,NIA
JA=IA(I)
L1=LIA(I)+1
L2=LIA(I+1)
J1=LPI(JA)+1
LLG=LPI(JA+1)-LPI(JA)
N1=10*MAX0(1-LLG,0)+LLG
IF(L2-L1)54,55,55
55 DO 56 L=L1,L2
JB=IB(L)
LU=LB(L)+1
C LU=COLUMN IN IC WHICH IS FIRST COLUMN FOR THIS TABLE
J2=LPI(JB)+1
LLG=LPI(JB+1)-LPI(JB)
N2=10*MAX0(1-LLG,0)+LLG
C N2=NO.OF COLUMNSIN TABLE
C PERFORM CALCULATIONS FOR ONE TABLE
56 CALL CHI(IC(1,LU),Q(JA),Q(JB),N1,N2,P(J1),P(J2),VN(JA),VN(JB),MD)
54 CONTINUE
IF(LSW)2,2,58
58 L=0
C RESTORE UNUSED TABLE SPECIFIERS THEN READ MORE TABLE CARDS
DO 59 J=LS1,LS2
L=L+1
59 IB(L)=IB(J)
L0=1
LL=0
I=1
IA(1)=IA(NIA)
GO TO 41
7777 WRITE (30,7778) NV,MV1
7778 FORMAT (1X,I4,' VARIABLES IS ILLEGAL. MUST BE LESS THAN '
2 ,I4)
GO TO 103
7779 WRITE (30,7780) I1
7780 FORMAT (1X,I5,' IS ILLEGAL INDEX OF VARIABLE TO BE CATEGORIZED AS
2SPECIFIED ON INTERVL CARD')
GO TO 103
C
END
C SUBROUTINE CHI FOR BMD02S JUNE 15, 1966
C PERFORM CALCULATION ON TABLE MC AND PRINT RESULTS
SUBROUTINE CHI(MC,VB,VA,NR,NC,TB,TA,BN,AN,MD)
DOUBLE PRECISION QQ1,QQ3,A2,Q1,A,BN,AN
DIMENSION MC(MD,MD),TA(2),TB(2),JT(23,23),KC(23),KR(23),IR(23),IC(
123),R(23),C(23),PA(23),PB(23),CC(23),RR(23)
COMMON ICOL , RP , CP , TP
ITT=1000
IF(TP)300,301,300
300 WRITE (30,101)
101 FORMAT(20H1ORIGINAL CATEGORIES )
301 NN=0
DATA QQ1,QQ2,QQ3,QQ4,A1,A2,YES/6H ORIGI,3HNAL,6H COLLA,4HPSED,
13HROW,6HCOLUMN,3HYES/
Q1=QQ1
Q2=QQ2
DO 6 I=1,MD
PA(I)=TA(I)
6 PB(I)=TB(I)
NR1=NR-1
NC1=NC-1
DO 10 I=1,NR1
10 IR(I)=0
DO 11 J=1,NC1
11 IC(J)=0
IT=0
DO 12 J=1,NC1
DO 13 I=1,NR1
MC(I,J)=MIN0(MC(I,J),9999)
IC(J)=IC(J)+MC(I,J)
13 IR(I)=IR(I)+MC(I,J)
12 IT=IT+IC(J)
TT=IT
8 NNC=NC1+1
NEMT=0
NNR=NR1+1
A=A1
IF(TP)302,303,302
302 WRITE (30,1)A,VB,BN,(PB(I),I=1,NNR)
1 FORMAT(/1X,A6,9H VARIABLE F6.1,1X,A6,16H BOUNDARY POINTS 10F7.2/(4
15X,10F7.2))
A=A2
WRITE (30,1)A,VA,AN,(PA(I),I=1,NNC)
83 WRITE (30,5)
5 FORMAT(16H0FREQUENCY TABLE )
IF (IT.GT.0) GO TO 82
WRITE (30,8001)
8001 FORMAT (/ 6X, 'THERE ARE NO CASES IN THE TABLE')
GO TO 93
82 WRITE (30,33)(I,I=1,NC1)
33 FORMAT(/3X,23I5)
WRITE (30,30)
30 FORMAT(/)
DO 3 I=1,NR1
3 WRITE (30,31)I,(MC(I,J),J=1,NC1),IR(I)
WRITE (30,32)(IC(I),I=1,NC1),IT
32 FORMAT(3X,23I5)
31 FORMAT(I3,23I5)
303 RML = 0.0
IF(TT .GT. 0.0) RML = -TT * ALOG(TT)
IF (TT.GT.0.0) TT = TT
DO 80 J=1,NC1
CC(J)=IC(J)
IF(CC(J))80,80,40
40 RML=RML+ALOG(CC(J))*CC(J)
80 CONTINUE
CH=-1.0
DO 14 I=1,NR1
RR(I)=IR(I)
IF(RR(I))85,85,41
41 RML=RML+ALOG(RR(I))*RR(I)
85 DO 14 J=1,NC1
B=MC(I,J)
IF(B)400,400,43
400 NEMT=NEMT+1
GO TO 14
43 RML=RML-ALOG(B)*B
CH=CH+B*B/(CC(J)*RR(I))
14 CONTINUE
81 CHIS=CH*TT
IDF=(NR1-1)*(NC1-1)
IF (IDF.LE.0) GO TO 810
IF(IDF.LE.0 .OR. CHIS.LT.0.0) GO TO 810
CH=CHIS/FLOAT(IDF)
CT=SQRT(CHIS/(CHIS+TT))
811 RML=-2.0*RML
IF(TP)305,306,305
810 CHIS = 0.0
CH = 0.0
CT=0.0
GO TO 811
305 WRITE (30,15)CHIS,IDF,CH,CT
15 FORMAT(10X,'CHI-SQUARE',11X,F10.4/6X,'DEGREES OF FREEDOM',
19X,I9/8X,'CHI-SQUARE/D.F.',7X,F10.4/5X,'CONTINGENCY COEFFICIENT',
34X,F10.4)
WRITE (30,46)RML
46 FORMAT(30H -2 LOG(MLR) F10.4,4X,
134H (MLR = MAXIMUM LIKELIHOOD RATIO) )
C TEST FOR PRINTING OF PERCENTAGES
777 IF(TP.EQ.YES)GO TO 71
GO TO 50
306 WRITE (30,307)Q1,Q2,VB,BN,VA,AN,NEMT,CHIS,IDF,CH,CT,RML
307 FORMAT(A6,A4,1X,2(F6.1,1X,A6),I6,F15.4,6X,I6,F16.4,F17.4,F18.4)
Q1=QQ3
Q2=QQ4
GO TO 70
50 IF(CP.EQ.YES)GO TO 61
60 IF(RP.EQ.YES)GO TO 51
70 IF(-NN)93,91,91
90 CONTINUE
93 RETURN
91 NN=1
IJM=ICOL*IT
CALL COLAP(MC,NR1,NC1,IR,IC,IJM,PA,PB,MD)
IF(NR1) 90,90,92
92 IF(TP)310,8,310
310 WRITE (30,95)ICOL
95 FORMAT(52H1COLLAPSED CATEGORIES (EXPECTED NUMBERS GREATER THAN
1 I3,1H))
C GO RECALCULATE FOR COLAPSED TABLE
GO TO 8
51 DO 53 I=1,NR1
DO 52 J=1,NC1
IF (CC(J)) 600,601,600
601 JT(I,J)=0
GO TO 52
600 JT(I,J)=FLOAT(MC(I,J))/CC(J)*1000.0+.5
52 CONTINUE
IF (TT) 602,603,602
603 KR(I)=0
GO TO 53
602 KR(I)=RR(I)/TT*1000.0+.5
53 CONTINUE
DO 54 J=1,NC1
IF(CC(J)) 604,605,604
605 KC(J)=0
GO TO 54
604 KC(J)=CC(J)/CC(J)*1000.0+.5
54 CONTINUE
WRITE (30,56)
M=1
GO TO 100
61 DO 63 I=1,NR1
DO 62 J=1,NC1
IF (RR(I)) 606,607,606
607 JT(I,J)=0
GO TO 62
606 JT(I,J)=FLOAT(MC(I,J))/RR(I)*1000.0+.5
62 CONTINUE
IF (RR(I)) 608,609,608
609 KR(I)=0
GO TO 63
608 KR(I)=RR(I)/RR(I)*1000.0+.5
63 CONTINUE
DO 64 J=1,NC1
IF (TT) 610,611,610
611 KC(J)=0
GO TO 64
610 KC(J)=CC(J)/TT*1000.0+.5
64 CONTINUE
M=2
WRITE (30,66)
GO TO 100
71 DO 72 I=1,NR1
DO 73 J=1,NC1
IF (TT) 612,613,612
613 JT(I,J)=0
GO TO 73
612 JT(I,J)=FLOAT(MC(I,J))/TT*1000.0+.5
73 CONTINUE
IF(TT) 614,615,614
615 KR(I)=0
GO TO 72
614 KR(I)=RR(I)/TT*1000.0+.5
72 CONTINUE
DO 74 J=1,NC1
IF (TT) 616,617,616
617 KC(J)=0
GO TO 74
616 KC(J)=CC(J)/TT*1000.0+.5
74 CONTINUE
WRITE (30,76)
M=3
100 WRITE (30,33)(I,I=1,NC1)
WRITE (30,30)
DO 35 I=1,NR1
35 WRITE (30,31)I,(JT(I,J),J=1,NC1),KR(I)
WRITE (30,32)(KC(I),I=1,NC1),ITT
GO TO (70,60,50),M
76 FORMAT(//40H TABLE PERCENTAGES (TENTHS OF A PERCENT))
56 FORMAT(//41H COLUMN PERCENTAGES (TENTHS OF A PERCENT))
66 FORMAT(//38H ROW PERCENTAGES (TENTHS OF A PERCENT))
RETURN
END
C SUBROUTINE COLAP FOR BMD02S JUNE 15, 1966
SUBROUTINE COLAP(MC,NR,NC,IR,IC,IJM,PA,PB,MD)
DIMENSION MC(MD,MD),PA(2),PB(2),IR(2),IC(2),MM(23)
L1=1
33 II = 1000
IK=0
DO 1 I=1,NR
IF(IR(I))1,1,2
2 IF(IR(I)-II)3,3,1
3 II=IR(I)
IK=I
1 CONTINUE
GO TO (31,32),L1
31 JJ = 1000
JK=0
DO 4 I=1,NC
IF(IC(I))4,4,5
5 IF(IC(I)-JJ)6,6,4
6 JJ=IC(I)
JK=I
4 CONTINUE
32 IF(II*JJ-IJM)7,7,50
7 IF(II-JJ)9,9,10
9 N=IK
11 N=N+1
IF(IR(N))12,11,12
12 M=IK
13 M=M-1
IF(IR(M))14,13,14
14 IF(M)16,16,17
16 IF(N-NR)21,21,50
17 IF(N-NR)19,19,20
19 IF(IR(M)-IR(N))20,21,21
20 IR(IK)=IR(IK)+IR(M)
IR(M)=0
GO TO 22
21 IR(N)=IR(N)+IR(IK)
IR(IK)=0
22 L1=2
GO TO 33
10 N=JK
41 N=N+1
IF(IC(N))42,41,42
42 M=JK
43 M=M-1
IF(IC(M))44,43,44
44 IF(M)46,46,47
46 IF(N-NC)61,61,50
47 IF(N-NC)49,49,70
49 IF(IC(M)-IC(N))70,61,61
70 IC(JK)=IC(JK)+IC(M)
IC(M)=0
GO TO 31
61 IC(N)=IC(N)+IC(JK)
IC(JK)=0
GO TO 31
50 DO 80 I=1,NR
80 MM(I)=0
N=0
DO 81 J=1,NC
DO 82 I=1,NR
82 MM(I)=MM(I)+MC(I,J)
IF(IC(J))81,81,83
83 N=N+1
IC(N)=IC(J)
PA(N+1)=PA(J+1)
DO 84 I=1,NR
MC(I,N)=MM(I)
84 MM(I)=0
81 CONTINUE
NC1=N
DO 90 I=1,NC1
90 MM(I)=0
N=0
DO 91 J=1,NR
DO 92 I=1,NC1
92 MM(I)=MM(I)+MC(J,I)
IF(IR(J))91,91,93
93 N=N+1
IR(N)=IR(J)
PB(N+1)=PB(J+1)
DO 94 I=1,NC1
MC(N,I)=MM(I)
94 MM(I)=0
91 CONTINUE
IF(NR-N+NC-NC1)96,95,96
95 NR=0
RETURN
96 NR=N
NC=NC1
RETURN
END