Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/stp/stp6.for
There is 1 other file named stp6.for in the archive. Click here to see a list.
C *** STAT PACK ***
C SUBROUTINE FOR FACTOR ANALYSIS.
C CALLING SEQUENCE: CALL STFACT(NV,NC,MV,MC,DATA,STD,VMN,COR,SP,H,NAMES)
C WHERE NV - IS THE NUMBER OF COLUMNS ACTUALLY FILLED (VARIABLES)
C NC - IS THE NUMBER OF ROWS ACTUALLY FILLED (CASES)
C MV - IS THE MAXIMUM NUMBER OF COLUMNS AS SPECIFIED IN THE MAIN
C MC - IS THE MAXIMUM NUMBER OF ROWS, AS SPECIFIED IN THE MAIN
C DATA - STORAGE FOR DATA DIMENSIONED FOR MAXIMUM.
C STD - IS A VECTOR CONTAINING STANDARD DEVIATIONS.
C VMN - IS A VECTOR CONTAINING VARIABLE MEANS.
C COR - IS THE CORRELATION MATRIX, DIMENSIONED FOR MAXIMUM
C SP - IS AN EXTRA VECTOR DIMENSIONED AT LEAST NV.
C H - IS AN EXTRA VECTOR DIMENSIONED AT LEAST NV.
C NAMES - IS A VECTOR CONTIANING VARIABLE NAMES
C
C FACTOR ANALYSIS SET FOR A MAXIMUM OF 40 VARIABLES. OPTIONS
C INCLUDE PRINTING OF EIGEN VECTORS AND COMMUNALITIES.
C METHOD USED HERE IS VARMAX ROTATION, PROGRAM IS A MODERATE
C REWRITE OF THE EXISTING IBM FACTOR ANALYSIS ROUTINES. CALLED
C IN ADDITION ARE: EIGEN - CALCULATE EIGEN VECTORS
C LOAD - LOAD FACTOR MATRIX
C TRACE -
C VARMX - VARMAX ROTATION METHODE
C ALL ROUTINES ARE PART OF THE SAME OVERLAY.
C
C SAVING FACTOR SCORES AS VARIABLES WAS PROGRAMMED INTO THE EXISTING
C FACTOR ANALYSIS BY LOUISANA STATE UNIVERSITY - CLYDE POOLE.
C
C FACTOR SCORE CALCULATION USED IS NOW A REGRESSION METHOD AS
C IS USED IN OSIRIS AND BMD PACKAGES. INDIANA UNIVERSITY WAS
C ORIGINALLY RESPONSIBLE FOR CHANGE, AFTER TALKING TO SEVERAL WMU
C USERS THEY AGREED THAT THE MORE STANDARD FACTOR SCORE METHOD WAS
C DESIRABLE. REFERENCE ON FACTOR SCORE METHOD IS
C
SUBROUTINE STFACT(NV,NC,MV,MC,DATA,STD,VMN,COR,SP,H,NAMES)
DIMENSION DATA(MC,MV),COR(MV,MV),DCOR(40,40),DDC(40,40),SP(1)
DIMENSION H(1),TV(40),F(40),IV(40),OPT(7),NAME(7),VMN(1),STD(1)
DIMENSION NAMES(1),IVA(40)
COMMON /DEV/ICC,IDATA,IOUT,IDLG,IDSK
COMMON /PRNT/ LINPP,ICOPS,RUNPRG
COMMON/EXTRA/HEDR(70),NSZ
CON=1.0
9 DO 8 I=1,6
8 OPT(I)=0
18 IF(ICC.NE.2) WRITE(IDLG,16)
16 FORMAT('0LIST THE OPTIONS SEPARATED BY COMMAS'/)
READ(ICC,17) NAME
17 FORMAT(7(A5,1X))
IF(NAME(1).EQ.'!') RETURN
DO 20 I=1,7
IF(NAME(I).EQ.' ') GO TO 36
30 IF(NAME(I).NE.'HELP') GO TO 22
WRITE(IDLG,21)
21 FORMAT('0POSSIBLE OPTIONS:'//' "EVECT" - PRINT EIGEN VECTORS'/
1' "FACTM" - PRINT FACTOR MATRIX'/' "CONST" - SPECIFY ',
2'CONSTANT TO DECIDE HOW MANY EIGENVALUES TO RETAIN'/
3' "VARIA" - PRINT VARIANCES'/' "COMMU" - PRINT COMMUNALITIES'/
4' "FSCOR" - PRINT FACTOR SCORES'/' "SFSCR" - SAVE FACTOR SCORES'
5/' "ALL" - ALL ABOVE OPTIONS'/
6'0IF NO OPTIONS ARE DESIRED TYPE A CARRIAGE RETURN')
GO TO 18
22 IF(NAME(I).NE.'EVECT') GO TO 23
OPT(1)=1
GO TO 20
23 IF(NAME(I).NE.'FACTM') GO TO 24
OPT(2)=1
GO TO 20
24 IF(NAME(I).NE.'VARIA') GO TO 25
OPT(3)=1
GO TO 20
25 IF(NAME(I).NE.'COMMU') GO TO 26
OPT(4)=1
GO TO 20
26 IF(NAME(I).NE.'FSCOR') GO TO 27
OPT(5)=1
GO TO 20
27 IF(NAME(I).NE.'CONST') GO TO 9004
OPT(6)=1
GO TO 20
9004 IF(NAME(I).NE.'SFSCR') GO TO 28
OPT(7)=1
GO TO 20
28 IF(NAME(I).NE.'ALL') GO TO19
OPT(1)=1
OPT(2)=1
OPT(3)=1
OPT(4)=1
OPT(5)=1
OPT(6)=1
OPT(7)=1
GO TO 20
19 WRITE(IDLG,29) NAME(I)
29 FORMAT('0OPTION ',A5,' DOES NOT EXIST')
GO TO 9
20 CONTINUE
36 IF(OPT(6).NE.1) GO TO 1
IF(ICC.NE.2) WRITE(IDLG,33)
33 FORMAT('0WHAT CONSTANT VALUE WOULD YOU LIKE? ',$)
READ(ICC,34) CON
34 FORMAT(G)
1 IF(ICC.NE.2) WRITE(IDLG,2)
2 FORMAT('0TYPE IN THE VARIABLES, SEPARATED BY COMMAS'/)
CALL ALPHA(IVA,40,NVV,IRET,IHELP,IERR,NAMES,NV)
IF(IRET.EQ.1) RETURN
IF(IERR.EQ.1) GO TO 1
IF(IHELP.EQ.1) GO TO 1
K=1
DO 3 I=1,NVV
IV(I)=IVA(I)
IF(IVA(I).GT.0) GO TO 3
IV(I)=K
K=K+1
3 CONTINUE
DO 4 I=1,NVV-1
IF(IVA(I).LT.0) GO TO 4
DO 5 J=I+1,NVV
IF(IVA(J).LT.0) GO TO 5
IF(IVA(I).NE.IVA(J)) GO TO 5
WRITE(IDLG,6)
6 FORMAT(' THE SAME VARIABLE HAS BEEN USED TWICE')
GO TO 1
5 CONTINUE
4 CONTINUE
IF(NVV.LE.NV) GO TO 208
WRITE(IDLG,7)
7 FORMAT(' YOU HAVE LISTED MORE VARIABLES THAN IS POSSIBLE'/
1' WITH THE DATA SET AVAILABLE')
GO TO 1
200 J=NVV
204 IF(IVA(J).GT.0) GO TO 205
IV(J)=IV(J)+1
IF(IV(J).LE.NV) GO TO 206
205 J=J-1
IF(J.GE.1) GOTO 204
RETURN
206 K=IV(J)
IF(J.EQ.NVV) GO TO 208
DO 207 I=J+1,NVV
IF(IVA(I).GT.0) GO TO 207
K=K+1
IF(K.GT.NV) GO TO 205
IV(I)=K
207 CONTINUE
208 DO 203 I=1,NVV-1
DO 203 K=I+1,NVV
IF(IV(I).EQ.IV(K)) GO TO 200
203 CONTINUE
DO 15 I=1,NVV
II=IV(I)
DO 15 J=1,NVV
15 DCOR(I,J)=COR(II,IV(J))
35 IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566 FORMAT('1',70A1)
IF(IOUT.EQ.21) CALL PRNTHD
WRITE(IOUT,69)(NAMES(IV(I)),I=1,NVV)
LINES=5+(NVV+8)/9
CALL EIGEN(DCOR,DDC,NVV)
CALL TRACE(NVV,DCOR,CON,K,SP)
C PRINT EIGEN VALUES
IF(IOUT.NE.21) GO TO 400
LINES=LINES+2+(K+3)/4
IF(LINES.LE.LINPP) GO TO 400
CALL PRNTHD
LINES=4+(K+3)/4
400 WRITE(IOUT,40)(DCOR(I,I),I=1,K)
C PRINT OUT CUMMULATIVE PERCENTAGES OF EIGEN VALUES
IF(IOUT.NE.21) GO TO 401
LINES=LINES+(K+11)/4
IF(LINES.LE.LINPP) GO TO 401
CALL PRNTHD
LINES=4+(K+3)/4
401 WRITE(IOUT,41)(SP(J),J=1,K)
C PRINT EIGEN VECTORS (OPTIONAL)
IF(OPT(1).NE.1) GO TO 45
IF(IOUT.NE.21) GO TO 402
LINES=LINES+3
IF(LINES.LE.(LINPP-3*((NVV+7)/4))) GO TO 402
CALL PRNTHD
LINES=5
402 WRITE(IOUT,42)
DO 43 J=1,K
IF(IOUT.NE.21) GO TO 43
LINES=LINES+(NVV+7)/4
IF(LINES.LE.LINPP) GO TO 43
CALL PRNTHD
WRITE(IOUT,42)
LINES=4+(NVV+7)/4
43 WRITE(IOUT,44) J,(DDC(I,J),I=1,NVV)
45 CALL LOAD (NVV,K,DCOR,DDC)
C PRINT FACTOR MATRIC(OPTIONAL)
IF(OPT(2).NE.1) GO TO 49
IF(IOUT.NE.21) GO TO 403
LINES=LINES+3
IF(LINES.LE.(LINPP-((K+5)/3)*2)) GO TO 403
CALL PRNTHD
LINES=5
403 WRITE(IOUT,46) K
DO 47 I=1,NVV
IF(IOUT.NE.21) GO TO 47
LINES=LINES+(K+5)/3
IF(LINES.LE.LINPP) GO TO 47
CALL PRNTHD
WRITE(IOUT,46) K
LINES=5+(K+5)/3
47 WRITE(IOUT,48) NAMES(IV(I)),(DDC(I,J),J=1,K)
49 IF(K.GT.1) GO TO 52
IF(IOUT.NE.21) GO TO 404
LINES=LINES+2
IF(LINES.LE.LINPP) GO TO 404
CALL PRNTHD
LINES=4
404 WRITE(IOUT,50) K
GO TO 200
52 CALL VARMX(NVV,K,DDC,INC,TV,H,F,SP)
C PRINT VARIANCES (OPTIONAL)
IF(OPT(3).NE.1) GO TO 55
INV=INC+1
IF(IOUT.NE.21) GO TO 405
LINES=LINES+4
IF(LINES.LE.(LINPP-5)) GO TO 405
CALL PRNTHD
LINES=6
405 WRITE(IOUT,154)
DO 53 I=1,INV
INC=I-1
IF(IOUT.NE.21) GO TO 53
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 53
CALL PRNTHD
WRITE(IOUT,154)
LINES=7
53 WRITE(IOUT,54) INC,TV(I)
C PRINT ROTATED FACTOR MATRIX
55 IF(IOUT.NE.21) GO TO 406
LINES=LINES+3
IF(LINES.LE.(LINPP-((K+5)/3)*2)) GO TO 406
CALL PRNTHD
LINES=5
406 WRITE(IOUT,56) K
DO 57 I=1,NVV
IF(IOUT.NE.21) GO TO 57
LINES=LINES+(K+5)/3
IF(LINES.LE.LINPP) GO TO 57
CALL PRNTHD
LINES=3+(K+2)/3
57 WRITE(IOUT,48) NAMES(IV(I)),(DDC(I,J),J=1,K)
C PRINT COMMUNALITIES (OPTIONAL)
IF(OPT(4).NE.1) GO TO 62
IF(IOUT.NE.21) GO TO 407
LINES=LINES+5
IF(LINES.LE.(LINPP-5)) GO TO 407
CALL PRNTHD
LINES=7
407 WRITE(IOUT,59)
DO 60 I=1,NVV
IF(IOUT.NE.21) GO TO 60
LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 60
CALL PRNTHD
WRITE(IOUT,59)
LINES=8
60 WRITE(IOUT,61)NAMES(IV(I)),H(I),F(I),SP(I)
C PRINT OUT FACTOR SCORES (OPTIONAL)
62 IF((OPT(5).NE.1).AND.(OPT(7).NE.1)) GO TO 200
IF(OPT(7).NE.1) GO TO 9001
IF(NV+K.LE.MV) GO TO 9001
WRITE(IDLG,9002)
9002 FORMAT(//,' TOO MANY FACTOR SCORES TO SAVE')
OPT(7)=0
GO TO 62
9001 IF(OPT(5).NE.1) GO TO 408
IF(IOUT.NE.21) GO TO 409
LINES=LINES+3
IF(LINES.LE.(LINPP-((K+3)/4)*3)) GO TO 409
CALL PRNTHD
LINES=5
409 WRITE(IOUT,63)
63 FORMAT(1H-,'FACTOR SCORES')
408 DO 80 I=1,K
DO 80 J=1,K
DCOR(I,J)=0
DO 80 L=1,NVV
80 DCOR(I,J)=DDC(L,I)*DDC(L,J)+DCOR(I,J)
CALL INVT(DCOR,SP,TV,K)
DO 84 J=1,NVV
DO 85 I=1,K
H(I)=0
DO 85 L=1,K
85 H(I)=H(I)+DDC(J,L)*DCOR(L,I)
DO 86 L=1,K
86 DDC(J,L)=H(L)
84 CONTINUE
DO 65 L=1,NC
DO 9003 J=1,K
SP(J)=0
DO 66 I=1,NVV
66 SP(J)=SP(J)+((DATA(L,IV(I))-VMN(IV(I)))/STD(IV(I)))*DDC(I,J)
IF(OPT(7).EQ.1)DATA(L,NV+J)=SP(J)
9003 CONTINUE
IF(OPT(5).NE.1) GO TO 65
IF(IOUT.NE.21) GO TO 410
LINES=LINES+(K+3)/4
IF(LINES.LE.LINPP) GO TO 410
CALL PRNTHD
WRITE(IOUT,63)
LINES=5+(K+3)/4
410 WRITE(IOUT,67) L,(SP(I),I=1,K)
65 CONTINUE
67 FORMAT(1X,I4,4G/(5X,4G))
69 FORMAT('0',20X,'***** FACTOR ANALYSIS *****'/'0VARIABLES: '
1,9(A5,1X)/(12X,9(A5,1X)))
IF(OPT(7).NE.1) GO TO 200
DO 9006 I = NV+1,NV+K
9006 ENCODE(5,9005,NAMES(I))I
9005 FORMAT(I3,2X)
J=NV+1
NV=NV+K
IF(K.EQ.1)WRITE(IDLG,9008)NV
IF(K.NE.1)WRITE(IDLG,9007)J,NV
9007 FORMAT(//,' FACTOR SCORES SAVED AS VARIABLES 'I2,' THRU 'I2,/)
9008 FORMAT(//,' FACTOR SCORE SAVED AS VARIABLE 'I2,/)
70 DO 90 I=1,NV
VMN(I)=0.0
STD(I)=0.0
DO 90 J=1,NV
90 COR(I,J)=0.0
DO 71 I=1,NC
DO 71 J=1,NV
VMN(J)=VMN(J)+DATA(I,J)
DO 71 K=J,NV
71 COR(J,K)=COR(J,K)+DATA(I,J)*DATA(I,K)
DO 72 I=1,NV
DO 72 J=I,NV
72 COR(J,I)=NC*COR(I,J)-VMN(I)*VMN(J)
DO 73 I=1,NV
STD(I)=SQRT(COR(I,I)/(NC*(NC-1)))
73 VMN(I)=VMN(I)/NC
DO 74 I=1,NV
DO 74 J=I,NV
IF(I.EQ.J) GO TO 74
IF(COR(I,I)*COR(J,J).EQ.0) GO TO 75
COR(I,J)=COR(J,I)/SQRT(COR(I,I)*COR(J,J))
COR(J,I)=COR(I,J)
GO TO 74
75 COR(I,J)=0.0
COR(J,I)=0.0
74 CONTINUE
DO 76 I=1,NV
76 COR(I,I)=1.0
GO TO 200
40 FORMAT('0EIGEN VALUES'/(4G))
41 FORMAT('0CUMULATIVE PERCENTAGE OF EIGENVALUES'/(4G))
42 FORMAT(1H0/' EIGEN VECTORS')
44 FORMAT('0VECTOR',I3,4G/(10X,4G))
46 FORMAT(1H0/' FACTOR MATRIX (',I3,' FACTORS)')
48 FORMAT('0VARIABLE: ',A5,1X,3G/(17X,3G))
50 FORMAT('0ONLY ',I2,' FACTORS RETAINED. NO ROTATION')
154 FORMAT('0'/' ITERATION',7X,'VARIANCES'/' CYCLE')
54 FORMAT(I6,G20.6)
56 FORMAT(1H0/' ROTATED FACTOR MATRIX (',I3,' FACTORS)')
59 FORMAT(1H0/' CHECK ON COMMUNALITIES'//' VARIABLE',7X,
1'ORIGINAL',12X,'FINAL',10X,'DIFFERENCE')
61 FORMAT(1X,A5,1X,3F18.5)
END
C *** STAT PACK ***
C SUBROUTINE IS PART OF FACTOR ANALYSIS ROUTINE.
C CALLING SEQUENCE: CALL EIGEN(DCOR,DDC,NNV)
C WHERE DCOR - DIMENSIONED 40,40 CONTAINS CORRELATION MATRIIX.
C DDC - DIMENSIONED 40 BY 40
C NVV - NUMBER OF VARIABLES ACTUALLY USED.
C
SUBROUTINE EIGEN(DCOR,DDC,NVV)
DIMENSION DCOR(40,40),DDC(40,40)
REAL NORMIN,NORMAX
RANGE=1.0E-6
FNVV=NVV
DO 1 I=1,NVV
DO 1 J=1,NVV
DDC(I,J)=0
IF(I.EQ.J) DDC(I,J)=1.0
1 CONTINUE
C COMPUTE INITIAL AND FINAL NORMS(NORMIN,NORMAX)
2 NORMIN=0
DO 3 I=1,NVV
DO 3 J=I,NVV
IF(I.EQ.J) GO TO 3
NORMIN=NORMIN+DCOR(I,J)**2
3 CONTINUE
IF(NORMIN.LE.0) GO TO 14
NORMIN=SQRT(2.*NORMIN)
NORMAX=NORMIN*RANGE/FNVV
C INITIALIZE INDICATORS AND COMPUTE THRESHOLD
IND=0
THR=NORMIN
7 THR=THR/FNVV
4 L=1
5 M=L+1
C COMPUTE SIN AND COS
6 IF(ABS(DCOR(L,M)).LT.THR) GO TO 10
IND=1
X=.5*(DCOR(L,L)-DCOR(M,M))
Y=-DCOR(L,M)/SQRT(DCOR(L,M)**2+X**2)
IF(X.LT.0) Y=-Y
SINX=Y/SQRT(2.0*(1.0+SQRT(1.0-Y**2)))
SINX2=SINX**2
COSX=SQRT(1.0-SINX2)
COSX2=COSX**2
SINCS=SINX*COSX
C ROTATE L AND M COLUMNS
DO 8 I=1,NVV
IF(I.EQ.L) GO TO 9
IF(I.EQ.M) GO TO 9
X=DCOR(I,L)*COSX-DCOR(I,M)*SINX
DCOR(I,M)=DCOR(I,L)*SINX+DCOR(I,M)*COSX
DCOR(M,I)=DCOR(I,M)
DCOR(I,L)=X
DCOR(L,I)=X
9 X=DDC(I,L)*COSX-DDC(I,M)*SINX
DDC(I,M)=DDC(I,L)*SINX+DDC(I,M)*COSX
DDC(I,L)=X
8 CONTINUE
X=2.0*DCOR(L,M)*SINCS
Y=DCOR(L,L)*COSX2+DCOR(M,M)*SINX2-X
X=DCOR(L,L)*SINX2+DCOR(M,M)*COSX2+X
DCOR(L,M)=(DCOR(L,L)-DCOR(M,M))*SINCS+DCOR(L,M)*(COSX2-SINX2)
DCOR(M,L)=DCOR(L,M)
DCOR(L,L)=Y
DCOR(M,M)=X
C TEST FOR COMPLETION
C TEST FOR M=LAST ROW
10 IF(M.EQ.NVV) GO TO 11
M=M+1
GO TO 6
C TEST FOR L=SECOND FROM LAST COLUMN
11 IF(L.EQ.(NVV-1)) GO TO 12
L=L+1
GO TO 5
12 IF(IND.NE.1) GO TO 13
IND=0
GO TO 4
C COMPARE THRESHOLD WITH FINAL NORM
13 IF(THR.GT.NORMAX) GO TO 7
14 DO 15 I=1,NVV
DO 15 J=I+1,NVV
IF(DCOR(I,I).GE.DCOR(J,J)) GO TO 15
X=DCOR(I,I)
DCOR(I,I)=DCOR(J,J)
DCOR(J,J)=X
DO 16 K=1,NVV
S=DDC(K,I)
DDC(K,I)=DDC(K,J)
16 DDC(K,J)=S
15 CONTINUE
RETURN
END
C *** STAT PACK ***
C SUBROUTINE IS PART OF FACTOR ANALYSIS ROUTINE.
C CALLING SEQUENCE: CALL LOAD(NVV,K,DCOR,DDC)
C WHERE NVV - IS THE NUMBER OF VARIABLES
C K - NUMBER OF FACTORS LOADED
C DCOR - DIMENSIONED 40 BY 40
C DDC - DIMENSIONED 40 BY 40
C
SUBROUTINE LOAD(NVV,K,DCOR,DDC)
DIMENSION DCOR(40,40),DDC(40,40)
C COMPUTE A FACTOR MATRIX(LOADING) FROM EIGEN VALUES AND
C ASSOCIATED EIGEN VALUES. THIS SUBROUTINE NORMALLY OCCURS
DO 160 J=1,K
SQ=SQRT(DCOR(J,J))
DO 160 I=1,NVV
160 DDC(I,J)=DDC(I,J)*SQ
RETURN
END
C *** STAT PACK ***
C SUBROUTINE IS PART OF FACTOR ANALYSIS ROUTINE.
C CALLING SEQUENCE: CALL TRACE(NVV,DCOR,CON,K,SP)
C WHERE NVV - NUMBER OF VARIABLES.
C CON - CONSTANT VALUE CONCERNED WITH DETERMINING WHETHER
C A EIGEN VECTOR IS SAVED.
C K - NUMBER SAVED.
C SP - VECTOR DIMENSIONED AT LEAST NVV
C
SUBROUTINE TRACE(NVV,DCOR,CON,K,SP)
DIMENSION DCOR(40,40),SP(1)
FNVV=NVV
DO 100 I=1,NVV
100 SP(I)=DCOR(I,I)
K=0
C TEST WHETHER I-TH EIGEN VALUE IS GREATER THAN OR EQUAL TO THE
C CONSTANT
DO 110 I=1,NVV
IF(SP(I).LT.CON) GO TO 120
K=K+1
110 SP(I)=SP(I)/FNVV
C COMPUTE CUMULATIVE PERCENTAGE OF EIGENVALUES
120 DO 130 I=2,K
130 SP(I)=SP(I)+SP(I-1)
RETURN
END
C *** STAT PACK ***
C SUBROUTINE IS PART OF FACTOR ANALYSIS ROUTINES.
C CALLING SEQUENCE: CALL VARMX(NVV,K,DDC,INC,TV,H,F,SP)
C WHERE NVV - IS NUMBER OF VARIABLES.
C K - IS THE NUMBER SAVED.
C DDC - IS A MATRIX DIMENSIONED 40 BY 40
C INC - 1 SINGLE VARIABLE
C TV - A VECTOR DIMENSIONED AT 40
C H - A VECTOR DIMENSIONED AT 40
C F - A VECTOR DIMENSIONED AT 40
C SP - A VECTOR AT LEAST NVV
C
SUBROUTINE VARMX(NVV,K,DDC,INC,TV,H,F,SP)
DIMENSION DDC(40,40),TV(1),SP(1),H(1),F(1)
C INITIALIZATION
EPS=.00116
TVLT=0
LL=K-1
INV=1
INC=0
FNVV=NVV
FFN=FNVV*FNVV
CONS=.7071066
C CALCULATE ORIGINGL COMMUNALITIES
DO 1 I=1,NVV
H(I)=0
DO 1 J=1,K
1 H(I)=H(I)+DDC(I,J)*DDC(I,J)
C CALCULATE NORMALIZED FACTOR MATRIX
DO 2 I=1,NVV
H(I)=SQRT(H(I))
DO 2 J=1,K
2 DDC(I,J)=DDC(I,J)/H(I)
GO TO 4
C CALCULATE VARIANCE FOR FACTOR MATREX
3 INV=INV+1
TVLT=TV(INV-1)
4 TV(INV)=0
DO 6 J=1,K
AA=0
BB=0
DO 5 I=1,NVV
CC=DDC(I,J)*DDC(I,J)
AA=AA+CC
5 BB=BB+CC*CC
6 TV(INV)=TV(INV)+(FNVV*BB-AA**2)/FFN
IF(INV.GE.51) GO TO 21
C PERFORM CONVERGENCE TEST
IF((TV(INV)-TVLT).GT.(1.E-7)) GO TO 7
INC=INC+1
IF(INC.GT.3) GO TO 21
C ROTATION OF TWO FACTORS CONTINUES
7 DO 20 J=1,LL
C CALUCULATE NUM AND DEM
DO 20 K1=J+1,K
AA=0
BB=0
CC=0
DD=0
DO 8 I=1,NVV
U=(DDC(I,J)+DDC(I,K1))*(DDC(I,J)-DDC(I,K1))
T=DDC(I,J)*DDC(I,K1)
T=T+T
CC=CC+(U+T)*(U-T)
DD=DD+2.0*U*T
AA=AA+U
8 BB=BB+T
T=DD-2.0*AA*BB/FNVV
B=CC-(AA*AA-BB*BB)/FNVV
C COMPARISON OF NUM AND DEN
IF(T-B) 10,9,12
9 IF((T+B).LT.EPS) GO TO 20
C NUM+DEN IS GREATER THAN OR EQUAL TO THE
COS4T=CONS
C TOLERANCE FACTOR
SIN4T=CONS
GO TO 15
C NUM IS LESS THA DEN
10 TAN4T=ABS(T)/ABS(B)
IF(TAN4T.LT.EPS) GO TO 11
COS4T=1.0/SQRT(1.+TAN4T**2)
SIN4T=TAN4T*COS4T
GO TO 15
11 IF(B.GE.0) GO TO 20
SINP=CONS
COSP=CONS
GO TO 18
C NUM IS GREATER THAN DEN
12 CTN4T=ABS(T/B)
IF(CTN4T.LT.EPS) GO TO 13
SIN4T=1.0/SQRT(1.0+CTN4T**2)
COS4T=CTN4T*SIN4T
GO TO 15
13 COS4T=0
SIN4T=0
C DETERMINE COS THETA AND SINTHETA
15 COS2T=SQRT((1.+COS4T)/2.0)
SIN2T=SIN4T/(2.0*COS2T)
COST=SQRT((1.0+COS2T)/2.0)
SINT=SIN2T/(2.0*COST)
C DETERMINE COS PHI AND SIN PHI
IF(B.LE.0) GO TO 16
COSP=COST
SINP=SINT
GO TO 17
16 COSP=CONS*COST+CONS*SINT
SINP=ABS(CONS*COST-CONS*SINT)
17 IF(T.GT.0) GO TO 18
SINP=-SINP
C PERFORM ROTATION
18 DO 19 I=1,NVV
AA=DDC(I,J)*COSP+DDC(I,K1)*SINP
DDC(I,K1)=-DDC(I,J)*SINP+DDC(I,K1)*COSP
19 DDC(I,J)=AA
20 CONTINUE
GO TO 3
C DENORMALIZE VARIMAX LOADINGS
21 DO 22 I=1,NVV
DO 22 J=1,K
22 DDC(I,J)=DDC(I,J)*H(I)
C CHECK ON COMMUNALITIES
INC=INV-1
DO 23 I=1,NVV
23 H(I)=H(I)*H(I)
DO 24 I=1,NVV
F(I)=0
DO 25 J=1,K
25 F(I)=F(I)+DDC(I,J)**2
24 SP(I)=H(I)-F(I)
RETURN
END
C *** STAT PACK ***
C SUBROUTINE IS PART OF FACTOR ANALYSIS USED FOR FACTOR SCORES
C CALLING SEQUENCE: CALL INVT(A,L,M,N)
C WHERE: A - IS A MATRIX CONTAINING ROTATED FACTOR MATRIX TIMES THE
C TRANSPOSE. DIMENSIONED 40 X 40.
C : L - IS A VECTOR AT LEAST N LONG.
C : M - IS A VECTOR AT LEAST N LONG
C : N - NUMBER OF FACTORS AND SIZE OF MATRIX TO BE INVERTED
C
C METHOD USE TO INVERT MATRIX IN PLACE IS A TAKE OFF ON OSIRIS II.
C
SUBROUTINE INVT(A,L,M,N)
DIMENSION A(40,40),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