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





















































