Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/cross/cross.for
There are 3 other files named cross.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C CROSS.F4 (FILE NAME ON LIBRARY DECTAPE)
C CROSS, 1.1.6 (CALLING NAME, SUBLST NO.)
C CROSS TABULATION WITH CHI-SQUARE
C THIS PROGRAM WAS TAKEN FROM THE BOOK "DATA PROCESSING"
C BY K. JANDA. SEVERAL MODIFICATIONS WERE MADE TO MAKE
C IT COMPATIBLE AND SOME ADDITIONAL FEATURES WERE ADDED BY
C SAM ANEMA.
C LIBRARY DECTAPE PROGS. USED: USAGE.MAC
C FORWMU PROGS. USED: TTYPTY, EXISTS, PRINTS, DEVCHG
C APLIB PROGS. USED: IO, GETFOR
C INTERNAL SUBR. USED: COEFT
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
DIMENSION LAD(200,40), MAX(40), KIND(40),KONST(40),MBRAK(16)
DIMENSION FMT(48), NTAB(101,21), NRTOT(101), NCLT(21),LOB(16)
DIMENSIONTAB(101,21),RTOT(101), CLT(21), PCENT(21)
DIMENSION NDEEP(72), NDIM(4,72), MXTEMP(4)
DIMENSION STORE(40,12)
DIMENSION CKK(16)
C
C
EQUIVALENCE (NTAB,TAB)
NOUTD=-1
IND=-4
INP=2
IOUT=3
C---------------TTYPTY RETURNS 0 - TTY JOB, MINUS ONE - BATCH JOB
CALL TTYPTY(ICODE)
WRITE(NOUTD,8745)
8745 FORMAT(20X,'PDP-10 VERSION OF NUCROS')
7 FORMAT(40I2)
C CALL USAGE('CROSS')
C
C READ AND PRINTTITLES
C---------------NDV1, NDV2 ARE RETURNED. OTHER ARGS. ARE INPUT.
C--------------- 1 MEANS OUTPUT? PRINTS, 0 - INPUT? PRINTS
CALL IO(IOUT,NDV1,NOUTD,IND,1,ICODE)
500 CALL IO(INP,NDV2,NOUTD,IND,0,ICODE)
C---------------FMT, ISTD ARE RETURNED. OTHER ARGS. ARE INPUT.
C--------------- 48=NO. OF OBJ. TIME FORMAT WORDS (3 LINES).
C--------------- 1 MEANS I TYPE FORMAT ONLY.
CALL GETFOR(NOUTD,IND,FMT,ISTD,48,1)
REWIND 20
REWIND 21
REWIND 22
800 FORMAT(1H1)
WRITE(NOUTD,1002)
1002 FORMAT(' ENTER IDENTIFICATION'/)
READ(IND,10)(CKK(I),I=1,16)
WRITE(IOUT,5)(CKK(I),I=1,16)
5 FORMAT(1X,16A5)
C
C READ CONTROL CARDS 4-8
230 WRITE(NOUTD,1003)
1003 FORMAT(' ENTER PROBLEM CONTROL PARAMETERS'/)
READ(IND,1001)NOBS,NVAR,NTABLE,IDVAR,NCODE,NCHI,NPCT,
1NTAU
1001 FORMAT(10I)
IF(NVAR.LT.2.OR.NVAR.GT.40)GO TO 217
IF(NTABLE .LT.1.OR.NTABLE.GT.72)GO TO 218
GO TO 219
217 WRITE(NOUTD,220)
220 FORMAT(' ILLEGAL NUMBER OF VARIABLES SPECIFIED.'/)
IF(ICODE.EQ.-1)CALL EXIT
GO TO 230
218 WRITE(NOUTD,221)
221 FORMAT(' ILLEGAL NUMBER OF TABLES SPECIFIED'/)
IF(ICODE.EQ.-1)CALL EXIT
GO TO 230
219 LCHI=0
IF(NCHI.GT.0.OR.NTAU.GT.0)LCHI=1
10 FORMAT(16A5)
IF(NDV1.NE.'TTY')WRITE(IOUT,690)(FMT(I),I=1,48)
690 FORMAT(' FORMAT OF DATA'/(' ',16A5))
LPCT = MIN0(NPCT,1)
WRITE (IOUT,185) NOBS
WRITE(NOUTD,1006)
1006 FORMAT(//' ENTER MAXIMUM VALUES'/)
READ (IND,1012) (MAX(J),J=1,NVAR)
WRITE (IOUT,9)(J,J=1,NVAR)
WRITE(IOUT,423) (MAX(J),J=1,NVAR)
6 FORMAT(I4,9I2)
185 FORMAT(/20H TOTAL CASES USED = I4)
780 FORMAT (1H0)
423 FORMAT(11H MAX.VALUE ,40I3)
9 FORMAT(//11H VARIABLE , 40I3)
C
C READ + PRINT CONTOL CDS 9,10,11,12.. DETERMINE ALL TABLE DEPTHS
WRITE(NOUTD,1007)
1007 FORMAT(//' ENTER TABLE CONTROL PARAMETERS'/)
1012 FORMAT(30I2)
DO 1000 I=1,4
1000 READ (IND,1012) (NDIM(I,K),K=1,NTABLE)
DO472 K=1,NTABLE
NDEEP(K) =2
DO 472 I=3,4
IF (NDIM(I,K))472,472,470
470 NDEEP(K)=I
472 CONTINUE
WRITE (IOUT,424)
IF(NTABLE-36)9001,9001,9002
9001 WRITE(IOUT,781)
GO TO 9003
9002 WRITE (IOUT,426)
9003 CONTINUE
LTABLE = MIN0(NTABLE,36)
DO 4 L=1,LTABLE
K=L
NDEPTH = NDEEP(K)
WRITE (IOUT,8) K,(NDIM(I,K),I=1,NDEPTH)
DO 550 I = 1,NDEPTH
J = NDIM(I,K)
550 MXTEMP(I) = MAX(J)
WRITE(IOUT,425) (MXTEMP(I),I=1,NDEPTH)
K=L+36
IF (K .GT. NTABLE) GO TO 427
NDEPTH = NDEEP(K)
WRITE (IOUT,428) K,(NDIM(I,K),I=1,NDEPTH)
DO 429I=1,NDEPTH
J=NDIM(I,K)
429 MXTEMP(I) = MAX(J)
WRITE (IOUT,430) (MXTEMP(I),I=1,NDEPTH)
GO TO 4
427 CONTINUE
4 CONTINUE
8 FORMAT (1X,9XI6,7X,4I4)
424 FORMAT (//1H+,9X9HTABLE NO.,6X9HVARIABLES,7X14HMAXIMUM VALUES)
425 FORMAT (1H+,38X,4I4)
426 FORMAT ('+'75X9HTABLE NO.,6X9HVARIABLES,7X14HMAXIMUM VALUES)
428 FORMAT (1H+,74X,I6,7X,4I4)
430 FORMAT ('+'104X,4I4)
781 FORMAT (' ')
C
C READ AND PRINT RECODE-CONTROL-CARDS, IF RECODING OPTED
IF (NCODE) 20,20,19
19 WRITE(NOUTD,1014)
1014 FORMAT(' ENTER TWO RECODING LINES'/)
READ (IND,1012) (KIND(J),J=1,NVAR)
READ (IND,1012) (KONST(J),J=1,NVAR)
DO 400 J=1,NVAR
NOPER = KIND(J)
IF (NOPER-3) 400,400,403
403 WRITE(NOUTD,1015)
1015 FORMAT(' ENTER GROUP BOUNDARIES'/)
READ(IND,1019) MBRAK
1019 FORMAT(16I)
IF(NOPER.EQ.8)GO TO 650
DO 651 I=0,15
651 LOB(I+1)=I
GO TO 652
650 WRITE(NOUTD,653)
653 FORMAT(' ENTER CATEGORIZATION'/)
READ(IND,1019)LOB
652 WRITE (21) MBRAK,LOB
400 CONTINUE
WRITE(IOUT,186)
REWIND 21
DO 2200 J=1,NVAR
NOPER=KIND(J)
IF(NOPER)2200,2200,2202
2202 IF(NOPER-3)2201,2201,2203
2201 WRITE(IOUT,188)J,NOPER,KONST(J)
GO TO 2200
2203 READ(21)MBRAK,LOB
IF(NOPER.EQ.4)WRITE(IOUT,187)J,NOPER,MBRAK
IF(NOPER.GT.4.AND.NOPER.LT.8)WRITE(IOUT,9004)J,NOPER,KONST(J),MBRA
1K
IF(NOPER.EQ.8)WRITE(IOUT,660)J,NOPER,MBRAK,LOB
2200 CONTINUE
186 FORMAT (27H1RECODING OPTION CALLED FOR // 4X,67H VAR.NO. CODE
X CONSTANT RECODED CATEGORIES (UPPER LIMITS) )
188 FORMAT (8XI2,5XI1,4XI3,10X4HNONE)
187 FORMAT ( 8XI2,5XI1,4X4HNONE,16I3)
9004 FORMAT(8XI2,5XI1,4XI3,1X,16I3)
660 FORMAT(8X,I2,5X,I1,4X,'NONE',16I3/24X,16I3)
REWIND 21
C
C READ AND PRINT VAR NAMES IF OPTED
20 IF (IDVAR) 480,480,495
495 WRITE(NOUTD,1008)
1008 FORMAT(' ENTER VARIABLE NAMES'/)
DO 490 I=1,NVAR
490 READ(IND,82)(STORE(I,J),J=1,12)
WRITE(IOUT,83)
DO 1009 I=1,NVAR
1009 WRITE(IOUT,57)I,(STORE(I,J),J=1,12)
83 FORMAT (8H1VAR.NO.//)
82 FORMAT(12A5)
57 FORMAT(I5,2X,12A5)
C
C INITIALIZE FOR FIRST BATCH
480 NBATCH = MIN0 (NOBS,200)
NOBS = NOBS- NBATCH
IBATCH = 1
IF(NDV2.EQ.'TTY')GO TO 214
WRITE(NOUTD,3576)
GO TO 215
214 IF(ISTD.EQ.1)WRITE(NOUTD,216)
216 FORMAT(/' ENTER DATA (AT MOST 10 PER LINE)'/)
IF(ISTD.EQ.0)WRITE(NOUTD,3577)
3576 FORMAT(/' YOUR DATA IS BEING READ. PLEASE WAIT.'/)
3577 FORMAT(/' ENTER DATA'/)
215 IOSW = 24
LINE = 60
C
C
C READ A BATCH OF DATA CARDS
299 DO 3817 N=1,NBATCH
3817 READ(INP,FMT)(LAD(N,J),J=1,NVAR)
C
C RECODE VARIABLES IF OPTED
302 IF (NCODE) 485,485,483
483 DO 32 J = 1,NVAR
NOPER = KIND(J)
IF(NOPER) 32,32,176
176 GO TO (23,24,25,52,52,52,52,52 ), NOPER
C
C DO TYPE 1,2,OR 3 RECODING
23 DO 1201 I=1,NBATCH
1201 LAD(I,J) = LAD(I,J)+KONST (J)
GO TO 32
24 DO1203 I = 1,NBATCH
1203 LAD(I,J) = MAX0(LAD(I,J)-KONST(J),0)
GO TO 32
25 DO1202 I= 1,NBATCH
1202 LAD(I,J)= LAD(I,J) /KONST(J)
GO TO 32
C
C DO TYPE 4 REC DING
52 READ (21) MBRAK,LOB
DO 29 I = 1,NBATCH
DO 28 K = 1,16
IF(LAD(I,J)-MBRAK(K)) 27,27,28
28 CONTINUE
LAD(I,J) = 20
GO TO 29
27 LAD(I,J) =LOB(K)
29 CONTINUE
IF(NOPER.EQ.8)GO TO 32
IF(NOPER-4) 32,32,7638
7638 NOPER= NOPER-4
GO TO 176
32 CONTINUE
REWIND 21
C
C DEFINE TABLE DEPTH,MAX VAL, VARS TO BE USED...FOR EACH TABLE
485 DO 115 KT=1,NTABLE
NDEPTH = NDEEP(KT)
J1 = NDIM(1,KT)
J2 = NDIM(2,KT)
MAXA = MAX(J1) +1
MAXB = MAX(J2) +1
MAXC = 0
MAXD = 0
LREQ = 10 + MAXB + 2*LCHI
LREQPC = LPCT * (MAXB-1+6)
GO TO (36,36,33,34),NDEPTH
34 J4 = NDIM(4,KT)
MAXD = MAX(J4)
33 J3 = NDIM(3,KT)
MAXC = MAX(J3)
LREQ = LREQ +2
C
36 DO 115 L3 = 0,MAXC
DO 115 L4 = 0,MAXD
IF (IBATCH-1) 41,41,42
C ZERO OUT TABLE AREA IF FIRST BATCH
41 DO 43 I=1,MAXB
DO 43 J=1,MAXA
43 NTAB(I,J) = 0
GO TO 44
C OTHERWISE READ INTO TABLE-AREA STORED-RESULTS FROM PREVIOUS BATCH(ES)
42 IF (IOSW-24) 46,46,47
46 READ (20) ((NTAB(I,J),J=1,MAXA),I=1,MAXB)
GO TO 44
47 READ (22) ((NTAB(I,J),J=1,MAXA),I=1,MAXB)
C
C DO A TABLE. TEST 3RD AND 4TH-DIMENSION VARIABLES FOR APPROP VALUE
C BEFORE COUNTING THE OBSERVATION (WHERE 3RD OR 4TH DIM IS USED)
44 DO 45 N=1,NBATCH
IF (NDEPTH-3) 40,35,39
39 IF (LAD(N,J4) -L4) 45,35,45
35 IF (LAD(N,J3) - L3) 45,40,45
40 L1 = LAD(N,J1) +1
L2 = LAD(N,J2) +1
IF (L1 .LT. 1 .OR. L2 .LT. 1) GO TO 45
IF (L1.GT.MAXA .OR. L2.GT.MAXB) GO TO 45
NTAB(L2,L1) = NTAB(L2,L1) +1
45 CONTINUE
C
C IF MORE OBSERVATIONS REMAIN, STORE THIS TABLE (RESULTS TO DATE)
IF (NOBS) 75,75,48
48 IF (IOSW-24) 50,50,60
50 WRITE (22) ((NTAB(I,J),J=1,MAXA),I=1,MAXB)
GO TO 115
60 WRITE (20) ((NTAB(I,J),J=1,MAXA),I=1,MAXB)
GO TO 115
C
C
C IF NO MORE OBSERVATIONS WRITE OUT THIS TABLE
75 LINE= LINE + LREQ
IF(LINE+2*LREQPC-58)431,431,811
431 WRITE(IOUT,780)
GO TO 76
811 WRITE (IOUT,800)
LINE=LREQ-1
76 GO TO (822,822,823,824),NDEPTH
C TWO-DIMENSION TITLE
822 IF(IDVAR)522,522,523
523 WRITE(IOUT,524)KT,J1,(STORE(J1,J),J=1,12),J2,(STORE(J2,J),J=1,12)
1 ,MAXA,MAXB
524 FORMAT(7X,11H TABLE NO. ,I3,10X,13HVARIABLE NO. ,I2,2X, 12A5,/
1 25X,6HVRS. ,13HVARIABLE NO. ,I2,2X,12A5, / 7X,13HTABLE SIZE =
2 I3,5H BY ,I3)
GO TO 86
522 WRITE (IOUT,80) KT,J1,J2,MAXA,MAXB
80 FORMAT ( 10X,11H TABLE NO. I3,5X,13HVARIABLE NO. I2,3X,
1 9HVRS. NO. , I2//13X,13HTABLE SIZE = I3,5H BY I3)
GO TO 86
C THREE DIMENSIOANL TITLE
823 IF (IDVAR)525, 525,526
526 WRITE(IOUT,527) KT,J1,(STORE (J1,J),J=1,12),J2,(STORE(J2,J),J=1,12
1),J3,(STORE(J3,J),J=1,12), MAXA,MAXB,L3
527 FORMAT(7X,11H TABLE NO. , I3,10X,13HVARIABLE NO. I2,2X,12A5,
1 2(/25X,6HVRS ,13HVARIABLE NO. I2,2X,12A5) //7X,
213HTABLE SIZE = I3,5H BY I3,10X,16H 3RD DIMENSION = I3)
GO TO 86
525 WRITE (IOUT,77) KT,J1,J2,J3, MAXA,MAXB,L3
77 FORMAT(10X,11HTABLE NO. I3,5X,13HVARIABLE NO I2,3X,9HVRS NO.
1I2,3X,9HVRS. NO. I2,///13X,13HTABLE SIZE = I3,5H BY I3//13X,
216H3RD DIMENSION = I3)
GO TO 86
C FOUR DIMENSION TITLE
824 IF(IDVAR)528,528,529
529 WRITE(IOUT,530)KT,J1,(STORE(J1,J),J=1,12),J2,(STORE(J2,J),J=1,12),
1J3,(STORE(J3,J),J=1,12),J4,(STORE(J4,J),J=1,12),MAXA,MAXB,L3,L4
530 FORMAT(7X,11H TABLE NO. I3,10X,13HVARIABLE NO. I2,2X, 12A5,
13(/25X,6HVRS. ,13HVARIABLE NO. I2,2X, 12A5) //7X,
213HTABLE SIZE = I3,5H BY I3,10X,16H3RD DIMENSION = I3,5X,
316H4TH DIMENSION = I3)
GO TO 86
528 WRITE(IOUT,72) KT,J1,J2,J3,J4,MAXA,MAXB,L3,L4
72 FORMAT (10X,11H TABLE NO. I3,5X,13HVARIABLE NO. I2,3X,9HVRS.
1NO. I2,3X,9HVRS NO. I2,3X,9HVRS NO. I2//13X,13HTABLE SIZE =
2I3,5H BY I3//13X,16H3RD DIMENSION = I3,5X,16H4TH DIMENSION = I3)
86 CONTINUE
C
C
C PRINT BODY OF MAP,COMPUTE AND PRINT COLUMN AND ROW TOTALS
84 MSUB = MAXA-1
WRITE (IOUT,91) (IH,IH=0,MSUB)
WRITE (IOUT,781)
NGRTO=0
DO 85 J=1,MAXA
85 NCLT(J)=0
DO 100 I=1,MAXB
ISUB = I-1
NRTOT(I)=0
DO 90 J=1,MAXA
NRTOT(I)= NRTOT(I) + NTAB(I,J)
90 NCLT(J)= NCLT(J)+NTAB (I,J)
WRITE(IOUT,105) ISUB,NRTOT(I),(NTAB(I,J),J=1,MAXA)
100 NGRTO = NGRTO+NRTOT(I)
WRITE (IOUT,106) NGRTO,(NCLT(J),J=1,MAXA)
91 FORMAT (/17X,3HTOT,3X,21I5)
105 FORMAT(9X,I3,I8,3X,21I5)
106 FORMAT (/9X5HTOTAL,I6,3X,21I5)
C
C
C DO CHISQARES +/OR PERCENTS +/OR COEFFICIENTS IF OPTED
C GET FLOATING _PT EQUIVALENTS..DEDUCT ZEROROW,COLUMN FROM TOTALS
IF(NCHI +NPCT + NTAU) 115,115,404
404 DO 407 I=1,MAXB
RTOT(I)=NRTOT(I)
DO 407 J=1,MAXA
407 TAB(I,J) = NTAB (I,J)
DO 405 I=1,MAXA
405 CLT(I)=NCLT(I)
GRTO=NGRTO
C
IF(NCHI+NTAU) 410,410,605
605 IF (NCHI) 606,606,408
606 WRITE (IOUT,781)
GO TO 304
408 CHISQ = 0.0
IF(NGRTO.NE.0)GO TO 1071
CHISQ=0.
CRXC=0.
GO TO 1070
1071 DO 409 I=1,MAXB
DO409 J=1,MAXA
A= (CLT(J)*RTOT(I)) /GRTO
B= (TAB(I,J)-A )**2
IF(A.NE.0.)GO TO 409
A=1.
B=0.
409 CHISQ = CHISQ + B/A
CRXC=( CHISQ/(CHISQ+GRTO)) **.5
1070 WRITE(IOUT,608) CHISQ,CRXC
608 FORMAT(/1H+,6X,12HCHI SQUARE = ,F8.3,3X,4HC = ,F5.3)
IF(MAXB.NE.2.OR.MAXA.NE.2)GO TO 631
A=TAB(1,1)
B=TAB(1,2)
C=TAB(2,1)
D=TAB(2,2)
C22=GRTO*(ABS(A*D-B*C)-GRTO/2)**2/((A+B)*(C+D)*(A+C)*(B+D))
WRITE(IOUT,632)C22
632 FORMAT(4X,'2X2 CORRECTED CHI-SQUARE = ',F8.3)
631 IF(NTAU) 604,604,304
604 WRITE(IOUT,781)
GO TO 410
C
C
C DO ALL CORRELATION COEFFICIENTS
C NTAU IS NOT USED BY THE SUBROUTINE BUT MAY BE USED TO
C EXERCISE OPTIONS IN ALTERNATE SUBROUTINES
304 CALL COEFT (TAB,MAXB,MAXA,NTAU)
C
C
C DO COLUMN PERCENTAGES
410 IF (NPCT) 115,115,411
411 LINE = LINE +LREQPC
IF (LINE-58) 412,412,413
412 WRITE (IOUT,108)
GO TO 414
413 WRITE (IOUT,800)
WRITE (IOUT,108)
WRITE (IOUT,113) (IH,IH=0,MSUB)
WRITE (IOUT,781)
LINE = LREQPC + 3
108 FORMAT(//7X,40HPERCENTS BY COLUMN FROM THE ABOVE MATRIX //)
113 FORMAT (17X,3HTOT,8X,20I5 )
C
414 DO 416 I = 1,MAXB
ISUB = I-1
DO 415 J = 1,MAXA
IF(CLT(J).NE.0.)GO TO 1080
PCENT(J)=0.
GO TO 415
1080 PCENT(J) = TAB(I,J) / CLT(J)
415 CONTINUE
416 WRITE(IOUT,109) ISUB,NRTOT(I), (PCENT(J),J = 1,MAXA)
WRITE (IOUT,111) NGRTO,(NCLT(J), J=1,MAXA )
109 FORMAT (9X,I3,I8,9X,2P20F5.1)
111 FORMAT(/9X,5HTOTAL, I6,8X,20I5 )
C
C DO ROW PERCENTS
LINE = LINE+ LREQPC
IF(LINE-58) 417,417,418
417 WRITE( IOUT,114)
GO TO 419
418 WRITE( IOUT,800)
WRITE (IOUT,114)
WRITE (IOUT,113) (IH,IH=0,MSUB)
WRITE (IOUT,781)
LINE = LREQPC + 3
114 FORMAT(//7X,38HPERCENTS BY ROW FROM THE ABOVE MATRIX //)
C
419 DO 421 I= 1,MAXB
ISUB = I-1
DO 420 J= 1,MAXA
IF(RTOT(I).NE.0.)GO TO 1081
PCENT(J)=0.
GO TO 420
1081 PCENT(J) = TAB (I,J) / RTOT(I)
420 CONTINUE
421 WRITE (IOUT,109) ISUB,NRTOT(I),(PCENT(J),J=1,MAXA)
WRITE (IOUT,111) NGRTO, (NCLT(J),J=1,MAXA)
115 CONTINUE
C
C
C PREPARE PROCESS NEXT BATCH IF ANY
IF (NOBS) 500,500,130
130 NBATCH = MIN0 (NOBS,200)
NOBS = NOBS-NBATCH
IBATCH = IBATCH + 1
IF(IOSW-24) 12,12,13
12 IOSW = 42
GO TO 15
13 IOSW = 24
15 REWIND 20
REWIND 22
GO TO 299
END
C---------------NTAU APPARENTLY NOT USED. OTHER ARGS. ARE INPUT.
SUBROUTINE COEFT(A,II,JJ,NTAU)
DIMENSION A(101,21), TI(101), TJ(21)
IOUT=3
SUM = 0.
DO 2 I = 1,II
DO 2 J = 1,JJ
2 SUM=SUM+A(I,J)
P=0.
IX = II-1
JX = JJ-1
DO 10 I=1,IX
DO 10 J=1,JX
IY=I+1
JY=J+1
DO 10 K=IY,II
DO 10 L=JY,JJ
10 P=P+A(I,J)*A(K,L)
Q=0.
DO 11 I=1,IX
DO 11J=2,JJ
IN=I+1
JY=J-1
DO 11 K=IN,II
DO 11 L=1,JY
11 Q = Q+A(I,J)*A(K,L)
S=P-Q
GAMMA=(P-Q)/(P+Q)
100 DO 21 I=1,II
TI(I)=0.
DO21 J=1,JJ
21 TI(I)=TI(I)+A( I,J)
DO 22 J=1,JJ
TJ(J)=0.
DO22 I=1,II
22 TJ(J)=TJ(J)+A(I,J)
I1=0
DO 105 I=1,II
105 IF(TI(I).NE.0.0)I1=I1+1
J1=0
DO 106 J=1,JJ
106 IF(TJ(J).NE.0.0)J1=J1+1
82 IF (I1-J1) 86,85,86
85 IPT = 2
GO TO 84
86 IPT = 3
84 CONTINUE
XU=0.
DO55 J=1,JJ
55 XU=XU + TJ(J)**2
XU=.5*(SUM**2-XU)
YU=0.
DO56 I=1,II
56 YU=YU+TI(I)**2
YU=.5*(SUM**2-YU)
DYX = S/XU
DXY=S/YU
IF(IPT-2) 57,57,102
57 CONTINUE
T=0.
DO 23 I=1,II
23 T=T+.5*TI(I)*(TI(I)-1.)
U=0.
DO24 J=1,JJ
24 U=U+.5*TJ(J)*(TJ(J)-1.)
BI=(SUM*(SUM-1.) *.5-T)**(.5)
BJ=(SUM*(SUM-1.) *.5-U)**(.5)
TAUB=S/(BI*BJ)
104 WRITE (IOUT,301) TAUB,GAMMA,DXY,DYX
RETURN
102 IF(I1-J1) 5,5,6
5 AM=I1
GO TO 12
6 AM=J1
12 TAUC=S/((.5*SUM**2*(AM-1.)/AM))
WRITE(IOUT,302) TAUC,GAMMA,DXY,DYX
RETURN
301 FORMAT(1H+,41X,8HTAU-B = ,F6.3,3X,8HGAMMA = ,F6.3,3X,6HDXY = ,
1F6.3,3X, 6HDYX = ,F6.3)
302 FORMAT(1H+,41X,8HTAU-C = ,F6.3,3X,8HGAMMA = ,F6.3,3X,6HDXY = ,
1F6.3,3X, 6HDYX = ,F6.3)
END