Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/bmd/bmd05m.for
There is 1 other file named bmd05m.for in the archive. Click here to see a list.
C PROGRAM WAS CONVERTED FROM FORTRAN 2 TO 7090 FORTRAN 4
C IT WAS THEN CONVERTED TO 360 FORTRAN IV (H-LEVEL)
C DISCRIMINANT ANALYSIS - SEVERAL GROUPS JUNE 9, 1966
C THIS IS A SIFTED VERSION OF BMD05M ORIGINALLY WRITTEN IN
C FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
DIMENSION X(4,15,150),XMEAN(5,25),COV(20,20),SUMCOV(20,20),
1 DENT(20,20),LB(25),M(25),N(5),C(20,20),NDENT(5,6)
COMMON K , NN , N , LA , KP , X
COMMON ONN , XMEAN , NTO , M , K1 , K2
946 FORMAT(60H1BMD05M - DISCRIMINANT ANALYSIS-SEVERAL GROUPS - VERSION
1 OF ,18HJUNE 9, 1966 ,/
241H HEALTH SCIENCES COMPUTING FACILITY, UCLA)
DOUBLE PRECISION A1,A2,XZ
DATA A2,XZ/6HPROBLM,6HFINISH/
C
NTAPE=5
CALL USAGEB('BMD05M')
100 READ (5,900)A1,PROB,K,KP,NVG,NADD,MTAPE,(N(I),I=1,5), KVR
IF(A1 .EQ. A2) GO TO 15
5 IF(A1 .EQ. XZ) GO TO 80
10 WRITE (6,948) A1
GO TO 80
15 CALL TPWD(MTAPE,NTAPE)
21 JESSE=KP+NADD
IF(KVR.GT.0.AND.KVR.LE.10)GO TO 210
KVR=1
WRITE (6,4000)
210 WRITE (6,946)
WRITE (6,929)PROB,JESSE
IF((K-1)*(K-5))20,949,949
20 IF((KP-K+1)*(KP-16))22,951,951
22 DO 23 I=1,K
IF(N(I)-150)23,23,953
23 CONTINUE
IF ((JESSE-K+1)*(JESSE-16))25,955,955
25 CALL READ(KVR,NVG,NTAPE,NADD)
KP=JESSE
IF(-NVG)81,81,100
81 DO 24 I=1,KP
DO 24 J=1,KP
24 SUMCOV(I,J)=0.0
DO 30 I=1,K
NN=N(I)
DO 26 L=1,KP
DO 26 LL=1,KP
COV(L,LL)=0.0
DO 26 J=1,NN
26 COV(L,LL)=COV(L,LL)+((X(I,L,J)-XMEAN(I,L))*(X(I,LL,J)-XMEAN(I,LL))
1)
DO 28 L=1,KP
DO 28 LL=1,KP
28 SUMCOV(L,LL)=SUMCOV(L,LL)+COV(L,LL)
WRITE (6,930)
WRITE (6,905)I
DO 30 J=1,KP
WRITE (6,906)J
30 WRITE (6,907)(COV(J,L),L=1,KP)
DF=0.0
DO 37 I=1,K
DFF=N(I)
37 DF=DF+DFF
OK=K
DF=DF-OK
DO 38 I=1,KP
DO 38 J=1,KP
SUMCOV(I,J)=SUMCOV(I,J)/DF
38 COV(I,J)=SUMCOV(I,J)
WRITE (6,908)
DO 39 I=1,KP
WRITE (6,906)I
39 WRITE (6,907)(SUMCOV(I,J),J=1,KP)
CALL INVERT(COV,KP,D,LB,M)
WRITE (6,930)
WRITE (6,909)
DO 40 I=1,KP
WRITE (6,906)I
40 WRITE (6,907)(COV(I,J),J=1,KP)
DO 41 I=1,KP
DO 41 J=1,KP
C(I,J)=0.0
DO 41 II=1,KP
41 C(I,J)=C(I,J)+SUMCOV(I,II)*COV(II,J)
WRITE (6,930)
WRITE (6,910)
DO 42 I=1,KP
WRITE (6,906)I
42 WRITE (6,907)(C(I,J),J=1,KP)
CALL CHITES(N,COV,XMEAN,KP,K)
DO 43 I=1,K
DO 43 J=1,KP
SUMCOV(I,J)=0.0
DO 43 II=1,KP
43 SUMCOV(I,J)=SUMCOV(I,J)+XMEAN(I,II)*COV(II,J)
DO 44 I=1,K
DO 44 J=1,KP
44 DENT(J,I)=XMEAN(I,J)*(-0.5)
DO 47 I=1,K
DO 47 J=1,K
C(I,J)=0.0
DO 47 II=1,KP
47 C(I,J)=C(I,J)+SUMCOV(I,II)*DENT(II,J)
WRITE (6,930)
DO 45 I=1,K
45 M(I)=I
WRITE (6,911)(M(I),I=1,K)
WRITE (6,912)
DO 50 I=1,KP
50 WRITE (6,913)I,(SUMCOV(J,I),J=1,K)
WRITE (6,930)
DO 48 I=1,K
J=I
48 C(15,I)=C(I,J)
WRITE (6,914) (C(15,I),I=1,K)
DO 52 II=1,K
NN=N(II)
DO 52 I=1,NN
DO 49 J=1,K
DENT(1,J)=0.0
DO 49 I1=1,KP
49 DENT(1,J)=DENT(1,J)+X(II,I1,I)*SUMCOV(J,I1)
DO 52 J=1,K
52 X(II,J,I)=DENT(1,J)+C(15,J)
SMALL=-(10.0**36.0)
K1=K+1
K2=K+2
DO 55 I=1,K
NN=N(I)
DO 55 J=1,NN
X(I,K1,J)=SMALL
DO 55 JJ=1,K
IF(X(I,K1,J)-X(I,JJ,J))54,55,55
54 X(I,K1,J)=X(I,JJ,J)
X(I,K2,J)=JJ
55 CONTINUE
DO 57 I=1,K
NN=N(I)
DO 57 J=1,NN
SUM=0.0
DO 56 JJ=1,K
X(I,JJ,J)=EXP(X(I,JJ,J)-X(I,K1,J))
56 SUM=SUM+X(I,JJ,J)
X(I,K1,J)=1.0
DO 57 JJ=1,K1
57 X(I,JJ,J)=X(I,JJ,J)/SUM
WRITE (6,930)
CALL CLASSI
DO 70 I=1,K
NN=N(I)
DO 70 II=1,K
NDENT(I,II)=0
DO 70 J=1,NN
NSUM=0
DO 68 JJ=1,K
IF(X(I,II,J)-X(I,JJ,J))68,68,67
67 NSUM=NSUM+1
68 CONTINUE
IF(NSUM-(K-1)) 70, 69, 70
69 NDENT(I,II)=NDENT(I,II)+1
70 CONTINUE
WRITE (6,930)
WRITE (6,924)
GO TO (91,92,93,94,95),K
91 WRITE (6,941)(M(I),I=1,K)
GO TO 96
92 WRITE (6,942)(M(I),I=1,K)
GO TO 96
93 WRITE (6,943)(M(I),I=1,K)
GO TO 96
94 WRITE (6,944)(M(I),I=1,K)
GO TO 96
95 WRITE (6,945)(M(I),I=1,K)
96 WRITE (6,926)
DO 75 I=1,K
NDENT(I,K1)=0
DO 72 J=1,K
72 NDENT(I,K1)=NDENT(I,K1)+NDENT(I,J)
WRITE (6,927)I,(NDENT(I,J), J=1,K1)
75 WRITE (6,931)
GO TO 100
900 FORMAT (A6 ,A2,5I2,5I3,37X,I2)
905 FORMAT(52H0MATRIX OF CROSS PRODUCTS OF DEV. FROM MEANS - GROUPI2)
906 FORMAT(4H ROWI3)
907 FORMAT(1H F14.5,7F15.5)
908 FORMAT(18H0DISPERSION MATRIX)
909 FORMAT(29H0INVERSE OF DISPERSION MATRIX)
910 FORMAT(46H0TEST OF ACCURACY OF DISPERSION MATRIX INVERSE)
911 FORMAT(9H0FUNCTIONI13,4I16)
912 FORMAT(12H0COEFFICIENT)
913 FORMAT(5X,I2,F20.5,4F16.5)
914 FORMAT(9H0CONSTANTF18.5,4F16.5)
924 FORMAT(22H0CLASSIFICATION MATRIX)
926 FORMAT(6H0GROUP)
927 FORMAT(1H I3,8X,6I11)
929 FORMAT(14H0PROBLEM NO. A2/20H NUMBER OF VARIABLESI4///)
930 FORMAT(1H0)
931 FORMAT(1H )
941 FORMAT(9H0FUNCTIONI14,12H TOTAL)
942 FORMAT(9H0FUNCTIONI14, I11,12H TOTAL)
943 FORMAT(9H0FUNCTIONI14,2I11,12H TOTAL)
944 FORMAT(9H0FUNCTIONI14,3I11,12H TOTAL)
945 FORMAT(9H0FUNCTIONI14,4I11,12H TOTAL)
948 FORMAT (' PROBLM OR FINISH CARD EXPECTED. CARD READ WAS ',A6)
949 WRITE (6,950) K
950 FORMAT (1X,I3,' NUMBER OF GROUPS IS ILLEGAL')
GO TO 80
951 WRITE (6,952) KP
952 FORMAT (1X,I3,' NUMBER OF ORIGINAL VARIABLES IS ILLEGAL')
GO TO 80
953 WRITE (6,954) N(I),I
954 FORMAT (' SAMPLE SIZE OF ',I4,' FOR GROUP ',I2,' IS ILLEGAL')
GO TO 80
955 WRITE (6,956) JESSE
956 FORMAT (1X,I5,' NUMBER OF VARIABLES AFTER TRANSGENERATION IS ILLEG
1AL')
GO TO 80
4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
1IED, ASSUMED TO BE 1.)
80 IF(NTAPE-5)83,83,82
82 REWIND NTAPE
83 STOP
END
C SUBROUTINE CHITES FOR BMD05M JUNE 9, 1966
SUBROUTINE CHITES(N,COV,XMEAN,KP,K)
DIMENSION XBAR(25),XMEAN(5,25),COV(20,20),N(5)
DFSUM=0.0
CHI=0.0
DO 2 I=1,KP
2 XBAR(I)=0.0
DO 3 I=1,K
XN=N(I)
3 DFSUM=DFSUM+XN
DO 4 I=1,KP
DO 5 J=1,K
TEMP=N(J)
5 XBAR(I)=XBAR(I)+TEMP*XMEAN(J,I)
4 XBAR(I)=XBAR(I)/DFSUM
DO 6 I=1,KP
DO 6 J=1,KP
TEMP=0.0
DO 7 L=1,K
DNN=N(L)
7 TEMP=TEMP+DNN*(XMEAN(L,I)-XBAR(I))*(XMEAN(L,J)-XBAR(J))
6 CHI=CHI+TEMP*COV(I,J)
NDF=KP*(K-1)
WRITE (6,9)CHI
WRITE (6,8)CHI,NDF,K,KP
8 FORMAT(11H0THE VALUE F10.5,31H CAN BE USED AS CHI-SQUARE WITHI3,/5
16H DEGREES OF FREEDOM TO TEST THE HYPOTHESIS THAT THE MEAN/32H VAL
2UES ARE THE SAME IN ALL THE I3,18H GROUPS FOR THESE I2,/11H VARIAB
3LES.//)
9 FORMAT(1H0/33H GENERALIZED MAHALANOBIS D-SQUAREF12.5)
RETURN
END
C SUBROUTINE CLASSI FOR BMD05M JUNE 9, 1966
SUBROUTINE CLASSI
DIMENSION N(5),X(4,15,150),XMEAN(5,25),M(25)
COMMON K , NN , N , LA , KP , X
COMMON ONN , XMEAN , NTO , M , K1 , K2
WRITE (6,915)
GO TO (81,82,83,84,85),K
81 WRITE (6,931)(M(II),II=1,K)
GO TO 86
82 WRITE (6,932)(M(II),II=1,K)
GO TO 86
83 WRITE (6,933)(M(II),II=1,K)
GO TO 86
84 WRITE (6,934)(M(II),II=1,K)
GO TO 86
85 WRITE (6,935)(M(II),II=1,K)
86 DO 66 I=1,K
WRITE (6,917)I
WRITE (6,918)
NN=N(I)
DO 66 J=1,NN
NX=X(I,K2,J)
GO TO (61,62,63,64,65),K
61 WRITE (6,919)J,(X(I,JJ,J),JJ=1,K1),NX
GO TO 66
62 WRITE (6,920)J,(X(I,JJ,J),JJ=1,K1),NX
GO TO 66
63 WRITE (6,921)J,(X(I,JJ,J),JJ=1,K1),NX
GO TO 66
64 WRITE (6,922)J,(X(I,JJ,J),JJ=1,K1),NX
GO TO 66
65 WRITE (6,923)J,(X(I,JJ,J),JJ=1,K1),NX
66 CONTINUE
915 FORMAT(53H0EVALUATION OF CLASSIFICATION FUNCTIONS FOR EACH CASE)
917 FORMAT(6H GROUPI3)
918 FORMAT(6H CASE)
919 FORMAT(I5,F21.5,F15.5,I12)
920 FORMAT(I5,F21.5,2F15.5,I12)
921 FORMAT(I5,F21.5,3F15.5,I12)
922 FORMAT(I5,F21.5,4F15.5,I12)
923 FORMAT(I5,F21.5,5F15.5,I12)
931 FORMAT(9H0FUNCTIONI12, 41H LARGEST FN.NO.FOR LARGES
1T/32X,27HPROBABILITY PROBABILITY)
932 FORMAT(9H0FUNCTIONI12, I15,41H LARGEST FN.NO.FOR LA
1RGEST/47X,27HPROBABILITY PROBABILITY)
933 FORMAT(9H0FUNCTIONI12,2I15,41H LARGEST FN.NO.FOR LA
1RGEST/62X,27HPROBABILITY PROBABILITY)
934 FORMAT(9H0FUNCTIONI12,3I15,41H LARGEST FN.NO.FOR LA
1RGEST/77X,27HPROBABILITY PROBABILITY)
935 FORMAT(9H0FUNCTIONI12,4I15,39H LARGEST FN.NO.FOR LARG
1EST/92X,27HPROBABILITY PROBABILITY)
RETURN
END
C SUBROUTINE INVERT FOR BMD05M JUNE 9, 1966
SUBROUTINE INVERT (A,N,D,L,M)
C PROGRAM FOR FINDING THE INVERSE OF A NXN MATRIX
DIMENSION A(20,20),L(25),M(25)
C SEARCH FOR LARGEST ELEMENT
D=1.0
DO80 K=1,N
L(K)=K
M(K)=K
BIGA=A(K,K)
DO20 I=K,N
DO20 J=K,N
IF(ABS(BIGA)-ABS(A(I,J))) 10,20,20
10 BIGA=A(I,J)
L(K)=I
M(K)=J
20 CONTINUE
C INTERCHANGE ROWS
J=L(K)
IF(L(K)-K) 35,35,25
25 DO30 I=1,N
HOLD=-A(K,I)
A(K,I)=A(J,I)
30 A(J,I)=HOLD
C INTERCHANGE COLUMNS
35 I=M(K)
IF(M(K)-K) 45,45,37
37 DO40 J=1,N
HOLD=-A(J,K)
A(J,K)=A(J,I)
40 A(J,I)=HOLD
C DIVIDE COLUMN BY MINUS PIVOT
45 DO55 I=1,N
46 IF(I-K)50,55,50
50 A(I,K)=A(I,K)/(-A(K,K))
55 CONTINUE
C REDUCE MATRIX
DO65 I=1,N
DO65 J=1,N
56 IF(I-K) 57,65,57
57 IF(J-K) 60,65,60
60 A(I,J)=A(I,K)*A(K,J)+A(I,J)
65 CONTINUE
C DIVIDE ROW BY PIVOT
DO75 J=1,N
68 IF(J-K)70,75,70
70 A(K,J)=A(K,J)/A(K,K)
75 CONTINUE
C CONTINUED PRODUCT OF PIVOTS
D=D*A(K,K)
C REPLACE PIVOT BY RECIPROCAL
A(K,K)=1.0/A(K,K)
80 CONTINUE
C FINAL ROW AND COLUMN INTERCHANGE
K=N
100 K=(K-1)
IF(K) 150,150,103
103 I=L(K)
IF(I-K) 120,120,105
105 DO110 J=1,N
HOLD=A(J,K)
A(J,K)=-A(J,I)
110 A(J,I)=HOLD
120 J=M(K)
IF(J-K) 100,100,125
125 DO130 I=1,N
HOLD=A(K,I)
A(K,I)=-A(J,I)
130 A(J,I)=HOLD
GO TO 100
150 RETURN
END
C SUBROUTINE READ FOR BMD05M JUNE 9, 1966
SUBROUTINE READ (KVR,NVG,NTAPE,NADD)
DIMENSION N(5),X(4,15,150),XMEAN(5,25),M(25),FMT(180)
COMMON K , NN , N , LA , KP , X
COMMON ONN , XMEAN , NTO , M , K1 , K2
11 KVR=KVR*18
IERROR=0
READ (5,937)(FMT(I), I=1,KVR)
WRITE (6,938) (FMT(I),I=1,KVR)
938 FORMAT (' VARIABLE FORMAT IS',/(5X,18A4))
DO 10 I=1,K
NN=N(I)
DO 10 J=1,NN
10 READ (NTAPE,FMT)(X(I,L,J), L=1,KP)
IF(NVG)60,60,70
70 CALL TRNGEN (X,KP,N,IERROR,K,NVG)
IF(-IERROR)400,400,500
400 KP=KP+NADD
60 DO21I=1,K
NN=N(I)
ONN=NN
DO 21 L=1,KP
XMEAN(I,L)=0.0
DO 19 J=1,NN
19 XMEAN(I,L)=XMEAN(I,L)+X(I,L,J)
21 XMEAN(I,L)=XMEAN(I,L)/ONN
NTO=0
DO 22 I=1,K
NTO=NTO+N(I)
22 M(I)=I
GO TO (121,122,123,124,125), K
121 WRITE (6,931)(M(I), I=1,K)
GO TO 129
122 WRITE (6,932)(M(I), I=1,K)
GO TO 129
123 WRITE (6,933)(M(I), I=1,K)
GO TO 129
124 WRITE (6,934)(M(I), I=1,K)
GO TO 129
125 WRITE (6,935)(M(I), I=1,K)
129 WRITE (6,902)(N(I),I=1,K),NTO
WRITE (6,903)
DO 23 I=1,KP
23 WRITE (6,904)I,(XMEAN(J,I),J=1,K)
902 FORMAT(7H0SAMPLEI17,5I16)
903 FORMAT(12H0MEAN SCORES)
904 FORMAT(I7,5H 5F16.5)
931 FORMAT(6H0GROUPI18,17H TOTAL)
932 FORMAT(6H0GROUPI18,I16,17H TOTAL)
933 FORMAT(6H0GROUPI18,2I16,17H TOTAL)
934 FORMAT(6H0GROUPI18,3I16,17H TOTAL)
935 FORMAT(6H0GROUPI18,4I16,17H TOTAL)
937 FORMAT(18A4)
100 RETURN
500 NVG=-175
GO TO 100
END
C SUBROUTINE TPWD FOR BMD05M JUNE 9, 1966
SUBROUTINE TPWD(NT1,NT2)
IF(NT1)40,10,12
10 NT1=5
12 IF(NT1-NT2)14,19,14
14 IF(NT2.EQ.5)GO TO 18
17 REWIND NT2
19 IF(NT1-5)18,24,18
18 IF(NT1-6)22,40,22
22 REWIND NT1
24 NT2=NT1
28 RETURN
40 WRITE (6,49)
49 FORMAT(25H ERROR ON TAPE ASSIGNMENT)
STOP
END
C SUBROUTINE TRNGEN FOR BMD05M JUNE 9, 1966
SUBROUTINE TRNGEN (DATA,NVAR,NSAM,IERROR,NGRP,NVG)
DIMENSION DATA(4,15,150),NSAM(5)
ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
DATA A2/4HTRNG/
MARY=0
WRITE (6,1403)
WRITE (6,1400)
DO1000J=1,NVG
READ (5,1100)A1,NEWA,LCODE,LVA,BNEW
IF(A1 .EQ. A2) GO TO 2
1 WRITE (6,1408)J
IERROR=-999
GO TO 1000
2 WRITE (6,1402)J,NEWA,LCODE,LVA,BNEW
IF(LCODE*(LCODE-15))250,1,1
250 IF(LCODE-10)4,5,5
5 NEWB=BNEW
4 DO3000JJ=1,NGRP
NN=NSAM(JJ)
DO 300 I=1,NN
D=DATA(JJ,LVA,I)
GOTO(10,20,30,40,50,60,70,80,90,100,110,120,130,140),LCODE
10 IF(D)99,7,8
7 D2=0.0
GOTO3
8 D2=SQRT(D)
GOTO3
20 IF(D)99,11,12
11 D2=1.0
GOTO3
12 D2=SQRT(D)+SQRT(D+1.0)
GOTO3
30 IF(D)99,99,14
14 D2=ALOG10(D)
GOTO3
40 D2=EXP(D)
GOTO3
50 IF(-D)17,7,99
17 IF(D-1.0)18,19,99
19 D2=3.14159265/2.0
GOTO3
18 D2=ASN(SQRT(D))
GOTO3
60 FN=NN
A=D/(FN+1.0)
B=A+1.0/(FN+1.0)
IF(A)99,23,24
23 IF(-B)27,7,99
27 D2=ASN(SQRT(B))
GOTO3
24 IF(B)99,28,29
28 D2=ASN(SQRT(A))
GOTO3
29 A=SQRT(A)
B=SQRT(B)
D2=ASN(A)+ASN(B)
GOTO3
70 IF(D)31,99,31
31 D2=1.0/D
GOTO3
80 D2=D+BNEW
GOTO3
90 D2=D*BNEW
GOTO3
100 IF(-D)33,7,99
33 D2=D**NEWB
GOTO3
110 D2=D+DATA(JJ,NEWB,I)
GOTO3
120 D2=D-DATA(JJ,NEWB,I)
GOTO3
130 D2=D*DATA(JJ,NEWB,I)
GOTO3
140 IF(DATA(JJ,NEWB,I))34,99,34
34 D2=D/DATA(JJ,NEWB,I)
GOTO3
99 IF(MARY)43,44,44
44 MARY=-999
IERROR=-999
WRITE (6,1404)J
43 WRITE (6,1405)I
GO TO 300
3 DATA(JJ,NEWA,I)=D2
300 CONTINUE
IF(IERROR)45,3000,3000
45 WRITE (6,1406)JJ
3000 CONTINUE
IF(IERROR)42,1000,1000
1000 CONTINUE
IF(IERROR) 42,1111,1111
42 WRITE (6,1407)
1100 FORMAT(A4,2X,I3,I2,I3,F6.0)
1400 FORMAT(46H0CARD NEW TRANS ORIG. ORIG. VAR(B)/45H NO.
1VARIABLE CODE VAR(A) OR CONSTANT)
1402 FORMAT(2H I2,I8,2I9,F15.5)
1403 FORMAT(1H06X,23HTRANS GENERATOR CARD(S))
1404 FORMAT(30H0THE INSTRUCTIONS INDICATED ON/25H TRANS GENERATOR CARD
1NO.I2,4H RE-/29H SULTED IN THE VIOLATION OF A/31H RESTRICTION FOR
2THIS TRANSFOR-/31H MATION. THE VIOLATION OCCURRED/28H FOR THE ITEM
3S LISTED BELOW./)
1405 FORMAT(10H ITEM NO. I3)
1406 FORMAT(32H0ITEMS ABOVE ARE FROM GROUP NO. I2//)
1407 FORMAT(41H0PROGRAM CANNOT CONTINUE FOR THIS PROBLEM)
1408 FORMAT(30H0ERROR ON TRANSGENERATION CARDI4)
1111 RETURN
END