Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/facto/facto.for
There is 1 other file named facto.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C FACTO.F4 (FILENAME ON LIBRARY DECTAPE)
C FACTO, 1.6.1 (CALLING NAME, SUBLST #)
C FACTOR ANALYSIS
C ADAPTED FROM IBM FSS (FORTRAN SCIENTIFIC SUBROUTINES)
C MODIFIED BY B. HOUCHARD, N. GRANT
C LIBRARY DECTAPE PROGS. USED: USAGE.MAC
C FORWMU PROGS. USED: TTYPTY, ALLCOR, DEVCHG, EXISTS, PRINTS,
C ZEROP
C APLIB PROGS. USED: GETFOR, IO
C INTERNAL SUBR. USED: FACANL, CORRN, TRACE, LOAD, VARMX,
C INVT, EIGEN
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C
C SUBROUTINES USED:
C
C TTYPTY (*) DETERMINE IF JOB IS FROM TELETYPE OR
C PSEUDO-TELETYPE
C
C USAGE (*) COUNTER FOR LIBRARY PROGRAMS USAGE
C
C IO INPUT/OUTPUT SUBROUTINE
C
C GETFOR FORMAT SUBROUTINE
C
C ALLCOR (*) ALLOCATES CORE DYNAMICALLY
C
C FACANL MAIN SUBROUTINE FOR THE PROGRAM WHICH
C IN TURN CALLS IN OTHER SUBROUTINES
C
C (*) MACRO SUBROUTINES
C
C
C LIMITATIONS:
C
C (1) NUMBER OF VARIABLES IS AT MOST 95 IF FACTOR SCORE IS ASKED;
C OTHERWISE MAXIMUM IS 123
C (2) AT MOST 3 LINES IF OBJECT TIME FORMAT IS USED
C (3) ONLY F-TYPE FORMAT ALLOWED
C
C***********************************************************************
C***********************************************************************
C
DIMENSION SPACE(1),NOTF(48),ID(16)
C
COMMON/BLOCK1/IDLG,ICC,INP,IOUT,IO2,ISTD,ICODE,NOTF,ID
COMMON/BLOCK2/ITMP,FSCORE
C
C***********************************************************************
C DEVICES USED:
C
C IDLG--DEVICE USED TO COMMUNICATE WITH USER
C IT IS ALWAYS SET TO -1
C ICC---DEVICE USED TO ACCEPT USER'S RESPONSES
C IT IS ALWAYS SET TO -4
C INP---DEVICE USED TO READ DATA
C ITS LOGICAL NUMBER IS DETERMINED BY SUBROUTINE IO
C IOUT--DEVICE USED TO WRITE OUT THE REPORT
C ITS LOGICAL NUMBER IS DETERMINED BY SUBROUTINE IO
C***********************************************************************
C
ITMP=1
IDLG=-1
ICC=-4
INP=2
IOUT=3
1 FORMAT(' HOW MANY VARIABLES? ',$)
3 FORMAT(I)
4 FORMAT(' HOW MANY OBSERVATIONS? ',$)
5 FORMAT(' ENTER VALUE USED TO SELECT EIGENVALUES(USUALLY 1.0):
1 ',$)
6 FORMAT(F)
7 FORMAT(A3)
C
C**********************************************************************
C DETERMINE IF JOB IS ON TELETYPE OR PSEUDO-TELETYPE
C IF ICODE = 0 JOB IS ON TELETYPE
C = -1 JOB IS ON PSEUDO-TELETYPE
C***********************************************************************
C
C---------------ICODE RETURNED
CALL TTYPTY(ICODE)
C
C************************************************************************
C CALL SUBROUTINE USAGE AND ADD 1 TO USAGE
C***********************************************************************
C
C CALL USAGE('FACTO')
C
C**********************************************************************
C GATHER INPUT/OUTPUT INFORMATION, OUTPUT OPTION IS AVAILABLE
C ONLY ONCE IN THE PROGRM
C***********************************************************************
C
C---------------IO3, IO2 ARE RETURNED. OTHER ARGS. ARE INPUT.
C--------------- 1 MEANS OUTPUT? PRINTS, 0 - INPUT? PRINTS.
CALL IO(IOUT,IO3,ICC,IDLG,1,ICODE)
203 CALL IO(INP,IO2,ICC,IDLG,0,ICODE)
C
C***********************************************************************
C FORMAT OPTION, ITYPE = 2 MEANS ONLY F-TYPE FORMAT ALLOWED
C***********************************************************************
C
ITYPE=2
C---------------NOTF, ISTD ARE RETURNED. OTHER ARGS ARE INPUT.
C--------------- 48=NO. OF OBJ. TIME FORMAT WORDS (3 LINES).
CALL GETFOR(IDLG,ICC,NOTF,ISTD,48,ITYPE)
C
C***********************************************************************
C GATHER OTHER INPUT INFORMATION, CALL SUBROUTINE ALLCOR AND
C MAIN SUBROUTINE FACANL
C***********************************************************************
C
WRITE(IDLG,200)
200 FORMAT(' ENTER HEADER'/)
READ(ICC,201) ID
201 FORMAT(16A5)
100 WRITE(IDLG,1)
C---------------M=NO. VARS.
READ(ICC,3,END=1000)M
IF (M.LE.0) GO TO 110
102 WRITE(IDLG,4)
C---------------N=NO. OF OBS.
READ(ICC,3,END=1000) N
IF (N.GT.1) GO TO 103
WRITE(IDLG,300)
300 FORMAT(' SAMPLE SIZE TOO SMALL'/)
IF (ICODE.GE.0) GO TO 102
CALL EXIT
103 WRITE(IDLG,5)
C---------------CON=NO. USED TO SELECT EIGENVALUES
READ(ICC,6,END=1000) CON
WRITE(IDLG,30)
30 FORMAT(' FACTOR SCORE? (YES OR NO)--',$)
READ(ICC,7) FSCORE
MAX=M*(M+6)+M*(M+1)/2
IF (FSCORE.EQ.'YES') MAX=MAX+M*M
CALL ALLCOR(MAX,IERR,IBASE,SPACE(1))
IF(IERR.EQ.0)GO TO 101
110 WRITE(IDLG,31)
31 FORMAT(' NUMBER OF VARIABLES OUTSIDE ALLOWABLE RANGE'/)
IF (ICODE.GE.0) GO TO 100
CALL EXIT
101 I2=IBASE+M*M
I3=I2+M*(M+1)/2
I4=I3+M
I5=I4+M
I6=I5+M
I7=I6+M
I8=I7+M
I9=IBASE
IF (FSCORE.EQ.'YES') I9=I8+M
CALL FACANL(M,N,CON,SPACE(IBASE),SPACE(I2),SPACE(I3),
1SPACE(I4),SPACE(I5),SPACE(I6),SPACE(I7),SPACE(I8),SPACE(I9))
C
C***********************************************************************
C END OF ONE DATA SET, BRANCH BACK TO DETERMINE IF MORE DATA
C IS TO BE PROCESSED
C***********************************************************************
C
WRITE(IDLG,202)
202 FORMAT('-')
GO TO 203
1000 STOP
END
C
C
C SUBROUTINE FACANL
C
C MAIN SUBROUTINE OF THE FACTOR ANALYSIS PROGRAM
C
C SUBROUTINES USED:
C
C CORRN
C EIGEN
C TRACE
C LOAD
C VARMX
C INVT
C***********************************************************************
C
C---------------M, N, CON ARE INPUT. OTHER ARGS. ARE SPACES
C--------------- RESERVED BY DYN. ALLOC.
SUBROUTINE FACANL(M,N,CON,V,R,B,D,S,T,XBAR,SD,DCOR)
C
DIMENSION B(1),D(1),S(1),T(1),XBAR(1),SD(1),V(1),R(1),NOTF(48),
1 ID(16),DCOR(M,1)
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO 51..
DIMENSION TV(56)
C---------------SUBR. CORRN, MAIN PROG. SHARE COMMON /BLOCK1/
C--------------- AND /BLOCK2/
COMMON/BLOCK1/IDLG,ICC,INP,IOUT,IO2,ISTD,ICODE,NOTF,ID
COMMON/BLOCK2/ITMP,FSCORE
1 FORMAT(' NO. OF OBSERVATIONS',4X,I6/' NO. OF VARIABLES',I6)
2 FORMAT('0MEANS'/(1X,5F12.4))
3 FORMAT('0STANDARD DEVIATIONS'/(1X,5F12.4))
4 FORMAT('0CORRELATION COEFFICIENTS')
5 FORMAT(4H0ROWI4/(1X,5F12.5))
6 FORMAT('-EIGENVALUES'/(1X,5F12.5))
7 FORMAT('0CUMULATIVE PERCENTAGE OF EIGENVALUES'/(1X,5F12.5))
8 FORMAT('-EIGENVECTORS')
9 FORMAT('0VECTOR'I4/(1X,5F12.5))
10 FORMAT('-FACTOR MATRIX (',I4,' FACTORS)')
11 FORMAT('0VARIABLE'I4/(1X,5F12.5))
12 FORMAT('-ITERATION',7X,'VARIANCES'/' CYCLE')
13 FORMAT(I6,F20.6)
14 FORMAT('-ROTATED FACTOR MATRIX ('I4,' FACTORS)')
15 FORMAT('0VARIABLE'I4/(1X,5F12.5))
16 FORMAT('-CHECK ON COMMUNALITIES'/'0VARIABLE',7X,'ORIGINAL',1
12X,'FINAL',10X,'DIFFERENCE')
17 FORMAT(I6,3F18.5)
19 FORMAT('0ONLY',I2,' FACTOR RETAINED. NO ROTATION')
39 FORMAT('1FACTOR ANALYSIS'/)
40 FORMAT(1X,16A5)
46 FORMAT('-FACTOR SCORES')
48 FORMAT(1X,I4,4F/(5X,4G))
C
C ..................................................................
C
IF (ISTD.NE.1) GO TO 201
NOTF(1)='(20F)'
DO 202 I=2,48
202 NOTF(I)=' '
201 IF (IO2.EQ.'DSK') GO TO 204
WRITE(IDLG,203)
203 FORMAT(' ENTER DATA'/)
IF (ISTD.EQ.1) WRITE(IDLG,2100)
2100 FORMAT('+(AT MOST 20 NUMBERS PER LINE,SEPARATED BY COMMAS)'/)
GO TO 206
204 WRITE(IDLG,205)
205 FORMAT(' PLEASE WAIT, YOUR DATA IS BEING PROCESSED'/)
206 X=0
CALL CORRN(N,M,X,XBAR,S,V,R,D,B,T)
WRITE(IOUT,39)
DO 207 I=1,16
IF (ID(I).NE.' ') GO TO 208
207 CONTINUE
GO TO 209
208 WRITE(IOUT,40) ID
209 WRITE(IOUT,1)N,M
C PRINT MEANS
WRITE (IOUT,2) (XBAR(J),J=1,M)
C COPY STANDARD DEVIATIONS
DO 420 J=1,M
420 SD(J)=S(J)
C PRINT STANDARD DEVIATIONS
WRITE (IOUT,3) (S(J),J=1,M)
C PRINT CORRELATION COEFFICIENTS
WRITE (IOUT,4)
DO 120 I=1,M
DO 110 J=1,M
IF(I-J) 102,104,104
102 L=I+(J*J-J)/2
GO TO 110
104 L=J+(I*I-I)/2
110 D(J)=R(L)
120 WRITE (IOUT,5) I,(D(J),J=1,M)
MV=0
CALL EIGEN (R,V,M,MV)
CALL TRACE (M,R,CON,K,D)
C PRINT EIGENVALUES
DO 130 I=1,K
L=I+(I*I-I)/2
130 S(I)=R(L)
WRITE (IOUT,6) (S(J),J=1,K)
C PRINT CUMULATIVE PERCENTAGE OF EIGENVALUES
WRITE (IOUT,7) (D(J),J=1,K)
C PRINT EIGENVECTORS
WRITE (IOUT,8)
L=0
DO 150 J=1,K
DO 140 I=1,M
L=L+1
140 D(I)=V(L)
150 WRITE(IOUT,9) J,(D(I),I=1,M)
CALL LOAD (M,K,R,V)
C PRINT FACTOR MATRIX
WRITE (IOUT,10) K
DO 180 I=1,M
DO 170 J=1,K
L=M*(J-1)+I
170 D(J)=V(L)
180 WRITE (IOUT,11) I,(D(J),J=1,K)
IF(K-1) 185,185,188
185 WRITE(IOUT,19) K
GO TO 100
188 CALL VARMX (M,K,V,NC,TV,B,T,D)
C PRINT VARIANCES
NV=NC+1
WRITE (IOUT,12)
DO 190 I=1,NV
NC=I-1
190 WRITE (IOUT,13) NC,TV(I)
C PRINT ROTATED FACTOR MATRIX
WRITE (IOUT,14) K
DO 220 I=1,M
DO 210 J=1,K
L=M*(J-1)+I
210 S(J)=V(L)
220 WRITE (IOUT,15) I,(S(J),J=1,K)
C PRINT COMMUNALITIES
WRITE (IOUT,16)
DO 230 I=1,M
230 WRITE (IOUT,17) I,B(I),T(I),D(I)
100 IF(FSCORE.NE.'YES')RETURN
CALL RELEAS (ITMP)
DO 400 I=1,K
L1=M*(I-1)
DO 400 J=1,K
L2=M*(J-1)
DCOR(I,J)=0
DO 400 L=1,M
400 DCOR(I,J)=DCOR(I,J)+V(L1+L)*V(L2+L)
CALL INVT(M,DCOR,S,T,K)
DO 410 J=1,M
DO 411 I=1,K
T(I)=0
DO 411 L=1,K
L1=M*(L-1)+J
C---------------SEE FACTOR SCORE FURMULA PAGE 4 OF WRITE UP
411 T(I)=T(I)+DCOR(L,I)*V(L1)
DO 412 L=1,K
L1=M*(L-1)+J
412 V(L1)=T(L)
410 CONTINUE
WRITE(IOUT,46)
DO 370 NOBS=1,N
READ(ITMP)(D(IJB),IJB=1,M)
DO 371 J=1,K
T(J)=0.
L1=M*(J-1)
DO 371 I=1,M
L=L1+I
C---------------SEE ST. 412 AND FACTOR SCORE FORMULA IN WRITE UP
371 T(J)=T(J)+((D(I)-XBAR(I))/SD(I))*V(L)
370 WRITE(IOUT,48)NOBS,(T(I),I=1,K)
CALL RELEASE(ITMP)
RETURN
1000 STOP
END
C---------------INVERT MATRIX A. MVAR, A, L, M, N, ARE INPUT.
C--------------- INVERSE OF A RETURNED IN A
SUBROUTINE INVT(MVAR,A,L,M,N)
DIMENSION A(MVAR,1),L(1),M(1)
D=1.0
DO 9 K=1,N
L(K)=K
M(K)=K
B=A(K,K)
DO 1 I=K,N
DO 1 J=K,N
IF (ABS(B).GE.ABS(A(I,J))) GO TO 1
B=A(I,J)
L(K)=I
M(K)=J
1 CONTINUE
J=L(K)
IF (L(K).LE.K) GO TO 3
DO 2 I=1,N
C=-A(K,I)
A(K,I)=A(J,I)
2 A(J,I)=C
3 I=M(K)
IF (M(K).LE.K) GO TO 5
DO 4 J=1,N
C=-A(J,K)
A(J,K)=A(J,I)
4 A(J,I)=C
5 DO 6 I=1,N
IF (I.NE.K) A(I,K)=A(I,K)/(-A(K,K))
6 CONTINUE
DO 7 I=1,N
DO 7 J=1,N
IF ((I.NE.K).AND.(J.NE.K)) A(I,J)=A(I,K)*A(K,J)+A(I,J)
7 CONTINUE
DO 8 J=1,N
IF (J.NE.K) A(K,J)=A(K,J)/A(K,K)
8 CONTINUE
D=D*A(K,K)
A(K,K)=1.0/A(K,K)
9 CONTINUE
K=N
10 K=K-1
IF (K.LE.0) GO TO 14
I=L(K)
IF (I.LE.K) GO TO 12
DO 11 J=1,N
C=A(J,K)
A(J,K)=-A(J,I)
11 A(J,I)=C
12 J=M(K)
IF (J.LE.K) GO TO 10
DO 13 I=1,N
C=A(K,I)
A(K,I)=-A(J,I)
13 A(J,I)=C
GO TO 10
14 RETURN
END
C
C
C
C SUBROUTINE CORRN
C
C MACRO SUBROUTINES USED:
C
C ZEROP (LISTED IN FORLIB.MAN - WG)
C DEVICE (APPARENTLY NOT USED - WG)
C
C***********************************************************************
C
C---------------N, M ARE INPUT; XBAR, STD, RX, R, B, D, T ARE RETURNED
C---------------X(1) APPARENTLY NOT USED IN CORRN
SUBROUTINE CORRN(N,M,X,XBAR,STD,RX,R,B,D,T)
DIMENSION X(1),XBAR(1),STD(1),RX(1),R(1),B(1),D(1),T(1)
DIMENSION NOTF(48),ID(16)
COMMON/BLOCK1/IDLG,ICC,INP,IOUT,IO2,ISTD,ICODE,NOTF,ID
COMMON/BLOCK2/ITMP,FSCORE
C INITIALIZATION
K=(M*M+M)/2
C---------------INITALIZE B, T, R TO ZERO
CALL ZEROP(B,M)
CALL ZEROP(T,M)
CALL ZEROP(R,K)
FN=N
L=0
C READ OBSERVATIONS AND CALCULATE TEMPORARY
C MEANS FROM THESE DATA IN T(J)
127 IF(N-M) 130,130,135
130 KK=N
GO TO 137
135 KK=M
137 DO 140 I=1,KK
352 READ (INP,NOTF,END=1000) (D(IJB),IJB=1,M)
351 IF(FSCORE.EQ.'YES')WRITE(ITMP)(D(IJB),IJB=1,M)
DO 140 J=1,M
T(J)=T(J)+D(J)
L=L+1
140 RX(L)=D(J)
FKK=KK
DO 150 J=1,M
XBAR(J)=T(J)
150 T(J)=T(J)/FKK
C CALCULATE SUMS OF CROSS-PRODUCTS OF DEVIATIONS
C FROM TEMPORARY MEANS FOR M OBSERVATIONS
L=0
DO 180 I=1,KK
JK=0
DO 170 J=1,M
L=L+1
170 D(J)=RX(L)-T(J)
DO 180 J=1,M
B(J)=B(J)+D(J)
DO 180 K=1,J
JK=JK+1
180 R(JK)=R(JK)+D(J)*D(K)
IF(N-KK)205,205,185
C READ THE REST OF OBSERVATIONS ONE AT A TIME, SUM
C THE OBSERVATIONS, AND CALCULATE SUMS OF CROSS-
C PRODUCTS OF DEVIATIONS FROM TEMPORARY MEANS
185 KK=N-KK
DO 200 I=1,KK
JK=0
362 READ (INP,NOTF,END=1000) (D(IJB),IJB=1,M)
361 IF(FSCORE.EQ.'YES')WRITE(ITMP)(D(IJB),IJB=1,M)
DO 190 J=1,M
XBAR(J)=XBAR(J)+D(J)
D(J)=D(J)-T(J)
190 B(J)=B(J)+D(J)
DO 200 J=1,M
DO 200 K=1,J
JK=JK+1
200 R(JK)=R(JK)+D(J)*D(K)
C CALCULATE MEANS
205 JK=0
DO 210 J=1,M
XBAR(J)=XBAR(J)/FN
C ADJUST SUMS OF CROSS-PRODUCTS OF DEVIATIONS
C FROM TEMPORARY MEANS
DO 210 K=1,J
JK=JK+1
210 R(JK)=R(JK)-B(J)*B(K)/FN
C CALCULATE CORRELATION COEFFICIENTS
JK=0
DO 220 J=1,M
JK=JK+J
220 STD(J)=SQRT(ABS(R(JK)))
DO 230 J=1,M
DO 230 K=J,M
JK=J+(K*K-K)/2
L=M*(J-1)+K
RX(L)=R(JK)
L=M*(K-1)+J
RX(L)=R(JK)
IF(STD(J)*STD(K)) 225,222,225
222 R(JK)=0.0
GO TO 230
225 R(JK)=R(JK)/(STD(J)*STD(K))
230 CONTINUE
C CALCULATE STANDARD DEVIATIONS
FN=SQRT(FN-1.0)
DO 240 J=1,M
240 STD(J)=STD(J)/FN
C COPY THE DIAGONAL OF THE MATRIX OF SUMS OF CROSS-PRODUCTS OF
C DEVIATIONS FROM MEANS.
L=-M
DO 250 I=1,M
L=L+M+1
250 B(I)=RX(L)
RETURN
1000 STOP
END
C---------------A, N, MV ARE INPUT. R IS OUTPUT.
SUBROUTINE EIGEN(A,R,N,MV)
DIMENSION A(1),R(1)
C GENERATE IDENTITY MATRIX
5 RANGE=1.0E-6
IF(MV-1) 10,25,10
10 IQ=-N
DO 20 J=1,N
IQ=IQ+N
DO 20 I=1,N
IJ=IQ+I
R(IJ)=0.0
IF(I-J) 20,15,20
15 R(IJ)=1.0
20 CONTINUE
C COMPUTE INITIAL AND FINAL NORMS(ANORM AND ANORMX)
25 ANORM=0.0
DO 35 I=1,N
DO 35 J=I,N
IF(I-J) 30,35,30
30 IA=I+(J*J-J)/2
ANORM=ANORM+A(IA)*A(IA)
35 CONTINUE
IF(ANORM) 165,165,40
40 ANORM=1.414*SQRT(ANORM)
ANRMX=ANORM*RANGE/FLOAT(N)
C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR
IND=0
THR=ANORM
45 THR=THR/FLOAT(N)
50 L=1
55 M=L+1
C COMPUTE SIN AND COS
60 MQ=(M*M-M)/2
LQ=(L*L-L)/2
LM=L+MQ
62 IF(ABS(A(LM))-THR) 130,65,65
65 IND=1
LL=L+LQ
MM=M+MQ
X=0.5*(A(LL)-A(MM))
68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X)
IF(X) 70,75,75
70 Y=-Y
75 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y))))
SINX2=SINX*SINX
78 COSX=SQRT(1.0-SINX2)
COSX2=COSX*COSX
SINCS=SINX*COSX
C ROTATE L AND M COLUMNS
ILQ=N*(L-1)
IMQ=N*(M-1)
DO 125 I=1,N
IQ=(I*I-I)/2
IF(I-L) 80,115,80
80 IF(I-M) 85,115,90
85 IM=I+MQ
GO TO 95
90 IM=M+IQ
95 IF(I-L) 100,105,105
100 IL=I+LQ
GO TO 110
105 IL=L+IQ
110 X=A(IL)*COSX-A(IM)*SINX
A(IM)=A(IL)*SINX+A(IM)*COSX
A(IL)=X
115 IF(MV-1) 120,125,120
120 ILR=ILQ+I
IMR=IMQ+I
X=R(ILR)*COSX-R(IMR)*SINX
R(IMR)=R(ILR)*SINX+R(IMR)*COSX
R(ILR)=X
125 CONTINUE
X=2.0*A(LM)*SINCS
Y=A(LL)*COSX2+A(MM)*SINX2-X
X=A(LL)*SINX2+A(MM)*COSX2+X
A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
A(LL)=Y
A(MM)=X
C TESTS FOR COMPLETION
C TEST FOR M=LAST COLUMN
130 IF(M-N) 135,140,135
135 M=M+1
GO TO 60
C TEST FOR L=SECOND FROM LAST COLUMN
140 IF(L-(N-1)) 145,150,145
145 L=L+1
GO TO 55
150 IF(IND-1) 160,155,160
155 IND=0
GO TO 50
C COMPARE THRESHOLD WITH FINAL NORM
160 IF(THR-ANRMX) 165,165,45
C SORT EIGENVALUES AND EIGENVECTORS
165 IQ=-N
DO 185 I=1,N
IQ=IQ+N
LL=I+(I*I-I)/2
JQ=N*(I-2)
DO 185 J=I,N
JQ=JQ+N
MM=J+(J*J-J)/2
IF(A(LL)-A(MM)) 170,185,185
170 X=A(LL)
A(LL)=A(MM)
A(MM)=X
IF(MV-1) 175,185,175
175 DO 180 K=1,N
ILR=IQ+K
IMR=JQ+K
X=R(ILR)
R(ILR)=R(IMR)
180 R(IMR)=X
185 CONTINUE
RETURN
END
C PURPOSE
C COMPUTE CUMULATIVE PERCENTAGE OF EIGENVALUES GREATER THAN
C OR EQUAL TO A CONSTANT SPECIFIED BY THE USER. THIS SUB-
C ROUTINE NORMALLY OCCURS IN A SEQUENCE OF CALLS TO SUB-
C ROUTINES CORRE,EIGEN,TRACE,LOAD,AND VARMX IN THE PER-
C FORMANCE OF A FACTOR ANALYSIS.
C USAGE
C CALL TRACE (M,R,CON,K,O)
C DESCRIPTION OF PARAMETERS
C M - NUMBER OF VARIABLES.
C R - INPUT MATRIX (SYMETRIC AND STORED IN COMPRESSED
C FORM WITH ONLY UPPER TRIANGLE BY COLUMN IN CORE)
C CONTAINING EIGENVALUES IN DIAGONAL. EIGENVALUES ARE
C ARRANGED IN DESCENDING ORDER. THE ORDER OF MATRIX R
C IS M BY M. ONLY M*(M+1)/2 ELEMENTS ARE IN STORAGE.
C (STORAGE MODE OF 1)
C CON - A CONSTANT USED TO DECIDE HOW MANY EIGENVALUES TO
C RETAIN. CUMULATIVE PERCENTAGE OF EIGENVALUES
C WHICH ARE GREATER THAN OR EQUAL TO THIS VALUE IS
C CALCULATED.
C K - OUTPUT VARIABLE CONTAINING THE NUMBER OF EIGENVALUES
C GREATER THAN OR EQUAL TO CON. (K IS THE NUMBER OF
C FACTORS.)
C D - OUTPUT VECTOR OF LENGTH M CONTAINING CUMULATIVE
C PERCENTAGE OF EIGENVALUES WHICH ARE GREATER THAN
C OR EQUAL TO CON.
C REMARKS
C NONE
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C METHOD
C EACH EIGENVALUE GREATER THAN OR EQUAL TO CON IS DIVIDED BY M
C AND THE RESULT IS ADDED TO THE PREVIOUS TOTAL TO OBTAIN
C THE CUMULATIVE PERCENTAGE FOR EACH EIGENVALUE.
C---------------M, R, CON ARE INPUT. K, D ARE OUTPUT.
SUBROUTINE TRACE (M,R,CON,K,D)
DIMENSION R(1),D(1)
C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C STATEMENT WHICH FOLLOWS.
C DOUBLE PRECISION R,D
C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C ROUTINE.
FM=M
L=0
DO 100 I=1,M
L=L+I
100 D(I)=R(L)
K=0
C TEST WHETHER 1-TH EIGENVALUE IS GREATER
C THAN OR EQUAL TO THE CONSTANT
DO 110 I=1,M
IF(D(I)-CON) 120, 105, 105
105 K=K+1
110 D(I)=D(I)/FM
C COMPUTE CUMULATIVE PERCENTAGE OF EIGENVALUES
120 DO 130 I=2,K
130 D(I)=D(I)+D(I-1)
RETURN
END
C PURPOSE
C COMPUTE A FACTOR MATRIX (LOADING) FROM EIGENVALUES AND
C ASSOCIATED EIGENVECTORS. THIS SUBROUTINE NORMALLY OCCURS
C IN A SEQUENCE OF CALLS TO SUBROUTINES CORRE, EIGEN, TRACE,
C LOAD, AND VARMX IN THE PERFORMANCE OF A FACTOR ANALYSIS.
C USAGE
C CALL LOAD (M,K,R,V)
C DESCRIPTION OF PARAMETERS
C M -NUMBER OF VARIABLES.
C K -NUMBER OF FACTORS.
C R -A MATRIX (SYMMETRIC AND STORED IN COMPRESSED FORM
C WITH ONLY UPPER TRIANGLE BY COLUMN IN CORE) CON-
C TAINING EIGENVALUES IN DIAGONAL. EIGENVALUES ARE
C ARRANGED IN DESCENDING ORDER, AND FIRST K
C EIGENVALUES ARE USED BY THIS SUBROUTINE. THE ORDER
C OF MATRIX R IS M BY M. ONLY M*(M+1)/2 ELEMENTS ARE
C IN STORAGE. (STORAGE MODE OF 1)
C V -WHEN THIS SUBROUTINE IS CALLED, MATRIX V (M X M)
C CONTAINS EIGENVECTORS COLUMNVISE. UPON RETURNING TO
C THE CALLING PROGRAM, MATRIX V CONTAINS A FACTOR
C MATRIX (M X K).
C REMARKS
C NONE
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C METHOD
C NORMALIZED EIGENVECTORS ARE CONVERTED TO THE FACTOR PATTERN
C BY MULTIPLYING THE ELEMENTS OF EACH VECTOR BY THE SQUARE
C ROOT OF THE CORRESPONDING EIGENVALUE.
C---------------M, K, R, V ARE INPUT. V IS MODIFIED.
SUBROUTINE LOAD (M,K,R,V)
DIMENSION R(1),V(1)
C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C STATEMENT WHICH FOLLOWS.
C DOUBLE PRECISION R,V,SQ
C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C ROUTINE.
C THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
C CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. SQRT IN STATEMENT
C 150 MUST BE CHANGED TO DSQRT.
L=0
JJ=0
DO 160 J=1,K
JJ=JJ+J
150 SQ= SQRT(R(JJ))
DO 160 I=1,M
L=L+1
160 V(L)=SQ*V(L)
RETURN
END
C---------------M, K, A, F, D ARE INUT. NC, TV, H ARE OUTPUT.
SUBROUTINE VARMX (M,K,A,NC,TV,H,F,D)
DIMENSION A(1),TV(1),H(1),F(1),D(1)
C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C DOUBLE PRECISION A,TV,H,F,D,TVLT,CONS,AA,BB,CC,DD,U,T,B,COS4T,
C 1 SIN4T,TAN4T,SINP,COSP,CTN4T,COS2T,SIN2T,COST,SINT
C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C ROUTINE.
C THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
C CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. SQRT IN STATEMENTS
C 115,290,330,350, AND 335 MUST BE CHANGED TO DSQRT. ABS IN
C STATEMENTS 280,320,AND 375 MUST BE CHANGED TO DABS.
C INITIALIZATION
EPS=0.00116
TVLT=0.0
LL=K-1
NV=1
NC=0
FN=M
FFN=FN*FN
CONS=0.7071066
C CALCULATE ORIGINAL COMMUNALITIES
DO 110 I=1,M
H(I)=0.0
DO 110 J=1,K
L=M*(J-1)+I
110 H(I)=H(I)+A(L)*A(L)
C CALCULATE NORMALIZED FACTOR MATRIX
DO 120 I=1,M
115 H(I)=SQRT(H(I))
DO 120 J=1,K
L=M*(J-1)+I
120 A(L)=A(L)/H(I)
GO TO 132
C CALCULATE VARIANCE FOR FACTOR MATRIX
130 NV=NV+1
TVLT=TV(NV-1)
132 TV(NV)=0.0
DO 150 J=1,K
AA=0.0
BB=0.0
LB=M*(J-1)
DO 140 I=1,M
L=LB+I
CC=A(L)*A(L)
AA=AA+CC
140 BB=BB+CC*CC
150 TV(NV)=TV(NV)+(FN*BB-AA*AA)/FFN
IF(NV-51) 160,430,430
C PERFORM CONVERGENCE TEST
160 IF((TV(NV)-TVLT)-(1.E-7))170, 170, 190
170 NC=NC+1
IF(NC-3) 190, 190, 430
C ROTATION OF TWO FACTORS CONTINUES UP TO
C THE STATEMENT 120.
190 DO 420 J=1,LL
L1=M*(J-1)
II=J+1
C CALCULATE NUM AND DEM PAGE 59
DO 420 K1=II,K
L2=M*(K1-1)
AA=0.0
BB=0.0
CC=0.0
DD=0.0
DO 230 I=1,M
L3=L1+I
L4=L2+I
U=(A(L3)+A(L4))*(A(L3)-A(L4))
T=A(L3)*A(L4)
T=T+T
CC=CC+(U+T)*(U-T)
DD=DD+2.0*U*T
AA=AA+U
230 BB=BB+T
T=DD-2.0*AA*BB/FN
B=CC-(AA*AA-BB*BB)/FN
C COMPARISON OF NUM AND DEN
IF(T-B) 280, 240, 320
240 IF((T+B)-EPS) 420,250,250
C NUM+DEN IS GREATER THAN OR EQUAL TO THE
250 COS4T=CONS
C TOLERANCE FACTOR
SIN4T=CONS
GO TO 350
C NUM IS LESS THAN DEN
280 TAN4T=ABS(T)/ ABS(B)
IF(TAN4T-EPS) 300, 290, 290
290 COS4T=1.0/SQRT(1.0+TAN4T*TAN4T)
SIN4T=TAN4T*COS4T
GO TO 350
300 IF(B) 310, 420, 420
310 SINP=CONS
COSP=CONS
GO TO 400
C NUM IS GREATER THAN DEN
320 CTN4T=ABS(T/B)
IF(CTN4T-EPS) 340, 330, 330
330 SIN4T=1.0/SQRT(1.0+CTN4T*CTN4T)
COS4T=CTN4T*SIN4T
GO TO 350
340 COS4T=0.0
SIN4T=1.0
C DETERMINE COS THETA AND SIN THETA
350 COS2T=SQRT((1.0+COS4T)/2.0)
SIN2T=SIN4T/(2.0*COS2T)
355 COST=SQRT((1.0+COS2T)/2.0)
SINT=SIN2T/(2.0*COST)
C DETERMINE COS PHI AND SIN PHI
IF(B) 370, 370, 360
360 COSP=COST
SINP=SINT
GO TO 380
370 COSP=CONS*COST+CONS*SINT
375 SINP=ABS(CONS*COST-CONS*SINT)
380 IF(T) 390,390,400
390 SINP=-SINP
C PERFORM ROTATION
400 DO 410 I=1,M
L3=L1+I
L4=L2+I
AA=A(L3)*COSP+A(L4)*SINP
A(L4)=-A(L3)*SINP+A(L4)*COSP
410 A(L3)=AA
420 CONTINUE
GO TO 130
C DENORMALIZE VARIMAX LOADINGS
430 DO 440 I=1,M
DO 440 J=1,K
L=M*(J-1)+I
440 A(L)=A(L)*H(I)
C CHECK ON COMMUNALITIES
NC=NV-1
DO 450 I=1,M
450 H(I)=H(I)*H(I)
DO 470 I=1,M
F(I)=0.0
DO 460 J=1,K
L=M*(J-1)+I
460 F(I)=F(I)+A(L)*A(L)
470 D(I)=H(I)-F(I)
RETURN
END