File: VARMX.FT of Tape: Various/ETH/eth11-2
(Source file text)
C
C ..................................................................
C
C SUBROUTINE VARMX
C
C PURPOSE
C PERFORM ORTHOGONAL ROTATIONS OF A FACTOR MATRIX. THIS
C SUBROUTINE NORMALLY OCCURS IN A SEQUENCE OF CALLS TO SUB-
C ROUTINES CORRE, EIGEN, TRACE, LOAD, VARMX IN THE PERFORMANCE
C OF A FACTOR ANALYSIS.
C
C USAGE
C CALL VARMX (M,K,A,NC,TV,H,F,D,IER)
C
C DESCRIPTION OF PARAMETERS
C M - NUMBER OF VARIABLES AND NUMBER OF ROWS OF MATRIX A.
C K - NUMBER OF FACTORS.
C A - INPUT IS THE ORIGINAL FACTOR MATRIX, AND OUTPUT IS
C THE ROTATED FACTOR MATRIX. THE ORDER OF MATRIX A
C IS M X K.
C NC - OUTPUT VARIABLE CONTAINING THE NUMBER OF ITERATION
C CYCLES PERFORMED.
C TV - OUTPUT VECTOR CONTAINING THE VARIANCE OF THE FACTOR
C MATRIX FOR EACH ITERATION CYCLE. THE VARIANCE PRIOR
C TO THE FIRST ITERATION CYCLE IS ALSO CALCULATED.
C THIS MEANS THAT NC+1 VARIANCES ARE STORED IN VECTOR
C TV. MAXIMUM NUMBER OF ITERATION CYCLES ALLOWED IN
C THIS SUBROUTINE IS 50. THEREFORE, THE LENGTH OF
C VECTOR TV IS 51.
C H - OUTPUT VECTOR OF LENGTH M CONTAINING THE ORIGINAL
C COMMUNALITIES.
C F - OUTPUT VECTOR OF LENGTH M CONTAINING THE FINAL
C COMMUNALITIES.
C D - OUTPUT VECTOR OF LENGTH M CONTAINING THE DIFFERENCES
C BETWEEN THE ORIGINAL AND FINAL COMMUNALITIES.
C IER - ERROR INDICATOR
C IER=0 - NO ERROR
C IER=1 - CONVERGENCE WAS NOT ACHIEVED IN 50 CYCLES
C OF ROTATION
C
C REMARKS
C IF VARIANCE COMPUTED AFTER EACH ITERATION CYCLE DOES NOT
C INCREASE FOR FOUR SUCCESSIVE TIMES, THE SUBROUTINE STOPS
C ROTATION.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C KAISER'S VARIMAX ROTATION AS DESCRIBED IN 'COMPUTER PROGRAM
C FOR VARIMAX ROTATION IN FACTOR ANALYSIS' BY THE SAME AUTHOR,
C EDUCATIONAL AND PSYCHOLOGICAL MEASUREMENT, VOL XIX, NO. 3,
C 1959.
C
C ..................................................................
C
SUBROUTINE VARMX (M,K,A,NC,TV,H,F,D,IER)
DIMENSION A(1),TV(1),H(1),F(1),D(1)
C
C ...............................................................
C
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
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 2, DABS, DSQRT
C
C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C ROUTINE.
C
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 355 MUST BE CHANGED TO DSQRT. ABS IN
C STATEMENTS 280, 320, AND 375 MUST BE CHANGED TO DABS.
C
C ...............................................................
C
C INITIALIZATION
C
IER=0
EPS=0.00116
TVLT=0.0
LL=K-1
NV=1
NC=0
FN=M
FFN=FN*FN
CONS=0.7071066
C
C CALCULATE ORIGINAL COMMUNALITIES
C
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
C CALCULATE NORMALIZED FACTOR MATRIX
C
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
C CALCULATE VARIANCE FOR FACTOR MATRIX
C
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,155,155
155 IER=1
GO TO 430
C
C PERFORM CONVERGENCE TEST
C
160 IF((TV(NV)-TVLT)-(1.E-7)) 170, 170, 190
170 NC=NC+1
IF(NC-3) 190, 190, 430
C
C ROTATION OF TWO FACTORS CONTINUES UP TO
C THE STATEMENT 120.
C
190 DO 420 J=1,LL
L1=M*(J-1)
II=J+1
C
C CALCULATE NUM AND DEN
C
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
C COMPARISON OF NUM AND DEN
C
IF(T-B) 280, 240, 320
240 IF((T+B)-EPS) 420, 250, 250
C
C NUM + DEN IS GREATER THAN OR EQUAL TO THE
C TOLERANCE FACTOR
C
250 COS4T=CONS
SIN4T=CONS
GO TO 350
C
C NUM IS LESS THAN DEN
C
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
C NUM IS GREATER THAN DEN
C
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
C DETERMINE COS THETA AND SIN THETA
C
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
C DETERMINE COS PHI AND SIN PHI
C
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
C PERFORM ROTATION
C
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
C DENORMALIZE VARIMAX LOADINGS
C
430 DO 440 I=1,M
DO 440 J=1,K
L=M*(J-1)+I
440 A(L)=A(L)*H(I)
C
C CHECK ON COMMUNALITIES
C
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