home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / statfctr / varmx.for < prev   
Encoding:
Text File  |  1985-11-29  |  6.5 KB  |  244 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE VARMX
  5. C
  6. C        PURPOSE
  7. C           PERFORM ORTHOGONAL ROTATIONS OF A FACTOR MATRIX.  THIS
  8. C           SUBROUTINE NORMALLY OCCURS IN A SEQUENCE OF CALLS TO SUB-
  9. C           ROUTINES CORRE, EIGEN, TRACE, LOAD, VARMX IN THE PERFORMANCE
  10. C           OF A FACTOR ANALYSIS.
  11. C
  12. C        USAGE
  13. C           CALL VARMX (M,K,A,NC,TV,H,F,D,IER)
  14. C
  15. C        DESCRIPTION OF PARAMETERS
  16. C           M     - NUMBER OF VARIABLES AND NUMBER OF ROWS OF MATRIX A.
  17. C           K     - NUMBER OF FACTORS.
  18. C           A     - INPUT IS THE ORIGINAL FACTOR MATRIX, AND OUTPUT IS
  19. C                   THE ROTATED FACTOR MATRIX.  THE ORDER OF MATRIX A
  20. C                   IS M X K.
  21. C           NC    - OUTPUT VARIABLE CONTAINING THE NUMBER OF ITERATION
  22. C                   CYCLES PERFORMED.
  23. C           TV    - OUTPUT VECTOR CONTAINING THE VARIANCE OF THE FACTOR
  24. C                   MATRIX FOR EACH ITERATION CYCLE.  THE VARIANCE PRIOR
  25. C                   TO THE FIRST ITERATION CYCLE IS ALSO CALCULATED.
  26. C                   THIS MEANS THAT NC+1 VARIANCES ARE STORED IN VECTOR
  27. C                   TV.  MAXIMUM NUMBER OF ITERATION CYCLES ALLOWED IN
  28. C                   THIS SUBROUTINE IS 50.  THEREFORE, THE LENGTH OF
  29. C                   VECTOR TV IS 51.
  30. C           H     - OUTPUT VECTOR OF LENGTH M CONTAINING THE ORIGINAL
  31. C                   COMMUNALITIES.
  32. C           F     - OUTPUT VECTOR OF LENGTH M CONTAINING THE FINAL
  33. C                   COMMUNALITIES.
  34. C           D     - OUTPUT VECTOR OF LENGTH M CONTAINING THE DIFFERENCES
  35. C                   BETWEEN THE ORIGINAL AND FINAL COMMUNALITIES.
  36. C           IER   - ERROR INDICATOR
  37. C                   IER=0 - NO ERROR
  38. C                   IER=1 - CONVERGENCE WAS NOT ACHIEVED IN 50 CYCLES
  39. C                           OF ROTATION
  40. C
  41. C        REMARKS
  42. C           IF VARIANCE COMPUTED AFTER EACH ITERATION CYCLE DOES NOT
  43. C           INCREASE FOR FOUR SUCCESSIVE TIMES, THE SUBROUTINE STOPS
  44. C           ROTATION.
  45. C
  46. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  47. C           NONE
  48. C
  49. C        METHOD
  50. C           KAISER'S VARIMAX ROTATION AS DESCRIBED IN 'COMPUTER PROGRAM
  51. C           FOR VARIMAX ROTATION IN FACTOR ANALYSIS' BY THE SAME AUTHOR,
  52. C           EDUCATIONAL AND PSYCHOLOGICAL MEASUREMENT, VOL XIX, NO. 3,
  53. C           1959.
  54. C
  55. C     ..................................................................
  56. C
  57.       SUBROUTINE VARMX (M,K,A,NC,TV,H,F,D,IER)
  58.       DIMENSION A(1),TV(1),H(1),F(1),D(1)
  59. C
  60. C        ...............................................................
  61. C
  62. C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
  63. C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
  64. C
  65. C     DOUBLE PRECISION A,TV,H,F,D,TVLT,CONS,AA,BB,CC,DD,U,T,B,COS4T,
  66. C    1                 SIN4T,TAN4T,SINP,COSP,CTN4T,COS2T,SIN2T,COST,SINT
  67. C
  68. C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
  69. C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
  70. C        ROUTINE.
  71. C
  72. C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
  73. C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  SQRT IN STATEMENTS
  74. C        115, 290, 330, 350, AND 355 MUST BE CHANGED TO DSQRT.  ABS IN
  75. C        STATEMENTS 280, 320, AND 375 MUST BE CHANGED TO DABS.
  76. C
  77. C        ...............................................................
  78. C
  79. C     INITIALIZATION
  80. C
  81.       IER=0
  82.       EPS=0.00116
  83.       TVLT=0.0
  84.       LL=K-1
  85.       NV=1
  86.       NC=0
  87.       FN=M
  88.       FFN=FN*FN
  89.       CONS=0.7071066
  90. C
  91. C     CALCULATE ORIGINAL COMMUNALITIES
  92. C
  93.       DO 110 I=1,M
  94.       H(I)=0.0
  95.       DO 110 J=1,K
  96.       L=M*(J-1)+I
  97.   110 H(I)=H(I)+A(L)*A(L)
  98. C
  99. C     CALCULATE NORMALIZED FACTOR MATRIX
  100. C
  101.       DO 120 I=1,M
  102.   115 H(I)= SQRT(H(I))
  103.       DO 120 J=1,K
  104.       L=M*(J-1)+I
  105.   120 A(L)=A(L)/H(I)
  106.       GO TO 132
  107. C
  108. C     CALCULATE VARIANCE FOR FACTOR MATRIX
  109. C
  110.   130 NV=NV+1
  111.       TVLT=TV(NV-1)
  112.   132 TV(NV)=0.0
  113.       DO 150 J=1,K
  114.       AA=0.0
  115.       BB=0.0
  116.       LB=M*(J-1)
  117.       DO 140 I=1,M
  118.       L=LB+I
  119.       CC=A(L)*A(L)
  120.       AA=AA+CC
  121.   140 BB=BB+CC*CC
  122.   150 TV(NV)=TV(NV)+(FN*BB-AA*AA)/FFN
  123.       IF(NV-51)160,155,155
  124.   155 IER=1
  125.       GO TO 430
  126. C
  127. C     PERFORM CONVERGENCE TEST
  128. C
  129.   160 IF((TV(NV)-TVLT)-(1.E-7)) 170, 170, 190
  130.   170 NC=NC+1
  131.       IF(NC-3) 190, 190, 430
  132. C
  133. C     ROTATION OF TWO FACTORS CONTINUES UP TO
  134. C     THE STATEMENT 120.
  135. C
  136.   190 DO 420 J=1,LL
  137.       L1=M*(J-1)
  138.       II=J+1
  139. C
  140. C        CALCULATE NUM AND DEN
  141. C
  142.       DO 420 K1=II,K
  143.       L2=M*(K1-1)
  144.       AA=0.0
  145.       BB=0.0
  146.       CC=0.0
  147.       DD=0.0
  148.       DO 230 I=1,M
  149.       L3=L1+I
  150.       L4=L2+I
  151.       U=(A(L3)+A(L4))*(A(L3)-A(L4))
  152.       T=A(L3)*A(L4)
  153.       T=T+T
  154.       CC=CC+(U+T)*(U-T)
  155.       DD=DD+2.0*U*T
  156.       AA=AA+U
  157.   230 BB=BB+T
  158.       T=DD-2.0*AA*BB/FN
  159.       B=CC-(AA*AA-BB*BB)/FN
  160. C
  161. C        COMPARISON OF NUM AND DEN
  162. C
  163.       IF(T-B) 280, 240, 320
  164.   240 IF((T+B)-EPS) 420, 250, 250
  165. C
  166. C        NUM + DEN IS GREATER THAN OR EQUAL TO THE
  167. C        TOLERANCE FACTOR
  168. C
  169.   250 COS4T=CONS
  170.       SIN4T=CONS
  171.       GO TO 350
  172. C
  173. C        NUM IS LESS THAN DEN
  174. C
  175.   280 TAN4T= ABS(T)/ ABS(B)
  176.       IF(TAN4T-EPS) 300, 290, 290
  177.   290 COS4T=1.0/ SQRT(1.0+TAN4T*TAN4T)
  178.       SIN4T=TAN4T*COS4T
  179.       GO TO 350
  180.   300 IF(B) 310, 420, 420
  181.   310 SINP=CONS
  182.       COSP=CONS
  183.       GO TO 400
  184. C
  185. C        NUM IS GREATER THAN DEN
  186. C
  187.   320 CTN4T= ABS(T/B)
  188.       IF(CTN4T-EPS) 340, 330, 330
  189.   330 SIN4T=1.0/ SQRT(1.0+CTN4T*CTN4T)
  190.       COS4T=CTN4T*SIN4T
  191.       GO TO 350
  192.   340 COS4T=0.0
  193.       SIN4T=1.0
  194. C
  195. C        DETERMINE COS THETA AND SIN THETA
  196. C
  197.   350 COS2T= SQRT((1.0+COS4T)/2.0)
  198.       SIN2T=SIN4T/(2.0*COS2T)
  199.   355 COST= SQRT((1.0+COS2T)/2.0)
  200.       SINT=SIN2T/(2.0*COST)
  201. C
  202. C        DETERMINE COS PHI AND SIN PHI
  203. C
  204.       IF(B) 370, 370, 360
  205.   360 COSP=COST
  206.       SINP=SINT
  207.       GO TO 380
  208.   370 COSP=CONS*COST+CONS*SINT
  209.   375 SINP= ABS(CONS*COST-CONS*SINT)
  210.   380 IF(T) 390, 390, 400
  211.   390 SINP=-SINP
  212. C
  213. C        PERFORM ROTATION
  214. C
  215.   400 DO 410 I=1,M
  216.       L3=L1+I
  217.       L4=L2+I
  218.       AA=A(L3)*COSP+A(L4)*SINP
  219.       A(L4)=-A(L3)*SINP+A(L4)*COSP
  220.   410 A(L3)=AA
  221.   420 CONTINUE
  222.       GO TO 130
  223. C
  224. C     DENORMALIZE VARIMAX LOADINGS
  225. C
  226.   430 DO 440 I=1,M
  227.       DO 440 J=1,K
  228.       L=M*(J-1)+I
  229.   440 A(L)=A(L)*H(I)
  230. C
  231. C     CHECK ON COMMUNALITIES
  232. C
  233.       NC=NV-1
  234.       DO 450 I=1,M
  235.   450 H(I)=H(I)*H(I)
  236.       DO 470 I=1,M
  237.       F(I)=0.0
  238.       DO 460 J=1,K
  239.       L=M*(J-1)+I
  240.   460 F(I)=F(I)+A(L)*A(L)
  241.   470 D(I)=H(I)-F(I)
  242.       RETURN
  243.       END
  244.