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

