Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/bmd/bmd01m.for
There is 1 other file named bmd01m.for in the archive. Click here to see a list.
C PRINCIPAL COMPONENT ANALYSIS MAY 10, 1966
C THIS IS A SIFTED VERSION OF BMD01M 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 A123,B123, TODE,FINISH,PROBLM,NPROB,NAMES
DIMENSION X(250,20),XMEAN(25),COV(20,20),VALU(25),SCALE(25),
1C(250,20),Z(20,20),FMT(180),NAMES(25)
COMMON X , COV,C,Z
DATA A123,B123,C123/6HFINISH,6HPROBLM,4HYES /
C
209 FORMAT('1BMD01M - COMPONENT ANALYSIS - REVISED ',
1'APRIL 30, 1969'/
240H HEALTH SCIENCES COMPUTING FACILITY,UCLA//
314H PROBLEM CODE A6,/
421H NUMBER OF VARIABLES I3,/
517H NUMBER OF CASES I6,/
627H NUMBER OF VARIABLES ADDED I4,/
735H NUMBER OF TRANSGENERATION CARD(S) I4,/
8 35H NUMBER OF VARIABLE FORMAT CARD(S) I3,///)
C
NTAPE=5
CALL USAGEB('BMD01M')
10 READ (5,901)TODE,NPROB,NV,N,RNCR,GCK,NADD,NVG ,NLV,MTAPE,KVR
IERROR=0
IF(A123 .EQ. TODE) GO TO 201
GO TO 200
202 WRITE (6,204) TODE
201 IF(NTAPE-5)12,12,11
11 REWIND NTAPE
12 STOP
200 IF(B123 .NE. TODE) GO TO 202
203 CALL TPWD(MTAPE,NTAPE)
9 IF((NV-1)*(NV-21)) 205,300,300
205 IF((N-2)*(N-251)) 206,301,301
206 IF((NV+NADD-1)*(NV+NADD-21)) 207,302,302
207 IF(KVR.GT.0.AND.KVR.LE.10)GO TO 208
KVR=1
WRITE(6,4000)
208 WRITE (6,209)NPROB,NV,N,NADD,NVG,KVR
211 NV1=NV+NADD
CALL RDLBL(NLV,NV1,NAMES)
17 KVR=KVR*18
READ (5,942)(FMT(I),I=1,KVR)
WRITE (6,310) (FMT(I),I=1,KVR)
310 FORMAT (' VARIABLE FORMAT IS'/(1X,18A4))
DO 13 I=1,N
13 READ (NTAPE,FMT)(X(I,J),J=1,NV)
19 ON=N
IF(NVG) 303,600,601
601 CALL TRANS (X,NV,N,IERROR,NVG)
IF(IERROR)10,600,600
600 IF(NV) 300,300,212
212 NV=NV1
DO 21 J=1,NV
XMEAN(J)=0.0
DO 20 I=1,N
20 XMEAN(J)=XMEAN(J)+X(I,J)
21 XMEAN(J)=XMEAN(J)/ON
DO 22 I=1,NV
DO 22 J=1,NV
COV(I,J)=0.0
DO 22 K=1,N
22 COV(I,J)=COV(I,J)+(X(K,I)-XMEAN(I))*(X(K,J)-XMEAN(J))
DO 23 I=1,NV
23 SCALE(I)=SQRT(COV(I,I))
DO 24 I=1,NV
DO 24 J=1,NV
24 Z(I,J)=COV(I,J)/(SCALE(I)*SCALE(J))
WRITE (6,923)
WRITE (6,904)
CALL PATTY2(Z,NV,NAMES,1)
CALL EIGEN(VALU,NV,NV)
WRITE (6,923)
WRITE (6,907)
WRITE (6,906)(VALU(I),I=1,NV)
RANK=0.0
DO 26 I=1,NV
26 RANK=RANK+VALU(I)
SMALL=0.0
DO 18 I=1,NV
SMALL=SMALL+VALU(I)
18 VALU(I)=SMALL/RANK
WRITE (6,937)
WRITE (6,938)(VALU(I),I=1,NV)
WRITE (6,923)
WRITE (6,908)
CALL PATTY2(Z,NV,NAMES,0)
DO 29 J=1,NV
DO 29 I=1,N
29 X(I,J)=(X(I,J)-XMEAN(J))/SCALE(J)
ONN=N-1
SQ=SQRT(ONN)
DO 43 I=1,N
DO 43 J=1,NV
C(I,J)=0.0
DO 42 K=1,NV
42 C(I,J)=C(I,J)+X(I,K)*Z(K,J)
43 C(I,J)=C(I,J)*SQ
IF(GCK .NE. C123) GO TO 57
41 DO 51 J=1,NV
XMEAN(J)=0.0
DO 50 I=1,N
50 XMEAN(J)=XMEAN(J)+C(I,J)
51 XMEAN(J)=XMEAN(J)/ON
DO 52 I=1,NV
DO 52 J=1,NV
COV(I,J)=0.0
DO 525 K=1,N
525 COV(I,J)=COV(I,J)+(C(K,I)-XMEAN(I))*(C(K,J)-XMEAN(J))
52 COV(I,J)=COV(I,J)/ONN
WRITE (6,923)
WRITE (6,922)
CALL PATTY2(COV,NV,NAMES,1)
57 IF(RNCR .NE. C123) GO TO 40
30 WRITE (6,923)
WRITE (6,909)
WRITE (6,910)
SMALL=-(10.0**36.0)
DO 39 II=1,NV
WRITE (6,912)II
DO 32 I=1,N
C(I,1)=0.0
C(I,2)=0.0
DO 31 K=1,NV
31 C(I,1)=C(I,1)+X(I,K)*Z(K,II)
32 C(I,1)=C(I,1)*SQ
DO 39 I=1,N
RANK=SMALL
DO 38 J=1,N
IF(C(J,1)-RANK)38,38,36
36 IF(C(J,2)-999.0)37,38,38
37 RANK=C(J,1)
NJ=J
38 CONTINUE
C(NJ,2)=999.0
WRITE (6,911)RANK,NJ
39 CONTINUE
40 GO TO 10
204 FORMAT (' PROGRAM EXPECTED PROBLM OR FINISH CARD BUT READ ',A6)
300 WRITE (6,312) NV
312 FORMAT (1X,I3,' ORIGINAL VARIABLES IS ILLEGAL')
GO TO 201
301 WRITE (6,313) N
313 FORMAT (1X,I4,' CASES IS ILLEGAL')
GO TO 201
302 NDOPE=NV+NADD
WRITE (6,314) NDOPE
314 FORMAT (1X,I4,' TOTAL VARIABLES IS ILLEGAL')
GO TO 201
303 WRITE (6,315) NVG
315 FORMAT (1X,I3,' IS ILLEGAL NUMBER OF TRANSGENERATION CARDS')
GO TO 201
901 FORMAT(2A6,I2,I3,2A3,I3,2I2,38X,2I2)
902 FORMAT(19H0COMPONENT ANALYSIS)
903 FORMAT(12H PROBLEM NO.I4)
904 FORMAT(31H0CORRELATION COEFFICIENT MATRIX)
907 FORMAT(12H0EIGENVALUES)
908 FORMAT(13H0EIGENVECTORS)
909 FORMAT(48H0RANK ORDER OF EACH STANDARDIZED CASE ORDERED BY)
910 FORMAT(44H SIZE OF EACH PRINCIPAL COMPONENT SEPARATELY)
911 FORMAT(F18.6,I10)
912 FORMAT(16H0 COMPONENT NO.I3,12H CASE NO.)
922 FORMAT(25H0EIGEN VALUE CHECK MATRIX)
923 FORMAT(1H0)
937 FORMAT(40H0CUMULATIVE PROPORTION OF TOTAL VARIANCE)
938 FORMAT(1H F11.2,7F15.2)
942 FORMAT(18A4)
906 FORMAT(6F16.7)
4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
1IED, ASSUMED TO BE 1.)
END
C FUNCTION ANUMB FOR BMD 02M AUGUST 17, 1966
C
C THE FUNCTION 'ANUMB' CONVERTS THE INTEGER 'I' TO RIGHT JUSTIFIED
C ALPHANUMERIC CHARACTERS WHICH ARE RETURNED AS THE HIGH ORDER FOUR
C BYTES OF THE DOUBLE PRECISION VARIABLE 'ANUMB'.
C
DOUBLE PRECISION FUNCTION ANUMB(II)
DIMENSION N(2)
DOUBLE PRECISION B
EQUIVALENCE (N(1),B)
ENCODE(10,11,N)II
11 FORMAT(4X,I4,2X)
ANUMB=B
RETURN
END
C SUBROUTINE EIGEN FOR BMD01M MAY 10, 1966
SUBROUTINE EIGEN(VALU,N,M)
C
C EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX
C
DIMENSION A(20,20), B(20,20), VALU(25), DIAG(25), SUPERD(24),
1 Q(24), VALL(25), S(24), C(24), D(25), IND(25), U(25),
2 DUMMY(5000)
DIMENSION X123(250,20),Y123(250,20),DUMY1(250,20),DUMY2(400)
COMMON X123 , DUMY2 , Y123
COMMON A
EQUIVALENCE (X123,DUMY1),(B,DUMY2),(SUPERD,DUMMY(26)),(TAU,BETA),
1 (VALL,D,DUMMY(50)),(Q,S,DUMMY(75)),(IND,U),(II,MATCH),
2 (DIAG,DUMMY,Y123),(ANORM,ANORM2),(P,PRODS),(T,SMALLD)
C
C CALCULATE NORM OF MATRIX
C
3 DPA=0.0
4 DO 6 I=1,N
5 DO 6 J=1,N
6 DPA=DPA+A(I,J)**2
7 ANORM=DPA
ANORM=SQRT(ANORM)
C
C GENERATE IDENTITY MATRIX
C
9 IF (M) 10, 45, 10
10 DO 40 I=1,N
12 DO 40 J=1,N
20 IF(I-J) 35, 25, 35
25 B(I,J)=1.0
30 GO TO 40
35 B(I,J)=0.0
40 CONTINUE
C
C PERFORM ROTATIONS TO REDUCE MATRIX TO JACOBI FORM
C
45 IEXIT=1
50 NN=N-2
52 IF (NN) 890, 170, 55
55 DO 160 I=1,NN
60 II=I+2
65 DO 160 J=II,N
70 T1=A(I,I+1)
75 T2=A(I,J)
80 GO TO 900
90 DO 105 K=I,N
95 T2=COS*A(K,I+1)+SIN*A(K,J)
100 A(K,J)=COS*A(K,J)-SIN*A(K,I+1)
105 A(K,I+1)=T2
110 DO 125 K=I,N
115 T2=COS*A(I+1,K)+SIN*A(J,K)
120 A(J,K)=COS*A(J,K)-SIN*A(I+1,K)
125 A(I+1,K)=T2
128 IF (M) 130, 160, 130
130 DO 150 K=1,N
135 T2=COS*B(K,I+1)+SIN*B(K,J)
140 B(K,J)=COS*B(K,J)-SIN*B(K,I+1)
150 B(K,I+1)=T2
160 CONTINUE
C
C MOVE JACOBI FORM ELEMENTS AND INITIALIZE EIGENVALUE BOUNDS
C
170 DO 200 I=1,N
180 DIAG(I)=A(I,I)
190 VALU(I)=ANORM
200 VALL(I)=-ANORM
210 DO 230 I=2,N
220 SUPERD(I-1)=A(I-1,I)
230 Q(I-1)=(SUPERD(I-1))**2
C
C DETERMINE SIGNS OF PRINCIPAL MINORS
C
235 TAU=0.0
240 I=1
260 MATCH=0
270 T2=0.0
275 T1=1.0
277 DO 450 J=1,N
280 P=DIAG(J)-TAU
290 IF(T2)300,330,300
300 IF(T1)310,370,310
310 T=P*T1-Q(J-1)*T2
320 GO TO 410
330 IF(T1)335,350,350
335 T1=-1.0
340 T=-P
345 GO TO 410
350 T1=1.0
355 T=P
360 GO TO 410
370 IF(Q(J-1))380,350,380
380 IF(T2)400,390,390
390 T=-1.0
395 GO TO 410
400 T=1.0
C
C COUNT AGREEMENTS IN SIGN
C
410 IF(T1)425,420,420
420 IF(T)440,430,430
425 IF(T)430,440,440
430 MATCH=MATCH+1
440 T2=T1
450 T1=T
C
C ESTABLISH TIGHTER BOUNDS
C
460 DO 530 K=1,N
465 IF(K-MATCH)470,470,520
470 IF(TAU-VALL(K))530,530,480
480 VALL(K)=TAU
490 GO TO 530
520 IF(TAU-VALU(K))525,530,530
525 VALU(K)=TAU
530 CONTINUE
540 IF(ABS(VALU(I)-VALL(I)).GT.ABS(VALU(I)+VALL(I))/1.E6)GO TO 580
570 I=I+1
575 IF(I-N)540,540,590
580 TAU=(VALL(I)+VALU(I))/2.0
585 GO TO 260
C
C JACOBI EIGENVECTORS BY ROTATIONAL TRIANGULARIZATION
C
590 IF (M) 593, 890, 593
593 IEXIT=2
595 DO 610 I=1,N
600 DO 610 J=1,N
610 A(I,J)=0.0
615 DO 850 I=1,N
620 IF (I-1) 625, 625, 621
621 IF(ABS(VALU(I-1)-VALU(I)).LT.ABS(VALU(I-1)+VALU(I))/1.E+6) GOTO730
625 COS=1.0
628 SIN=0.0
630 DO 700 J=1,N
635 IF(J-1) 680, 680, 640
640 GO TO 900
650 S(J-1)=SIN
660 C(J-1)=COS
670 D(J-1)=T1*COS+T2*SIN
680 T1=(DIAG(J)-VALU(I))*COS-BETA*SIN
690 T2=SUPERD(J)
700 BETA=SUPERD(J)*COS
710 D(N)=T1
720 DO 725 J=1,N
725 IND(J)=0
730 SMALLD=ANORM
735 DO 780 J=1,N
740 IF (IND(J)-1) 750, 780, 780
750 IF (ABS(SMALLD)-ABS(D(J)))780, 780, 760
760 SMALLD=D(J)
770 NN=J
780 CONTINUE
790 IND(NN)=1
800 PRODS=1.0
805 IF (NN-1) 810, 850, 810
810 DO 840 K=2,NN
820 II=NN+1-K
830 A(II+1,I)=C(II)*PRODS
840 PRODS=-PRODS*S(II)
850 A(1,I)=PRODS
C
C FORM MATRIX PRODUCT OF ROTATION MATRIX WITH JACOBI VECTOR MATRIX
C
855 DO 885 J=1,N
860 DO 865 K=1,N
865 U(K)=A(K,J)
870 DO 885 I=1,N
875 A(I,J)=0.0
880 DO 885 K=1,N
885 A(I,J)=B(I,K)*U(K)+A(I,J)
890 GO TO 941
C
C CALCULATE SINE AND COSINE OF ANGLE OF ROTATION
C
900 IF (T2) 910, 940, 910
910 T=SQRT(T1**2+T2**2)
920 COS=T1/T
925 SIN=T2/T
930 GO TO (90,650), IEXIT
940 GO TO (160,910), IEXIT
941 RETURN
END
C SUBROUTINE PATTY2 FOR BMD01M MAY 10, 1966
SUBROUTINE PATTY2(A,N,NAMES,JK)
DOUBLE PRECISION NN,NAMES
DIMENSION A(20,20),NAMES(25),NN(8)
IT=1
KK=0
K1=IT
K2=MIN0(8,N)
5 KK=KK+8
IF(N-KK)3,3,4
4 IT=IT+1
GO TO 5
3 DO 50 JX=1,IT
LLL=K2-K1+1
LL=0
IF(JK)35,35,37
35 WRITE (6,350)(IG,IG=1,LLL)
GO TO 45
37 DO 40 JJ=K1,K2
LL=LL+1
40 NN(LL)=NAMES(JJ)
WRITE (6,300)(NN(II),II=1,LLL)
45 DO 10 I=1,N
10 WRITE (6,20)NAMES(I),(A(I,J),J=K1,K2)
K1=K2+1
K2=K1+7
K2=MIN0(K2,N)
300 FORMAT(1H013X,A6,7(8X,A6)/)
20 FORMAT(1H A6,1X,8F14.4)
350 FORMAT(1H017X,I2,7(12X,I2)/)
50 CONTINUE
RETURN
END
C SUBROUTINE RDLBL FOR BMD01M MAY 10, 1966
C SUBROUTINE TO READ IN LABELS CARDS, STORE THEM IN ARRAY,
C AND SUBSTITUTE NUMBERS FOR UNLABELED VARIABLES
C NVAR IS TOTAL NUMBER OF VARIABLES
C NLBVAR IS NUMBER OF LABELED VARIABLES EXPECTED
C
SUBROUTINE RDLBL(NLBVAR,NVAR,ARRAY)
C EQUIVALENCE INTEGER AND FLOATING NAMES SO THAT INTEGER SUBTRACTION
C MAY BE USED TO TEST ALPHABETIC EQUALITY
DOUBLE PRECISION ARRAY,ANUMB,DUMY,TEST,ALABEL
DIMENSION ARRAY(1),IDUM(7),DUMY(7)
DATA ALABEL/6HLABELS/
C NUMBER VARIABLES
DO 1 I=1,NVAR
1 ARRAY(I)=ANUMB(I)
C IF NO LABELS, RETURN
IF(NLBVAR) 9,9,2
2 N=0
C READ 1 LABELS CARD
20 READ (5,3) TEST,(IDUM(J),DUMY(J),J=1,7)
3 FORMAT(A6 ,7(I4,A6))
C TEST FOR 'LAB' IN FIRST 3 COLS.
IF(TEST.EQ.ALABEL)GO TO 6
C ERROR--PRINT MESSAGE AND QUIT
4 WRITE (6,5) TEST
5 FORMAT (' LABELS CARD EXPECTED BUT READ ',A6)
STOP
C EXAMINE 7 FIELDS
6 DO 8 J=1,7
K=IDUM(J)
C TEST INDEX. IF 0, IGNORE. IF ILLEGAL, PRINT MESSAGE AND
C IGNORE EXCEPT TO COUNT
IF(K) 11,8,10
10 IF(K-NVAR) 7,7,11
11 WRITE (6,12)K,DUMY(J)
12 FORMAT(18H0LABELS CARD INDEX,I7,18H INCORRECT. LABEL ,A6,9H IGNORE
1D.)
GO TO 13
C MOVE LABEL TO ARRAY
7 ARRAY(K)=DUMY(J)
C STEP NUMBER OF VARIABLES
13 N=N+1
C TEST FOR END. IF END, RETURN. IF NOT, SCAN OTHER FIELDS.
IF(N-NLBVAR) 8,9,9
8 CONTINUE
GO TO 20
9 RETURN
END
C SUBROUTINE TPWD FOR BMD01M MAY 10, 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 TRANS FOR BMD01M MAY 10, 1966
SUBROUTINE TRANS(DATA,NVAR,NSAM,IERROR,NVG)
DOUBLE PRECISION TODE,D123
DIMENSION COV(20,20),C123(250,20),Z123(20,20)
DIMENSION DATA(250,20)
COMMON XDATA(250,20),COV,C123,Z123
ASN(XX)=ATAN(XX /SQRT(1.0-XX**2))
DATA D123/6HTRNGEN/
MARY=0
FN=NSAM
WRITE (6,1403)
WRITE (6,1400)
DO1000 J=1,NVG
READ (5,1100)TODE,NEWA,LCODE,LVA,BNEW
IF (D123 .EQ. TODE) GO TO 201
200 NVAR=-NVAR
GO TO 1111
201 WRITE (6,1402)J,NEWA,LCODE,LVA,BNEW
IF(LCODE*(15-LCODE)) 710,710,714
710 WRITE (6,712)
714 IF(LCODE-10)4,5,5
5 NEWB=BNEW
4 DO 35 I=1,NSAM
D=DATA(I,LVA)
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
GO TO 3
8 D2=SQRT(D)
GO TO 3
20 IF(D)99,11,12
11 D2=1.0
GO TO 3
12 D2=SQRT(D)+SQRT(D+1.0)
GO TO 3
30 IF(D)99,99,14
14 D2=ALOG10(D)
GO TO 3
40 D2=EXP(D)
GO TO 3
50 IF(D)99, 7,17
17 IF(D-1.0)18,19,99
19 D2=3.14159265/2.0
GO TO 3
18 A=SQRT(D)
DATA(I,NEWA)=ASN(A)
GO TO 3
60 A=D/(FN+1.0)
B=A+1.0/(FN+1.0)
IF(A)99,23,24
23 IF(B)99, 7,27
27 D2=ASN(SQRT(B))
GO TO 3
24 IF(B)99,28,29
28 D2=ASN(SQRT(A))
GO TO 3
29 A=SQRT(A)
B=SQRT(B)
D2=ASN(A)+ASN(B)
GO TO 3
70 IF(D)31,99,31
31 D2=1.0/D
GO TO 3
80 D2=D+BNEW
GO TO 3
90 D2=D*BNEW
GO TO 3
100 IF(D)99, 7,33
33 D2=D**NEWB
GO TO 3
110 D2=D+DATA(I,NEWB)
GO TO 3
120 D2=D-DATA(I,NEWB)
GO TO 3
130 D2=D*DATA(I,NEWB)
GO TO 3
140 IF(DATA(I,NEWB))34,99,34
34 D2=D/DATA(I,NEWB)
GO TO 3
99 IF(MARY)43,44,44
44 MARY=-999
IERROR=-999
WRITE (6,1404)J
43 WRITE (6,1405)I
GO TO 35
3 DATA(I,NEWA)=D2
35 CONTINUE
IF(IERROR)42,1000,1000
1000 CONTINUE
712 FORMAT(29HERROR ON TRANSGENERATION CODE)
1100 FORMAT(A6,I3,I2,I3,F6.0)
1400 FORMAT(46H0CARD K TRANS ORIG. ORIG. VAR(J)/45H NO.
1VARIABLE CODE VAR(I) OR CONSTANT)
1401 FORMAT(41H0PROGRAM CANNOT CONTINUE FOR THIS PROBLEM)
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)
IF(IERROR)42,1111,1111
42 WRITE (6,1401)
1111 RETURN
END