Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/bmd/bmdx75.for
There is 1 other file named bmdx75.for in the archive. Click here to see a list.
C CANONICAL ANALYSIS - MAIN PROGRAM JANUARY 15, 1966
DIMENSION DATE(2)
DOUBLE PRECISION DATE,PR,FN,SET1,SET2,SET3,P,PC,A,ONO,FMO
DATA DATE/'MAY 26, ','1969 '/
DIMENSION C(50,50),X(50),KK(50),NN(3),LN(3,50),S(50),
1D(50),LX(50),LY(50),IN(50)
DIMENSION F(162),FMO(162)
DIMENSION NTAPE(16)
COMMON/BIG/NX,LX,LY,KK1,NC,ITX,D,NV,A2,FMO,X
COMMON/NYPASS/NY
DATA PR,FN,SET1,SET2,ONO/6HPROBLM,6HFINISH,6HSETONE,6HSETTWO,
16HNO /
DATA SET3/6HSETTHR/
LON=0
CALL USAGEB('BMDX75')
P10005 = 0
100 N=0
4 READ (5,1) P,PC,NV,NC,MMM,NF,IT,A1,A2,NFX,ITX
11 IF(P.EQ.PR) GO TO 7
IF(P.EQ.FN) GO TO 5
WRITE (6,1010) P
1010 FORMAT (' PROBLM OR FINISH CARD EXPECTED BUT READ ',A6)
N=N+1
GO TO 4
5 IF(N.NE.0) WRITE (6,10) N
10 FORMAT(1H0I5,17H RECORDS NOT READ)
IF(LON.EQ.0) STOP
DO 201 I=1,LON
IF (NTAPE(I).EQ.6) GO TO 201
II=NTAPE(I)
ENDFILE II
REWIND II
201 CONTINUE
STOP
C 7 CALL USEBUF
C READ (5,1) P,PC,NV,NC,MMM,NF,IT,A1,A2,NFX,ITX
7 IF (MMM.EQ.0) MMM=4
1 FORMAT(2A6,I2,I6,2I1,I2,A2,A3,I1,I2)
IF(ITX.EQ.0) GO TO 200
LON=LON+1
NTAPE(LON)=ITX
200 IF(IT.EQ.0) IT=5
IF (A1 .NE. ONO .AND. IT .NE. 5) REWIND IT
REWIND 2
REWIND 1
IF(N.NE.0) WRITE (6,10) N
WRITE (6,2) DATE,PC,NV,NC
NN(1)=0
NN(2)=0
NN(3)=0
2 FORMAT('1BMDX75 - CANONICAL ANALYSIS - REVISED ',2A8/
X20H0PROBLEM CODE ,2X,A6/
120H0NUMBER OF VARIABLESI8,/20H0NUMBER OF CASES I8)
56 READ (5,3) A,(KK(I),I=1,19)
3 FORMAT (A6,18A4,A2)
REWIND 1
WRITE (1,3) A,(KK(I),I=1,19)
IF(A.EQ.SET1) GO TO 6
IF(A.EQ.SET2) GO TO 31
C
C SET3 WAS INTENDED FOR USE IN A FUTURE MODIFICATION OF THIS
C PROGRAM TO ALLOW FOR THE ANALYSIS OF THREE SETS OF DATA
C
IF(A.EQ.SET3) GO TO 33
C
GO TO 53
33 K=3
GO TO 32
6 K=1
GO TO 32
31 K=2
32 REWIND 1
READ (1,61) A,(KK(I),I=1,22)
61 FORMAT(A6,22I3)
DO 62 I=1,22
IF(KK(I))63,62,65
63 L0=KK(I-1)+1
L1=-KK(I)
IF((L0-1)*(NV-L0).LT.0 .OR. (L1-1)*(NV-L1).LT.0) GO TO 10000
DO 66 L=L0,L1
NN(K)=NN(K)+1
N=NN(K)
66 LN(K,N)=L
GO TO 62
65 NN(K)=NN(K)+1
N=NN(K)
LN(K,N)=KK(I)
IF((KK(I)-1)*(NV-KK(I)))10000,62,62
62 CONTINUE
GO TO 56
10000 WRITE (6,10001) A
10001 FORMAT(' INDEX NUMBER EXCEEDS NUMBER OF VARIABLES FOR ',A6)
GO TO 80
C
C
10006 WRITE(6,10007)
GO TO 80
10008 WRITE(6,10009)
C
10005 FORMAT ('NUMBER OF VARIABLE FORMAT CARDS INCORRECT. ASSUMED TO BE
11')
10007 FORMAT(' TOO MANY VARIABLES')
10009 FORMAT(' NUMBER OF VARIABLES CANNOT BE NEGATIVE')
80 WRITE (6,81)
81 FORMAT(45H0CONTROL CARDS INCORRECTLY ORDERED OR PUNCHED)
GO TO 100
53 IF (NF.GT.0.AND.NF.LE.9) GO TO 54
P10005 = 1
NF=1
54 IF (ITX .EQ . 0) GO TO 55
IF (NFX .GT. 0 .AND. NFX .LE. 9) GO TO 55
P10005 = 1
NFX=1
55 NF=NF*18
NFX=NFX*18
IF(NN(1).EQ.0 .OR. NN(2).EQ.0)GO TO 10000
IF(NV.GT.50)GO TO 10006
IF(NV.LT.0)GO TO 10008
REWIND 1
READ (1,41) (F(I),I=1,NF)
REWIND 1
41 FORMAT(18A4)
IF (P10005 .EQ. 0) GO TO 58
WRITE(6,10005)
P10005 = 0
58 WRITE (6,1011) (F(I),I=1,NF)
IF(ITX .EQ. 0) GO TO 59
1011 FORMAT (' INPUT FORMAT'/(5X,18A4))
READ (5,41) (FMO(I),I=1,NFX)
WRITE (6,1012) (FMO(I),I=1,NFX)
1012 FORMAT (' OUTPUT FORMAT'/(5X,18A4))
59 NNT =21
IF(NN(3).GT.0) NNT=22
REWIND NNT
DO 20 I=1,NV
S(I)=0.
DO 20 J=1,NV
20 C(I,J)=0.
DO 21 L=1,NC
H=L
H1=H*(H-1.)
READ (IT,F) (X(I),I=1,NV)
WRITE (NNT) (X(I),I=1,NV)
23 DO 21 I=1,NV
D(I)=(X(I)-S(I))/H
S(I)=S(I)+D(I)
PQ=D(I)*H1
DO 21 J=1,I
21 C(I,J)=C(I,J)+PQ*D(J)
REWIND NNT
H=NC-1
DO 110 I=1,NV
D(I)=1.
110 X(I)=SQRT(C(I,I)/H)
GO TO (101,102,101,102) ,MMM
105 DO 120 I=1,NV
DO 120 J=1,I
120 C(I,J)=C(I,J)/H
GO TO 108
101 H=NC
DO 104 I=1,NV
PQ=H*S(I)
DO 104 J=1,I
104 C(I,J)=C(I,J)+PQ*S(J)
102 GO TO (105,105,106,106),MMM
106 DO 107 I=1,NV
D(I)=C(I,I)
IF(D(I).LE.0.) D(I)=1.
D(I)=SQRT(D(I))
DO 107 J=1,I
107 C(I,J)=C(I,J)/(D(I)*D(J))
H=SQRT(H)
DO 109 I=1,NV
109 D(I)=D(I)/H
108 DO 24 I=1,NV
DO 24 J=1,I
24 C(J,I)=C(I,J)
WRITE (6,85)
85 FORMAT(1H-,4(32H VARIABLE MEAN STANDARD )/1X,4(23X,9HDEVIA
1TION))
K=(NV-1)/4+1
DO 87 I=1,K
87 WRITE (6,86) (J,S(J),X(J),J=I,NV,K)
86 FORMAT(1X,4(I7,F13.6,F12.6))
GO TO (821,822,823,824),MMM
821 WRITE (6,921)
GO TO 960
822 WRITE (6,922)
GO TO 960
823 WRITE (6,923)
GO TO 960
824 WRITE (6,924)
921 FORMAT(76H1THE COVARIANCE MATRIX ABOUT THE ORIGIN IS USED IN THE F
XOLLOWING CALCULATION)
922 FORMAT(59H1THE COVARIANCE MATRIX IS USED IN THE FOLLOWING CALCULAT
XION)
923 FORMAT(77H1THE CORRELATION MATRIX ABOUT THE ORIGIN IS USED IN THE
XFOLLOWING CALCULATION)
924 FORMAT(60H1THE CORRELATION MATRIX IS USED IN THE FOLLOWING CALCULA
XTION)
960 K=1
IF(NN(1).LE.NN(2)) K=2
KK1=K
L=3-K
NX=NN(K)
DO 70 I=1,NX
70 LX(I)=LN(K,I)
NY=NN(L)
J=NX
DO 71 I=1,NY
LY(I)=LN(L,I)
J=J+1
71 LX(J)=LY(I)
NZ=NN(3)
IF(NZ.EQ.0) GO TO 35
WRITE (6,38) (LN(3,I),I=1,NZ)
38 FORMAT(22H0SET THREE VARIABLES -/(5X,40I3))
DO 34 I=1,NV
34 IN(I)=1
DO 36 I=1,NZ
J=LN(3,I)
36 IN(J)=0
CALL MATINV(C,NV,IN,1.E-6)
DO 305 I=1,NV
DO 305 J=1,I
305 C(J,I)=C(I,J)
DO 303 L=1,NC
READ(22) (X(I),I=1,NV)
DO 302 I=1,NZ
J=LN(3,I)
XJ=X(J)
DO 302 K=1,NV
302 X(K)=X(K)-C(J,K)*XJ
DO 301 I=1,NZ
J=LN(3,I)
301 X(J)=0.
303 WRITE (21) (X(I),I=1,NV)
1000 FORMAT (20A4)
REWIND 21
DO 37 I=1,NZ
J=LN(3,I)
DO 37 K=1,NV
C(K,J)=0.
37 C(J,K)=0.
35 NV1=NX+NY
DO 72 I=1,NV
72 IN(I)=I
DO 73 I=1,NV1
J=LX(I)
76 IF(LX(I).EQ.IN(J)) GO TO 75
J=IN(J)
GO TO 76
75 IF(I.EQ.J) GO TO 73
M=IN(J)
IN(J)=IN(I)
IN(I)=M
DO 78 K=1,NV
T=C(I,K)
C(I,K)=C(J,K)
78 C(J,K)=T
DO 77 K=1,NV
T=C(K,I)
C(K,I)=C(K,J)
77 C(K,J)=T
73 CONTINUE
NX1=NX+1
CALL CANON(C,C(NX1,1),C(NX1,NX1),C(1,NX1))
C CALL CANON(C,C(NX1,1),C(NX1,NX1),C(1,NX1),NX,NY,LX,LY,KK1,NC,ITX,D
C 1,NX,A2,FMO)
GO TO 100
END
C SUBROUTINE CANON FOR BMDX75 JANUARY 15, 1966
SUBROUTINE CANON(A,B,C,E)
DOUBLE PRECISION FMO
DATA YES/3HYES/
DIMENSION FMO(162)
DIMENSION D(50),X(50)
DIMENSION A(50,50),B(50,50),C(50,50),E(50,50),R(50),IN(5
10),DX(50),DY(50),LX(50),LY(50)
COMMON/BIG/NX,LX,LY,KK1,NC,ITS,D,MV,A2,FMO,X
COMMON/NYPASS/NY
COMMON/NRPASS/NR
T1=1.E-4
T2=T1
NV=NX+NY
DO 2 J=1,NX
2 IN(J)=0
DO 1 I=1,NY
J=I+NX
IN(J)=1
DO 1 J=1,I
E(I,J)=C(I,J)
1 C(I,J)=0.
CALL MATINV(A,NV,IN,T1)
DO 3 I=1,NY
DO 3 J=1,I
C(I,J)=-C(I,J)
C(J,I)=C(I,J)
3 E(J,I)=E(I,J)
CALL JACOB2 (C,E,1.E-7,1,T2)
DO 4 I=1,NY
4 R(I)=SIGN(SQRT(ABS(C(I,I))),C(I,I))
DO 5 I=1,NX
DO 5 J=1,NY
T=0.
DO 8 K=1,NY
IF(R(J).LE.0.) E(K,J)=0.
8 T=T+E(K,J)*B(K,I)
IF(R(J).LE.0.) A(I,J)=0.
5 IF(R(J).GT.0.) A(I,J)=T/R(J)
L1=0
30 L0=L1+1
L1= MIN0(L1+10,NY)
LL=0
WRITE (6,31) (I,I=L0,L1)
31 FORMAT(1H-15X,22HCANONICAL CORRELATIONS//(2X,10I12))
WRITE (6,32) (R(I),I=L0,L1)
32 FORMAT(6X,10F12.5)
WRITE (6,33)
GO TO (35,34),KK1
33 FORMAT(66H0VARIABLE COEFFICIENTS FOR CANONICAL VARIABLES OF THE
1 FIRST SET//)
35 DO 48 I=1,NX
48 WRITE (6,39) LX(I),(A(I,L),L=L0,L1)
39 FORMAT(I6,10F12.5)
49 IF(LL)40,40,41
40 LL=1
WRITE (6,43)
GO TO (34,35) ,KK1
43 FORMAT(67H0VARIABLE COEFFICIENTS FOR CANONICAL VARIABLES OF THE
1 SECOND SET//)
34 DO 38 I=1,NY
38 WRITE (6,39) LY(I),(E(I,L),L=L0,L1)
GO TO 49
41 IF(L1-NY)30,47,30
47 IF(A2.NE.YES) RETURN
DO 304 I=1,NY
DO 302 J=1,NX
K=LX(J)
302 A(J,I)=A(J,I)/D(K)
DO 304 J=1,NY
K=LY(J)
304 E(J,I)=E(J,I)/D(K)
WRITE (6,180) (I,I=1,NY)
180 FORMAT(1H1,20X,19HCANONICAL VARIABLES//6(10X,I7,4X))
WRITE (6,181)
181 FORMAT(6H CASE/6H NO.)
REWIND 21
DO 100 L=1,NC
LL=0
DO 101 IX=1,50
101 X(IX)=0.0
8765 FORMAT(7F10.5)
READ (21) (X(I),I=1,MV)
1000 FORMAT (20A4)
DO 102 I=1,NY
DX(I)=0.
DO 102 J=1,NX
K=LX(J)
102 DX(I)=DX(I)+A(J,I)*X(K)
DO 103 I=1,NY
DY(I)=0.
DO 103 J=1,NY
K=LY(J)
103 DY(I)=DY(I)+E(J,I)*X(K)
GO TO (160,161),KK1
160 IF(ITS.NE.0) WRITE (ITS,FMO) (X(I),I=1,NV),(DX(I),DY(I),I=1,NY)
WRITE (6,170) L,(DX(I),DY(I),I=1,NY)
170 FORMAT(1H0,I4,6(2F10.5,1X)/(5X,2F10.5,1X,2F10.5,1X,2F10.5,1X,2F10.
X5,1X,2F10.5,1X,2F10.5,1X))
GO TO 100
161 IF(ITS.NE.0) WRITE (ITS,FMO) (X(I),I=1,NV),(DY(I),DX(I),I=1,NY)
WRITE (6,170) L,(DY(I),DX(I),I=1,NY)
100 CONTINUE
REWIND 21
RETURN
END
C SUBROUTINE JACOB2 FOR BMDX75 JANUARY 15, 1966
SUBROUTINE JACOB2 (A,B,ACC,IV,TOL)
DIMENSION A(50,2),B(50,2),U(50),IN(50),V(50)
COMMON/NYPASS/N
COMMON/NRPASS/NR
Q=1.
DO 1 I=1,N
V(I)=B(I,I)
1 IN(I)=0
K=1
DO 13 L=1,N
DO 2 I=1,N
U(I)=B(I,K)
2 B(K,I)=0.
P=U(K)
C
U(K)=1.
S=Q-TOL
IF(S.GT.0.) RP=SQRT(P)
Q=-1.
IN(K)=-1
DO 11 I=1,N
Y=U(I)/P
IF(IN(I))6,3,11
3 IF(S)16,16,17
17 DO 4 J=1,N
4 B(J,I)=B(J,I)-U(J)*Y
16 IF(B(I,I)/V(I)-Q)9,9,5
5 Q=B(I,I)/V(I)
KK=I
9 IF(S)11,11,18
18 T=A(I,I)-Y*A(I,K)
DO 10 J=1,N
A(I,J)=A(I,J)-Y*A(K,J)
10 A(J,I)=A(I,J)
A(I,I)=T-Y*A(I,K)
GO TO 11
6 IF(S)11,11,23
23 DO 8 J=1,N
IF(IN(J))7,14,7
14 B(J,I)=0.
GO TO 8
7 B(J,I)=U(J)/RP
8 CONTINUE
11 CONTINUE
IF(S)19,19,20
19 DO 21 I=1,N
B(I,K)=0.
A(I,K)=0.
21 A(K,I)=0.
GO TO 22
20 DO 12 I=1,N
A(I,K)=A(I,K)/RP
12 A(K,I)=A(I,K)
A(K,K)=A(K,K)/RP
22 IN(K)=1
13 K=KK
CALL JACOBI(A,B,ACC,2*IV,NR)
RETURN
END
C SUBROUTINE JACOBI FOR BMDX75 JANUARY 15, 1966
SUBROUTINE JACOBI (A,B,ACC,IV)
DIMENSION A(50,50),B(50,50),LK(50),Q(50)
COMMON/NYPASS/N
COMMON/NRPASS/NR
IF(IV-1)200,201,200
201 DO 202 I=1,N
DO 203 J=1,N
203 B(I,J)=0.
202 B(I,I)=1.
200 NR=0
Q(1)=0.0
W=0.
H=0.5*A(1,1)*A(1,1)
DO 1 I=2,N
H=H+.5*A(I,I)*A(I,I)
Q(I)=0.
I1=I-1
DO 2 J=1,I1
Z=ABS(A(I,J))
H=H+Z*Z
IF(Z-Q(I))2,2,3
3 Q(I)=Z
LK(I)=J
2 CONTINUE
IF(Q(I)-W)1,1,4
4 W=Q(I)
III=I
1 CONTINUE
H=ACC*SQRT(2.*H)/FLOAT(N)
30 II=LK(III)
JJ=III
X=A(II,II)
Y=A(JJ,II)
Z=A(JJ,JJ)
W=X-Z
T=.5*(W+SQRT(W*W+4.*Y*Y))/Y
W=SQRT(1.+T*T)
S=T/W
C=1./W
CC=C*C
SS=S*S
SC=S*C*2.
Q1=0.
Q2=0.
W=0.
NR=NR+1
LK(1)=0
DO 27 I=1,N
LKI=LK(I)
IF(I-II)10,11,12
10 U=A(II,I)
V=A(JJ,I)
E=U*S+V*C
A(II,I)=E
IF(ABS(E)-Q1)15,15,14
14 Q1=ABS(E)
I1=I
15 F=V*S-U*C
A(JJ,I)=F
IF(LKI-II)100,101,100
101 Q(I)=-1.E20
I3=I-1
DO 103 J=1,I3
IF(ABS(A(I,J))-Q(I))103,103,104
104 LK(I)=J
Q(I)=ABS(A(I,J))
103 CONTINUE
100 IF(ABS(F)-Q2)9,9,16
16 Q2=ABS(F)
I2=I
GO TO 9
11 A(II,I)=SS*X+SC*Y+CC*Z
Q(I)=Q1
LK(I)=I1
GO TO 9
12 IF(I-JJ)17,18,19
17 U=A(I,II)
V=A(JJ,I)
E=S*U+C*V
A(I,II)=E
IF(ABS(E)-Q(I))15,15,21
21 LK(I)=II
Q(I)=ABS(E)
GO TO 15
18 A(JJ,I)=CC*X-SC*Y+SS*Z
A(I,II)=0.
Q(I)=Q2
LK(I)=I2
GO TO 9
19 U=A(I,II)
V=A(I,JJ)
E=U*S+V*C
F=V*S-U*C
A(I,II)=E
A(I,JJ)=F
IF((LK(I)-II)*(LK(I)-JJ)) 105,101,105
105 G=AMAX1(ABS(E),ABS(F))
IF(G-Q(I))9,9,13
13 Q(I)=G
IF(ABS(E)-ABS(F))23,24,24
24 LK(I)=II
GO TO 9
23 LK(I)=JJ
9 IF(Q(I)-W)40,25,25
25 W=Q(I)
III=I
40 IF(IV)27,27,33
33 U=B(I,II)
V=B(I,JJ)
B(I,II)=U*S+V*C
B(I,JJ)=V*S-U*C
27 CONTINUE
IF(W-H)31,31,30
31 DO 50 I=1,N
U=-1.E20
DO 51 J=I,N
IF(A(J,J)-U) 51,51,52
52 U=A(J,J)
K=J
51 CONTINUE
IF(K-I)53,50,53
53 A(K,K)=A(I,I)
A(I,I)=U
IF(IV)50,50,54
54 DO 55 J=1,N
T=B(J,I)
B(J,I)=B(J,K)
55 B(J,K)=T
50 CONTINUE
RETURN
END
C SUBROUTINE MATINV FOR BMDX75 JANUARY 15, 1966
SUBROUTINE MATINV(A,N,IN,T)
DIMENSION A(50,50),IN(50),U(50),V(50)
H=1.
DO 1 I=1,N
IF(IN(I))1,10,1
10 K=I
1 V(I)=A(I,I)
8 DO 2 I=1,K
U(I)=A(K,I)
2 A(K,I)=0.
P=U(K)
DO 3 I=K,N
U(I)=A(I,K)
3 A(I,K)=0.
S=H-T
H=-1.
IN(K)=1
U(K)=-1.
DO 7 I=1,N
IF(S)11,11,12
12 Y=U(I)/P
DO 4 J=1,I
4 A(I,J)=A(I,J)-Y*U(J)
11 IF(IN(I))5,5,7
5 IF(H*V(I).GE.A(I,I))GO TO 7
6 H=A(I,I)/V(I)
K=I
7 CONTINUE
IF(H+1.) 8,9,8
9 RETURN
END