C        DESCRIPTION OF STRATA WITH HISTOGRAMS         MAY 20, 1966
C
C             UCLA MEDICAL CENTER - HEALTH SCIENCES COMPUTING FACILITY
C        THIS IS A SIFTED VERSION OF BMD07D ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
C
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3 XMN(100), SX(10), DX(10), XL(10), SXX(100,10)
      DIMENSION FNG(10)
       DIMENSION OUT(17)
      DOUBLE PRECISION     OCC,FRMT(120),PROBLM,FINISH,SPVALU,GRPDIV,COB,JBCD
       DOUBLE PRECISION TRYME(10), BLANK
       DATA OUT(1),OUT(2),OUT(4),OUT(5),OUT(6),OUT(7),OUT(8),OUT(10),
     * OUT(11),OUT(12),OUT(14),OUT(15),OUT(16),OUT(17) / 4H(16X, 4H,   ,
     * 4H('GR, 4HOUP', 4H,I3,, 4H3X)/, 4H16X,,4H(A8,, 4H3X)/, 4H15X,,
     * 4H('.., 4H...., 4H...., 4H+')) / , YES / 3HYES /
      DATA FNG/'1','2','3','4','5','6','7','8','9','10'/
       DATA BLANK / '        ' /
      DATA PROBLM,FINISH,SPVALU,GRPDIV,SPCVAL/6HPROBLM,6HFINISH,6HSPCVAL
     1,6HGRPDIV,4H-7CV/
       COMMON / LABS / TRYME
      COMMON  ZZ     , X
      COMMON  NVAR   , NCAS   , NTR    , CP     , IGR    , NCP
      COMMON  NCI    , KCICP  , NERR   , NG     , NADD   , IPNT
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
C
C           CP       GROUP CUT POINTS
C           IPNT     PRINT CORRELATION MATRIX
C           IAO      PRINT ALL OF ORDERED DATA
C           IPA      PRINT ALL OF INPUT DATA
C           IPO      PRINT PART OF ORDERED INPUT
C           IPP      PRINT PART OF INPUT DATA
C           ISV      NO. OF VARIABLE WITH SPECIAL VALUE
C           IVF      NO. OF VARIABLE FORMAT CARDS
C           JBCD     PROBLEM NAME
C           NADD     NO. OF VARIABLES ADDED IN TRANSGENERATION
C           NCAS     NO. OF CASES
C           NCI      NO. OF CLASS INTERVALS
C           NCP      NO. OF GROUP CUT POINTS
C           NNC      NO. OF NAME CARDS OR VARIABLES WITH SPECIFIED NAMES
C           NSVC     NO. OF SPECIAL VALUE CARDS
C           NTR      NO. OF TRANS GENERATOR CARDS
C           NVAB     NO. OF BASE VARIABLE FOR GROUP DIVISION
C           NVAR     NO. OF VARIABLES
C
C
  17  FORMAT ('1BMD07D - DESCRIPTION OF STRATA WITH HISTOGRAMS',
     * ' - REVISED MAY 10, 1968' /
     241H HEALTH SCIENCES COMPUTING FACILITY, UCLA//)
      NTAPE=5
	CALL USAGEB('BMD07D')
C
  998 READ (5,10)COB,JBCD,NVAR,NSVC,NNC,NCAS,NG,NCP,NCI,NTR, NADD, PA, P
     1P,AO,PO,PNT,SNGLE,NULABS,MTAPE,IVF
      IF(COB.EQ.FINISH)GO TO 999
      IF(COB.EQ.PROBLM)GO TO 50
      WRITE(6,10012)
10012 FORMAT(49H PROBLEM CARD OUT OF ORDER OR INCORRECTLY PUNCHED)
 999  IF(5-NTAPE)9998,9999,9999
 9998 REWIND NTAPE
 9999 STOP
   50 CALL TPWD(MTAPE,NTAPE)
    3 IF(IVF .GT. 0.AND. IVF .LE. 10) GO TO 1152
      WRITE (6,888)
      IVF=1
 1152 IVF = IVF*12
      READ (5,12)(FRMT(I),I=1,IVF)
 53   L1=1
 1153 IF(NG-10)1155,1155,10000
10000 WRITE(6,10001)
10001 FORMAT(28H NUMBER OF GROUPS EXCEEDS 10)
      GO TO 1151
 1155 IF(NCI-30)1156,1156,10002
10002 WRITE(6,10003)
10003 FORMAT(47H NUMBER OF CLASSES SPECIFIED IS GREATER THAN 30)
      GO TO 1151
 1156 IF(NCI-5)1157,1158,1158
 1157 IF(NCI)10004,1158,10004
10004 WRITE(6,10005)
10005 FORMAT(43H NUMBER OF CLASSES SPECIFIED IS LESS THAN 5)
      GO TO 1151
 1158 L2 = NVAR
      IF(-NADD)2158,2158,1150
 2158 L3=1+NADD
      NVAR=NVAR+NADD
      IF(NCAS*NVAR-4000)1159,1159,11006
11006 WRITE(6,11007)
11007 FORMAT(' DATA MATRIX IS TOO LARGE, CANNOT EXCEED 16,000')
      GO TO 999
10006 WRITE(6,10007)
10007 FORMAT(19H TOO MANY VARIABLES)
      GO TO 999
 1159 IF(NVAR-100)1154,1154,10006
 1154 DO 531 J=1,NCAS
      READ (NTAPE,FRMT)(X(I),I=L1,L2)
      L1=L2+L3
      L2=L2+NVAR
 531  CONTINUE
      DO 533 I=1,NVAR
      OCC(I) =0 .0
      BIGX(I) =0.0
      SMLX(I) =0.0
      NS(I)=0
      DO 532 J=1,5
      SPV(I,J)=0.0
 532  CONTINUE
 533  CONTINUE
      IF (NSVC) 5305, 5305, 54
   54 IF(NSVC-NVAR)541,541,10010
10010 WRITE(6,10011)
10011 FORMAT(' NUMBER OF SPECIAL VALUES INCORRECTLY SPECIFIED')
      GO TO 1151
C
C     READ IN SPECIAL VALUE CARDS
C
 541  DO 55 I=1,NSVC
      READ (5,1300)COB,KM,NVIV,(SPV(KM,K), K=1,NVIV)
      IF(COB.NE.SPVALU)GO TO 10013
  555 IF(NVIV-5)556,556,20014
20014 WRITE(6,10011)
      GO TO 1151
10013 WRITE(6,10014)
10014 FORMAT(55H SPECIAL VALUE CARD OUT OF ORDER OR INCORRECTLY PUNCHED
     1)
      GO TO 999
  556 ISKIP=1
      IF(SNGLE.NE.YES)GO TO 5567
 5563 ISKIP=2
 5567 IF(KM-100)557,557,10016
10016 WRITE(6,10017)
10017 FORMAT(50H INDEX OF VARIABLE WITH SPECIAL VALUES EXCEEDS 100)
      GO TO 1151
  557 NS(KM)=NVIV
      DO 5460 II = 1, NVIV
      DO 5450 J = II, NVIV
      IF(SPV(KM,II).GE.SPV(KM,J))GO TO 5450
 5425 HOLD = SPV(KM,II)
      SPV(KM,II) = SPV(KM,J)
      SPV(KM,J) = HOLD
 5450 CONTINUE
 5460 CONTINUE
 55   CONTINUE
      GO TO (5305,5500),ISKIP
 5500 DO 5501 I=1,NVAR
      NS(I)=NVIV
      DO 5501 J=1,NVIV
 5501 SPV(I,J)=SPV(KM,J)
 5305 IF(NG)10018,5315,5306
10018 WRITE(6,10019)
10019 FORMAT(36H NUMBER OF GROUPS CANNOT BE NEGATIVE)
      GO TO 1151
 5315 IF(-NCP)560,275,10020
10020 WRITE(6,10021)
10021 FORMAT(46H NUMBER OF GROUP CUT POINTS CANNOT BE NEGATIVE)
      GO TO 1151
C
C     READ IN GROUP DIVISION CARD
C
  560 READ (5,14)COB,NVAB,(CP(I),I=1,NCP)
      IF(COB.NE.GRPDIV)GO TO 20020
  565 IF(NVAB-100)570,570,10022
10022 WRITE(6,10023)
10023 FORMAT(41H NO. OF BASE VARIABLE IS GREATER THAN 100)
      GO TO 1151
20020 WRITE(6,20021)
20021 FORMAT(47H GROUP DIVISION CARD OUT OF ORDER OR MISPUNCHED)
      GO TO 999
  570 GO TO (572,571),ISKIP
  571 NS(NVAB)=0
  572 IF(NS(NVAB))5306,5306,270
 5306 WRITE (6,17)
      IF (NG) 101, 100, 101
 101  WRITE (6,19)NCAS
      GO TO 105
C
 100  WRITE (6,18)NVAB, NCAS, NCP
      DO 1001 I=1,NCP
      WRITE (6,45)I, CP(I)
 1001 CONTINUE
 105  WRITE (6,20)
      WRITE (6,21)JBCD,NCAS,NVAR,NG,NCP,NSVC,NCI,NNC, PA,  PP, AO, PO, P
     1NT
      WRITE(6,35000)(FRMT(I),I=1,IVF)
35000 FORMAT(' VARIABLE FORMAT CARD(S)'/10(1X,12A6/))
      IF (NTR) 56,56,5300
 5300 IF(NTR-99)5301,5301,10024
10024 WRITE(6,10025)
10025 FORMAT(48H NO. OF TRANSGENERATION CARDS IS GREATER THAN 99)
      GO TO 999
 5301 CALL TVA
      IF(-NERR)56,56,999
 56   NVNC=NVAR*NCAS
  205 IF(PA.NE.YES)GO TO 225
  210 WRITE (6,41)
      ASSIGN 5302 TO KBRNCH
  214 L4=NCAS
  215 WRITE (6,47)
      L3=1
      L1=1
      L2=NVAR
  220 DO 2110 J=L3,L4
      WRITE (6,40)J,(X(I),I=L1,L2)
      L1=L2+1
      L2=L2+NVAR
 2110 CONTINUE
      GO TO KBRNCH,(5302,240,71,249)
C
C     PRINT FIRST AND LAST TEN ROWS OF INPUT MATRIX AFTER
C     TRANSGENERATION, IF DESIRED
C
  225 IF(PP.NE.YES)GO TO 5302
 230  IF(20-NCAS)236,210,210
  236 WRITE (6,41)
      ASSIGN 240 TO KBRNCH
  238 WRITE (6,43)
      L4=10
      GO TO 215
C
  240 ASSIGN 5302 TO KBRNCH
  242 WRITE (6,42)
      L1=NVNC-(10*NVAR)+1
      L2=L1+NVAR-1
      L3=NCAS-9
      L4=NCAS
      GO TO 220
C
 5302 IF(-NG)60,63,63
 60   NC=NCAS
      NCP=NG-1
      NGG=NG
      DO 61 I=1,NG
      IGR(I)=NC/(NGG)
      NC=NC-IGR(I)
      NGG=NGG-1
 61   CONTINUE
      GO TO 70
C
 63   CALL DIVIDE (NVAB)
C
C     PRINT ORDERED MATRIX, IF DESIRED
C
   70 IF(AO.NE.YES)GO TO 325
  212 WRITE (6,46)
      ASSIGN 71 TO KBRNCH
      GO TO 214
C
C     PRINT FIRST AND LAST TEN ROWS OF ORDERED MATRIX, IF DESIRED
C
  325 IF(PO.NE.YES)GO TO 71
 330  IF(20-NCAS)248,212,212
 248  WRITE (6,46)
      ASSIGN 249 TO KBRNCH
      GO TO 238
C
  249 ASSIGN 71 TO KBRNCH
      GO TO 242
C
 71    DO 3681 IZQ=1,10
 3681  TRYME(IZQ) = BLANK
       CALL RDLBL(NNC,NVAR,OCC)
       IF (NULABS.LE.0) GO TO 74
       READ (5,8912) TRYME
 8912  FORMAT (10A8)
 74   NG=NCP+1
      DO 997 I=1,10
      INTEG(I)=I
 997  CONTINUE
      DO 900 NV=1, NVAR
      NCICP(NV)=0
      CALL MAMIN(NV)
      KCICP=0
      CALL CLINT(NV)
      IF(KCICP)999,106,999
  106 IF(OCC(NV)-0.0)75,76,75
 75   WRITE (6,22)NV, OCC(NV)
      GO TO 78
C
 76   WRITE (6,23)NV
  78  WRITE (6,16)
      L=NG
       OUT(3) = FNG(L)
       OUT(9) = FNG(L)
       OUT(13) = FNG(L)
 79    WRITE (6,OUT) (I,I=1,NG), (TRYME(I), I=1,NG)
 90   JY=NS(NV)
      J=NCICP(NV)
      IF(-JY)91,911,911
C
C            NTR IS USED TO GIVE RITE THE NUMBER OF THE VARIABLE
C
   91 NTR = NV
      WRITE (6,1112)
      NABC = 1
      CALL PLUG (NV,NABC)
      DO 910 JJ = 1, JY
      CALL RITE (JJ,NG,NABC)
  910 CONTINUE
 911  WRITE (6,1111)
      NABC=2
      CALL PLUG(NV,NABC)
      GO TO 93
   92 J=J-1
   93 CALL RITE(J,NG,NABC)
      IF (J-1) 95, 95, 92
 95   R=0.0
      S=0.0
      CAS=NCAS
      CAS1=CAS-1.0
      YMN=0.0
      YSD=0.0
      DO 9100 I=1,NG
      SX(I)=0.0
      DX(I)=0.0
      XL(I)=0.0
 9100 CONTINUE
      J1=1
      SSW=0.0
      DO 9200 I=1,NG
      II=I-1
      IG=IGR(I)
      GI=IGR(I)
      IF(-II)9106,9107,9107
 9106 J1=J1+IGR(II)
 9107 IG1=J1+IG-1
      IF(-IG)9104,9200,9200
 9104 DO 9150 J=J1,IG1
      K=(J -1)*NVAR+NV
      IF(-JY)9110,9125,9125
 9110 IF(X(K).EQ.SPCVAL)GO TO 9119
 9125 SX(I)=SX(I)+X(K)
      GO TO 9150
C
 9119 XL(I)=XL(I)+1.0
 9150 CONTINUE
       IF (GI.EQ.XL(I)) SX(I) = 0.0
       IF (GI.NE.XL(I)) SX(I) = SX(I)/(GI-XL(I))
C
C     MEAN OF GROUP I IN SX(I)---SPECIAL VALUES EXCLUDED
C
      DO 9175 J=J1,IG1
      K=(J -1)*NVAR+NV
      IF (JY) 9170, 9170, 9160
 9160 IF(X(K).EQ.SPCVAL)GO TO 9175
 9170 DX(I)=DX(I)+(X(K)-SX(I))**2
 9175 CONTINUE
      SSW=SSW+DX(I)
       IF (GI.LE.XL(I)+1.0) DX(I) = 0.0
       IF (GI.GT.XL(I)+1.0) DX(I) = SQRT(DX(I)/(GI-XL(I)-1.0))
C
C     STANDARD DEVIATION OF GROUP I IN DX(I)
C
 9200 CONTINUE
      DO 125 J=1,NCAS
      K=(J-1)*NVAR+NV
      IF (JY) 109, 121, 109
  109 IF(X(K).EQ.SPCVAL)GO TO 110
 121  YMN=YMN+X(K)
      GO TO 125
C
 110  R=R+1.0
 125  CONTINUE
      YMN=YMN/(CAS-R)
      XMN(NV)=YMN
      DO 145 J=1,NCAS
      K=(J-1)*NVAR+NV
      IF(-JY)128,141,141
  128 IF(X(K).EQ.SPCVAL)GO TO 130
 141  YSD=YSD+(X(K)-YMN)**2
      GO TO 145
C
 130  S=S+1.0
 145  CONTINUE
      IDFT=CAS1-S
      IDFG=NG-1
      IDFW=IDFT-IDFG
      SST=YSD
      SSG=SST-SSW
      IF(IDFG.NE.0) SGM=SSG/IDFG
      IF(IDFG.EQ.0) SGM=0.0
      IF(IDFW.NE.0) SWM=SSW/IDFW
      IF(IDFW.EQ.0) SWM=0.0
       IF (SWM.EQ.0) FST = 0.0
       IF (SWM.NE.0) FST = SGM/SWM
      IF(CAS1.NE.S) YSD=SQRT(YSD/(CAS1-S))
      IF(CAS1.EQ.S) YSD=0.0
      WRITE (6,36)( SX(I), I=1,NG)
      WRITE (6,37)( DX(I), I=1,NG)
      DO 150 I=1,NG
      LX=XL(I)
      XL(I)=IGR(I)-LX
 150  CONTINUE
      WRITE (6,38)( XL(I), I=1,NG)
      WRITE (6,1010)YMN,SSG,IDFG,SGM,FST,YSD,SSW,IDFW,SWM,BIGX(NV),SST,I
     1DFT,SMLX(NV)
 900  CONTINUE
  250 IF(PNT.NE.YES)GO TO 998
  255 CALL CORR
      GO TO 998
C
 270  WRITE (6,15)NS(NVAB)
      GO TO 999
C
 275  WRITE (6,24)
      GO TO 999
C
 1150 WRITE (6,1501)
 1151 WRITE (6,1500)
      GO TO 999
C
 10    FORMAT(2A6,3I3,I4,I2,I1,3I2,6A3,15X,I1,2I2)
 11   FORMAT (1H1)
  12  FORMAT(12A6)
  14  FORMAT(A6,I3,9F7.0)
 15   FORMAT(1H0,9HTHERE ARE,I4,86H SPECIAL VALUES FOR THE VARIABLE ON W
     1HICH DATA IS STRATIFIED. PROGRAM CANNOT CONTINUE.)
  16  FORMAT(75X32HLOWER LIMITS OF CLASS INTERVALS)//)
  18  FORMAT(1H034HDATA STRATIFIED ON VARIABLE NUMBERI4,1H,I5,6H CASES/I
     16,21H CUT POINTS SPECIFIED//)
  19  FORMAT(1H0,34HDATA STRATIFIED BY ORDER OF ENTRY.I5,6H CASES)
  20  FORMAT(13H0PROBLEM CARD)
  21  FORMAT(1H0,8X,12HPROBLEM NAME25XA6/9X15HNUMBER OF CASES21XI7/9X19H
     1NUMBER OF VARIABLES17XI7/9X26HNUMBER OF GROUPS SPECIFIED10XI7/9X26
     2 HNUMBER OF GROUP CUT POINTS,I17/9X29HNUMBER OF SPECIAL VALUE CARD
     3S7XI7/9X35HNUMBER OF CLASS INTERVALS SPECIFIED3XI5/9X30HNUMBER OF
     4VARIABLES WITH NAMES,6X,I7,//9X25HPRINT ENTIRE INPUT MATRIX15XA3/9
     5X37HPRINT FIRST AND LAST 10 ROWS OF INPUT3XA3/9X27HPRINT ENTIRE OR
     6DERED MATRIX13XA3/9X40HPRINT FIRST AND LAST 10 ROWS OF ORDERED A3/
     79X24HPRINT CORRELATION MATRIX16XA3)
 22   FORMAT(1H1/40X8HVARIABLEI4,2X2H( A6,2H )10X34H(PRINTED INTERVAL DE
     1SIGNATIONS ARE)
 23   FORMAT(1H1/40X8HVARIABLEI4,22X34H(PRINTED INTERVAL DESIGNATIONS AR
     1E)
 24   FORMAT(1H0,11X,94HBOTH THE NUMBER OF GROUPS AND THE NUMBER OF GROU
     1P CUT POINTS ARE ZERO. PROGRAM CANNOT PROCEED.)
  36  FORMAT(1H05X5HMEAN 10F11.3)
  37  FORMAT(6X5HS DEV10F11.3)
  38  FORMAT(7X1HN10F11.0,//)
  40  FORMAT(I5,8F14.5,/(5X,8F14.5))
  41  FORMAT(1H136X37HINPUT MATRIX  (AFTER TRANSGENERATION))
  42  FORMAT(3(3X1H./)/)
  43  FORMAT(41X29H-FIRST AND LAST 10 ROWS ONLY-)
  45  FORMAT(1H I4,1H.F14.4)
  46  FORMAT(1H136X47HINPUT MATRIX  (ORDERED - AFTER TRANSGENERATION))
  47  FORMAT(3X3HROW//)
 1010 FORMAT(46H-ALL GROUPS COMBINED (SPECIAL VALUES EXCLUDED) 20X,46HSU
     1M OF SQUARES    DF    MEAN SQUARE    F RATIO/1H0,13X,7HMEAN   F12.
     24,24X,7HBETWEENF13.4,I9,F14.4,F11.4/13X,8H S DEV  F12.4,24X,7HWITH
     3IN F13.4,I9,F14.4/13X,8H MAXIMUMF12.4,24X,7HTOTAL  F13.4,I9/13X,8H
     4 MINIMUMF12.4)
 1111 FORMAT (1H013X,64HTABULATIONS AND COMPUTATIONS WHICH FOLLOW EXCLUD
     1E SPECIAL VALUES /6X8HINTERVAL//)
 1112 FORMAT(1H 5X7HSPECIAL/7X6HVALUES//)
 1300 FORMAT(A6,I3,I1,5F7.0)
 1500 FORMAT(1H0,12X75HCONTROL CARDS OUT OF ORDER OR NOT PROPERLY PUNCHE
     1D. PROGRAM CANNOT PROCEED.)
 1501 FORMAT(1H036X45HNUMBER OF VARIABLES ADDED CANNOT BE NEGATIVE.)
  888 FORMAT(1H0,23X,71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPEC
     1IFIED, ASSUMED TO BE 1.)
C
      END
C             SUBROUTINE CLINT FOR BMD07D              MAY 20, 1966
      SUBROUTINE CLINT (NV)
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3 XMN(100), SX(10), DX(10), XL(10), SXX(100,10)
      COMMON  ZZ
      COMMON  X      , NVAR   , NCAS   , NTR    , CP     , IGR
      COMMON  NCP    , NCI    , KCICP  , NERR   , NG     , NADD
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
C
C
      IF(NCI)10,10,11
 10   FNCA=NCAS
      G=NG
      CIN=10.0*ALOG10(FNCA/G)
      NCI=CIN
   11 IF(NCI-5)12,15,15
 12   NCI=5
   15 FNCI=NCI-1
      IF(30-NCI)16,18,18
 16   NCI=30
      FNCI=29.0
   18 IF(SMLX(NV)-BIGX(NV))185,22,185
  185 CALL SCALE1(SMLX(NV),BIGX(NV),FNCI,XMIN,XIJ)
      ICI=NCI
   19 NCICP(NV)=ICI
      CICP(1)=XMIN
      ICI=ICI+1
      DO 20 J=2,ICI
      CICP(J)=CICP(J-1)+XIJ
   20 CONTINUE
      RETURN
   22 ICI=5
      XIJ=0.5
      XMIN=SMLX(NV)-1.0
C
C     PRINT ALL OF INPUT MATRIX AFTER TRANSGENERATION, IF DESIRED
C
      GO TO 19
      END
      SUBROUTINE SCALE1(YMIN,YMAX,YINT,TYMIN,YIJ)
      DIMENSION C(10)
      DATA C/1.0,1.5,2.0,3.0,4.0,5.0,7.5,10.0,15.0,20.0/
      TEST=1.0/(2.0**20)
   50 YR=YMAX-YMIN
      TT=YR/YINT
      J=ALOG10(TT)+TEST
      E=10.0**J
      TT=TT/E
      I=0
      IF(TT-1.0+TEST)205,201,201
  205 TT=TT*10.0
      E=E/10.0
  201 I=I+1
      IF(10-I)1,2,2
    1 E=E*10.0
      TT=TT/10.0
      I=1
    2 IF(TT-C(I))233,202,201
  233 YIJ=C(I)*E
      GO TO 203
  202 Y=YMIN/C(I)
      J=Y
      T=J
      IF(0.0001-ABS(T-Y))204,233,233
  204 YIJ=C(I+1)*E
  203 X=((YMAX+YMIN)/YIJ-YINT )/2.0+.00001
      K=X
      IF(K)235,234,240
  235 Y=K
      IF(X-Y)236,240,236
  236 K=K-1
  240 TYMIN=K
      TYMIN=YIJ*TYMIN
      TYMAX=TYMIN + YINT*YIJ
      IF(YMAX-TYMAX-TEST)10,10,201
  234 IF(YMIN)235,240,240
   10 RETURN
      END
C             SUBROUTINE CORR  FOR BMD07D              MAY 20, 1966
      SUBROUTINE CORR
C
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3XMN(100),SX(10),DX(10),XL(10),SXX(100,10)
       DOUBLE PRECISION TRYME(10), BLANK
       COMMON / LABS / TRYME
       DATA BLANK / '        ' /
      COMMON  ZZ     , X
      COMMON  NVAR   , NCAS   , NTR    , CP     , IGR    , NCP
      COMMON  NCI    , KCICP  , NERR   , NG     , NADD   , IPNT
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
C
C     SUBROUTINE CORR COMPUTES CORRELATION COEFFICIENTS FOR EACH
C      GROUP AND PRINTS THEM OUT IN THE FORM OF THE LOWER TRIANGLE
C      INCLUDING THE PRINCIPAL DIAGONAL.  SPECIAL VALUES ARE EXCLUDED
C     FROM THE COMPUTATIONS.
C
      DATA SPCVAL/4H-7CV/
C
  85  IGF = 1
      LAST = (NVAR*(1+NVAR))/2
      WRITE (6,9101)NG
      ASSIGN 100 TO KSKIP
      DO 40 NV=1,NG
      DO 810 I = 1,LAST
      ZZ(I) = 0.0
  810 CONTINUE
 811  IG = IGR(NV)
      KKK = (IGF - 1)*NVAR
      IKK = KKK
      IGF = IG + IGF
 255  IF (IG) 256,256,7
 256   IF (TRYME(NV).EQ.BLANK     ) WRITE(6,1000) NV
       IF (TRYME(NV).NE.BLANK     ) WRITE(6,1001) NV,TRYME(NV)
      ASSIGN 105 TO KSKIP
      GO TO 40
C
 7    KK=(IGF-2)*NVAR
      JKK = KK
 8    M=0
      DO 1 I = 1, NVAR
      SUMX = 0.0
      SUMSQX = 0.0
      XIG=0.0
      KK = KK + 1
      KKK = KKK + 1
      DO 6 K = KKK,KK,NVAR
      IF(X(K).EQ.SPCVAL)GO TO 6
 15   SUMX = X(K) + SUMX
      SUMSQX = X(K)**2 + SUMSQX
  17  XIG=XIG+1.0
    6  CONTINUE
      DO 2 J = 1, I
      SUMX1=SUMX
 60   SUMY1=0.0
      SUMY = 0.0
      SUMSQY = 0.0
      SUMXY = 0.0
      YIG=0.0
      XYIG=0.0
      IK = I + IKK
      JKIG = J + JKK
      KKJ = J + IKK
      DO 5 K = KKJ,JKIG,NVAR
      IF(X(K).EQ.SPCVAL)GO TO 705
 65   SUMY = X(K) + SUMY
      SUMSQY = X(K)**2 + SUMSQY
  68  YIG=YIG+1.0
      IF(X(IK).EQ.SPCVAL)GO TO 715
  71  SUMXY = X(IK)*X(K) + SUMXY
  72  XYIG=XYIG+1.0
      GO TO 75
C
  705 IF(X(IK).EQ.SPCVAL)GO TO 75
 710  SUMX1=SUMX1-X(IK)
      GO TO 75
C
 715  SUMY1=SUMY1+X(K)
  75  IK = IK + NVAR
    5  CONTINUE
      SUMY1=SUMY-SUMY1
 50   M = M + 1
      IF(XIG*YIG*XYIG)51,505,51
 51   A=SQRT((SUMSQX-(SUMX**2)/XIG)*(SUMSQY-(SUMY**2)/YIG))
      IF(A)515,505,515
 515  IF(XYIG-1.0)550,505,550
 550  ZZ(M)=((SUMXY-(SUMX1*SUMY1)/XYIG)* (SQRT((XIG-1.0)*(YIG-1.0))))/(A
     1*(XYIG-1.0))
      GO TO 2
C
 505  ZZ(M)=99.0
    2 CONTINUE
    1 CONTINUE
       NPAGE = 0
      KVAR = (NVAR+7)/8
      ISTART = 1
      LLL = 0
      GO TO KSKIP,(100,105)
 100  ASSIGN 105 TO KSKIP
      GO TO 110
C
 105  WRITE (6,9105)
 110  DO912 II = 1, KVAR
      IF(II-1)90,900,90
 90   WRITE (6,9105)
 900  KOUNT = LLL
      LL = LLL + 9
       NPAGE = NPAGE + 1
       WRITE (6,9103)  TRYME(NV), NV, NPAGE
      KK = II-1
      KK = KK-5*(KK/5)
      IF(-KK)99,913,913
   99 DO 97 III= 1, KK
      WRITE (6,9104)
   97 CONTINUE
  913 I = ISTART
       IBEGN = KOUNT+1
       IFINSH = MIN0(NVAR,IBEGN+7)
       WRITE (6,9136) (MXL, MXL=IBEGN,IFINSH)
 9136  FORMAT (I16,7I14)
      DO91 IJ = 1, 8
      K = IJ + I- 1
       KOUNT = KOUNT+1
       WRITE (6,9102) KOUNT, (ZZ(N),N=I, K)
      I = KOUNT + I
      IF (KOUNT - NVAR) 91,40,40
   91 CONTINUE
      ISTART = K + LL
      GO TO 92
C
  910 DO 93 L = 1, 8
      JK = I + 7
       KOUNT = KOUNT+1
       WRITE (6,9102) KOUNT, (ZZ(N),N=I,JK)
      I = I + KOUNT
      IF ( KOUNT - NVAR ) 93,911,911
   93 CONTINUE
   92 IF (40 - KOUNT) 914,94,910
  914 IF ( KOUNT - 80) 910,94,910
   94 NPAGE = NPAGE + 1
      WRITE (6,9105)
       WRITE (6,9103)  TRYME(NV), NV, NPAGE
       IBEGN = KOUNT+1
       IFINSH = MIN0(NVAR,IBEGN+7)
       WRITE (6,9136)(MXL,MXL=IBEGN,IFINSH)
      GO TO910
C
  911 LLL   = LLL+8
  912 CONTINUE
  40  CONTINUE
      RETURN
C
 1000 FORMAT(1H0,5HGROUPI3,41H IS EMPTY, THEREFORE, THERE IS NO MATRIX.)
 1001  FORMAT(1X,'GROUP',I3,', ',A8,', IS EMPTY, THEREFORE, THERE IS NO
     *MATRIX.')
 9101 FORMAT (1H1 40X, 22HCORRELATION MATRICES ( I3,1H)//10X,102HTHESE A
     1RE SYMMETRIC MATRICES AND THEREFORE ARE PRINTED IN HALF FORM, INCL
     2UDING THE PRINCIPAL DIAGONAL./10X86HA VALUE OF 99.0 INDICATES THE
     3COEFFICIENT IS NOT COMPUTABLE DUE TO A ZERO DENOMINATOR.//)
 9102  FORMAT(I4,8F14.4)
 9103  FORMAT (7X,A8,'  GROUP',I3,6X,'(PAGE',I3,')'//)
 9104 FORMAT (1H ////////)
 9105 FORMAT(1H1/////)
      END
C             SUBROUTINE DESQ  FOR BMD07D              MAY 20, 1966
      SUBROUTINE DESQ (L1, DEVSQ)
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3XMN(100),SX(10),DX(10),XL(10),SXX(100,10)
      COMMON  ZZ     , X
      COMMON  NVAR   , NCAS   , NTR    , CP     , IGR    , NCP
      COMMON  NCI    , KCICP  , NERR   , NG     , NADD
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
      EXTERNAL SIGN
      CAS=NCAS-1
      CALL YMEAN (L1, YMEA )
      L=NS(L1)
      SUM=0.0
      DO 100 J=1,NCAS
      I=(J-1)*NVAR+L1
      IF (L) 52, 52, 45
 45   DO 50 K=1,L
      IF(X(I)-SPV(L1,K))50,46,50
 46   B=SIGN(1.0,SPV(L1,K))
      C=SIGN(1.0, X(I))
      IF(B+C)55,50,55
 50   CONTINUE
   52 SUM=SUM+(X(I)-YMEA)**2
      GO TO 100
C
 55   CAS=CAS-1.0
 100  CONTINUE
      DEVSQ=SQRT((1.0/CAS)*SUM)
      RETURN
      END
C             SUBROUTINE DIVIDE FOR BMD07D             MAY 20, 1966
      SUBROUTINE DIVIDE (NVAB)
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3 XMN(100), SX(10), DX(10), XL(10), SXX(100,10)
      COMMON  ZZ
      COMMON  X      , NVAR   , NCAS   , NTR    , CP     , IGR
      COMMON  NCP    , NCI    , KCICP  , NERR   , NG     , NADD
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
C
C         MATRIX ORDERING BASED ON ONE VARIABLE - NVAB
C
      NVNC=NVAR*NCAS
      NVN=NVNC-NVAR-NVAR+NVAB
      DO 1000 K=NVAB,NVN,NVAR
      KK=K+NVAR
      L=K-NVAB+1
      M=L
      CHECK=X(K)
 75   IF(X(KK)-CHECK)100,115,115
 100  M=KK-NVAB+1
      CHECK=X(KK)
 115  KK=KK+NVAR
 110  IF(NVNC-KK)200,75,75
 200  IF(M-L)201,1000,201
 201  DO 205 J=1,NVAR
      DUM=X(L)
      X(L)=X(M)
      X(M)=DUM
      L=L+1
      M=M+1
 205  CONTINUE
 1000 CONTINUE
      K=NCP+1
      M=K
      CP(K)=10.0**6
      DO 63 I=1,K
      IGP(I)=0
 63   CONTINUE
      I=NVN+NVAR
 64   DO 66 J=NVAB,I,NVAR
      IF (X(J)-CP(K)) 65, 66, 66
 65   IGP(K)=IGP(K)+1
 66   CONTINUE
      K=K-1
      IF (K) 67, 67, 64
 67   DO 68 K=2,M
      L=K-1
      IGR(K)=IGP(K)-IGP(L)
 68   CONTINUE
      IGR(1)=IGP(1)
      RETURN
      END
C             SUBROUTINE MAMIN  FOR BMD07D             MAY 20, 1966
      SUBROUTINE MAMIN (NV)
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3 XMN(100), SX(10), DX(10), XL(10), SXX(100,10)
      COMMON  ZZ
      COMMON  X      , NVAR   , NCAS   , NTR    , CP     , IGR
      COMMON  NCP    , NCI    , KCICP  , NERR   , NG     , NADD
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
      EXTERNAL SIGN
      TENSIX=10.0**6
      SMLX(NV)=TENSIX
      BIGX(NV)=-SMLX(NV)
      JJ=NS(NV)
      IF(JJ)10,10,15
 10   ASSIGN 52 TO JTIMES
      GO TO 20
C
 15   ASSIGN 50 TO JTIMES
 20   J=((NCAS-1)*NVAR)+NV
      DO 60 K=NV,J,NVAR
      GO TO JTIMES,(52,50)
 50   DO 51 I=1,JJ
      IF (X(K)-SPV(NV,I)) 51, 501, 51
 501  B=SIGN(1.0, SPV(NV,I))
      C=SIGN(1.0, X(K))
      IF (B+C) 60, 51, 60
 51   CONTINUE
 52   BIGX(NV)=AMAX1(BIGX(NV),X(K))
      SMLX(NV)=AMIN1(SMLX(NV), X(K))
 60   CONTINUE
      IF(SMLX(NV)-TENSIX)100,70,100
 70   IF(BIGX(NV)+TENSIX)100,80,100
 80   SMLX(NV)=0.0
      BIGX(NV)=0.0
 100  RETURN
      END
C          SUBROUTINE PLUG FOR BMD07D                  MAY 20, 1966
      SUBROUTINE PLUG (NV,NABC)
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3XMN(100),SX(10),DX(10),XL(10),SXX(100,10)
      COMMON  ZZ     , X
      COMMON  NVAR   , NCAS   , NTR    , CP     , IGR    , NCP
      COMMON  NCI    , KCICP  , NERR   , NG     , NADD
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500),SXX)
C
      DATA SPCVAL/4H-7CV/
       EXTERNAL SIGN
C
      ASSIGN 86 TO KSKIP
      JR=NS(NV)
      IF (NABC - 1) 10,10,15
  10  ASSIGN 88 TO KSKIP
      IF(-JR)15,998,998
   15 DO 50 I=1,10
      DO 50 J=1,30
      K000FX(I,J)=0
 50   CONTINUE
      KKP=NCICP(NV)+1
      KOUNT=0
      KOU=0
      DO 200 K=1,NG
      MN=IGR(K)
      IF(-MN)80,200,200
 80   DO 150 J=1,MN
      L=1
      I=NVAR*(KOUNT+J-1)+NV
      IF(-JR)85,100,100
  85  GO TO KSKIP,(86,88)
   86 IF(X(I).EQ.SPCVAL)GO TO 140
 100   IF ((CICP(L)-X(I))-.00001) 101,101,110
 101  L=L+1
      IF (L-KKP) 100, 100, 200
C
  88  DO 90 JS=1,JR
      IF(X(I).NE.SPV(NV,JS))GO TO 90
 89   B=SIGN(1.0, SPV(NV,JS))
      C=SIGN(1.0, X(I))
 897  IF(B+C) 899,90,899
  899 K000FX(K,JS) = K000FX(K,JS) + 1
      X(I)=SPCVAL
 90   CONTINUE
      GO TO 140
C
 110  LL=L-1
      K000FX(K,LL)=K000FX(K,LL)+1
 140  KOU=KOU+1
 150  CONTINUE
      KOUNT=KOU
 200  CONTINUE
 998  RETURN
      END
C             SUBROUTINE RDLBL FOR BMD07D              MAY 20, 1966
C     SUBROUTINE TO READ IN LABELS CARDS, STORE THEM IN ARRAY,
C     AND SUBSTITUTE ZEROS FOR UNLABELED VARIABLES
C     NVAR IS TOTAL NUMBER OF VARIABLES
C     NLBVAR IS NUMBER OF LABELED VARIABLES EXPECTED
C
      SUBROUTINE RDLBL(NLBVAR,NVAR,ARRAY)
      DOUBLE PRECISION DUMY,ARRAY
      DIMENSION ARRAY(1),IDUM(7),DUMY(7)
      DATA ALABEL/'LAB'/
C     CLEAR VARIABLES
      DO 1 I=1,NVAR
   1  ARRAY(I)=0.0
C     IF NO LABELS, RETURN
      IF(-NLBVAR)2,9,9
   2  N=0
C     READ 1 LABELS CARD
  20  READ (5,3) TEST,(IDUM(J),DUMY(J),J=1,7)
   3  FORMAT(A3,3X,7(I4,A6))
C     TEST FOR 'LAB' IN FIRST 3 COLS.
      IF(TEST.EQ.ALABEL)GO TO 6
C     ERROR--PRINT MESSAGE AND QUIT
   4  WRITE (6,5)
   5  FORMAT(36H0LABELS CARD NOT FOUND WHEN EXPECTED)
      STOP
C     EXAMINE 7 FIELDS
   6  DO 8 J=1,7
      K=IDUM(J)
C     TEST INDEX.  IF 0, IGNORE.  IF ILLEGAL, PRINT MESSAGE AND
C     IGNORE EXCEPT TO COUNT
      IF(-K)10,8,11
  10  IF(K-NVAR) 7,7,11
  11  WRITE (6,12)K,DUMY(J)
   12 FORMAT('0LABELS CARD INDEX',I7,' INCORRECT. LABEL ',A8,'IGNORED.')
C
      GO TO 13
C     MOVE LABEL TO ARRAY
   7  ARRAY(K)=DUMY(J)
C     STEP NUMBER OF VARIABLES
  13  N=N+1
C     TEST FOR END. IF END, RETURN. IF NOT, SCAN OTHER FIELDS.
      IF(N-NLBVAR) 8,9,9
   8  CONTINUE
      GO TO 20
   9  RETURN
      END
      SUBROUTINE RITE (L,NG,NABC)
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3 XMN(100), SX(10), DX(10), XL(10), SXX(100,10),OC(10)
      DATA OC/1H ,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/
       DATA ZERO / '0' /
      COMMON  ZZ
      COMMON  X      , NVAR   , NCAS   , NTR    , CP     , IGR
      COMMON  NCP    , NCI    , KCICP  , NERR   , QZ     , NADD
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
C
      J=L
      NABC = NABC
      DATA STAR  /1H*/
      DO 700 K=2,113
      VEC(K) = OC(1)
 700  CONTINUE
      DO 799 I=1,NG
      NN=K000FX(I,J)
      IF(-NN)711,799,799
 711  MA=((I-1)*11)+2
      MB=MA+(NN-1)
      MC=MA+8
      MD=MA+9
      ME=MA+10
      MF=MA+7
      MG=MA+6
      IF (NN-10) 712, 712, 714
 712  DO 713 K=MA,MB
      VEC(K)=STAR
 713  CONTINUE
      GO TO 799
C
 714  IF (NN-99) 715, 715, 717
 715  NNN=1
      DO 716 K=MA,MC
      VEC(K)=STAR
 716  CONTINUE
      NN1=NN/10
      NN2=NN-(10*NN1)
      GO TO 750
C
 717  IF (NN-999) 718, 718, 720
 718   NNN=2
      DO 719 K=MA,MF
      VEC(K)=STAR
 719  CONTINUE
      KNN=NN/10
      NN1=NN/100
      NN2=KNN-(10*NN1)
      NN3=NN-(10*KNN)
      GO TO 750
C
 720  IF (NN-9999) 721, 721, 799
 721  NNN=3
      DO 722 K=MA, MG
      VEC(K)=STAR
 722  CONTINUE
      NN1=NN/1000
      INN=NN/100
      JNN=NN/10
      NN2=INN-(10*NN1)
      NN3=JNN-(10*INN)
      NN4=NN-(10*JNN)
 750  IF (NNN-2) 751, 756, 771
 751  MZ=MD
      IJ=2
 7511 IF(-NN1)752,760,760
 752  GO TO (761, 762, 763, 764, 765, 766, 767, 768, 769), NN1
 753  IJ=IJ-1
      IF(-IJ)754,799,799
 754  MZ=ME
 7541 IF(-NN2)755,760,760
 755  GO TO (761, 762, 763, 764, 765, 766, 767, 768, 769), NN2
 756  IJ=2
      MZ=MC
      IF (NN1) 760, 760, 752
C
 757  IJ=IJ-1
      IF(-IJ)758,759,799
 758  MZ=MD
      IF (NN2) 760, 760, 755
C
 759  MZ=ME
 7591 IF(-NN3)770,760,760
 770  GO TO (761, 762, 763, 764, 765, 766, 767, 768, 769), NN3
C
 771  IJ=2
      IJK=2
      MZ=MF
      GO TO 7511
C
 772  IJ=IJ-1
      IF(-IJ)773,774,775
 773  MZ=MC
      GO TO 7541
C
 774  MZ=MD
      GO TO 7591
C
 775  IJK=IJK-1
      IF(-IJK)776,799,799
 776  MZ=ME
      IF(-NN4)777,760,760
 777  GO TO (761, 762, 763, 764, 765, 766, 767, 768, 769), NN4
C
 760   VEC(MZ) = ZERO
      GO TO (753, 757, 772), NNN
C
  761 VEC(MZ) = OC(2)
      GO TO (753, 757, 772), NNN
C
  762 VEC(MZ) = OC(3)
      GO TO (753, 757, 772), NNN
C
  763 VEC(MZ) = OC(4)
      GO TO (753, 757, 772), NNN
C
  764 VEC(MZ) = OC(5)
      GO TO (753, 757, 772), NNN
C
  765 VEC(MZ) = OC(6)
      GO TO (753, 757, 772), NNN
C
  766 VEC(MZ) = OC(7)
      GO TO (753, 757, 772), NNN
C
  767 VEC(MZ) = OC(8)
      GO TO (753, 757, 772), NNN
C
  768 VEC(MZ) = OC(9)
      GO TO (753, 757, 772), NNN
C
  769 VEC(MZ) = OC(10)
      GO TO (753, 757, 772), NNN
C
 799  CONTINUE
      HOLD = CICP(J)
       GO TO (4001,4002), NABC
 4001 CICP(J) =  SPV(NTR,J)
      GO TO 4004
C
 4002 IF(CICP(J))4004,4003,4004
 4003 CICP(J)=0.0
 4004 WRITE (6,4000)CICP(J), (VEC(MM), MM = 2, 113)
      CICP(J) = HOLD
      RETURN
 4000 FORMAT(1H ,F12.3,2H )113A1)
      END
C             SUBROUTINE  TEST  FOR BMD07D             MAY 20, 1966
      SUBROUTINE TEST(N1,N2,N3,C)
C
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3XMN(100),SX(10),DX(10),XL(10),SXX(100,10)
      COMMON  ZZ     , X
      COMMON  NVAR   , NCAS   , NTR    , CP     , IGR    , NCP
      COMMON  NCI    , KCICP  , NERR   , NG     , NADD
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
      EXTERNAL SIGN
      IF(N3)600,600,100
 100  DO 500 L=1,N3
      IF(X(N1).EQ.SPV(N2,L))GO TO 1000
 500  CONTINUE
 600  C=1.0
      RETURN
C
 1000 IF(X(N1))1050,1025,1050
 1025 B=SIGN(1.0,SPV(N2,L))
      C=SIGN(1.0,X(N1))
      IF(B+C)1050,500,1050
 1050 C=0.0
 2000 RETURN
      END
C        SUBROUTINE TPWD FOR BMD07D                    MAY 20, 1966
      SUBROUTINE TPWD(NT1,NT2)
      IF(NT1)40,10,12
 10   NT1=5
   12 IF(NT1 .EQ. NT2) GO TO 19
      IF(NT2.EQ.5)GO TO 18
   17 REWIND NT2
   19 IF(NT1 .EQ. 5) GO TO 24
   18 IF(NT1 .EQ. 6) GO TO 40
      REWIND NT1
 24   NT2=NT1
      RETURN
 40   WRITE (6,49)
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      STOP
      END
C        SUBROUTINE TVA FOR BMD07D                     MAY 20, 1966
      SUBROUTINE TVA
C
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3XMN(100),SX(10),DX(10),XL(10),SXX(100,10)
      COMMON  ZZ     , X
      COMMON  NVAR   , NCAS   , NTR    , CP     , IGR    , NCP
      COMMON  NCI    , KCICP  , NERR   , NG     , NADD
      ASN(XX) = ATAN(XX/SQRT(1.0-XX**2))
C
      DOUBLE PRECISION TRNGEN,COB
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
      DATA TRNGEN/6HTRNGEN/
C
      ON=NCAS+1
      NORM=0
      WRITE (6,10)
      WRITE (6,11)
      NERR=0
      DO 1000 I=1,NTR
      READ (5,12)COB,NE,NC,NV,CO
      IF(COB.NE.TRNGEN)GO TO 250
 45   WRITE (6,13)I, NE, NC, NV, CO
      IF(NE-NVAR)50,50,290
 50   NSN=NS(NV)
      NTIMES=1
      LL=CO
      IF(NV-NVAR)5005,5005,290
 5005 IF(10-NC)51,59,59
 51   IF(15-NC)52,59,55
 52   IF(16-NC)53,55,59
 53   IF(NC-23)59,55,275
 55   NSM=NS(LL)
      NTIMES=2
      IF(LL-NVAR)59,59,290
 59   IF(-NC)60,1000,275
   60 IF(-NERR)605,605,1000
  605 NVNC=NVAR*(NCAS-1)+NV
      KK=NE-NVAR
      KKK=LL-NVAR
      DO 910 K=NV,NVNC,NVAR
      D1=X(K)
      KK=KK+NVAR
      KKK=KKK+NVAR
      IF (NSN) 61, 61, 64
 61   GO TO (75,62),NTIMES
 62   IF(-NSM)63,75,75
 63   CALL TEST(KKK,LL,NSM,CHECK)
      IF(CHECK)75,70,75
 64   CALL TEST(K,NV,NSN,CHECK)
      IF(CHECK)61,70,61
 70   IF (NS(NE)) 71, 71, 72
 71   NS(NE)=1
      SPV(NE, 1) = -9909.9909
   72 D2=SPV(NE,1)
      GO TO 900
C
 75   GO TO (100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155,
     1160, 165, 170, 175, 180, 185, 190, 195, 200, 205, 210 ), NC
  100 IF(-D1)102,101,300
  102 D2=SQRT(D1)
      GO TO 900
  101 D2=0.0
      GO TO 900
  105 IF(-D1)107,106,300
  107 D2=SQRT(D1)+SQRT(D1+1.0)
      GO TO 900
  106 D2=1.0
      GO TO 900
  110 IF(-D1)111,300,300
  111 D2=ALOG10(D1)
      GO TO 900
  115 D2=EXP(D1)
      GO TO 900
  120 IF(-D1)121,101,300
  121 IF(D1-1.0)123,122,300
  123 D2=ASN(SQRT(D1))
      GO TO 900
  122 D2=1.57079632
      GO TO 900
  125 A=D1/ON
      B=A+1.0/ON
      IF(-A)128,126,300
 128  IF(-B)1290,129,300
 1290 A=SQRT(A)
      B=SQRT(B)
      D2=ASN(A)+ASN(B)
      GO TO 900
 126  IF(-B)127,101,300
  127 D2=ASN(SQRT(B))
      GO TO 900
  129 D2=ASN(SQRT(A))
      GO TO 900
  130 IF(D1)131,300,131
  131 D2=1.0/D1
      GO TO 900
  135 D2=D1+CO
      GO TO 900
  140 D2=D1*CO
      GO TO 900
 145  DO 146 JJJ=1,15
      AB=ABS(CO*2.0**JJJ)
      IF (AB-.5) 146, 147, 146
 146  CONTINUE
      GO TO 148
  147 IF(-D1)148,101,300
  148 D2=D1**CO
      GO TO 900
  150 D2=D1+X(KKK)
      GO TO 900
  155 D2=D1-X(KKK)
      GO TO 900
  160 D2=D1*X(KKK)
      GO TO 900
  165 IF(X(KKK))166,300,166
  166 D2=D1/X(KKK)
      GO TO 900
  170 IF(D1-CO)101,106,106
  175 IF(D1-X(KKK))101,106,106
  180 IF(-D1)181,300,300
  181 D2=ALOG(D1)
      GO TO 900
 185  GO TO(186,187),LTIMES
 186  LTIMES=2
      CALL YMEAN (NV, YMEA)
  187 D2=D1-YMEA
      GO TO 900
 190  GO TO (191,192),LTIMES
 191  LTIMES=2
      CALL DESQ(NV,DEVSQ)
  192 D2=D1/DEVSQ
      GO TO 900
  195 D2=SIN(D1)
      GO TO 900
  200 D2=COS(D1)
      GO TO 900
  205 D2=ATAN(D1)
      GO TO 900
 210  DO 211 JJJ=1,15
      AC=ABS(X(KKK)*2.0**JJJ)
      IF (AC-.5) 211, 212, 211
 211  CONTINUE
      GO TO 213
  212 IF(-D1)213,101,300
  213 IF(X(KKK))214,214,215
  214 IF(D1)215,300,215
  215 D2=D1**X(KKK)
  900 X(KK)=D2
  910 CONTINUE
      IF (NORM) 2000, 1000, 1000
C
 250  NERR=-999
      WRITE (6,17)
      GO TO 2000
C
  275 WRITE (6,18)I
  280 NERR=-999
      GO TO 1000
  290 WRITE (6,19)
      GO TO 280
C
 300  NORM=-999
      NERR=-999
      WRITE (6,14)I
      WRITE (6,15)K
 1000 CONTINUE
      IF (NERR) 2000, 3000, 3000
 2000 WRITE (6,16)
 3000 RETURN
C
 10   FORMAT (1H06X, 23HTRANS GENERATOR CARD(S))
 11   FORMAT (46H0CARD    NEW     TRANS    ORIG.   ORIG. VAR(B)/45H  NO.
     1 VARIABLE   CODE    VAR(A)   OR CONSTANT  )
  12  FORMAT(A6,I3,I2,I3,F6.0)
 13   FORMAT(I4,I8,2I9,4XF10.5)
 14   FORMAT (55H0THE INSTRUCTIONS INDICATED ON TRANS GENERATOR CARD NO.
     1 I2, 1X, 3HRE-/60H SULTED IN THE VIOLATION OF A RESTRICTION FOR TH
     2IS TRANSFOR-/58H MATION. THE VIOLATION OCCURED FOR THE ITEMS LISTE
     3D BELOW.)
 15   FORMAT (10H0ITEM NO. I5)
 16   FORMAT (41H0PROGRAM CANNOT CONTINUE FOR THIS PROBLEM )
 17   FORMAT(33H0TRANSGENERATOR CARD OUT OF ORDER)
 18   FORMAT(1H039X19HTRANSGENERATOR CARDI3,18H HAS ILLEGAL CODE.)
   19 FORMAT(1H0,26X,65HVARIABLE INDEX ON ABOVE TRANSGENERATION CARD IS   
     1GREATER THAN P+Q.)
      END
C             SUBROUTINE YMEAN FOR BMD07D              MAY 20, 1966
      SUBROUTINE YMEAN (M1,YMEA)
C
      DIMENSION X(4000),NS(100),SPV(100,5),CP(10),
     1 OCC(100), INTEG(10), ZZ(5050), IGR(10), BIGX(100), SMLX(100),
     2 CICP(31), K000FX(10,30), NCICP(100), VEC(113),IGP(10),NCN(100),
     3XMN(100),SX(10),DX(10),XL(10),SXX(100,10)
      COMMON  ZZ
      COMMON  X      , NVAR   , NCAS   , NTR    , CP     , IGR
      COMMON  NCP    , NCI    , KCICP  , NERR   , NG     , NADD
      EQUIVALENCE (ZZ(101),BIGX),   (ZZ(202),SMLX),               (ZZ(30
     13), NCICP), (ZZ(520), VEC), (ZZ(634), INTEG), (ZZ(755), K000
     2FX), (ZZ(2551), XMN), (ZZ(2751), SPV), (ZZ(2651), NS), (ZZ(1056),
     3NCN), (ZZ(1289), CICP), (ZZ(1320), SX), (ZZ(1331),DX), (ZZ(1342),
     4XL), (ZZ(1353), IGP), (ZZ(3500), SXX)
       EXTERNAL SIGN
C
      L=NS(M1)
      CAS=NCAS
      XNN=0
      DO 100 J=1,NCAS
      I=(J-1)*NVAR+M1
      IF (L) 51, 51, 45
 45   DO 50 K=1,L
      IF(X(I)-SPV(M1,K))50,55,50
 55    IF (SIGN(1.0,X(I))-SIGN(1.0,SPV(M1,K))) 50,57,50
 50   CONTINUE
   51 XNN=XNN+X(I)
      GO TO 100
 57    CAS = CAS-1.0
 100  CONTINUE
      YMEA=XNN/CAS
      RETURN
      END


















































