C        CANONICAL ANALYSIS - MAIN PROGRAM         JANUARY 15, 1966
      DIMENSION DATE(2)
      DOUBLE PRECISION DATE,PR,FN,SET1,SET2,SET3,P,PC,A,ONO,FMO
      DATA DATE/'MAY 26, ','1969    '/
      DIMENSION C(50,50),X(50),KK(50),NN(3),LN(3,50),S(50),
     1D(50),LX(50),LY(50),IN(50)
      DIMENSION F(162),FMO(162)
      DIMENSION NTAPE(16)
      COMMON/BIG/NX,LX,LY,KK1,NC,ITX,D,NV,A2,FMO,X
      COMMON/NYPASS/NY
      DATA PR,FN,SET1,SET2,ONO/6HPROBLM,6HFINISH,6HSETONE,6HSETTWO,
     16HNO    /
      DATA SET3/6HSETTHR/
      LON=0
	CALL USAGEB('BMDX75')
      P10005 = 0
  100 N=0
    4 READ (5,1) P,PC,NV,NC,MMM,NF,IT,A1,A2,NFX,ITX
 11   IF(P.EQ.PR) GO TO 7
      IF(P.EQ.FN) GO TO 5
      WRITE (6,1010) P
 1010 FORMAT (' PROBLM OR FINISH CARD EXPECTED BUT READ ',A6)
      N=N+1
      GO TO 4
 5    IF(N.NE.0) WRITE (6,10) N
 10   FORMAT(1H0I5,17H RECORDS NOT READ)
      IF(LON.EQ.0) STOP
      DO 201 I=1,LON
      IF (NTAPE(I).EQ.6) GO TO 201
      II=NTAPE(I)
      ENDFILE II
      REWIND II
  201 CONTINUE
      STOP
C   7 CALL USEBUF
C     READ (5,1) P,PC,NV,NC,MMM,NF,IT,A1,A2,NFX,ITX
    7 IF (MMM.EQ.0) MMM=4
 1    FORMAT(2A6,I2,I6,2I1,I2,A2,A3,I1,I2)
      IF(ITX.EQ.0) GO TO 200
      LON=LON+1
      NTAPE(LON)=ITX
 200  IF(IT.EQ.0) IT=5
      IF (A1 .NE. ONO .AND. IT .NE. 5) REWIND IT
      REWIND 2
      REWIND 1
      IF(N.NE.0) WRITE (6,10) N
      WRITE (6,2) DATE,PC,NV,NC
      NN(1)=0
      NN(2)=0
      NN(3)=0
    2 FORMAT('1BMDX75 - CANONICAL ANALYSIS - REVISED ',2A8/
     X20H0PROBLEM CODE        ,2X,A6/
     120H0NUMBER OF VARIABLESI8,/20H0NUMBER OF CASES           I8)
   56 READ (5,3) A,(KK(I),I=1,19)
    3 FORMAT (A6,18A4,A2)
      REWIND 1
      WRITE (1,3) A,(KK(I),I=1,19)
      IF(A.EQ.SET1) GO TO 6
      IF(A.EQ.SET2) GO TO 31
C
C      SET3 WAS INTENDED FOR USE IN A FUTURE MODIFICATION OF THIS
C           PROGRAM TO ALLOW FOR THE ANALYSIS OF THREE SETS OF DATA
C
      IF(A.EQ.SET3) GO TO 33
C
      GO TO 53
 33   K=3
      GO TO 32
 6    K=1
      GO TO 32
 31   K=2
   32 REWIND 1
      READ (1,61) A,(KK(I),I=1,22)
 61   FORMAT(A6,22I3)
      DO 62 I=1,22
      IF(KK(I))63,62,65
 63   L0=KK(I-1)+1
      L1=-KK(I)
      IF((L0-1)*(NV-L0).LT.0 .OR. (L1-1)*(NV-L1).LT.0) GO TO 10000
      DO 66 L=L0,L1
      NN(K)=NN(K)+1
      N=NN(K)
 66   LN(K,N)=L
      GO TO 62
 65   NN(K)=NN(K)+1
      N=NN(K)
      LN(K,N)=KK(I)
      IF((KK(I)-1)*(NV-KK(I)))10000,62,62
 62   CONTINUE
      GO TO 56
10000 WRITE (6,10001) A
10001 FORMAT(' INDEX NUMBER EXCEEDS NUMBER OF VARIABLES FOR ',A6)
      GO TO 80
C
C
10006 WRITE(6,10007)
      GO TO 80
10008 WRITE(6,10009)
C
10005 FORMAT ('NUMBER OF VARIABLE FORMAT CARDS INCORRECT. ASSUMED TO BE
     11')
10007 FORMAT(' TOO MANY VARIABLES')
10009 FORMAT(' NUMBER OF VARIABLES CANNOT BE NEGATIVE')
 80   WRITE (6,81)
 81   FORMAT(45H0CONTROL CARDS INCORRECTLY ORDERED OR PUNCHED)
      GO TO 100
   53 IF (NF.GT.0.AND.NF.LE.9) GO TO 54
      P10005 = 1
      NF=1
   54 IF (ITX .EQ . 0) GO TO 55
      IF (NFX .GT. 0 .AND. NFX .LE. 9) GO TO 55
      P10005 = 1
      NFX=1
   55 NF=NF*18
      NFX=NFX*18
      IF(NN(1).EQ.0 .OR. NN(2).EQ.0)GO TO 10000
      IF(NV.GT.50)GO TO 10006
      IF(NV.LT.0)GO TO 10008
      REWIND 1
      READ (1,41) (F(I),I=1,NF)
      REWIND 1
 41   FORMAT(18A4)
      IF (P10005 .EQ. 0) GO TO 58
      WRITE(6,10005)
      P10005 = 0
   58 WRITE (6,1011) (F(I),I=1,NF)
      IF(ITX .EQ. 0) GO TO 59
 1011 FORMAT (' INPUT FORMAT'/(5X,18A4))
      READ (5,41) (FMO(I),I=1,NFX)
      WRITE (6,1012) (FMO(I),I=1,NFX)
 1012 FORMAT (' OUTPUT FORMAT'/(5X,18A4))
   59 NNT =21
      IF(NN(3).GT.0) NNT=22
      REWIND NNT
      DO 20 I=1,NV
      S(I)=0.
      DO 20 J=1,NV
 20   C(I,J)=0.
      DO 21 L=1,NC
      H=L
      H1=H*(H-1.)
      READ (IT,F) (X(I),I=1,NV)
      WRITE (NNT) (X(I),I=1,NV)
 23   DO 21 I=1,NV
      D(I)=(X(I)-S(I))/H
      S(I)=S(I)+D(I)
      PQ=D(I)*H1
      DO 21 J=1,I
 21   C(I,J)=C(I,J)+PQ*D(J)
      REWIND NNT
      H=NC-1
      DO 110 I=1,NV
      D(I)=1.
 110  X(I)=SQRT(C(I,I)/H)
      GO TO (101,102,101,102) ,MMM
 105  DO 120 I=1,NV
      DO 120 J=1,I
 120  C(I,J)=C(I,J)/H
      GO TO 108
 101  H=NC
      DO 104 I=1,NV
      PQ=H*S(I)
      DO 104 J=1,I
 104  C(I,J)=C(I,J)+PQ*S(J)
 102  GO TO (105,105,106,106),MMM
 106  DO 107 I=1,NV
      D(I)=C(I,I)
      IF(D(I).LE.0.) D(I)=1.
      D(I)=SQRT(D(I))
      DO 107 J=1,I
 107  C(I,J)=C(I,J)/(D(I)*D(J))
      H=SQRT(H)
      DO 109 I=1,NV
 109  D(I)=D(I)/H
 108  DO 24 I=1,NV
      DO 24 J=1,I
 24   C(J,I)=C(I,J)
      WRITE (6,85)
 85   FORMAT(1H-,4(32H  VARIABLE   MEAN      STANDARD )/1X,4(23X,9HDEVIA
     1TION))
      K=(NV-1)/4+1
      DO 87 I=1,K
 87   WRITE (6,86) (J,S(J),X(J),J=I,NV,K)
 86   FORMAT(1X,4(I7,F13.6,F12.6))
      GO TO (821,822,823,824),MMM
 821  WRITE (6,921)
      GO TO 960
 822   WRITE (6,922)
      GO TO 960
 823  WRITE (6,923)
      GO TO 960
 824  WRITE (6,924)
 921  FORMAT(76H1THE COVARIANCE MATRIX ABOUT THE ORIGIN IS USED IN THE F
     XOLLOWING CALCULATION)
 922  FORMAT(59H1THE COVARIANCE MATRIX IS USED IN THE FOLLOWING CALCULAT
     XION)
 923  FORMAT(77H1THE CORRELATION MATRIX ABOUT THE ORIGIN IS USED IN THE 
     XFOLLOWING CALCULATION)
 924  FORMAT(60H1THE CORRELATION MATRIX IS USED IN THE FOLLOWING CALCULA
     XTION)
 960  K=1
      IF(NN(1).LE.NN(2)) K=2
      KK1=K
      L=3-K
      NX=NN(K)
      DO 70 I=1,NX
 70   LX(I)=LN(K,I)
      NY=NN(L)
      J=NX
      DO 71 I=1,NY
      LY(I)=LN(L,I)
      J=J+1
 71   LX(J)=LY(I)
      NZ=NN(3)
      IF(NZ.EQ.0) GO TO 35
      WRITE (6,38) (LN(3,I),I=1,NZ)
 38   FORMAT(22H0SET THREE VARIABLES -/(5X,40I3))
      DO 34 I=1,NV
 34   IN(I)=1
      DO 36 I=1,NZ
      J=LN(3,I)
 36   IN(J)=0
      CALL MATINV(C,NV,IN,1.E-6)
      DO 305 I=1,NV
      DO 305 J=1,I
 305  C(J,I)=C(I,J)
      DO 303 L=1,NC
      READ(22) (X(I),I=1,NV)
      DO 302 I=1,NZ
      J=LN(3,I)
      XJ=X(J)
      DO 302 K=1,NV
 302  X(K)=X(K)-C(J,K)*XJ
      DO 301 I=1,NZ
      J=LN(3,I)
 301  X(J)=0.
  303 WRITE (21) (X(I),I=1,NV)
 1000 FORMAT (20A4)
      REWIND 21
      DO 37 I=1,NZ
      J=LN(3,I)
      DO 37 K=1,NV
      C(K,J)=0.
 37   C(J,K)=0.
 35   NV1=NX+NY
      DO 72 I=1,NV
 72   IN(I)=I
      DO 73 I=1,NV1
      J=LX(I)
 76   IF(LX(I).EQ.IN(J)) GO TO 75
      J=IN(J)
      GO TO 76
 75   IF(I.EQ.J) GO TO 73
      M=IN(J)
      IN(J)=IN(I)
      IN(I)=M
      DO 78 K=1,NV
      T=C(I,K)
      C(I,K)=C(J,K)
 78   C(J,K)=T
      DO 77 K=1,NV
      T=C(K,I)
      C(K,I)=C(K,J)
 77   C(K,J)=T
 73   CONTINUE
      NX1=NX+1
      CALL CANON(C,C(NX1,1),C(NX1,NX1),C(1,NX1))
C     CALL CANON(C,C(NX1,1),C(NX1,NX1),C(1,NX1),NX,NY,LX,LY,KK1,NC,ITX,D
C    1,NX,A2,FMO)
      GO TO 100
      END
C        SUBROUTINE CANON FOR BMDX75               JANUARY 15, 1966
      SUBROUTINE CANON(A,B,C,E)
      DOUBLE PRECISION FMO
      DATA YES/3HYES/
      DIMENSION FMO(162)
      DIMENSION D(50),X(50)
      DIMENSION A(50,50),B(50,50),C(50,50),E(50,50),R(50),IN(5
     10),DX(50),DY(50),LX(50),LY(50)
      COMMON/BIG/NX,LX,LY,KK1,NC,ITS,D,MV,A2,FMO,X
      COMMON/NYPASS/NY
      COMMON/NRPASS/NR
      T1=1.E-4
      T2=T1
      NV=NX+NY
      DO 2 J=1,NX
 2    IN(J)=0
      DO 1 I=1,NY
      J=I+NX
      IN(J)=1
      DO 1 J=1,I
      E(I,J)=C(I,J)
 1    C(I,J)=0.
      CALL MATINV(A,NV,IN,T1)
      DO 3 I=1,NY
      DO 3 J=1,I
      C(I,J)=-C(I,J)
      C(J,I)=C(I,J)
 3    E(J,I)=E(I,J)
      CALL JACOB2 (C,E,1.E-7,1,T2)
      DO 4 I=1,NY
 4    R(I)=SIGN(SQRT(ABS(C(I,I))),C(I,I))
      DO 5 I=1,NX
      DO 5 J=1,NY
      T=0.
      DO 8 K=1,NY
      IF(R(J).LE.0.) E(K,J)=0.
 8    T=T+E(K,J)*B(K,I)
      IF(R(J).LE.0.) A(I,J)=0.
 5    IF(R(J).GT.0.) A(I,J)=T/R(J)
      L1=0
 30   L0=L1+1
      L1=   MIN0(L1+10,NY)
      LL=0
      WRITE (6,31) (I,I=L0,L1)
 31   FORMAT(1H-15X,22HCANONICAL CORRELATIONS//(2X,10I12))
      WRITE (6,32) (R(I),I=L0,L1)
 32   FORMAT(6X,10F12.5)
      WRITE (6,33)
      GO TO (35,34),KK1
 33   FORMAT(66H0VARIABLE    COEFFICIENTS FOR CANONICAL VARIABLES OF THE
     1 FIRST SET//)
 35   DO 48 I=1,NX
 48   WRITE (6,39) LX(I),(A(I,L),L=L0,L1)
 39   FORMAT(I6,10F12.5)
 49   IF(LL)40,40,41
 40   LL=1
      WRITE (6,43)
      GO TO (34,35) ,KK1
 43   FORMAT(67H0VARIABLE    COEFFICIENTS FOR CANONICAL VARIABLES OF THE
     1 SECOND SET//)
 34   DO 38 I=1,NY
 38   WRITE (6,39) LY(I),(E(I,L),L=L0,L1)
      GO TO 49
 41   IF(L1-NY)30,47,30
 47   IF(A2.NE.YES) RETURN
      DO 304 I=1,NY
      DO 302 J=1,NX
      K=LX(J)
 302  A(J,I)=A(J,I)/D(K)
      DO 304 J=1,NY
      K=LY(J)
 304  E(J,I)=E(J,I)/D(K)
      WRITE (6,180) (I,I=1,NY)
 180  FORMAT(1H1,20X,19HCANONICAL VARIABLES//6(10X,I7,4X))
      WRITE (6,181)
 181  FORMAT(6H  CASE/6H   NO.)
      REWIND 21
      DO 100 L=1,NC
      LL=0
      DO 101 IX=1,50
  101 X(IX)=0.0
8765  FORMAT(7F10.5)
      READ (21) (X(I),I=1,MV)
 1000 FORMAT (20A4)
      DO 102 I=1,NY
      DX(I)=0.
      DO 102 J=1,NX
      K=LX(J)
 102  DX(I)=DX(I)+A(J,I)*X(K)
      DO 103 I=1,NY
      DY(I)=0.
      DO 103 J=1,NY
      K=LY(J)
 103  DY(I)=DY(I)+E(J,I)*X(K)
      GO TO (160,161),KK1
 160  IF(ITS.NE.0) WRITE (ITS,FMO) (X(I),I=1,NV),(DX(I),DY(I),I=1,NY)
      WRITE (6,170) L,(DX(I),DY(I),I=1,NY)
 170  FORMAT(1H0,I4,6(2F10.5,1X)/(5X,2F10.5,1X,2F10.5,1X,2F10.5,1X,2F10.
     X5,1X,2F10.5,1X,2F10.5,1X))
      GO TO 100
 161  IF(ITS.NE.0) WRITE (ITS,FMO) (X(I),I=1,NV),(DY(I),DX(I),I=1,NY)
      WRITE (6,170) L,(DY(I),DX(I),I=1,NY)
 100  CONTINUE
      REWIND 21
      RETURN
      END
C        SUBROUTINE JACOB2 FOR BMDX75                JANUARY 15, 1966
      SUBROUTINE JACOB2 (A,B,ACC,IV,TOL)
      DIMENSION A(50,2),B(50,2),U(50),IN(50),V(50)
      COMMON/NYPASS/N
      COMMON/NRPASS/NR
      Q=1.
      DO 1 I=1,N
      V(I)=B(I,I)
 1    IN(I)=0
      K=1
      DO 13 L=1,N
      DO 2 I=1,N
      U(I)=B(I,K)
 2    B(K,I)=0.
      P=U(K)
C
      U(K)=1.
      S=Q-TOL
      IF(S.GT.0.) RP=SQRT(P)
      Q=-1.
      IN(K)=-1
      DO 11 I=1,N
      Y=U(I)/P
      IF(IN(I))6,3,11
 3    IF(S)16,16,17
 17   DO 4 J=1,N
 4    B(J,I)=B(J,I)-U(J)*Y
 16   IF(B(I,I)/V(I)-Q)9,9,5
 5    Q=B(I,I)/V(I)
      KK=I
 9    IF(S)11,11,18
 18   T=A(I,I)-Y*A(I,K)
      DO 10 J=1,N
      A(I,J)=A(I,J)-Y*A(K,J)
 10   A(J,I)=A(I,J)
      A(I,I)=T-Y*A(I,K)
      GO TO 11
 6    IF(S)11,11,23
 23   DO 8 J=1,N
      IF(IN(J))7,14,7
 14   B(J,I)=0.
      GO TO 8
 7    B(J,I)=U(J)/RP
 8    CONTINUE
 11   CONTINUE
      IF(S)19,19,20
 19   DO 21 I=1,N
      B(I,K)=0.
      A(I,K)=0.
 21   A(K,I)=0.
      GO TO 22
 20   DO 12 I=1,N
      A(I,K)=A(I,K)/RP
 12   A(K,I)=A(I,K)
      A(K,K)=A(K,K)/RP
 22   IN(K)=1
 13   K=KK
      CALL JACOBI(A,B,ACC,2*IV,NR)
      RETURN
      END
C        SUBROUTINE JACOBI FOR BMDX75              JANUARY 15, 1966
      SUBROUTINE JACOBI (A,B,ACC,IV)
      DIMENSION A(50,50),B(50,50),LK(50),Q(50)
      COMMON/NYPASS/N
      COMMON/NRPASS/NR
      IF(IV-1)200,201,200
 201  DO 202 I=1,N
      DO 203 J=1,N
 203  B(I,J)=0.
 202  B(I,I)=1.
 200  NR=0
      Q(1)=0.0
      W=0.
      H=0.5*A(1,1)*A(1,1)
      DO 1 I=2,N
      H=H+.5*A(I,I)*A(I,I)
      Q(I)=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.*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.*Y*Y))/Y
      W=SQRT(1.+T*T)
      S=T/W
      C=1./W
      CC=C*C
      SS=S*S
      SC=S*C*2.
      Q1=0.
      Q2=0.
      W=0.
      NR=NR+1
      LK(1)=0
      DO 27 I=1,N
      LKI=LK(I)
      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(LKI-II)100,101,100
 101  Q(I)=-1.E20
      I3=I-1
      DO 103 J=1,I3
      IF(ABS(A(I,J))-Q(I))103,103,104
 104  LK(I)=J
      Q(I)=ABS(A(I,J))
 103  CONTINUE
 100  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.
      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
      IF((LK(I)-II)*(LK(I)-JJ)) 105,101,105
 105  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   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
      IF(IV)50,50,54
 54   DO 55 J=1,N
      T=B(J,I)
      B(J,I)=B(J,K)
 55   B(J,K)=T
 50   CONTINUE
      RETURN
      END
C        SUBROUTINE MATINV FOR BMDX75              JANUARY 15, 1966
      SUBROUTINE MATINV(A,N,IN,T)
      DIMENSION A(50,50),IN(50),U(50),V(50)
      H=1.
      DO 1 I=1,N
      IF(IN(I))1,10,1
 10   K=I
 1    V(I)=A(I,I)
 8    DO 2 I=1,K
      U(I)=A(K,I)
 2    A(K,I)=0.
      P=U(K)
      DO 3 I=K,N
      U(I)=A(I,K)
 3    A(I,K)=0.
      S=H-T
      H=-1.
      IN(K)=1
      U(K)=-1.
      DO 7 I=1,N
      IF(S)11,11,12
 12   Y=U(I)/P
      DO 4 J=1,I
 4    A(I,J)=A(I,J)-Y*U(J)
 11   IF(IN(I))5,5,7
 5    IF(H*V(I).GE.A(I,I))GO TO 7
 6    H=A(I,I)/V(I)
      K=I
 7    CONTINUE
      IF(H+1.) 8,9,8
 9    RETURN
      END