Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-05 - decus/20-0137/bmd/bmd03m.for
There is 1 other file named bmd03m.for in the archive. Click here to see a list.
C        FACTOR ANALYSIS - MAIN PROGRAM                MAY  2, 1966
C        THIS IS A SIFTED VERSION OF BMD03M ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
C        IT WAS THEN CONVERTED TO 360 FORTRAN IV (H-LEVEL)
      COMMON  SUMXX
      COMMON  A      , R      , MM     , SUMDV  , C
      DIMENSION A(60,60),P(60,60),R(60,60),MM(60),SUMXX(60,60),
     1 SUMDV(108),C(60)
      DIMENSION LL(60)
      DIMENSION TC(255),XM(60),SD(60),FFMT(18)
      EQUIVALENCE(SUMXX,P)
      DOUBLE PRECISION A123,B123,C123,TODE,RNN,COMMUN
C
  412 FORMAT(39H1BMD03M - FACTOR ANALYSIS - VERSION OF   
     X18HMAY  2, 1966      /
     135H HEALTH SCIENCES COMPUTING FACILITY//
     214H PROBLEM CODE A6,/
     318H NO. OF VARIABLES I3,/
     414H NO. OF CASES I5,//)
C
      MTAPE =5
	CALL USAGEB('BMD03M')
      DATA A123/6HFINISH/
      DATA B123/6HPROBLM/
      DATA C123/6HCOMMUN/
 10   READ (5,900)TODE,RNN,NV,N,INTYPE,IIDIAG,NRO,CNFLAM,ITF,NFP,NTAPE,N
     1N
      NTAPE2=NTAPE
      LO=0
      NJ=1
      IF(TODE .EQ. A123) GO TO 401
      GO TO 400
  402 WRITE (6,404) TODE
 401  IF(NTAPE-5)371,371,370
 370  REWIND NTAPE
 371  CALL EXIT
      STOP
 400  IF(TODE .NE. B123) GO TO 402
 403  CALL TPWD(NTAPE,MTAPE)
      IF(NTAPE)476,931,405
  476 N12=IABS(NTAPE)
      NJ=2
      GO TO 307
  405 IF(NN.GT.0.AND.NN.LE.6)GO TO 408
      WRITE (6,4000)
      NN=1
 408  NN=NN*18
      READ (5,901)(SUMDV(I),I=1,NN)
  307 IF((NV-1)*(NV-81))406,933,933
  406 IF(INTYPE*(INTYPE-4)) 407,935,935
  407 IF(IIDIAG-5)20,937,937
   20 IF (INTYPE.EQ.3.AND.IIDIAG.NE.0) GO TO 6732
      IF (INTYPE.NE.3.AND.IIDIAG.EQ.0) GO TO 6732
      ON=N
      IF(-ITF)7192,7191,7191
 7192 IF(INTYPE-1) 7193,7191,7193
 7193 WRITE (6,7194)
      ITF=0
      READ  (5,900) XAX
 7191 ONN=ON-1.
      WRITE (6,412)RNN,NV,N
      WRITE (6,930) (SUMDV(I),I=1,NN)
      IF(N*(N-NV))550,554,554
 550  IF(IIDIAG-2)554,552,554
 552  IIDIAG=1
 554  GO TO (30,35,40),INTYPE
   30 DO 25 I=1,NV
      XM(I)=0.
 553  DO 25 J=1,NV
 25   SUMXX(I,J)=0.0
      IF(ITF) 3240,3240,3241
 3241 READ (5,901)FFMT
      REWIND 3
 3240 IF(N-9999)49,49,939
 49   IF(N)939,939,50
 50   DO 100 I=1,N
      GO TO(55,60),NJ
   55 READ (NTAPE,SUMDV)(C(J),J=1,NV)
      GO TO 80
   60 READ (N12)(C(J),J=1,NV)
 80   H=I
      H1=H*(H-1.)
      DO 85 J=1,NV
      SD(J)=(C(J)-XM(J))/H
      XM(J)=XM(J)+SD(J)
      DD=SD(J)*H1
      DO 85 K=1,J
 85   SUMXX(K,J)=SUMXX(K,J)+DD*SD(K)
      IF(ITF) 100,100,3280
 3280 DO 3202 J=1,NV
      IF(LO-255) 3200,3201,3201
 3201 WRITE (3) TC
      LO=0
 3200 LO=LO+1
 3202 TC(LO)=C(J)
 100  CONTINUE
      IF(ITF) 3281,3281,3283
 3283 WRITE (3) TC
      ENDFILE 3
      REWIND 3
 3281 NJ=NV-1
      DO 103 I=1,NJ
      NK=I+1
      DO 103 J=NK,NV
 103  SUMXX(J,I)=SUMXX(I,J)
      WRITE (6,912)
      WRITE (6,913) (XM(I), I=1,NV)
      DO 135 I=1,NV
 135  SD(I) = SQRT(SUMXX(I,I)/ONN)
      WRITE (6,914)
      WRITE (6,913) (SD(I), I=1,NV)
      DO 140 I=1,NV
      DO 140 J=I,NV
 140  R(I,J)=SUMXX(I,J)
      DO 141 I=1,NJ
      NK=I+1
      DO 141 J=NK,NV
 141  R(J,I)=R(I,J)
      NN=NV-1
 215  DO 220 I=1,NV
 220  SUMDV(I)=SQRT(R(I,I))
      DO 225 I=1,NV
      DO 225 J=I,NV
 225  P(I,J)=R(I,J)/(SUMDV(I)*SUMDV(J))
      DO 226 I=1,NN
      NK=I+1
      DO 226 J=NK,NV
 226  P(J,I)=P(I,J)
      GO TO 227
   35 DO 142 I=1,NV
      GO TO(145,150),NJ
  145 READ (NTAPE,SUMDV)(P(I,J),J=1,NV)
      GO TO 142
  150 READ (N12)(P(I,J),J=1,NV)
  142 CONTINUE
 227  IF(ITF) 7123,7123,7234
 7123 IF(IIDIAG-2) 7345,7234,7345
 7234 DO 200 I=1,NV
      DO 200 J=1,NV
 200  R(I,J) = P(I,J)
      CALL MULTR(R,NV,SUMDV)
 7345 GO TO (160,180,185,190),IIDIAG
 180  DO 230 I=1,NV
 230  P(I,I)=SUMDV(I)
      GO TO 160
 185  DO 110 I=1,NV
      P(I,I)=0.0
      DO 110 J=1,NV
      IF(P(I,I)-ABS(P(I,J)))115,110,110
 115  P(I,I)=ABS(P(I,J))
 110  CONTINUE
      GO TO 160
  950 FORMAT(A6,11F6.0)
  190 DO 8888 JJ=1,NV,11
      KK=MIN0(JJ+10,NV)
      READ (5,950)COMMUN,(P(J,J),J=JJ,KK)
      IF(COMMUN .NE. C123) GO TO 402
 8888 CONTINUE
  160 WRITE (6,915)
      DO 238 I=1,NV
      WRITE (6,916)I
 238  WRITE (6,917)(P(I,J), J=1,NV)
 241  CALL JACOBI(P,A,NV,1.E-7,1,NROT)
      LLL=0
      NZ=0
      DO 243 I=1,NV
      SUMDV(I)=P(I,I)
      DO 1724 J=1,NV
 1724 P(I,J)=A(I,J)
      IF(SUMDV(I)) 2421,2422,2422
 2421 NZ=NZ+1
      GO TO 242
 2422 IF(SUMDV(I)-CNFLAM) 242,243,243
 242  LLL=LLL+1
 243  CONTINUE
      NZ=NV-NZ
  224 IF (NZ)10,10,500
  500 WRITE (6,918)
  918 FORMAT(12H0EIGENVALUES//)
  920 FORMAT(13H0EIGENVECTORS//)
      WRITE (6,917)(SUMDV(I),I=1,NV)
      IF(NZ-NV)1775,1774,1775
 1775 WRITE (6,1789)
 1789 FORMAT(95H0ONLY THE POSITIVE EIGENVALUES AND ASSOCIATED VECTORS WI
     1LL BE USED IN THE FOLLOWING COMPUTATION)
 1774 SUM=.0
      FNV=NV
      DO 250 I=1,NZ
      SUM=SUM+SUMDV(I)
  250 A(I,1)= SUM  /FNV
      WRITE (6,919)
      WRITE (6,917)(A(I,1), I=1,NZ)
      WRITE (6,920)
      DO 255 J=1,NZ
      WRITE (6,921)J
 255  WRITE (6,917)(P(I,J), I=1,NV)
      DO 260 I=1,NZ
      SUMDV(I)=SQRT(SUMDV(I))
      DO 260 J=1,NV
 260  A(J,I)=SUMDV(I)*P(J,I)
      WRITE (6,922)
      DO 256 I3=1,NZ
      NMINUS = 0
      DO 252 I1=1,NV
      IF(A(I1,I3)) 251,252,252
 251  NMINUS = NMINUS + 1
 252  CONTINUE
      IF(NV - 2 * NMINUS) 253,256,256
 253  DO 254 I2=1,NV
 254  A(I2,I3) = -A(I2,I3)
 256  CONTINUE
      DO 265 I=1,NV
      WRITE (6,923)I
 265  WRITE (6,917)(A(I,J), J=1,NZ)
      IF(NV-NZ) 269, 269, 266
 266  DO 267 I=1,NZ
      DO 267 J=1,NZ
      P(I,J)=0.0
      DO 267 K=1,NV
 267  P(I,J)=P(I,J)+A(K,I)*A(K,J)
      WRITE (6,911)
      WRITE (6,926)
      GO TO 278
 269  DO 270 I=1,NV
      DO 270 J=1,NV
      P(I,J)=0.0
      DO 270 K=1,NV
 270  P(I,J)=P(I,J)+A(I,K)*A(J,K)
      WRITE (6,911)
      WRITE (6,924)
 278  DO 275 I=1,NZ
      WRITE (6,916)I
 275  WRITE (6,917) (P(I,J), J=1,NZ)
      NZ=NV-LLL
      NRO=MIN0(NRO,NZ)
      IF(NRO*(1-NRO)) 325,320,320
 320  WRITE (6,911)
      WRITE (6,927)
      GO TO 10
   40 DO 282 I=1,NV
       GO TO (285,295),NJ
  285 READ (NTAPE,SUMDV)(A(I,J),J=1,NRO)
      GO TO 282
  295 READ (N12)(A(I,J),J=1,NRO)
  282 CONTINUE
      WRITE (6,922)
      DO 305 I=1,NV
      WRITE (6,923)I
  305 WRITE (6,917)(A(I,J),J=1,NRO)
      IF(NRO-NV)325,325,941
 325  CALL ROTATE(NV,NRO)
      IF(-ITF)3212,10,10
 3212 IF(NFP) 3310,3310,3311
 3310 NFP=NRO
 3311 DO 5110 I=1,NRO
      DO 5111 K=1,NV
      SUMDV(K)=A(K,I)
 5111 A(K,I)=0.0
      DO 5110 J=1,NV
      DO 5110 K=1,NV
 5110 A(J,I)=A(J,I)+R(J,K)*SUMDV(K)
      WRITE (6,3314) ITF,FFMT
 3314 FORMAT(47H0FACTOR SCORES ARE COMPUTED AND WRITTEN ON TAPEI3,31H UN
     1DER THE FORMAT PRINTED BELOW/1X,18A4)
      LO=255
      DO 3218 L=1,N
      DO 3216 I=1,NV
      IF(LO-255) 3213,2314,2314
 2314 READ (3) TC
  300 FORMAT (20G)
      LO=0
 3213 LO=LO+1
 3216 C(I)=(TC(LO)-XM(I))/SD(I)
      DO 3217 J=1,NRO
      SUMDV(J)=0.0
      DO 3217 I=1,NV
 3217 SUMDV(J)=SUMDV(J)+A(I,J)*C(I)
 3218 WRITE (ITF,FFMT)L,(SUMDV(J),J=1,NFP)
 3252 IF(ITF-6)7771,10,7771
 7771 ENDFILE ITF
      REWIND ITF
      GO TO 10
 900  FORMAT(2A6,I2,I4,2I1,I2,F6.0,2I2,35X,I3,I2)
  404 FORMAT (' PROBLM CARD EXPECTED. CARD READ WAS ',A6)
 901  FORMAT(18A4)
 902  FORMAT(12F6.0)
 903  FORMAT(36I2)
 904  FORMAT(23H1FACTOR ANALYSIS  (CASEI2,1H)/12H PROBLEM NO.A6)
 906  FORMAT(20H0NUMBER OF VARIABLESI4/12H SAMPLE SIZE6X,I6)
 908  FORMAT(20H0DATA TRANSFORMATION50I2)
 909  FORMAT(20X,50I2)
 910  FORMAT(23H0NO DATA TRANSFORMATION)
 911  FORMAT(1H0)
 912  FORMAT(1H0//6H MEANS)
 913  FORMAT(8F15.5)
 914  FORMAT(1H0//20H STANDARD DEVIATIONS)
 915  FORMAT(1H0//25H CORRELATION COEFFICIENTS)
 916  FORMAT(4H0ROWI3)
 917  FORMAT(10F12.5)
 919  FORMAT(40H0CUMULATIVE PROPORTION OF TOTAL VARIANCE)
 921  FORMAT(7H0VECTORI3)
 922  FORMAT(1H0//14H FACTOR MATRIX)
 923  FORMAT(9H0VARIABLEI3)
 924  FORMAT(47H0FACTOR CHECK MATRIX (CORRELATION COEFFICIENTS))
 925  FORMAT(1H09X,62HCOMPUTATIONAL ACCURACY NOT SUFFICIENT FOR ANY MORE
     1 EIGENVALUES)
 926  FORMAT(45H0FACTOR CHECK MATRIX (EIGENVALUES DIAGONALLY))
 927  FORMAT(61H0THE FACTOR MATRIX CANNOT BE ROTATED.  CHANGE THE CODE N
     1UMBER)
 929  FORMAT(61H0VALUE SPECIFIED TO LIMIT THE NUMBER OF FACTORS TO BE RO
     1TATEDF9.4)
  930 FORMAT (' VARIABLE FORMAT IS '/(5X,18A4))
  931 WRITE (6,932) NTAPE2
  932 FORMAT (1X,I4,' IS ILLEGAL TAPE UNIT')
      GO TO 401
  933 WRITE (6,934) NV
  934 FORMAT (1X,I3,' VARIABLES IS ILLEGAL')
      GO TO 401
  935 WRITE (6,936) INTYPE
  936 FORMAT (1X,I2,' IS ILLEGAL INPUT CODE (COL. 19)')
      GO TO 401
  937 WRITE (6,938) IDIAG
  938 FORMAT (1X,I2,' IS ILLEGAL DIAGONAL ELEMENT SPECIFICATION CODE (CO
     1L. 20)')
      GO TO 401
  939 WRITE (6,940) N
  940 FORMAT (1X,I5,' IS ILLEGAL NUMBER OF CASES')
      GO TO 401
  941 WRITE (6,942) NRO
  942 FORMAT (1X,I3,' IS ILLEGAL NUMBER OF FACTORS TO BE ROTATED')
      GO TO 401
 6732 WRITE (6,943)
  943 FORMAT (' COLUMNS 19 AND 20 OF PROBLM CARD ARE CONTRIDICTORY')
      GO TO 401
 7194 FORMAT(56H0FACTOR SCORES CAN BE COMPUTED ONLY IF INPUT IS RAW DATA
     X)
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
      END
      SUBROUTINE MULTR(A,N,U)
C        SUBROUTINE MULTR FOR BMD03M                   MAY  2, 1966
      DIMENSION V(60),U(60),A(60,60)
      DO 1 I=1,N
      IF(A(I,I)) 1,1,16
 16   K=I
 1    V(I)=A(I,I)
 2    DO 3 I=1,K
      U(I)=A(K,I)
 3    A(K,I)=0.
      P=-U(K)
      DO 4 I=K,N
      U(I)=A(I,K)
 4    A(I,K)=0.
      U(K)=-1.
      V(K)=-V(K)
      H=0.
      DO 8 I=1,N
      Y=U(I)/P
      DO 5 J=1,I
 5    A(I,J)=A(I,J)+U(J)*Y
      IF(V(I)) 8,8,6
 6    IF(A(I,I)-H) 8,7,7
 7    H=A(I,I)
      K=I
 8    CONTINUE
      IF(H- .00001) 9,9,2
 9    DO 10 I=1,N
      DO 10 J=1,I
      A(I,J)=-A(I,J)
 10   A(J,I)=A(I,J)
      DO 15 K=1,N
      IF(V(K)) 11,14,14
 11   DO 13 I=1,N
      IF(V(I)) 13,12,12
 12   IF(A(I,I)-A(I,K)*A(I,K)/A(K,K)+.00001) 14,14,13
 13   CONTINUE
      U(K)=1.-1./A(K,K)
      GO TO 15
 14   U(K)=-V(K)
 15   CONTINUE
      DO 17 I=1,N
      DO 17 J=1,I
      IF(V(I)) 18,19,19
 18   IF(V(J)) 17,19,19
 19   A(I,J)=0.
      A(J,I)=0.
 17   CONTINUE
      RETURN
      END
      SUBROUTINE JACOBI(A,B,N,ACC,IV,NR)
C        SUBROUTINE JACOBI FOR BMD03M                  MAY  2, 1966
      DIMENSION A(60,60),B(60,60),LK(60),Q(60)
      IF(IV-1)200,201,200
 201  DO 202 I=1,N
      DO 203 J=1,N
 203  B(I,J)=0.0
 202  B(I,I)=1.0
 200  NR=0
      Q(1)=0.0
      W=0.0
      H=.5*A(1,1)*A(1,1)
      DO 1 I=2,N
      H=H+.5*A(I,I)*A(I,I)
      Q(I)=0.0
      I1=I-1
      DO 2 J=1,I1
      Z=ABS(A(I,J))
      H=H+Z*Z
      IF(Z-Q(I))2,2,3
 3    Q(I)=Z
      LK(I)=J
 2    CONTINUE
      IF(Q(I)-W)1,1,4
 4    W=Q(I)
      III=I
 1    CONTINUE
      H=ACC*SQRT(2.0*H)/FLOAT(N)
 30   II=LK(III)
      JJ=III
      X=A(II,II)
      Y=A(JJ,II)
      Z=A(JJ,JJ)
      W=X-Z
      T=.5*(W+SQRT(W*W+4.0*Y*Y))/Y
      W=SQRT(1.0+T*T)
      S=T/W
      C=1.0/W
      CC=C*C
      SS=S*S
      SC=S*C*2.0
      Q1=0.0
      Q2=0.0
      W=0.0
      NR=NR+1
      DO 27 I=1,N
      IF(I-II)10,11,12
 10   U=A(II,I)
      V=A(JJ,I)
      E=U*S+V*C
      A(II,I)=E
      IF(ABS(E)-Q1)15,15,14
 14   Q1=ABS(E)
      I1=I
 15   F=V*S-U*C
      A(JJ,I)=F
      IF(ABS(F)-Q2)9,9,16
 16   Q2=ABS(F)
      I2=I
      GO TO 9
 11   A(II,I)=SS*X+SC*Y+CC*Z
      Q(I)=Q1
      LK(I)=I1
      GO TO 9
 12   IF(I-JJ)17,18,19
 17   U=A(I,II)
      V=A(JJ,I)
      E=S*U+C*V
      A(I,II)=E
      IF(ABS(E)-Q(I))15,15,21
 21   LK(I)=II
      Q(I)=ABS(E)
      GO TO 15
 18   A(JJ,I)=CC*X-SC*Y+SS*Z
      A(I,II)=0.0
      Q(I)=Q2
      LK(I)=I2
      GO TO 9
 19   U=A(I,II)
      V=A(I,JJ)
      E=U*S+V*C
      F=V*S-U*C
      A(I,II)=E
      A(I,JJ)=F
      G=AMAX1(ABS(E),ABS(F))
      IF(G-Q(I))9,9,13
 13   Q(I)=G
      IF(ABS(E)-ABS(F))23,24,24
 24   LK(I)=II
      GO TO 9
 23   LK(I)=JJ
 9    IF(Q(I)-W)40,25,25
 25   W=Q(I)
      III=I
 40   IF(IV)27,27,33
 33   U=B(I,II)
      V=B(I,JJ)
      B(I,II)=U*S+V*C
      B(I,JJ)=V*S-U*C
 27   CONTINUE
      IF(W-H)31,31,30
 31   IF(IV)55,55,56
 56   DO 50 I=1,N
      U=-1.E20
      DO 51 J=I,N
      IF(A(J,J)-U)51,51,52
 52   U=A(J,J)
      K=J
 51   CONTINUE
      IF(K-I)53,50,53
 53   A(K,K)=A(I,I)
      A(I,I)=U
      DO 54 J=1,N
      U=B(J,K)
      B(J,K)=B(J,I)
 54   B(J,I)=U
 50   CONTINUE
 55   RETURN
      END
      SUBROUTINE ROTATE(N,L)
C        SUBROUTINE ROTATE FOR BMD03M                  MAY  2, 1966
      COMMON  SUMXX
      COMMON  A      , R      , MM     , H      , TV
      DIMENSION A(60,60),TV(60),H(60)
      DIMENSION P(60,60),R(60,60),SUMXX(60,60),MM(60)
      EQUIVALENCE(SUMXX,P)
      EPS=0.00116
      P(1,1)=1.0
      NCOUNT=0
      MM(1)=(L*(L-1))/2
      TV(1)=0.0
      LL=L-1
      NV=1
      FN=N
      FFN=FN**2
      CONS=1.0/SQRT(2.0)
      ZERO=1.E-4
      DO 3 I=1,N
    3 H(I)=0.0
      DO 4 I=1,N
      DO 4 J=1,L
    4 H(I)=H(I)+A(I,J)*A(I,J)
      DO 5 I=1,N
      H(I)=SQRT(H(I))
      DO 5 J=1,L
    5 A(I,J)=A(I,J)/H(I)
      WRITE (6,10006)
10006 FORMAT ('1 ITERATION   CRITERION VALUE ',/)
  222 NV=NV+1
      TV(NV)=0.0
      LV=NV-1
      DO 88 J=1,L
      AA=0.0
      BB=0.0
      DO 77 I=1,N
      CC=A(I,J)*A(I,J)
      AA=AA+CC
   77 BB=BB+CC**2
   88 TV(NV)=TV(NV)+(FN*BB-AA**2)/FFN
      IVNQ=NV-2
      WRITE (6,10007) IVNQ,TV(NV)
10007 FORMAT (1X,I10,F20.10)
      IF(NV-50)9,999,999
 9    IF(ABS(TV(NV)-TV(LV))-1.E-6)999,999,13
   13 DO 500 J=1,LL
      II=J+1
      DO 500 K=II,L
      IF(NCOUNT-MM(1))32,500,500
   32 IF(ABS(P(1,1)-ZERO))14,14,10
   10 NCOUNT=0
   14 AA=0.0
      BB=0.0
      CC=0.0
      DD=0.0
      DO15I=1,N
      U=(A(I,J)+A(I,K))*(A(I,J)-A(I,K))
      T=A(I,J)*A(I,K)
      T=T+T
      CC=CC+(U+T)*(U-T)
      DD=DD+2.0*U*T
      AA=AA+U
   15 BB=BB+T
      T=DD-2.0*AA*BB/FN
      B=CC-(AA**2-BB**2)/FN
      P(1,1)=0.25*ATAN(T/B)
      IF (ABS(P(1,1))-ZERO) 7,7,69
    7 NCOUNT=NCOUNT+1
      GO TO 500
   69 TAN4P=T/B
      IF(T-B)1041,1433,1042
 1433 IF(T+B-EPS)500,1043,1043
 1043 COS4T=CONS
      SIN4T=CONS
      GOTO5000
 1041 TAN4T=ABS(T)/ABS(B)
      IF(TAN4T-EPS)8000,1100,1100
 1100 COS4T=1.0/SQRT(1.0+TAN4T**2)
      SIN4T=TAN4T*COS4T
      GOTO5000
 8000 IF(B)1150,500,500
 1150 SINP=CONS
      COSP=CONS
      GOTO1000
 1042 CTN4T=ABS(T/B)
      IF(CTN4T-EPS)9000,1200,1200
 1200 SIN4T=1.0/SQRT(1.0+CTN4T**2)
      COS4T=CTN4T*SIN4T
      GOTO5000
 9000 COS4T=0.0
      SIN4T=1.0
 5000 COS2T=SQRT((1.0+COS4T)/2.0)
      SIN2T=SIN4T/(2.0*COS2T)
      COST=SQRT((1.0+COS2T)/2.0)
      SINT=SIN2T/(2.0*COST)
      IF(B)1250,1250,1300
 1300 COSP=COST
      SINP=SINT
      GOTO7000
 1250 COSP=CONS*COST+CONS*SINT
      SINP=ABS(CONS*COST-CONS*SINT)
 7000 IF(T)1400,1400,1000
 1400 SINP=-SINP
 1000 DO 100 I=1,N
      AA=A(I,J)*COSP+A(I,K)*SINP
      BB=-A(I,J)*SINP+A(I,K)*COSP
      A(I,J)=AA
  100 A(I,K)=BB
  500 CONTINUE
      GO TO 222
  999 DO 6 I=1,N
      DO 6 J=1,L
    6 A(I,J)=A(I,J)*H(I)
      NC=NV-2
      WRITE (6,901)L
      WRITE (6,902)NC
      WRITE (6,903)
      DO 256 I3=1,L
      NMINUS = 0
      DO 252 I1=1,N
      IF(A(I1,I3)) 251,252,252
 251  NMINUS = NMINUS + 1
 252  CONTINUE
      IF(N - 2 * NMINUS) 253,256,256
 253  DO 254 I2=1,N
 254  A(I2,I3) = -A(I2,I3)
 256  CONTINUE
      DO 310 I=1,N
      WRITE (6,910)I
 310  WRITE (6,911)(A(I,J), J=1,L)
C
C
C
C
C
C
      WRITE (6,900)
      WRITE (6,907)
      WRITE (6,908)
      DO 330 I=1,N
      H(I)=H(I)*H(I)
      BB=0.0
      DO 322 J=1,L
 322  BB=BB+A(I,J)**2
      AA=BB-H(I)
 330  WRITE (6,909)I,H(I),BB,AA
 900  FORMAT(1H0)
 901  FORMAT(26H0NUMBER OF FACTORS ROTATEDI4)
 902  FORMAT(27H NUMBER OF ITERATION CYCLESI3//)
 903  FORMAT(22H0ROTATED FACTOR MATRIX)
C
C
C
 907  FORMAT(23H0CHECK ON COMMUNALITIES)
 908  FORMAT(12H0  VARIABLES8X,8HORIGINAL9X,5HFINAL8X,10HDIFFERENCE)
 909  FORMAT(I8,5X,3F15.5)
 910  FORMAT(9H0VARIABLEI3)
 911  FORMAT(10F12.5)
      RETURN
      END
      SUBROUTINE TPWD(LT1,LT2)
C        SUBROUTINE TPWD FOR BMD03M                    MAY  2, 1966
C     VERSION II ALLOWS ALTERNATE BINARY INPUT
      NT1=IABS(LT1)
      NT2=IABS(LT2)
      IF(NT1)12,10,12
 10   LT1=5
      NT1=5
 12   IF(NT1-NT2)14,19,14
   14 IF(NT2.EQ.5)GO TO 18
 17   REWIND NT2
   19 IF(NT1-5)18,24,18
 18   IF(NT1-6)22,40,22
 22   REWIND NT1
 24   LT2=LT1
 28   RETURN
 40   WRITE (6,49)
      STOP
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      END