Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - 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