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

