Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/bmd/bmdx72.for
There is 1 other file named bmdx72.for in the archive. Click here to see a list.
00100 C
00200 C***********************************************************************
00300 C
00400 C PROJECTED REVISIONS IN THE FORTRAN IV (7094) VERSION OF BMD72X
00500 C TO ALLOW ITS EXECUTION ON THE IBM SYSTEM/360. AUGUST 9, 1966.
00600 C
00700 C***********************************************************************
00800 C
00900 C FACTOR ANALYSIS - MAIN PROGRAM MARCH 28, 1966
01000 C
01100 DIMENSION DATE(5)
01200 DOUBLE PRECISION FINISH
01300 DOUBLE PRECISION P,PC,PROBLM
01400 DATA DATE/'APRI','L 25',', 19','69 ',' '/
01500 DOUBLE PRECISION CUNST,GAMO
01600 DATA MAXST/6700/
01700 DIMENSION A(6700)
01800 DATA BINARY,YES,PROBLM,FINISH,ONO/4HBNRY,3HYES,6HPROBLM,6HFINISH,
01900 .2HNO/
02000 C EXTERNAL SIGN STATEMENT WAS REMOVED FROM HERE
02100 LOGICAL FSC
02200 COMMON PC,NC,NV,MAT,ICOM,IROT,GAMA,NROT,CONST,COR,VEC,PFS,WFS,MT,N
02300 XF,RE,MFT,NFO,REO,AABBCC,FSC,MXC,MXR,AMIS,NVV
02400 C
02500 C
02600 AMIS=ONO
02700 C THIS VARIABLE (AMIS) IS A DUMMY VARIABLE AT THIS TIME. IT WAS
02800 C INTENDED THAT AT SOME FUTURE TIME IT WOULD BE USED IN
02900 C ALLOWING FOR MISSING VARIABLES. THIS HAS NOT YET BEEN DONE.
03000 C
03100 C
03200 CALL USAGEB('BMDX72')
03300 55 REWIND 1
03400 REWIND 2
03500 REWIND 3
03600 81 FORMAT(37H1BMDX72 - FACTOR ANALYSIS - REVISED ,5A4
03700 X/41H0HEALTH SCIENCES COMPUTING FACILITY, UCLA)
03800 READ(5,2) P,PC,NC,NV,MAT,ICOM,MXC,IROT,MXR,GAMO,NROT,CUNST,AKN,
03900 1COR,VEC,PFS,WFS,MT,NF,RE,MFT,NFO,REO,UPL
04000 *,NLF,KL,KRL
04100 IF(P.EQ.FINISH) STOP
04200 IF(P.NE.PROBLM) GO TO 121
04300 C PRIOR TO MAY 1,1968 THIS PROGRAM
04400 C PRINTED THE PROBLM CARD AS IT WAS READ IN
04500 C UNDER A-FORMAT, THEN USED 'USEBUF' TO
04600 C REREAD THIS CARD. THIS A-FORMAT INFORMATION
04700 C WAS STORED IN A REAL * 8 ' FF(14) ' ARRAY
04800 C USING ' FORMAT (13A6,A2) ' TYPE READ. THIS
04900 C SEEMED UNNECESSARY AND WAS DROPPED.
05000 2 FORMAT(2A6,I6,I3,2I1,I2,I1,I2,A6,I3,A6,A2,4A3,2(2I2,A2),F3.3,3I2)
05100 KN=1
05200 IF(AKN.EQ.ONO) KN=0
05300 WRITE(6,81) DATE
05400 IF(UPL.EQ.0.0) UPL=.95
05500 IF(ICOM.EQ.0) ICOM=1
05600 MAT=4-MAT
05700 IF(MAT.LE.0) MAT=5-MAT
05800 IF (NROT.EQ.0) NROT=NV*0.5
05900 IF(MT.EQ.0) MT=5
06000 IF(NFO.EQ.0 .AND. WFS.EQ.YES) NFO=1
06100 IF(NF.EQ.0) NF=1
06200 CONST=1.0
06300 CALL ATOF (CUNST,6,CONST)
06400 GAMA=100.00
06500 CALL ATOF (GAMO,6,GAMA)
06600 IF(MXC.EQ.0) MXC=1
06700 IF(MXR.EQ.0) MXR=50
06800 NF=18*NF
06900 NFO=18*NFO
07000 IF(NF.LT.0) NF=-1
07100 IF(NFO.LT.0) NFO=-1
07200 L00=1+IABS(NF)
07300 L0=L00+IABS(NFO)
07400 L8=L0+NV
07500 L9=L8+NV
07600 M3=L9+NV
07700 L1=M3+NV
07800 L2=L1+NV
07900 L3=L2+NV
08000 L4=L3+NV
08100 L5=L4+NV
08200 L6=L5+NV
08300 L7=MAX0(L6+NV,L4+255)
08400 LL=L7+NV*NROT
08500 NNR=0
08600 IF(IROT.GT.1) NNR=NROT*NROT
08700 NVV=(NV*(NV+1))/2
08800 NLF=NLF*20
08900 K8=L7+NVV
09000 K9=K8+NLF
09100 K91=K9-1
09200 IF(K9.GT.MAXST) GO TO 791
09300 LA=1
09400 A(1)=BINARY
09500 A(L00)=BINARY
09600 FSC=.FALSE.
09700 IF(PFS.EQ.YES .OR. WFS.EQ.YES) FSC=.TRUE.
09800 CALL DOIT1(A,A(L00),A(L0),A(L8),A(L9),A(L7),A(L4),UPL)
09900 IF(NLF.NE.0) READ (5,20) (A(I),I=K8,K91)
10000 20 FORMAT(20A4)
10100 17 CALL DOIT2(A,A(L4),A(L00),A(M3),A(L0),A(L8),A(L9),A(L7),A(L7),NV,
10200 XLA,LAL)
10300 GO TO (11,11,12,13,57),LA
10400 11 CALL MULTR(A(L7),A(M3),A(L0),A(L1),NV,A(L4))
10500 GO TO 17
10600 12 A(L0)=AABBCC
10700 CALL RAV(A(L7),A(M3),NV,A(L0),A(L1),A(L2),A(L3),A(L4),A(L5),A(L6),
10800 X.00001,CONST,NROT,MXC)
10900 LL=L7+NV*NROT
11000 NNR=0
11100 IF(IROT.GT.1) NNR=NROT*NROT
11200 IF(LL.GT.MAXST) GO TO 371
11300 GO TO 17
11400 13 IF(KL.NE.0) CALL LOUT(A(L7),NV,NROT,KL,A(K8))
11500 IF(IROT.EQ.0) GO TO 57
11600 IF(LL+NNR.GT.MAXST) GO TO 371
11700 GO TO (8,9,10),IROT
11800 8 CALL ROTAT1(A(L7),A(LL),NV,NROT,GAMA,.00001,MXR,1,KN,A(L0),A(L1)
11900 X,A(L2),A(L4),A(L3),A(L5),A(L6))
12000 IF(KRL.NE.0) CALL LOUT(A(L7),NV,NROT,KRL,A(K8))
12100 GO TO 17
12200 9 CALL ROTAT2(A(L7),A(LL),NV,NROT,GAMA,.00001,MXR,1,KN,A(L0),A(L1)
12300 X,A(L2),A(L4),A(L3),A(L5),A(L6),UPL)
12400 IF(KRL.NE.0) CALL LOUT(A(L7),NV,NROT,KRL,A(K8))
12500 GO TO 17
12600 10 CALL ROTAT3(A(L7),A(LL),NV,NROT,GAMA,.00001,MXR,1,KN,A(L0),A(L1)
12700 X,A(L2),A(L4),A(L3),A(L5),A(L6))
12800 IF(KRL.NE.0) CALL LOUT(A(L7),NV,NROT,KRL,A(K8))
12900 GO TO 17
13000 121 WRITE (6,122)
13100 122 FORMAT(28H0UNRECOGNIZABLE PROBLEM CARD)
13200 STOP
13300 371 WRITE (6,372)
13400 372 FORMAT(47H0INSUFFICIENT STORAGE IS AVAILABLE FOR ROTATION)
13500 STOP
13600 791 WRITE (6,792)
13700 792 FORMAT(36H0THIS PROBLEM HAS TOO MANY VARIABLES)
13800 STOP
13900 57 IF(.NOT.FSC) GO TO 55
14000 502 FORMAT(1H0,'HENCE FACTOR SCORES ARE NOT COMPUTED OR PRINTED, ALTHO
14100 XUGH REQUESTED')
14200 IF(MAT.EQ.5) WRITE(6,500)
14300 500 FORMAT(1H0,'A CORRELATION OR COVARIANCE MATRIX IS READ AS DATA')
14400 IF(MAT.EQ.6) WRITE(6,501)
14500 501 FORMAT(1H0,'A FACTOR LOADING MATRIX IS READ AS DATA ')
14600 IF(MAT.GE.5) WRITE(6,502)
14700 LL7=LL+NNR
14800 IF(MAT.GE.5) GO TO 55
14900 IF(LL7+NV*NROT.GT.MAXST) GO TO 356
15000 CALL FACSCO(A(L00),A(L7),A(LL),A(LL7),A(L1),A(L2),A(L3),A(L8),A(L9
15100 X),A(L4),NV,NROT)
15200 GO TO 55
15300 356 WRITE (6,558)
15400 558 FORMAT(67H0INSUFFICIENT STORAGE IS AVAILABLE FOR COMPUTATION OF FA
15500 XCTOR SCORES)
15600 GO TO 55
15700 END
15800 SUBROUTINE LOUT(A,NV,NR,KL,F)
15900 DIMENSION A(NV,NR),F(16)
16000 DO 1 I=1,NV
16100 1 WRITE(KL,F) I,(A(I,J),J=1,NR)
16200 RETURN
16300 END
16400 C FUNCTION NBL FOR BMDX72 MARCH 28, 1966
16500 LOGICAL FUNCTION NBL(X)
16600 EXTERNAL SIGN
16700 NBL=.TRUE.
16800 IF (X.EQ.0.0.AND.SIGN(1.0,X).LT.0.0) NBL=.FALSE.
16900 RETURN
17000 END
17100 C SUBROUTINE FACSCO FOR BMDX72 MARCH 28, 1966
17200 SUBROUTINE FACSCO(FO,A,T,H,X,Y,Z,S1,D1,C,NV,NROT)
17300 DATA ST,BL,YES/1H*,1H ,3HYES/
17400 DOUBLE PRECISION PC
17500 DIMENSION A(NV,NROT),T(NROT,NROT),H(NV,NROT),S1(2),D1(2),X(2),Y(2)
17600 X,Z(2),C(255),FO(180)
17700 COMMON PC,NC,DM,MAT,ICOM,IROT,GAMA,DUM2,CONST,COR,VEC,PFS,WFS,MT,N
17800 XF,RE,MFT,NFO,REO,AABBCC,FSC,MXC,MXR,AMIS,NVV
17900 LOGICAL NBL
18000 IF(IROT.LT.2) GO TO 37
18100 DO 1 I=1,NV
18200 DO 2 K=1,NROT
18300 X(K)=A(I,K)
18400 2 A(I,K)=0.
18500 DO 1 J=1,NROT
18600 DO 1 K=1,NROT
18700 1 A(I,J)=A(I,J)+T(J,K)*X(K)
18800 37 M=255
18900 DO 3 I=1,NV
19000 DO 21 K=1,NV
19100 IF(M.LT.255) GO TO 20
19200 M=0
19300 READ (1) C
19400 200 FORMAT (20A4)
19500 20 M=M+1
19600 21 X(K)=-C(M)
19700 DO 3 J=1,NROT
19800 H(I,J)=0.
19900 DO 3 K=1,NV
20000 3 H(I,J)=H(I,J)+X(K)*A(K,J)
20100 REWIND 1
20200 WRITE (6,300)
20300 300 FORMAT(26H1FACTOR SCORE COEFFICIENTS)
20400 L1=0
20500 301 L0=L1+1
20600 L1=MIN0(L1+10,NROT)
20700 WRITE (6,302) (L,L=L0,L1)
20800 302 FORMAT(//20X,6HFACTOR,/2X,10I12)
20900 WRITE (6,303)
21000 303 FORMAT(9H VARIABLE)
21100 DO 304 I=1,NV
21200 304 WRITE (6,305) I,(H(I,L),L=L0,L1)
21300 305 FORMAT(I5,10F12.5)
21400 IF(L1.NE.NROT) GO TO 301
21500 M=255
21600 IF(PFS.NE.YES) GO TO 100
21700 IF(AMIS.EQ.YES) WRITE (6,101)
21800 IF(AMIS.NE.YES) WRITE (6,102)
21900 101 FORMAT(1H1,10X,47HFACTOR SCORES (* COMPUTED FROM INCOMPLETE DATA)/
22000 X5H CASE)
22100 102 FORMAT(1H1,10X,13HFACTOR SCORES/5H CASE)
22200 100 DO 71 L=1,NC
22300 HH=BL
22400 DO 5 I=1,NV
22500 IF(M.LT.255) GO TO 6
22600 M=0
22700 READ (2) C
22800 6 M=M+1
22900 Z(I)=C(M)
23000 IF (MAT.EQ.4) X(I)=(C(M)-S1(I))/D1(I)
23100 IF (MAT.EQ.3) X(I)=C(M)/D1(I)
23200 IF (MAT.EQ.2) X(I)=C(M)-S1(I)
23300 IF (MAT.EQ.1) X(I)=C(M)
23400 4 IF(AMIS.NE.YES .OR. NBL(C(M))) GO TO 5
23500 X(I)=0.
23600 HH=ST
23700 5 CONTINUE
23800 DO 17 I=1,NROT
23900 Y(I)=0.
24000 DO 17 J=1,NV
24100 17 Y(I)=Y(I)+H(J,I)*X(J)
24200 IF(WFS.NE.YES) GO TO 71
24300 IF(NFO.GT.0) WRITE (MFT,FO) (Z(I),I=1,NV),(Y(I),I=1,NROT)
24400 GO TO 1000
24500 1000 IF(NFO.LT.0) WRITE (MFT) (Z(I),I=1,NV),(Y(I),I=1,NROT)
24600 GO TO 71
24700 71 IF(PFS.EQ.YES) WRITE (6,75) L,HH,(Y(J),J=1,NROT)
24800 75 FORMAT(1X,I4,A1,10F12.5/(6X,10F12.5))
24900 RETURN
25000 END
25100 C SUBROUTINE DOIT1 FOR BMDX72 MARCH 28, 1966
25200 SUBROUTINE DOIT1(F,FO,X,S1,D1,A,C,UPL)
25300 DIMENSION C(255)
25400 DIMENSION F(180),FO(180),X(2),S1(2),D1(2),A(2)
25500 DOUBLE PRECISION PC
25600 COMMON PC,NC,NV,MAT,ICOM,IROT,GAMA,NROT,CONST,COR,VEC,PFS,WFS,MT,N
25700 XF,RE,MFT,NFO,REO,AABBCC,FSC,MXC,MXR,AMIS,NVV
25800 LOGICAL NBL
25900 DATA ONO,YES/2HNO,3HYES/
26000 LOGICAL FSC
26100 IF(GAMA.NE.100.00) GO TO 70
26200 GAMA=0.0
26300 IF(IROT.EQ.2) GAMA=0.5
26400 IF(IROT.EQ.1) GAMA=1.0
26500 70 MOT=MOD(MAT,2)
26600 IF(RE.NE.ONO .AND. MT.NE.5) REWIND MT
26700 IF(MFT.EQ.0) REO=ONO
26800 IF(REO.NE.ONO .AND. MFT.NE.6) REWIND MFT
26900 IF(NF.GT.0) READ (5,87) (F(I),I=1,NF)
27000 IF(NFO.GT.0) READ (5,87) (FO(I),I=1,NFO)
27100 87 FORMAT (18A4)
27200 WRITE (6,50) PC,NC,NV,MT,MXC,MXR,NROT,CONST
27300 50 FORMAT(40H0PROBLEM CODE
27400 XA6/40H0NUMBER OF CASES
27500 XI6/40H0NUMBER OF VARIABLES
27600 XI6/40H0INPUT TAPE
27700 XI6,/,'0MAX. ITERATIONS FOR COMMUNALITIES',8X,I3,/,'0MAX. ITERATION
27800 XS FOR ROTATION',13X,I3,/,'0NUMBER OF FACTORS TO BE ROTATED',10X,I4
27900 X,/,'0CONSTANT',31X,F10.6)
28000 WRITE(6,2000) UPL
28100 2000 FORMAT(40H0UPPER LIMIT ON CORRELATION COEFFICIENT ,2X,F6.5)
28200 WRITE (6,51) (F(I),I=1,NF)
28300 51 FORMAT(40H0INPUT FORMAT
28400 X18A4/(40X,18A4))
28500 IF(AMIS.EQ.YES) WRITE (6,500)
28600 500 FORMAT(35H0BLANKS ARE TREATED AS MISSING DATA)
28700 IF(WFS.EQ.YES) WRITE (6,52) MFT,(FO(I),I=1,NFO)
28800 52 FORMAT(40H0OUTPUT TAPE
28900 XI6/40H0OUTPUT FORMAT
29000 X18A4/(40X,18A4))
29100 IF(MAT.EQ.1) WRITE (6,601)
29200 IF(MAT.EQ.2) WRITE (6,602)
29300 IF(MAT.EQ.3) WRITE (6,603)
29400 IF(MAT.EQ.4) WRITE (6,604)
29500 IF(MAT.EQ.5) WRITE (6,605)
29600 IF(MAT.EQ.6) WRITE (6,606)
29700 IF(MAT.GE.6) GO TO 666
29800 601 FORMAT(49H0THE COVARIANCE MATRIX ABOUT THE ORIGIN IS FORMED)
29900 602 FORMAT(32H0THE COVARIANCE MATRIX IS FORMED)
30000 603 FORMAT(50H0THE CORRELATION MATRIX ABOUT THE ORIGIN IS FORMED)
30100 604 FORMAT(33H0THE CORRELATION MATRIX IS FORMED)
30200 605 FORMAT(37H0THE MATRIX TO BE FACTORED IS READ IN)
30300 606 FORMAT(27H0A FACTOR MATRIX IS READ IN)
30400 IF(ICOM.EQ.1) WRITE (6,701)
30500 IF(ICOM.EQ.2) WRITE (6,702)
30600 IF(ICOM.EQ.3) WRITE (6,703)
30700 IF(ICOM.EQ.4) WRITE (6,704)
30800 701 FORMAT(32H0DIAGONAL ELEMENTS ARE UNALTERED)
30900 702 FORMAT(64H0INITIAL COMMUNALITY ESTIMATES ARE SQUARED MULTIPLE CORR
31000 XELATIONS)
31100 703 FORMAT(63H0INITIAL COMMUNALITY ESTIMATES ARE LARGEST ABSOLUTE ROW
31200 X VALUES)
31300 704 FORMAT(42H0INITIAL COMMUNALITY ESTIMATES ARE READ IN)
31400 666 IF(IROT.EQ.0) WRITE (6,800)
31500 IF(IROT.EQ.1 .AND. GAMA.EQ.0.0) WRITE (6,811)
31600 IF(IROT.EQ.1 .AND. GAMA.EQ.1.0) WRITE (6,821)
31700 IF(IROT.EQ.1 .AND.GAMA.NE.1.0 .AND. GAMA.NE.0.0) WRITE (6,831)
31800 XGAMA
31900 IF(IROT.EQ.2) WRITE (6,802) GAMA
32000 IF(IROT.EQ.3) WRITE (6,803) GAMA
32100 800 FORMAT(25H0NO ROTATION IS PERFORMED)
32200 811 FORMAT(32H0QUARTIMAX ROTATION IS PERFORMED)
32300 821 FORMAT(30H0VARIMAX ROTATION IS PERFORMED)
32400 831 FORMAT(53H0ORTHOGONAL ROTATION IS PERFORMED WITH GAMMA EQUAL TOF7.
32500 X4)
32600 801 FORMAT(33H0ORTHOGONAL ROTATION IS PERFORMED)
32700 802 FORMAT(50H0OBLIMIN ROTATION IS PERFORMED WITH GAMMA EQUAL TOF7.4)
32800 803 FORMAT(70H0OBLIQUE ROTATION FOR SIMPLE LOADINGS IS PERFORMED WITH
32900 XGAMMA EQUAL TOF7.4)
33000 IF(MAT.NE.5) GO TO 20
33100 L1=0
33200 DO 21 J=1,NV
33300 L0=L1+1
33400 L1=L0+NV-J
33500 IF(NF.GT.0) READ (MT,F) (X(I),I=1,NV)
33600 1000 IF(NF.LE.0) READ (MT) (X(I),I=1,NV)
33700 1001 I=J
33800 DO 21 L=L0,L1
33900 A(L)=X(I)
34000 21 I=I+1
34100 RETURN
34200 20 IF(MAT.EQ.6) RETURN
34300 22 L=0
34400 DO 1 I=1,NV
34500 D1(I)=0.
34600 S1(I)=0.
34700 DO 1 J=I,NV
34800 L=L+1
34900 1 A(L)=0.
35000 ANC=NC
35100 IF(AMIS.EQ.YES) RETURN
35200 30 M=0
35300 DO 4 LL=1,NC
35400 IF(NF.GT.0) READ (MT,F) (X(I),I=1,NV)
35500 GO TO 1004
35600 1004 IF(NF.LE.0) READ (MT) (X(I),I=1,NV)
35700 GO TO 1005
35800 1005 IF(.NOT.FSC) GO TO 35
35900 DO 36 I=1,NV
36000 IF(M.NE.255) GO TO 37
36100 M=0
36200 WRITE(2) C
36300 200 FORMAT (20A4)
36400 37 M=M+1
36500 36 C(M)=X(I)
36600 35 H=LL
36700 DO 3 I=1,NV
36800 D1(I)=(X(I)-S1(I))/H
36900 S1(I)=S1(I)+D1(I)
37000 3 IF(MOT.EQ.1) D1(I)=X(I)
37100 H=H*(H-1.)
37200 IF(MOT.EQ.1) H=1.
37300 L=0
37400 DO 4 I=1,NV
37500 DD=D1(I)*H
37600 DO 4 J=I,NV
37700 L=L+1
37800 4 A(L)=A(L)+DD*D1(J)
37900 IF(.NOT.FSC) GO TO 10
38000 WRITE (2) C
38100 ENDFILE 2
38200 REWIND 2
38300 10 L=1
38400 H=ANC
38500 M=255
38600 DO 60 I=1,NV
38700 AL=A(L)
38800 IF(MOT.EQ.1) AL=A(L)-H*S1(I)*S1(I)
38900 D12=SQRT(AL/(H-1.))
39000 LL=L
39100 DO 220 J=I,NV
39200 7 GO TO (400,6,400,6),MAT
39300 400 A(L)=A(L)/H
39400 GO TO 220
39500 6 A(L)=A(L)/(H-1.)
39600 220 L=L+1
39700 X(I)=A(LL)
39800 60 D1(I)=D12
39900 WRITE (6,80) (I,S1(I),D1(I),I=1,NV)
40000 80 FORMAT('1VARIABLE MEAN ST.DEV.'//(I6,2F15.6))
40100 DO 61 I=1,NV
40200 61 D1(I)=SQRT(X(I))
40300 IF(MAT.LT.3) RETURN
40400 L=0
40500 DO 11 I=1,NV
40600 DO 11 J=I,NV
40700 L=L+1
40800 11 A(L)=A(L)/(D1(I)*D1(J))
40900 RETURN
41000 END
41100 C SUBROUTINE DOIT2 FOR BMDX72 MARCH 28, 1966
41200 SUBROUTINE DOIT2(F,H,FO,R,X,S1,D1,A,AB,NV,LA,LAL)
41300 DOUBLE PRECISION COMM,COM
41400 DIMENSION FO(180),R(180),X(2),S1(2),D1(2),A(2),AB(NV,2)
41500 X,F(2),H(2)
41600 DOUBLE PRECISION PC
41700 COMMON PC,NC,DUMMY,MAT,ICOM,IROT,GAMA,NROT,CONST,COR,VEC,PFS,WFS,
41800 XMT,NF,RE,MFT,NFO,REO,AABBCC,FSC,MXC,MXR,AMIS,NVV
41900 LOGICAL FSC
42000 DATA COMM,YES/6HCOMMUN,3HYES/
42100 GO TO (571,31,573,574),LA
42200 571 IF(MAT.NE.6) GO TO 771
42300 DO 723 I=1,NV
42400 IF(NF.GT.0) READ (MT,F) (AB(I,J),J=1,NROT)
42500 723 IF(NF.LE.0) READ (MT) (AB(I,J),J=1,NROT)
42600 GO TO 100
42700 771 L=1
42800 K=NV
42900 AABBCC=0.
43000 DO 50 I=1,NV
43100 AABBCC=AABBCC+A(L)
43200 L=L+K
43300 50 K=K-1
43400 IF(.NOT.(FSC.OR.ICOM.EQ.2.OR.MXC.GT.1)) GO TO 310
43500 C
43600 C
43700 C
43800 WRITE (1) (A(L),L=1,NVV)
43900 200 FORMAT (20A4)
44000 LA=2
44100 IF(FSC.OR.ICOM.EQ.2) RETURN
44200 31 ENDFILE 1
44300 REWIND 1
44400 IF (FSC.OR.ICOM.EQ.2) GO TO 1798
44500 310 GO TO (36,42,32,33),ICOM
44600 32 DO 40 I=1,NV
44700 40 R(I)=0.
44800 L=0
44900 NV1=NV-1
45000 DO 41 I=1,NV1
45100 L=L+1
45200 I1=I+1
45300 DO 41 J=I1,NV
45400 L=L+1
45500 R(I)=AMAX1(ABS(A(L)),R(I))
45600 41 R(J)=AMAX1(ABS(A(L)),R(J))
45700 GO TO 42
45800 1798 CONTINUE
45900 C
46000 C
46100 READ (1) (A(L),L=1,NVV)
46200 C
46300 GO TO 310
46400 33 L1=0
46500 38 L0=L1+1
46600 L1=MIN0(L1+11,NV)
46700 37 FORMAT(A6,11F6.6)
46800 READ (5,37) COM,(R(L),L=L0,L1)
46900 IF(COM.NE.COMM) GO TO 89
47000 IF(L1.LT.NV) GO TO 38
47100 GO TO 42
47200 89 WRITE (6,98) COM,(R(L),L=L0,L1)
47300 98 FORMAT(42H0UNRECOGNIZABLE COMMUNALITY ESTIMATES CARD/1X,A6,11F7.6)
47400 STOP
47500 42 L=1
47600 K=NV
47700 DO 34 I=1,NV
47800 A(L)=R(I)
47900 L=L+K
48000 34 K=K-1
48100 36 IF(COR.NE.YES) GO TO 85
48200 IF(MAT.EQ.1) WRITE (6,831)
48300 IF(MAT.EQ.2) WRITE (6,832)
48400 IF(MAT.EQ.3) WRITE (6,833)
48500 IF(MAT.EQ.4) WRITE (6,834)
48600 831 FORMAT(35H1COVARIANCE MATRIX ABOUT THE ORIGIN)
48700 832 FORMAT(18H1COVARIANCE MATRIX)
48800 833 FORMAT(36H1CORRELATION MATRIX ABOUT THE ORIGIN)
48900 834 FORMAT(19H1CORRELATION MATRIX)
49000 L1=0
49100 81 L0=L1+1
49200 L1=MIN0(NV,L1+10)
49300 L11=L0
49400 K1=((L0-1)*(2*NV-L0+2))/2+1
49500 WRITE (6,82) (L,L=L0,L1)
49600 82 FORMAT(1H-,10I12)
49700 DO 84 I=L0,NV
49800 K=K1
49900 K1=K1+1
50000 K0=NV-L0
50100 DO 90 L=L0,L11
50200 X(L)=A(K)
50300 K=K+K0
50400 90 K0=K0-1
50500 WRITE (6,86) I,(X(L),L=L0,L11)
50600 86 FORMAT(I4,10F12.5)
50700 84 L11=MIN0(L11+1,L1)
50800 IF(L1.NE.NV) GO TO 81
50900 85 IF(MXC.GT.1) WRITE (6,928)
51000 928 FORMAT(28H-ITERATION FOR COMMUNALITIES/1H0,5X,14HMEAN ABS. DEV. 4X
51100 X,14HMAX. ABS. DEV.)
51200 LA=3
51300 RETURN
51400 573 ENDFILE 3
51500 REWIND 3
51600 IF (NROT.LE.0) GO TO 123
51700 DO 12 I=1,NROT
51800 H(I)=SQRT(H(I))
51900 12 READ (3) (AB(J,I),J=1,NV)
52000 REWIND 3
52100 123 DO 812 I=1,NV
52200 X(I)=0.
52300 DO 812 J=1,NROT
52400 AB(I,J)=AB(I,J)*H(J)
52500 812 X(I)=X(I)+AB(I,J)*AB(I,J)
52600 WRITE (6,813) (I,R(I),X(I),I=1,NV)
52700 813 FORMAT(37H-VARIABLE ESTIMATED FINAL/12X,28HCOMMUNALITY
52800 X COMMUNALITY/(I5,2F17.6))
52900 IF(VEC.NE.YES) GO TO 100
53000 LLL=1
53100 WRITE (6,61)
53200 61 FORMAT(30H1FACTOR MATRIX BEFORE ROTATION)
53300 68 L1=0
53400 62 L0=L1+1
53500 L1=MIN0(L1+10,NROT)
53600 WRITE (6,63) (L,L=L0,L1)
53700 63 FORMAT(//20X,6HFACTOR/2X,10I12)
53800 WRITE (6,633)
53900 633 FORMAT(1X,8HVARIABLE)
54000 DO 64 I=1,NV
54100 64 WRITE (6,65) I,(AB(I,J),J=L0,L1)
54200 IF(L1.NE.NROT) GO TO 62
54300 65 FORMAT(I5,10F12.5)
54400 LA=5
54500 IF(LLL.NE.1) RETURN
54600 100 LA=4
54700 IF(IROT.EQ.0) LA=5
54800 RETURN
54900 574 WRITE (6,69)
55000 69 FORMAT(22H1ROTATED FACTOR MATRIX)
55100 LLL=2
55200 GO TO 68
55300 END
55400 C SUBROUTINE MULTR FOR BMDX72 MARCH 28, 1966
55500 SUBROUTINE MULTR(A,U,V,W,N,C)
55600 DIMENSION A(2),U(2),V(2),W(2),C(255)
55700 LOGICAL W
55800 L=N
55900 M=1
56000 K=0
56100 DO 1 I=1,N
56200 V(I)=A(M)
56300 M=M+L
56400 L=L-1
56500 IF(V(I).NE.0.) K=I
56600 1 W(I)=.FALSE.
56700 IF (K.EQ.0) GO TO 50
56800 6 M=K-N
56900 L=N
57000 DO 2 I=1,N
57100 IF(I.GT.K) L=1
57200 M=M+L
57300 L=L-1
57400 U(I)=A(M)
57500 2 A(M)=0.
57600 P=U(K)
57700 W(K)=.TRUE.
57800 U(K)=-1.
57900 M=1
58000 T=0.
58100 K=0
58200 DO 5 I=1,N
58300 Y=-U(I)/P
58400 MM=M
58500 DO 4 J=I,N
58600 A(M)=A(M)+Y*U(J)
58700 4 M=M+1
58800 IF(W(I)) GO TO 5
58900 IF (V(I).EQ.0.0) GO TO 5
59000 H=A(MM)/V(I)
59100 IF(H.LT.T) GO TO 5
59200 T=H
59300 K=I
59400 5 CONTINUE
59500 IF(T.GT.1.E-5) GO TO 6
59600 50 L=N
59700 M=1
59800 DO 9 K=1,N
59900 IF(W(K)) GO TO 7
60000 U(K)=V(K)
60100 GO TO 11
60200 7 L1=N
60300 U(K)=V(K)+1./A(M)
60400 M1=1
60500 L2=N
60600 M2=K-N
60700 DO 8 I=1,N
60800 IF(I.GT.K) L2=1
60900 M2=M2+L2
61000 L2=L2-1
61100 IF(W(I)) GO TO 10
61200 IF(A(M1)-A(M2)*A(M2)/A(M) .GT. 1.E-5*V(I)) U(K)=V(K)
61300 10 M1=M1+L1
61400 L1=L1-1
61500 8 CONTINUE
61600 11 M=M+L
61700 9 L=L-1
61800 MM=0
61900 DO 12 K=1,N
62000 M=K-N
62100 L=N
62200 DO 12 I=1,N
62300 IF(I.GT.K) L=1
62400 M=M+L
62500 L=L-1
62600 IF(MM.NE.255) GO TO 20
62700 MM=0
62800 WRITE ( 1) C
62900 200 FORMAT (20A4)
63000 20 MM=MM+1
63100 C(MM)=A(M)
63200 12 IF(.NOT.(W(I).AND.W(K))) C(MM)=0.
63300 WRITE (1) C
63400 RETURN
63500 END
63600 C SUBROUTINE RAV FOR BMDX72 MARCH 28, 1966
63700 SUBROUTINE RAV(A,C,N,U,V,W,G,H,R,P,ACC,CONST,NNN,MXC)
63800 LOGICAL ITR
63900 DIMENSION A(N),C(N)
64000 DIMENSION U(N),V(N),W(N),G(N),H(N),R(N),P(N)
64100 DATA CON/"201400000400/
64200 EXTERNAL SIGN
64300 ITR=.TRUE.
64400 MX=0
64500 L=1
64600 M=N
64700 DO 167 I=1,N
64800 C(I)=A(L)
64900 L=L+M
65000 167 M=M-1
65100 AG=U(1)
65200 NNNN=NNN
65300 158 DO 98 I=1,N
65400 V(I)=0.
65500 98 U(I)=0.
65600 JJ=0
65700 MX=MX+1
65800 IF(MX.GE.MXC) ITR=.FALSE.
65900 IF(ITR) REWIND 1
66000 NM2=N-2
66100 DO 91 KK=1,NM2
66200 KP1=KK+1
66300 AN=0.
66400 II=JJ+1
66500 DO 92 I=KP1,N
66600 G(I)=0.
66700 II=II+1
66800 R(I)=A(II)+U(KK)*V(I)+U(I)*V(KK)
66900 92 AN=AN+R(I)*R(I)
67000 AN=SQRT(AN)
67100 A2=R(KP1)
67200 V(KK)=A(JJ+1)+2.*U(KK)*V(KK)
67300 X=SQRT(1.+ABS(A2)/AN)
67400 IF(X.GT.1.1) GO TO 45
67500 A2=-A2
67600 X=SQRT(1.-ABS(A2)/AN)
67700 45 U(KK)=-SIGN(AN,A2)
67800 Y=SIGN(1./(X*AN), A2)
67900 W(KP1)=X
68000 JJ=JJ+2
68100 A(JJ)=X
68200 KP2=KK+2
68300 DO 93 I=KP2,N
68400 JJ=JJ+1
68500 W(I)=R(I)*Y
68600 93 A(JJ)=W(I)
68700 DO 95 K=KP1,N
68800 II=II+1
68900 A(II)=A(II)+U(K)*V(K)*2.
69000 G(K)=G(K)+A(II)*W(K)
69100 IF(K.EQ.N) GO TO 99
69200 KPP1=K+1
69300 DO 95 I=KPP1,N
69400 II=II+1
69500 A(II)=A(II)+U(I)*V(K)+V(I)*U(K)
69600 G(I)=G(I)+A(II)*W(K)
69700 95 G(K)=G(K)+A(II)*W(I)
69800 99 UAU=0.
69900 DO 96 I=KP1,N
70000 UAU=UAU+G(I)*W(I)
70100 96 U(I)=W(I)
70200 DO 97 I=KP1,N
70300 97 V(I)=UAU*U(I)/2.-G(I)
70400 91 CONTINUE
70500 X=U(N)*V(N-1)+U(N-1)*V(N)
70600 Y=2.*U(N-1)*V(N-1)
70700 Z=2.*U(N)*V(N)
70800 U(N-1)=A(JJ+2)+X
70900 V(N-1)=A(JJ+1)+Y
71000 V(N)=A(JJ+3)+Z
71100 JJJ=JJ-2
71200 NM1=N-1
71300 DO 15 I=1,NM1
71400 U(I)=SIGN(AMAX1(1.E-9*ABS(V(I)),ABS(U(I))),U(I))
71500 15 W(I)=U(I)*U(I)
71600 D=AMAX1(ABS(V(1))+ABS(U(1)),ABS(V(N))+ABS(U(NM1)))
71700 DO 10 I=2,NM1
71800 10 D=AMAX1(ABS(U(I-1))+ABS(V(I))+ABS(U(I)),D)
71900 DO 20 I=1,N
72000 G(I)=-D
72100 20 H(I)=D
72200 DO 30 I=1,N
72300 70 X=(G(I)+H(I))/2.
72400 IF (X.EQ.G(I).OR.X.EQ.H(I)) GO TO 30
72500 K=0
72600 GG=V(1)-X
72700 IF(GG.GT.0.) K=K+1
72800 DO 40 J=2,N
72900 IF (GG.NE.0.0) GG=V(J)-X-W(J-1)/GG
73000 IF (GG.EQ.0.0) GG=V(J)-X
73100 40 IF(GG.GT.0.) K=K+1
73200 DO 50 J=I,N
73300 IF(K.LT.J) H(J)=AMIN1(H(J),X)
73400 50 IF(K.GE.J) G(J)=AMAX1(G(J),X)
73500 GO TO 70
73600 30 CONTINUE
73700 GO TO 207
73800 2077 WRITE (6,200) (H(I),I=1,N)
73900 200 FORMAT(13H1EIGENVALUES //(10F12.5))
74000 GG=0.
74100 DO 201 I=1,N
74200 GG=GG+H(I)
74300 IF(H(I).GT.0.) N1N=I
74400 IF (AG.NE.0.0) G(I)=GG/AG
74500 IF (AG.EQ.0.0) G(I)=0.0
74600 201 CONTINUE
74700 WRITE (6,202) (G(I),I=1,N1N)
74800 202 FORMAT(40H0CUMULATIVE PROPORTION OF TOTAL VARIANCE//(10F12.5))
74900 RETURN
75000 207 T=1.
75100 DO 203 I=1,NM1
75200 IF(H(I).NE.H(I+1) .OR. H(I).EQ.0.) GO TO 203
75300 T=0.
75400 WRITE(6,1) I
75500 1 FORMAT(' EIGEN VALUE OF I= 'I4,' AND I+1 EQUAL, CHECK ACCURACY ')
75600 H(I)=H(I)*CON
75700 203 CONTINUE
75800 IF(T.EQ.0.) GO TO 207
75900 L=1
76000 M=N
76100 DO 102 I=1,N
76200 A(L)=0.
76300 L=L+M
76400 102 M=M-1
76500 510 DO 107 II=1,NNNN
76600 IF(H(II).LE.CONST) GO TO 108
76700 DO 57 I=1,N
76800 57 R(I)=1.E-8
76900 LIP=0
77000 56 LIP=LIP+1
77100 DO 51 I=1,N
77200 G(I)=V(I)-H(II)
77300 P(I)=U(I-1)
77400 IF(ABS(R(I)).LT.1.E-17) R(I)=SIGN(1.E-17,R(I))
77500 51 W(I)=U(I)
77600 P(1)=0.
77700 W(N)=0.
77800 DO 52 I=1,NM1
77900 IF(ABS(G(I)).GT.ABS(P(I+1))) GO TO 53
78000 T=P(I+1)
78100 P(I+1)=G(I)
78200 G(I)=T
78300 T=G(I+1)
78400 G(I+1)=W(I)
78500 W(I)=T
78600 P(I)=W(I+1)
78700 W(I+1)=0.
78800 T=R(I+1)
78900 R(I+1)=R(I)
79000 R(I)=T
79100 53 IF(G(I).NE.0.0) D=P(I+1)/G(I)
79200 IF (G(I).EQ.0.0) D=0.0
79300 G(I+1)=G(I+1)-W(I)*D
79400 W(I+1)=W(I+1)-P(I)*D
79500 R(I+1)=R(I+1)-D*R(I)
79600 52 P(I+1)=0.
79700 L=N
79800 IF(G(L).EQ.0.0) G(L)=1.E-30
79900 DO 54 I=1,N
80000 IF (G(L).NE.0.0) R(L)=(R(L)-W(L)*R(L+1)-P(L)*R(L+2))/G(L)
80100 IF (G(L).EQ.0.0) R(L)=0.0
80200 IF(ABS(R(L)).LT.1.E8) GO TO 54
80300 DO 85 J=1,N
80400 IF(R(L).NE.0.0) R(J)=R(J)/R(L)
80500 85 IF(R(L).EQ.0.0) R(J)=0.0
80600 54 L=L-1
80700 T=0.
80800 DO 86 I=1,N
80900 86 T=T+R(I)*R(I)
81000 T=SQRT(T)
81100 IF(LIP.EQ.1) T=T*1.E8
81200 DO 83 I=1,N
81300 IF (T.NE.0.0) R(I)=R(I)/T
81400 IF (T.EQ.0.0) R(I)=0.0
81500 83 CONTINUE
81600 GO TO (56,55),LIP
81700 55 JJ=JJJ
81800 100 KQK=6
81900 DO 64 K=1,NM2
82000 NK=N-K
82100 T=0.
82200 DO 63 I=NK,N
82300 JJ=JJ+1
82400 G(I)=A(JJ)
82500 63 T=T+G(I)*R(I)
82600 JJ=JJ-KQK
82700 KQK=KQK+2
82800 DO 64 I=NK,N
82900 64 R(I)=R(I)-T*G(I)
83000 T=0.
83100 TM=0
83200 DO 65 I=1,N
83300 IF(ABS(R(I)).GT.ABS(TM)) TM=R(I)
83400 65 T=T+R(I)*R(I)
83500 T=SIGN(SQRT(T),TM)
83600 DO 66 I=1,N
83700 IF (T.NE.0.0) R(I)=R(I)/T
83800 IF (T.EQ.0.0) R(I)=0.0
83900 66 CONTINUE
84000 L=1
84100 M=N
84200 DO 101 I=1,N
84300 A(L)=A(L)+R(I)*R(I)*H(II)
84400 L=L+M
84500 101 M=M-1
84600 IF(.NOT.ITR) WRITE (3) (R(I),I=1,N)
84700 107 CONTINUE
84800 GO TO 109
84900 108 NNN=II-1
85000 109 IF (NNN.NE.0) GO TO 1099
85100 WRITE (6,1098) H(I)
85200 STOP
85300 1098 FORMAT (' PROGRAM TERMINATED BECAUSE THE MAX. EIGENVALUE IS ',
85400 1 F10.6/' OR THE FUNCTION SPECIFIED IS ZERO' )
85500 1099 CONTINUE
85600 X=0.
85700 Y=0.
85800 L=1
85900 M=N
86000 DO 110 I=1,N
86100 Z=ABS(A(L)-C(I))
86200 IF(ITR) C(I)=A(L)
86300 L=L+M
86400 M=M-1
86500 X=X+Z
86600 110 Y=AMAX1(Y,Z)
86700 IF (N.NE.0.0) X=X/FLOAT(N)
86800 IF (N.EQ.0.0) X=0.0
86900 IF(.NOT.ITR) GO TO 2077
87000 WRITE (6,916) MX,X,Y
87100 IF(Y.LT.0.001) ITR=.FALSE.
87200 916 FORMAT(I4,F11.4,F18.4)
87300 NVV=(N*(N+1))/2
87400 C
87500 C
87600 C
87700 READ (1) (A(L),L=1,NVV)
87800 300 FORMAT (20A4)
87900 C
88000 L=1
88100 M=N
88200 DO 152 I=1,N
88300 A(L)=C(I)
88400 L=L+M
88500 152 M=M-1
88600 GO TO 158
88700 END
88800 C SUBROUTINE ROTAT1 FOR BMDX72 MARCH 28, 1966
88900 SUBROUTINE ROTAT1(A,TT,LV,LR,GG,ACC,MR,INI,NOR,S,C,V,TL,FL,XL,YL)
89000 DIMENSION A(LV,LR),TT(LR,LR),C(2),S(2),V(2),TL(2),FL(2),XL(2),YL(2
89100 X)
89200 WRITE (6,882)
89300 882 FORMAT(30H1ORTHOGONAL ROTATION
89400 X/23H0ITERATION SIMPLICITY/13X,9HCRITERION)
89500 NV=LV
89600 NR=LR
89700 GA=GG/FLOAT(NV)
89800 IF(NOR.NE.1) GO TO 420
89900 DO 40 I=1,NV
90000 S(I)=0.
90100 DO 4 J=1,NR
90200 4 S(I)=S(I)+A(I,J)*A(I,J)
90300 FL(I)=SQRT(S(I))
90400 DO 2 J=1,NR
90500 2 A(I,J)=A(I,J)/FL(I)
90600 40 CONTINUE
90700 420 DO 43 I=1,NR
90800 41 C(I)=0.
90900 DO 43 J=1,NV
91000 43 C(I)=C(I)+A(J,I)*A(J,I)
91100 D=0.
91200 DO 44 J=1,NR
91300 44 D=D+C(J)
91400 F0=1.E20
91500 50 DO 32 L=1,MR
91600 G=0.
91700 H=0.
91800 DO 100 J=1,NR
91900 100 V(J)=0.
92000 DO 36 I=1,NV
92100 S(I)=0.
92200 DO 101 J=1,NR
92300 T=A(I,J)*A(I,J)
92400 S(I)=S(I)+T
92500 101 V(J)=V(J)+T*T
92600 36 H=H+S(I)*S(I)
92700 DO 102 J=1,NR
92800 102 G=G+V(J)-GA*C(J)*C(J)
92900 F=H-GA*D*D-G
93000 L1=L-1
93100 WRITE (6,78) L1,F
93200 78 FORMAT(I6,F16.6)
93300 IF(L1.EQ.0) FCC=ABS(F*ACC)
93400 IF(ABS(F-F0).LE.FCC) GO TO 38
93500 F0=F
93600 DO 32 IP=1,NR
93700 DO 32 IQ=1,IP
93800 IF(IP-IQ) 33,32,33
93900 33 Y=0.
94000 T=0.
94100 R=0.
94200 Z=0.
94300 DO 34 I=1,NV
94400 PQ=(A(I,IQ)+A(I,IP))*(A(I,IQ)-A(I,IP))
94500 AB=A(I,IP)*A(I,IQ)
94600 Z=Z+AB
94700 T=T+AB*AB
94800 Y=Y+AB*PQ
94900 34 R=R+PQ*PQ
95000 H=C(IQ)-C(IP)
95100 X=T-GA*Z*Z-.25*(R-GA*H*H)
95200 Y=Y-GA*Z*H
95300 PHI=ATAN2(-Y,-X)/4.
95400 S1=SIN(PHI)
95500 C1=COS(PHI)
95600 DO 35 I=1,NV
95700 X=C1*A(I,IP)+S1*A(I,IQ)
95800 A(I,IQ)=S1*A(I,IP)-C1*A(I,IQ)
95900 35 A(I,IP)=X
96000 Z=2.*S1*C1*Z
96100 S1=S1*S1
96200 C1=C1*C1
96300 T=C(IP)
96400 C(IP)=C1*T+S1*C(IQ)+Z
96500 C(IQ)=S1*T+C1*C(IQ)-Z
96600 32 CONTINUE
96700 38 IF(NOR) 88,88,711
96800 711 DO 46 I=1,NR
96900 DO 46 J=1,NV
97000 46 A(J,I)=A(J,I)*FL(J)
97100 88 RETURN
97200 END
97300 C SUBROUTINE PFC FOR BMDX72 MARCH 28, 1966
97400 SUBROUTINE PFC(A,N)
97500 DIMENSION A(N,N)
97600 WRITE (6,1)
97700 1 FORMAT(26H-FACTOR CORRELATION MATRIX)
97800 L1=0
97900 2 L0=L1+1
98000 L1=MIN0(N,L1+10)
98100 WRITE (6,3) (L,L=L0,L1)
98200 3 FORMAT(2H0 10I12)
98300 DO 4 I=1,N
98400 4 WRITE (6,5) I,(A(I,J),J=L0,L1)
98500 5 FORMAT(I5,10F12.5)
98600 IF(L1.NE.N) GO TO 2
98700 RETURN
98800 END
98900 SUBROUTINE ROTAT2 (A,TT,LV,LR,GG,ACC,MR,INI,NOR,S,C,V,TL,FL,XL,YL,
99000 XUPL)
99100 DIMENSION A(LV,LR),TT(LR,LR),C(2),S(2),V(2),TL(2),FL(2),XL(2),
99200 XYL(2)
99300 WRITE (6,882)
99400 882 FORMAT(30H1OBLIMIN ROTATION
99500 X/23H0ITERATION SIMPLICITY/13X,9HCRITERION)
99600 NR=LR
99700 NV=LV
99800 GA=GG/FLOAT(NV)
99900 DO 30 I=1,NV
00100 S(I)=0.0
00200 DO 4 J=1,NR
00300 4 S(I)=S(I)+A(I,J)*A(I,J)
00400 IF(NOR) 30,30,32
00500 32 FL(I)=SQRT(S(I))
00600 S(I)=1.
00700 DO 2 J=1,NR
00800 2 A(I,J)=A(I,J)/FL(I)
00900 30 CONTINUE
01000 D=0.
01100 DO 34 I=1,NR
01200 IF(INI) 31,31,3
01300 3 DO 6 J=1,NR
01400 6 TT(I,J)=0.
01500 TT(I,I)=1.
01600 31 C(I)=0.
01700 DO 33 J=1,NV
01800 33 C(I)=C(I)+A(J,I)*A(J,I)
01900 34 D=D+C(I)
02000 60 IF(INI) 38,38,39
02100 38 CALL MATINV(TT,NR,TOL,TL,XL,YL)
02200 DO 40 I=1,NR
02300 TL(I)=SQRT(TT(I,I))
02400 DO 40 J=1,I
02500 TT(I,J)=TT(I,J)/(TL(I)*TL(J))
02600 40 TT(J,I)=TT(I,J)
02700 F0=1.E20
02800 39 DO 11 L=1,MR
02900 F=0.
03000 DO 200 J=1,NR
03100 F=F-GA*C(J)*C(J)
03200 DO 200 I=1,NV
03300 T=A(I,J)*A(I,J)
03400 200 F=F+T*T
03500 F=-F-GA*D*D
03600 DO 201 I=1,NV
03700 201 F=F+S(I)*S(I)
03800 L1=L-1
03900 IF(L1.EQ.0) FCC=ABS(F*ACC)
04000 WRITE (6,52) L1,F
04100 52 FORMAT(I6,F16.6)
04200 IF(ABS(F-F0).LE.FCC) GO TO 17
04300 F0=F
04400 DO 11 IP=1,NR
04500 DO 11 IQ=1,NR
04600 IF(IP-IQ) 18,11,18
04700 45 A1=1.
04800 A2=0.
04900 GO TO 46
05000 18 X=0.
05100 Y=0.
05200 Z=0.
05300 D=D-C(IP)
05400 DO 12 I=1,NV
05500 AA=A(I,IP)*A(I,IP)
05600 S(I)=S(I)-AA
05700 W=S(I)-GA*D
05800 X=X+AA*W
05900 Y=Y+A(I,IP)*A(I,IQ)*W
06000 12 Z=Z+A(I,IQ)*A(I,IQ)*W
06100 O=TT(IP,IQ)
06200 A1=1.-O*O
06300 B1=Y*O-.5*(X+Z)
06400 C1=X*Z-Y*Y
06500 AL=B1*B1-A1*C1
06600 IF(AL) 45,49,49
06700 49 AL=-(B1+SQRT(AL))/A1
06800 A1=AL-Z
06900 A2=Y-AL*O
07000 IF(A1) 23,24,23
07100 24 A1=A2
07200 A2=AL-X
07300 23 W=A1*A1+A2*(A2+2.*A1*O)
07400 IF(ABS(A1).LE.ABS(.001*A2) .OR. W.LE.0.) GO TO 45
07500 W=SQRT(W)
07600 A1=A1/W
07700 A2=A2/W
07800 DO 14 I=1,NR
07900 IF(I.EQ.IP) GO TO 14
08000 V(I)=A1*TT(I,IP)+A2*TT(I,IQ)
08100 IF(ABS(V(I)).GT.UPL) GO TO 45
08200 14 CONTINUE
08300 DO 210 I=1,NR
08400 TT(I,IP)=V(I)
08500 210 TT(IP,I)=V(I)
08600 46 C(IP)=0.
08700 DO 13 I=1,NV
08800 A(I,IP)=A1*A(I,IP)+A2*A(I,IQ)
08900 AA=A(I,IP)*A(I,IP)
09000 S(I)=S(I)+AA
09100 13 C(IP)=C(IP)+AA
09200 D=D+C(IP)
09300 TT(IP,IP)=1.
09400 11 CONTINUE
09500 17 CALL MATINV(TT,NR,TOL,TL,XL,YL)
09600 DO 19 I=1,NR
09700 S(I)=SQRT(TT(I,I))
09800 DO 20 J=1,I
09900 TT(I,J)=TT(I,J)/(S(I)*S(J))
10000 20 TT(J,I)=TT(I,J)
10100 DO 19 J=1,NV
10200 19 A(J,I)=A(J,I)*S(I)
10300 89 CALL PFC(TT,NR)
10400 IF(NOR) 88,88,71
10500 71 DO 36 I=1,NR
10600 DO 36 J=1,NV
10700 36 A(J,I)=A(J,I)*FL(J)
10800 88 RETURN
10900 END
11000 SUBROUTINE MATINV(A,N,T,IN,V,U)
11100 DIMENSION A(N,N),IN(2),V(2),U(2)
11200 DO 1 I=1,N
11300 V(I)=A(I,I)
11400 1 IN(I)=0
11500 K=1
11600 DO 7 L=1,N
11700 DO 2 I=1,N
11800 U(I)=A(K,I)
11900 IF(I.GT.K) U(I)=A(I,K)
12000 A(I,K)=0.
12100 2 A(K,I)=0.
12200 P=U(K)
12300 T=H
12400 H=-1.E20
12500 IN(K)=1
12600 U(K)=-1.
12700 DO 7 I=1,N
12800 Y=U(I)/P
12900 DO 4 J=1,I
13000 4 A(I,J)=A(I,J)-Y*U(J)
13100 IF(IN(I)) 5,5,7
13200 5 IF(H-A(I,I)/V(I)) 6,7,7
13300 6 H=A(I,I)/V(I)
13400 K=I
13500 7 CONTINUE
13600 DO 8 I=1,N
13700 DO 8 J=1,I
13800 A(I,J)=-A(I,J)
13900 8 A(J,I)=A(I,J)
14000 RETURN
14100 END
14200 C SUBROUTINE ROTAT3 FOR BMDX72 MARCH 28, 1966
14300 SUBROUTINE ROTAT3(A,TT,LV,LR,GG,ACC,MR,INI,NOR,S,C,V,TL,FL,XL,YL)
14400 DIMENSION A(LV,LR),TT(LR,LR),C(2),S(2),V(2),TL(2),FL(2),XL(2),YL(2
14500 X)
14600 WRITE (6,882)
14700 882 FORMAT(40H1ROTATION FOR SIMPLE LOADINGS
14800 X/23H0ITERATION SIMPLICITY/13X,9HCRITERION)
14900 NV=LV
15000 NR=LR
15100 GA=GG/FLOAT(NV)
15200 DO 30 I=1,NV
15300 S(I)=0.
15400 DO 4 J=1,NR
15500 4 S(I)=S(I)+A(I,J)*A(I,J)
15600 IF(NOR) 30,30,32
15700 32 FL(I)=SQRT(S(I))
15800 S(I)=1.
15900 DO 2 J=1,NR
16000 2 A(I,J)=A(I,J)/FL(I)
16100 30 CONTINUE
16200 DO 33 I=1,NR
16300 IF(INI) 31,31,3
16400 3 DO 6 J=1,NR
16500 6 TT(I,J)=0.
16600 TT(I,I)=1.
16700 31 C(I)=0.
16800 V(I)=0.
16900 DO 33 J=1,NV
17000 AA=A(J,I)*A(J,I)
17100 C(I)=C(I)+AA
17200 33 V(I)=V(I)+AA*AA
17300 G=0.
17400 D=0.
17500 DO 34 J=1,NR
17600 D=D+C(J)
17700 V(J)=V(J)-GA*C(J)*C(J)
17800 34 G=G+V(J)
17900 H=0.
18000 DO 5 I=1,NV
18100 5 H=H+S(I)*S(I)
18200 H=H-GA*D*D
18300 F0=H-G
18400 FCC=F0*ACC
18500 L=0
18600 WRITE (6,52) L,F0
18700 52 FORMAT(I6,F16.6)
18800 70 DO 80 L=1,MR
18900 DO 81 IP=1,NR
19000 DO 81 IQ=1,NR
19100 IF(IP-IQ) 82,81,82
19200 82 D=D-C(IP)-C(IQ)
19300 G=G-V(IP)-V(IQ)
19400 P=0.
19500 R=0.
19600 T=0.
19700 O=0.
19800 Y=0.
19900 Z=0.
20000 DO 83 I=1,NV
20100 A1=A(I,IP)
20200 A2=A(I,IQ)
20300 AA=A1*A1
20400 BB=A2*A2
20500 AB=A1*A2
20600 S(I)=S(I)-AA-BB
20700 Z=Z+AB
20800 R=R+AA*BB
20900 P=P+AB*AA
21000 T=T+AA*S(I)
21100 83 O=O+AB*S(I)
21200 X=TT(IP,IQ)
21300 GAC=GA*C(IP)
21400 R=R-GAC*C(IQ)
21500 P=P-GAC*Z
21600 O=O-GA*Z*D
21700 U=U-GA*C(IQ)*D
21800 T=T-GAC*D
21900 100 P1=1.5*(X-P/V(IP))
22000 Q1=.5*(V(IP)-4.*X*P+R+2.*T)/V(IP)
22100 R1=.5*(X*(T+R)-P-O)/V(IP)
22200 CALL ROOT(P1,Q1,R1,A2)
22300 A22=A2*A2
22400 P=1.+2.*X*A2+A22
22500 IF(P)89,90,90
22600 89 WRITE (6,91) P
22700 91 FORMAT(2H $F20.7)
22800 P=-P
22900 90 A1=SQRT(P)
23000 A3=A2/A1
23100 A11=P*P
23200 C(IP)=P*C(IP)
23300 V(IP)=A11*V(IP)
23400 Z=0.
23500 Y=0.
23600 DO 84 I=1,NV
23700 A(I,IQ)=A(I,IQ)-A2*A(I,IP)
23800 A(I,IP)=A1*A(I,IP)
23900 BB=A(I,IQ)*A(I,IQ)
24000 Z=Z+BB*BB
24100 Y=Y+BB
24200 84 S(I)=S(I)+A(I,IP)*A(I,IP)+BB
24300 V(IQ)=Z-GA*Y*Y
24400 C(IQ)=Y
24500 D=D+Y+C(IP)
24600 G=G+V(IP)+V(IQ)
24700 A1=1./A1
24800 DO 85 I=1,NR
24900 TT(I,IP)=A1*TT(I,IP)+A3*TT(I,IQ)
25000 85 TT(IP,I)=TT(I,IP)
25100 TT(IP,IP)=1.
25200 81 CONTINUE
25300 H=0.
25400 DO 86 I=1,NV
25500 86 H=H+S(I)*S(I)
25600 F=H-GA*D*D-G
25700 WRITE (6,52) L,F
25800 157 IF(ABS(F-F0)-FCC) 38,38,80
25900 80 F0=F
26000 38 CONTINUE
26100 CALL PFC(TT,NR)
26200 IF(NOR) 88,88,71
26300 71 DO 36 I=1,NR
26400 DO 36 J=1,NV
26500 36 A(J,I)=A(J,I)*FL(J)
26600 88 RETURN
26700 END
26800 C SUBROUTINE ROOT FOR BMDX72 MARCH 28, 1966
26900 SUBROUTINE ROOT(P,Q,R,X)
27000 H=(ABS(P)+ABS(Q)+ABS(R))*1.E-5
27100 P2=2.*P
27200 A=P*P-3.*Q
27300 IF(A) 1,1,2
27400 1 X=0.
27500 GO TO 5
27600 2 A=SQRT(A)
27700 X=(A-P)/3.
27800 X1=-(P+A)/3.
27900 IF(X*(Q+X*(P+X))+X1*(Q+X1*(P+X1))+R+R)3,4,4
28000 3 X=X+1.
28100 GO TO5
28200 4 X=X1-1.
28300 5 DO 7 I=1,50
28400 F=R+X*(Q+X*(P+X))
28500 FP=Q+X*(P2+3.*X)
28600 DX=F/FP
28700 X=X-DX
28800 IF(ABS(F)-H) 6,6,7
28900 7 CONTINUE
29000 WRITE (6,8) P,Q,R,X
29100 8 FORMAT(2H *4E20.8)
29200 6 RETURN
29300 END
29400 SUBROUTINE ATOF(A,N,F)
29500 DIMENSION A(1)
29600 LOGICAL BLANK
29700 BLANK=.TRUE.
29800 S=1.0
29900 NUMB=0
30000 TEN=1.0
30100 DIV=1.0
30200 DO 10 I=1,N
30300 L=INTCHR(A,I)
30400 IF(L.EQ.36) GO TO 10
30500 BLANK=.FALSE.
30600 IF(L.NE.38) GO TO 2
30700 S=-1.0
30800 GO TO 10
30900 2 IF(L.NE.44) GO TO 4
31000 TEN=10.0
31100 GO TO 10
31200 4 IF(L.GT.9) GO TO 9
31300 NUMB=NUMB*10+L
31400 DIV=DIV*TEN
31500 9 CONTINUE
31600 10 CONTINUE
31700 IF(BLANK)RETURN
31800 F=S*FLOAT(NUMB)/DIV
31900 RETURN
32000 END
32100 FUNCTION INTCHR(STRING,N)
32200 DIMENSION SEQ(50),STRING(1),EBCD(5)
32300 DATA SEQ/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,
32400 X 1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,
32500 X 1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,
32600 X 1HU,1HV,1HW,1HX,1HY,1HZ,1H ,1H+,1H-,1H*,
32700 X 1H/,1H(,1H),1H,,1H.,1H',1H=,1H$,1H ,1H /
32800 DATA EBCD/1H+,1H(,1H),1H',1H=/
32900 CALL GETCHR(STRING,N,CHR)
33000 IF (CHR.NE.SEQ(37)) GO TO 2
33100 INTCHR = 36
33200 GO TO 10
33300 2 DO 1 I=1,48
33400 IF(SEQ(I).EQ.CHR) GO TO 9
33500 1 CONTINUE
33600 I=51
33700 IF(EBCD(1).EQ.CHR) I=38
33800 IF(EBCD(2).EQ.CHR) I=42
33900 IF(EBCD(3).EQ.CHR) I=43
34000 IF(EBCD(4).EQ.CHR) I=46
34100 IF(EBCD(5).EQ.CHR) I=47
34200 9 INTCHR=I-1
34300 10 RETURN
34400 END
34500
34600
34700
34800
34900
35000
35100
35200
35300
35400
35500
35600
35700
35800
35900
36000
36100
36200
36300
36400
36500
36600
36700
36800
36900
37000
37100
37200
37300
37400
37500
37600
37700
37800
37900
38000
38100
38200
38300
38400
38500
38600
38700
38800
38900
39000
39100
39200
39300
39400
39500
39600
39700
39800
39900
40000
40100