C             CANONICAL ANALYSIS                      JUNE     30, 1966
C        THIS IS A SIFTED VERSION OF BMD06M ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
      DIMENSION RHO(70,70),DATA(100),SX(70),SX2(70),D(35,35),
     1A(35,35),E(35,35),G(35,35),B(35),C(35),QQQ(70,35)
      DIMENSION L1(35),K1(35),FMT(180)
      DIMENSION NNEWA(99),LLCODE(99),LLVA(99),BBNEW(99)
      DIMENSION V(70)
      COMMON  QQQ
      COMMON  N      , A      , B      , C      , RHO    , D
      EQUIVALENCE (G,QQQ(1226)),(E,QQQ)
      EQUIVALENCE(DATA,V)
C
      DATA Q003HL/'YES '/
C
  405 FORMAT(42H1BMD06M - CANONICAL ANALYSIS - VERSION OF ,
     116HDECEMBER 8, 1964/
     2/40H HEALTH SCIENCES COMPUTING FACILITY,UCLA///)
C
C     THIS PROGRAM USES TWO SCRATCH TAPES.  THEY ARE DESIGNATED HERE
C     BY IT2 AND IT3.
C
      INFORM=5
	CALL USAGEB('BMD06M')
      IT2=2
      IT3=3
      REWIND IT2
      REWIND IT3
      DOUBLE PRECISION A123,B123,C123,XODE,KODE,TODE,PROB
      DATA A123/'FINISH'/,B123/'PROBLM'/,C123/'TRNGEN'/
  999 READ (5,5001)TODE,PROB,NCASES,NVAR,IP,IQ,CRITER, PRINTD,NVG,NTAPE,
     1KVR
 5001 FORMAT(2A6,I4,I3,2I2,F5.0,A3,35X,3I2)
      IF(TODE.EQ.A123) GO TO 401
      IF(TODE.EQ.B123) GO TO 403
  402 WRITE (6,404) TODE
  404 FORMAT (51H PROGRAM EXPECTED FINISH OR PROBLM CARD. IT READ A   ,
     1 A6,5H CARD  )
  401 IF(INFORM-5)27,27,26
   26 REWIND INFORM
   27 STOP
  403 IF((NTAPE-IT2)*(NTAPE-IT3))473,440,473
  473 CALL TPWD(NTAPE,INFORM)
      WRITE (6,405)
      WRITE (6,9002)PROB
 9002 FORMAT(10X,23HTHIS IS PROBLEM NUMBER A6//)
      WRITE (6,9003)NCASES,IP,IQ,CRITER,NVAR,NVG,KVR
 9003 FORMAT(10X,18HNUMBER OF CASES = I4//10X,35HNUMBER OF VARIABLES IN   
     1FIRST SET = I2//
     210X,36HNUMBER OF VARIABLES IN SECOND SET = I2//
     310X,75HONLY THE CANONICAL COEFFICIENTS WITH CORRESPONDING CORRELAT
     4ION LARGER THAN F6.5,13H ARE COMPUTED//
     510X,30HNUMBER OF VARIABLES READ IN = I3//
     610X,35HNUMBER OF TRANS-GENERATION CARDS = I2//
     710X, 39HNUMBER OF VARIABLE FORMAT CARDS USED = I2//)
      IF(INFORM-5)901,902,901
  902 WRITE (6,9005)
 9005 FORMAT(10X,22HINPUT DATA IS BY CARDS//)
      GO TO 531
C
  440 WRITE (6,441) NTAPE
  441 FORMAT (19H ILLEGAL TAPE UNIT  ,I3,22H LISTED ON PROBLM CARD  )
      GO TO 401
C
  901 WRITE (6,9006)INFORM
 9006 FORMAT(10X,47HINPUT DATA IS BY TAPE---LOGICAL TAPE NUMBER IS I2//)
  531 IPQ=IP+IQ
      IF(IQ-35) 532,532,910
 532  IF(IP-IQ)533,533,920
 533  IA=IP+1
      N=IP
      IF(NVAR*(NVAR-101)) 534,535,535
535   WRITE (6,536)NVAR
536   FORMAT(1H0,20X,65HNUMBER OF ENTERING VARIABLES (T) EXCEEDS PROBLEM
     1 LIMITS. YOU HAVE,I5,1H.)
      GO TO 401
  534 CASES=NCASES
      MERRY=0
      LCASE=0
      IF(KVR.GT.0.AND.KVR.LE.10) GO TO 540
      WRITE(6,6000)
 6000 FORMAT(1H0,23X,71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPEC
     XIFIED, ASSUMED TO BE 1.)
      KVR = 1
  540  CONTINUE
      KVR=KVR*18
      READ (5,5003)(FMT(I),I=1,KVR)
 5003 FORMAT(18A4)
      WRITE (6,5004) (FMT(I),I=1,KVR)
 5004 FORMAT (1X,18HVARIABLE FORMAT IS  /,(1H ,18A4))
C
   25 DO 4 I=1,IPQ
      SX(I)=0.0
      SX2(I)=0.0
      DO 4 J=1,IPQ
    4 RHO(I,J)=0.0
      IF(NVG) 29,29,28
   28 DO 24 I=1,NVG
      READ (5,5002)XODE,NNEWA(I),LLCODE(I),LLVA(I),BBNEW(I)
 5002 FORMAT(A6,I3,I2,I3,F6.0)
      IF(XODE.EQ.C123) GO TO 23
   22 NVG=-101
      GO TO 24
   23 IF(LLCODE(I)-41)232,24,22
  232 IF(LLCODE(I)*(LLCODE(I)-18))24,22,22
   24 CONTINUE
      IF(NVG)234,29,29
  234 J=I
      WRITE (6,5073) XODE,NNEWA(J),LLCODE(J),LLVA(J),BBNEW(J)
 5073 FORMAT (1X,A6,I4,I3,I4,F7.0)
      WRITE (6,4003)
 4003 FORMAT (49H ABOVE TRNGEN CARD ERROR. PROGRAM CANNOT PROCEED  )
      GO TO 401
   29 IF((+Q003HL).NE.PRINTD) GO TO 21
295   WRITE (6,1035)PROB
 1035 FORMAT (1H1,9X,22HINPUT DATA OF PROBLEM A6///2X,8HCASE NO./)
 21   DO 150 II=1,NCASES
      H=II
      READ (INFORM,FMT)(DATA(I),I=1,NVAR)
      IF((+Q003HL).NE.PRINTD) GO TO 33
 325  WRITE (6,1040)II,(DATA(I),I=1,NVAR)
 1040 FORMAT(1H0,I6,2X,10F11.6/(9X,10F11.6))
 33   IF(NVG) 702,702,34
   34 CALL TRNGEN (DATA,IPQ,NVG,NCASES,LCASE,NNEWA,LLCODE,LLVA,BBNEW,MER
     1RY,II)
      IF(LCASE) 35,36,36
   35 LCASE=0
      GO TO 150
C
   36 WRITE (IT2) (DATA(I),I=1,IPQ)
  702  DO 8 I=1,IPQ
      SX(I)=SX(I)+DATA(I)
      SX2(I)=DATA(I)-SX(I)/H
      IF(H-1.0) 718,728,718
  728 TT=0.
      GO TO 738
  718 TT=SX2(I)*H/(H-1.0)
      IF(H.EQ.1.) TT=0.
  738 DO 8 J=1,I
 8    RHO(I,J)=RHO(I,J)+TT*SX2(J)
  150 CONTINUE
 5759 DO  7589 I=1,IPQ
      DO  7589 J=1,I
 7589 RHO(J,I)=RHO(I,J)
      IF(NVG)706,706,707
 707  REWIND IT2
      WRITE (6,1045)IP, PROB,IQ
 1045 FORMAT(1H1,10X,38HDATA AFTER TRANSGENERATION, THE FIRST I2,
     1 42H BELONG TO THE FIRST SET IN THIS PROBLEM (A6,1H),/
     2 10X,10HTHE OTHER I2, 25H BELONG TO THE SECOND SET//)
      DO 40 II=1,NCASES
      READ (IT2) (DATA(I),I=1,IPQ)
   40 WRITE (6,1040)II,(DATA(I),I=1,IPQ)
      REWIND IT2
 706  WRITE (6,1100)PROB
 1100 FORMAT(1H1,10X,32HSTANDARD DEVIATIONS FOR PROBLEM ,A6/1H0,10X,
     1 8HVARIABLE,3X,8HST. DEV./1H0)
       DO 11 I=1,IPQ
      DATA(I)=SQRT(RHO(I,I)/(H-1.0))
 11   WRITE (6,1101)I,DATA(I)
 1101 FORMAT(10X,I6,F15.6)
      DO 12  I=1,IPQ
      DO 12 J=1,IPQ
 12   RHO(I,J)=RHO(I,J)/((H-1.0)*DATA(I)*DATA(J))
      WRITE (6,1018)PROB
 1018 FORMAT(1H1,10X,30HCORRELATION MATRIX OF PROBLEM A6///)
      DO 140 I=1,IPQ
      WRITE (6,1017)I,(RHO(I,J),J=1,IPQ)
 1017 FORMAT(4H0ROWI3,3X,15F8.4/(10X,15F8.4))
  140 WRITE (IT3) (RHO(I,J),J=1,IPQ)
      REWIND IT3
      DO 205 I=1,IP
      DO 205 J=1,IP
  205 D(I,J)=RHO(I,J)
      DO 210 II=IA,IPQ
      I=II-IP
      DO 210 JJ=IA,IPQ
      J=JJ-IP
  210 A(I,J)=RHO(II,JJ)
      CALL INV(A,A,IQ,1.E-7,0)
      DO 215 I=1,IP
      DO 215 JJ=IA,IPQ
      J=JJ-IP
  215 E(I,J)=RHO(I,JJ)
      DO 230 I=1,IQ
      DO 230 J=1,IP
      G(I,J)=0.0
      DO 230 K=1,IQ
  230 G(I,J)=G(I,J)+A(I,K)*E(J,K)
      DO 235 I=1,IP
      DO 235 J=1,IP
      A(I,J)=0.0
      DO 235 K=1,IQ
  235 A(I,J)=A(I,J)+E(I,K)*G(K,J)
      DO 240 I=1,IP
      DO 240 J=1,IP
  240 RHO(I,J)=A(I,J)
      CALL INV(D,A,IP,1.E-7,1)
      CALL HDIAG(A,IP,1,G,NNRR)
      DO 225 I=1,IP
 225  C(I)=A(I,I)
      WRITE (6,1021)PROB
 1021 FORMAT(1H1,10X,34HCANONICAL CORRELATIONS OF PROBLEM A6///)
      IPP=IP-1
      DO 250 I=1,IPP
      K=I+1
      DO 250 J=K,IP
      IF(C(I)-C(J))255,250,250
  255 TEMP=C(I)
      C(I)=C(J)
      C(J)=TEMP
  250 CONTINUE
      DO 260 I=1,IP
      IF(C(I))251,252,252
 251  WRITE (6,4002)I,C(I)
      SX(I)=CRITER-1.0
      GO TO 260
 252  SX(I)=SQRT(C(I))
  260 WRITE (6,1022)I,SX(I)
 1022 FORMAT(10X,I2,10X,F10.7)
      DO 265 I=1,IPQ
  265 READ (IT3) (RHO(I,J),J=1,IPQ)
      REWIND IT3
      WRITE (6,1023)CRITER
 1023 FORMAT(1H1,10X,77HCANONICAL COEFFICIENTS---ONLY THOSE WITH CORRESP
     1ONDING CANONICAL CORRELATIONS/
     2  38X, 13H LARGER THAN F6.4, 13H ARE COMPUTED/10X,80(1H.))
      WRITE (6,2023)
 2023 FORMAT(10X,80(1H.))
      DO 300 I=1,IP
      DO 300 J=1,IP
  300 A(I,J)=RHO(I,J)
      N=IPQ
      DO 777 L=1,IP
      IF(SX(L)-CRITER) 7777,7777,305
  305 WRITE (6,1024)SX(L)
 1024 FORMAT(1H0,5X,24HCANONICAL CORRELATION = F10.8)
      WRITE (6,2024)
 2024 FORMAT(5X,34(1H*)//)
      DO 307 I=1,IP
      DO 307 J=IA,IPQ
      RHO(I,J)=-RHO(I,J)/SX(L)
 307  RHO(J,I)=RHO(I,J)
      CALL VCTR (RHO,V,N,0.0)
      ZZZ=0.0
      DO 501 I=1,IP
      DO 501 J=1,IP
 501  ZZZ=ZZZ+V(I)*A(I,J)*V(J)
      ZZZ=1.0/SQRT(ZZZ)
      DO 502 I=1,IPQ
      V(I)=V(I)*ZZZ
 502  QQQ(I,L)=V(I)
      WRITE (6,1025)
 1025 FORMAT(10X,43HCOEFFICIENTS FOR THE FIRST SET OF VARIABLES//)
      WRITE (6,1026)(V(I),I=1,IP)
 1026 FORMAT(8(5X,F10.6))
      WRITE (6,1027)
 1027 FORMAT(1H0/1H0,9X,45HCOEFFICIENTS FOR THE SECOND SET OF VARIABLES
     1//)
      WRITE (6,1026)(V(I),I=IA,IPQ)
      WRITE (6,2000)
 2000 FORMAT(3(1H0/))
      DO 776 I=1,IPQ
  776 READ (IT3) (RHO(I,J),J=1,IPQ)
 8000 FORMAT (20A4)
      REWIND IT3
  777 CONTINUE
      L=IP+1
 7777 L2=L-1
      CALL ABC(RHO,QQQ,QQQ,D,IP,IP,L2,RHO)
      CALL ABC(RHO(IP+1,1),QQQ,QQQ(IP+1,1),D,IP,IQ,L2,RHO(L,1))
      CALL ABC(RHO(IP+1,IP+1),QQQ(IP+1,1),QQQ(IP+1,1),D,IQ,IQ,L2,RHO(L,L
     1))
      L1(1)=2*L2
      LL1 = L1(1)
      DO 7778 I=L,LL1
      DO 7778 J=1,L2
 7778 RHO(J,I)=RHO(I,J)
      WRITE(6,7779)
 7779 FORMAT(78H1CANONICAL CHECK MATRIX - ONLY 0., 1. AND CANONICAL CORR
     1ELATIONS SHOULD APPEAR)
      DO 7780 I=1,LL1
 7780 WRITE (6,1019) I,(RHO(I,J),J=1,LL1)
 1019 FORMAT (4H0ROW,I3,3X,15F7.3,/(10X,15F7.3))
      GO TO 999
 910  WRITE (6,4000)IQ
  915 WRITE (6,9004)
      IF(KVR.GT.0.AND.KVR.LE.10) GO TO 950
      WRITE(6,6000)
      KVR = 1
  950  CONTINUE
      KVR=KVR*18
      READ (5,5003)(FMT(I),I=1,KVR)
      WRITE (6,5004) (FMT(I),I=1,KVR)
      IF(-NVG)916,918,918
  916 DO 917 I=1,NVG
  917 READ (5,5002)KODE
  918 DO 919 I=1,NCASES
  919 READ (INFORM,FMT)(DATA(J),J=1,NVAR)
      GO TO 999
 4000 FORMAT(1H0,34X,50HMAXIMUM NUMBER OF VARIABLES (Q) EXCEEDED. YOU HA
     1VE,I5,1H.)
 920  WRITE (6,4001)IP,IQ
      GO TO 915
 4001 FORMAT(1H0,9X,41HTHE NUMBER OF VARIABLES IN THE FIRST SET,,I4,47H,
     1 IS GREATER THAN THE NUMBER IN THE SECOND SET,,I4,1H.)
 4002 FORMAT(1H028X28HCANONICAL CORRELATION NUMBERI3,20H SQUARED IS EQUA
     1L TO F10.8/1H019X80HTHE PROGRAM WILL CONTINUE AND DETERMINE COEFFI
     2CIENTS FOR THE VALID CORRELATIONS.)
 9004 FORMAT(1H ,31X,54HPROGRAM WILL GO TO NEXT PROBLEM (IF ANY) OR TERM
     1INATE.)
      END
C             SUBROUTINE    ABC FOR BMD06M            JUNE     30, 1966
      SUBROUTINE ABC(R,A,C,D,IP,IQ,L,S)
      DIMENSION R(70,70),A(70,35),C(70,35),D(35,35),S(70 ,70)
      DO 1 I=1,IQ
      DO 1 J=1,L
      D(I,J)=0.0
      DO 1 K=1,IP
 1    D(I,J)=D(I,J)+R(I,K)*A(K,J)
      DO 2 I=1,L
      DO 2 J=1,L
      S(I,J)=0.0
      DO 2 K=1,IQ
 2    S(I,J)=S(I,J)+C(K,J)*D(K,I)
      RETURN
      END
C             SUBROUTINE  HDIAG FOR BMD06M            JUNE     30, 1966
C     MIHDI3, FORTRAN II DIAGONALIZATION OF A REAL SYMMETRIC MATRIX BY
C        THE JACOBI METHOD.
C             SIFTED TO FORTRAN IV
C     CALLING SEQUENCE FOR DIAGONALIZATION
C           CALL  HDIAG( H, N, IEGEN, U, NR)
C            WHERE H IS THE ARRAY TO BE DIAGONALIZED.
C     N IS THE ORDER OF THE MATRIX, H.
C
C     IEGEN MUST BE SET UNEQUAL TO ZERO IF ONLY EIGENVALUES ARE
C            TO BE COMPUTED.
C     IEGEN MUST BE SET EQUAL TO ZERO IF EIGENVALUES AND EIGENVECTORS
C            ARE TO BE COMPUTED.
C
C     U IS THE UNITARY MATRIX USED FOR FORMATION OF THE EIGENVECTORS.
C
C     NR IS THE NUMBER OF ROTATIONS.
C
C     A DIMENSION STATEMENT MUST BE INSERTED IN THE SUBROUTINE.
C     DIMENSION H(N,N), U(N,N), X(N), IQ(N)
C
C     SUBROUTINE PLACES COMPUTER IN THE FLOATING TRAP MODE
C
C     THE SUBROUTINE OPERATES ONLY ON THE ELEMENTS OF H THAT ARE TO THE
C            RIGHT OF THE MAIN DIAGONAL.  THUS, ONLY A TRIANGULAR
C            SECTION NEED BE STORED IN THE ARRAY H.
C
C
      SUBROUTINE   HDIAG (H,N,IEGEN,U,NR)
      DIMENSION H(35,35),U(35,35),X(35),IQ(35)
C     EXTERNAL SIGN
C
      IF (IEGEN) 15,10,15
 10   DO 14 I=1,N
      DO 14 J=1,N
      IF(I-J)12,11,12
 11   U(I,J)=1.0
      GO TO 14
 12   U(I,J)=0.0
 14   CONTINUE
C
 15   NR = 0
      IF (N-1) 1000,1000,17
C
C     SCAN FOR LARGEST OFF DIAGONAL ELEMENT IN EACH ROW
C     X(I) CONTAINS LARGEST ELEMENT IN ITH ROW
C     IQ(I) HOLDS SECOND SUBSCRIPT DEFINING POSITION OF ELEMENT
C
 17   NMI1=N-1
      DO 30 I=1,NMI1
      X(I) = 0.0
      IPL1=I+1
      DO 30 J=IPL1,N
      IF ( X(I) - ABS( H(I,J))) 20,20,30
 20   X(I)=ABS(H(I,J))
      IQ(I)=J
 30   CONTINUE
C
C     SET INDICATOR FOR SHUT-OFF.RAP=2**-27,NR=NO. OF ROTATIONS
      RAP=.745058059E-08
      HDTEST=1.0E38
C
C     FIND MAXIMUM OF X(I) S FOR PIVOT ELEMENT AND
C     TEST FOR END OF PROBLEM
C
 40   DO  70  I=1,NMI1
      IF (I-1) 60,60,45
 45   IF ( XMAX- X(I)) 60,70,70
 60   XMAX=X(I)
      IPIV=I
      JPIV=IQ(I)
 70   CONTINUE
C
C     IS MAX. X(I) EQUAL TO ZERO, IF LESS THAN HDTEST, REVISE HDTEST
      IF ( XMAX) 1000,1000,80
 80   IF (HDTEST) 90,90,85
 85   IF (XMAX - HDTEST) 90,90,148
 90   HDIMIN = ABS( H(1,1) )
      DO 110  I= 2,N
      IF (HDIMIN- ABS( H(I,I))) 110,110,100
 100  HDIMIN=ABS(H(I,I))
 110  CONTINUE
C
      HDTEST=HDIMIN*RAP
C
C     RETURN IF MAX.H(I,J)LESS THAN(2**-27)ABSF(H(K,K)-MIN)
      IF (HDTEST- XMAX) 148,1000,1000
 148  NR = NR+1
C
C     COMPUTE TANGENT, SINE AND COSINE,H(I,I),H(J,J)
 150  TANG=SIGN(2.0,(H(IPIV,IPIV)-H(JPIV,JPIV)))*H(IPIV,JPIV)/(ABS(H(IPI
     1V,IPIV)-H(JPIV,JPIV))+SQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.0*H(IP
     2IV,JPIV)**2))
      COSINE=1.0/SQRT(1.0+TANG**2)
      SINE=TANG*COSINE
      HII=H(IPIV,IPIV)
      H(IPIV,IPIV)=COSINE**2*(HII+TANG*(2.0*H(IPIV,JPIV)+TANG*H(JPIV,JPI
     1V)))
      H(JPIV,JPIV)=COSINE**2*(H(JPIV,JPIV)-TANG*(2.0*H(IPIV,JPIV)-TANG*H
     1 II))
      H(IPIV,JPIV)=0.0
C
C      PSEUDO RANK THE EIGENVALUES
C      ADJUST SINE AND COS FOR COMPUTATION OF H(IK) AND U(IK)
      IF ( H(IPIV,IPIV) -  H(JPIV,JPIV)) 152,153,153
 152  HTEMP = H(IPIV,IPIV)
      H(IPIV,IPIV) = H(JPIV,JPIV)
      H(JPIV,JPIV) = HTEMP
C       RECOMPUTE SINE AND COS
      HTEMP = SIGN(1.0, -SINE) * COSINE
      COSINE = ABS(SINE)
      SINE = HTEMP
 153  CONTINUE
C
C     INSPECT THE IQS BETWEEN I+1 AND N-1 TO DETERMINE
C     WHETHER A NEW MAXIMUM VALUE SHOULD BE COMPUTED SINCE
C     THE PRESENT MAXIMUM IS IN THE I OR J ROW.
C
      DO 350 I=1,NMI1
      IF(I-IPIV)210,350,200
 200  IF(I-JPIV)210,350,210
 210  IF(IQ(I)-IPIV)230,240,230
 230  IF(IQ(I)-JPIV)350,240,350
 240  K=IQ(I)
 250  HTEMP=H(I,K)
      H(I,K)=0.0
      IPL1=I+1
      X(I) =0.0
C
C     SEARCH IN DEPLETED ROW FOR NEW MAXIMUM
C
      DO 320 J=IPL1,N
      IF ( X(I)- ABS( H(I,J)) ) 300,300,320
 300  X(I) = ABS(H(I,J))
      IQ(I)=J
 320  CONTINUE
      H(I,K)=HTEMP
 350  CONTINUE
C
      X(IPIV) =0.0
      X(JPIV) =0.0
C
C     CHANGE THE OTHER ELEMENTS OF H
C
      DO 530 I=1,N
C
      IF(I-IPIV)370,530,420
 370  HTEMP = H(I,IPIV)
      H(I,IPIV) = COSINE*HTEMP + SINE*H(I,JPIV)
      IF ( X(I) - ABS( H(I,IPIV)) )380,390,390
 380  X(I) = ABS(H(I,IPIV))
      IQ(I) = IPIV
 390  H(I,JPIV) = -SINE*HTEMP + COSINE*H(I,JPIV)
      IF ( X(I) - ABS( H(I,JPIV)) ) 400,530,530
 400  X(I) = ABS(H(I,JPIV))
      IQ(I) = JPIV
      GO TO 530
C
 420  IF(I-JPIV)430,530,480
 430  HTEMP = H(IPIV,I)
      H(IPIV,I) = COSINE*HTEMP + SINE*H(I,JPIV)
      IF ( X(IPIV) - ABS( H(IPIV,I)) ) 440,450,450
 440  X(IPIV) = ABS(H(IPIV,I))
      IQ(IPIV) = I
 450  H(I,JPIV) = -SINE*HTEMP + COSINE*H(I,JPIV)
      IF ( X(I) - ABS( H(I,JPIV)) ) 400,530,530
C
 480  HTEMP = H(IPIV,I)
      H(IPIV,I) = COSINE*HTEMP + SINE*H(JPIV,I)
      IF ( X(IPIV) - ABS( H(IPIV,I)) ) 490,500,500
 490  X(IPIV) = ABS(H(IPIV,I))
      IQ(IPIV) = I
 500  H(JPIV,I) = -SINE*HTEMP + COSINE*H(JPIV,I)
      IF ( X(JPIV) - ABS( H(JPIV,I)) ) 510,530,530
 510  X(JPIV) = ABS(H(JPIV,I))
      IQ(JPIV) = I
 530  CONTINUE
C
C     TEST FOR COMPUTATION OF EIGENVECTORS
C
      IF(IEGEN)40,540,40
 540  DO 550 I=1,N
      HTEMP=U(I,IPIV)
      U(I,IPIV)=COSINE*HTEMP+SINE*U(I,JPIV)
 550  U(I,JPIV)=-SINE*HTEMP+COSINE*U(I,JPIV)
      GO TO 40
 1000 RETURN
      END
C             SUBROUTINE    INV FOR BMD06M            JUNE     30, 1966
      SUBROUTINE INV(A,B,N,T,KD)
      DIMENSION U(35),A(35,35),B(35,35)
      DIMENSION IN(35)
      DO 1 I=1,N
 1    IN(I)=0
      K=1
 2    DO 3 I=K,N
      U(I)=A(I,K)
 3    A(I,K)=0.0
      P=U(K)
      DO 4 I=1,K
      U(I)=A(K,I)
 4    A(K,I)=0.0
      U(K)=-1.0
      IN(K)=1
      H=0.0
      DO 5 I=1,N
      Y=U(I)/P
      DO 6 J=1,I
 6    A(I,J)=A(I,J)-U(J)*Y
      IF(IN(I))5,7,5
 7    IF(KD)12,12,13
 13   Z=B(I,I)-Y*(2.0*B(I,K)-Y*B(K,K))
      DO 8 J=1,N
      B(I,J)=B(I,J)-Y*B(K,J)
 8    B(J,I)=B(I,J)
      B(I,I)=Z
 12   IF(A(I,I)-H)5,5,9
 9    H=A(I,I)
      KK=I
 5    CONTINUE
      IF(KD)20,20,21
 21   Z=B(K,K)/P
      P=SQRT(P)
      DO 14 I=1,N
      B(I,K)=B(I,K)/P
 14   B(K,I)=B(I,K)
      B(K,K)=Z
 20   K=KK
      IF(H-T)10,10,2
 10   DO 91 I=1,N
      IF(IN(I))90,90,81
 90   DO 92 J=1,I
      A(I,J)=0.0
      A(J,I)=0.0
      IF(KD)92,92,93
 93   B(I,J)=0.0
      B(J,I)=0.0
 92   CONTINUE
      GO TO 91
 81   DO 11 J=1,I
      A(I,J)=-A(I,J)
 11   A(J,I)=A(I,J)
 91   CONTINUE
      RETURN
      END
C             SUBROUTINE   TPWD FOR BMD06M            JUNE     30, 1966
      SUBROUTINE TPWD(NT1,NT2)
      IF(NT1)40,10,12
 10   NT1=5
 12   IF(NT1-NT2)14,19,14
   14 IF(NT2.EQ.5) GO TO 19
      REWIND NT2
   19 IF(NT1-5)18,24,18
 18   IF(NT1-6)22,40,22
 22   REWIND NT1
 24   NT2=NT1
 28   RETURN
 40   WRITE (6,49)
      STOP
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      END
C             SUBROUTINE TRNGEN FOR BMD06M            JUNE     30, 1966
      SUBROUTINE TRNGEN (DATA,IPQ,NVG,NCASES,LCASE,NNEWA,LLCODE,LLVA,BBN
     1EW,MERRY,II)
      DIMENSION DATA(100),NNEWA(99),LLCODE(99),LLVA(99),BBNEW(99)
C     EXTERNAL SIGN
      REAL LOG10
      ASN(Q000FL)=ATAN(Q000FL/SQRT(1.0-Q000FL**2))
      SAMP=NCASES
      ITEM=II
      DO 3 J=1,NVG
  305 NEWA=NNEWA(J)
      LCODE=LLCODE(J)
      LVA=LLVA(J)
      BNEW=BBNEW(J)
      IF(LCODE-10) 4,4,5
    5 NEWB=BNEW
    4 D=DATA(LVA)
      IF(LCODE-41)500,180,500
  500 GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,
     1 150,160,170),LCODE
   10 IF(D)99,7,8
    7 DATA(NEWA)=0.0
      GO TO 3
    8 DATA(NEWA)=SQRT(D)
      GO TO 3
   20 IF(D)99,11,12
   11 DATA(NEWA)=1.0
      GO TO 3
   12 DATA(NEWA)=SQRT(D)+SQRT(D+1.0)
      GO TO 3
   30 IF(D)99,99,14
   14 DATA(NEWA)=ALOG10(D)
      GO TO 3
   40 DATA(NEWA)=EXP(D)
      GO TO 3
   50 IF(D)99,16,17
   16 DATA(NEWA)=0.0
      GO TO 3
   17 IF(D-1.0)18,19,99
   19 DATA(NEWA)=3.14159265/2.0
      GO TO 3
   18 A=SQRT(D)
      DATA(NEWA)=ASN(A)
      GO TO 3
   60 A=D/(SAMP+1.0)
      B=A+1.0/(SAMP+1.0)
      IF(A)99,23,24
   23 IF(B)99,26,27
   26 DATA(NEWA)=0.0
      GO TO 3
   27 DATA(NEWA)=ASN(SQRT(B))
      GO TO 3
   24 IF(B)99,28,29
   28 DATA(NEWA)=ASN(SQRT(A))
      GO TO 3
   29 A=SQRT(A)
      B=SQRT(B)
      DATA(NEWA)=ASN(A)+ASN(B)
      GO TO 3
   70 IF(D)31,99,31
   31 DATA(NEWA)=1.0/D
      GO TO 3
   80 DATA(NEWA)=D+BNEW
      GO TO 3
   90 DATA(NEWA)=D*BNEW
      GO TO 3
  100 IF(D)33,32,33
   32 DATA(NEWA)=0.0
      GO TO 3
   33 DATA(NEWA)=D**BNEW
      GO TO 3
  110 DATA(NEWA)=D+DATA(NEWB)
      GO TO 3
  120 DATA(NEWA)=D-DATA(NEWB)
      GO TO 3
  130 DATA(NEWA)=D*DATA(NEWB)
      GO TO 3
  140 IF(DATA(NEWB))34,99,34
   34 DATA(NEWA)=D/DATA(NEWB)
  150 BNEW=NEWB
      IF(D-BNEW) 151,152,152
  152 DATA(NEWA)=1.0
      GO TO 3
  151 DATA(NEWA)=0.0
      GO TO 3
  160 IF(D-DATA(NEWB)) 161,162,162
  162 DATA(NEWA)=1.0
      GO TO 3
  161 DATA(NEWA)=0.0
  170 IF(-D)175,99,99
  175 DATA(NEWA)=ALOG(D)
      GO TO 3
  180 IF(D)3,503,3
  503 IF(SIGN(10.0,D)) 504,3,3
  504 DATA(NEWA)=BNEW
    3 CONTINUE
      GO TO 42
   99 LCASE=-999
      IF(MERRY-J) 402,401,402
  402 MERRY=J
      WRITE (6,1404)J
  401 WRITE (6,1405)ITEM
      WRITE (6,1408)
      NCASES=NCASES-1
 1408 FORMAT(45H0THIS CASE WILL BE DELETED FOR ALL VARIABLES )
 1404 FORMAT(30H0THE INSTRUCTIONS INDICATED ON/25H TRANS GENERATOR CARD   
     1NO.I2,4H RE-/29H SULTED IN THE VIOLATION OF A/31H RESTRICTION FOR   
     2THIS TRANSFOR-/31H MATION. THE VIOLATION OCCURRED/27H FOR THE CASE   
     3 LISTED BELOW./)
 1405 FORMAT( 9H CASE NO.I5)
   42 RETURN
      END
C             SUBROUTINE   VCTR FOR BMD06M            JUNE     30, 1966
      SUBROUTINE VCTR(A,V,N,ALPHA)
      DIMENSION A(70,70),V(70)
      A(1,1)=A(1,1)-ALPHA
 6    DO 15 I=2,N
      A(I,I)=A(I,I)-ALPHA
 70   II=I-1
 7    DO 15 J=1,II
 8    IF (A(I,J))9,15,9
 9    IF (ABS(A(J,J))-ABS(A(I,J)))11,10,10
 10   R=A(I,J)/A(J,J)
      GO TO 130
 11   R=A(J,J)/A(I,J)
      DO 12 K=1,N
      C=A(J,K)
      A(J,K)=A(I,K)
 12   A(I,K)=C
 130  JJ=J+1
 13   DO 14 K=JJ,N
 14   A(I,K)=A(I,K)-R*A(J,K)
 15   CONTINUE
      C=A(N,N)
      V(N)=1.0
      DO 29 I=2,N
      JJ=N-I+1
      R=0.0
      II=N-I+2
      DO 25 K=II,N
 25   R=R+A(JJ,K)*V(K)
      IF (ABS(A(JJ,JJ))-1.0E-10)27,27,28
 27   V(JJ)=1.0
      C=0.0
      DO 26 J=II,N
 26   V(J)=0.0
      GO TO 29
 28   V(JJ)=(C-R)/A(JJ,JJ)
 29   CONTINUE
       RETURN
      END




















