C        MULTIPLE REGRESSION WITH CASE COMBINATIONS    MAY  2, 1966
C        THIS PROGRAM IS A SIFTED VERSION OF BMD03R WITH SOME
C        MODIFICATIONS TO TAKE ADVANTAGE OF CERTAIN FEATURES OF FORTRAN
C        IV.
C        IT WAS THEN CONVERTED TO 360 FORTRAN IV (H-LEVEL)
C
      COMMON  BO     , SR     , SD     , ND1    , AA     , B
      COMMON  SB     , T      , ST     , AV     , ND     , P
      COMMON  DSR    , DSD    , NV     , NNEWA  , LLCODE , LLVA
      COMMON  RES    , DATA   , IND    , NOOB   , NDP    , FRMT(180)
      COMMON  NDEL   , OBNO   , COUNT  , JOB    , NOPR   , LN
      COMMON  CODET  , COR    , SSAR   , SSDR   , VARI   , SEE
      COMMON  N      , SS1    , F      , LLM    , SS2    , LLN
      COMMON  NSEL   , NVG    , ISAMP  , LCASE  , MERRY  , NADD
      COMMON  NGROUP , NTNOOB , KVR    , NTAPE  , NT2    , NT3
      DIMENSION DATA(50),SR(40,40),SD(40,40),ND1(28),AA(40,40),B(50),
     1SB(50),T(50),ST(50),AV(180),ND(28),P(40,40),DSR(40,40),DSD(40,40),
     2NV(50),NNEWA(50),LLCODE(50),LLVA(50),RES(50)
      DOUBLE PRECISION BO,SE1,SE2,SE3,SE4,SE5,SE6
      DATA SE5/'SELECT  '/,SE6/'FINISH  '/
 5010 FORMAT(55H1BMD03R - MULTIPLE REGRESSION WITH CASE COMBINATIONS -   
     1'REVISED JULY  7, 1969'/
     241H HEALTH SCIENCES COMPUTING FACILITY, UCLA,/
     330H-NO. OF SUBSAMPLES ...........,I2,/
     430H NO. OF ORIGINAL VARIABLES ...,I2,/
     527H TOT. COMB. SAMP. SIZE ....,I5,/
     630H ALTERNATE TAPE NO. ..........,I2,/
     730H NO. OF VAR. FORMAT CARDS ....,I2)
      READ(5,9)BO,NGROUP,IND,NTNOOB,MTAPE,KVR
 132  NTAPE = 5
      NT2=2
      NT3=3
      REWIND NT2
      REWIND NT3
C
      KVRH = KVR
      DATA SE1/6HPROBLM/
      IF(BO .EQ. SE1) GO TO 14
   13 WRITE (6,1013)
      WRITE(6,1017)
      GO TO 352
C
C
   14 IIND=IND
      IF((IND-1)*(IND-51)) 15,13,13
 15   IF(NTNOOB-99999)17,17,13
   17 IF(NGROUP*(NGROUP-29)) 18,13,13
 18   IF((MTAPE-NT2)*(MTAPE-NT3))20,13,20
 20   CALL TPWD(MTAPE,NTAPE)
      CALL DATAN
      READ (5,10)BO,JOB,NRE,NVG,NADD,NSEL,(ND(I),I=1,28)
      DATA SE2/6HSELECT/
      IF(BO .EQ. SE2) GO TO 16
 3900 WRITE (6,1014)
      WRITE(6,1017)
      GO TO 3915
 16   JOBD=JOB
      IF(NSEL*(NSEL-29)) 19,3900,3900
   19 IND=IIND
      NOPR=0
      NDEL=0
      ISAMP=0
      MERRY=0
      LCASE=0
      WRITE(6,5010) NGROUP,IND,NTNOOB,MTAPE,KVRH
      WRITE(6,451) (FRMT(I), I=1,KVR)
 451  FORMAT(' THE VARIABLE FORMAT IS  ',18A4/(1H ,24X,18A4))
      WRITE (6,499)JOB,NOPR
      WRITE (6,1008)(ND(I), I=1,NSEL)
      WRITE (6,510)
      IF(NVG)4000,4000,4001
 4001 IF(NVG-50)400,400,3900
 400  WRITE (6,5004)
      WRITE (6,5005)
      DATA SE3/6HTRNGEN/
      DO 4002 I=1,NVG
      READ (5,5006)BO,NNEWA(I),LLCODE(I),LLVA(I),RES(I)
      IF(BO .NE. SE3) GO TO 3910
  405 IF(LLCODE(I)*(18-LLCODE(I)))3910,3910,4002
 3910 WRITE(6,1015)
      WRITE(6,1017)
      IND=-IND
 4002 WRITE (6,5007)I,NNEWA(I),LLCODE(I),LLVA(I),RES(I)
      IF(IND)352,4000,4000
 4000 CALL READ
      IF(IND)3900,3900,410
 410  IF(NOOB-ISAMP)4003,4003,4004
 4004 OBNO=ISAMP
      NOOB=OBNO
 4003 IND=IND+NADD
      NDP=IND
      GO TO 4006
C
C
 3915 MERRY=-999
      GO TO 4007
 4006 WRITE (6,91)
      WRITE (6,92)(SB(I), I=1,IND)
      WRITE (6,510)
      WRITE (6,93)
      WRITE (6,92)(T(I), I=1,IND)
      WRITE (6,510)
      WRITE (6,1000)
      DO 2000 I=1,IND
      WRITE (6,173)I
 2000 WRITE (6,92)(DSR(I,J), J=1,IND)
      DO 112 I=1,IND
      AV(I)=SB(I)/OBNO
      DATA(I)=T(I)-((SB(I)*SB(I))/OBNO)
  112 ST(I)=SQRT(DATA(I)/(OBNO-1.0))
      DO171I=1,IND
      DO 171 J=1,IND
      DSD(I,J)=DSR(I,J)-((SB(I)*SB(J))/OBNO)
      DSD(I,J)=DSD(I,J)+1.E-25
  171 DSR(I,J)=DSD(I,J)/SQRT(DATA(I)*DATA(J))
      WRITE (6,510)
      WRITE (6,1011)
      DO 175 I=1,IND
      WRITE (6,173)I
  175 WRITE (6,92)(DSD(I,J), J=1,IND)
   56 LN=IND
      MM=LN-1
      DO 54 I=1,LN
   54 NV(I)=I
      IF(NDEL) 2987, 2987, 1118
 1118 DO 2985 I=1,NDEL
      KK=ND(I)
 2985 NV(KK)=0
      GO TO 2990
 2987 ND(1)=0
 2990 IF(NDP-LN) 2993, 3000, 2992
 2992 WRITE (6,1016)
      WRITE(6,1017)
      GO TO 4007
 2993 DO 2994 I=NDP,MM
 2994 NV(I)=NV(I+1)
      NV(LN)=NDP
 3000 IF(NDEL) 3050, 3050, 3030
 3030 LN=0
      DO 3045 I=1,IND
      IF(NV(I)) 3045, 3045, 3042
 3042 LN=LN+1
      NV(LN)=NV(I)
 3045 CONTINUE
 3050 DO 3060 J=1,LN
      MM=NV(J)
      DO 3060 I=1,LN
      KK=NV(I)
      SR(I,J)=DSR(KK,MM)
 3060 SD(I,J)=DSD(KK,MM)
   60 WRITE (6,500)JOB,NOPR
      WRITE (6,172)
      DO 181 I=1,LN
      WRITE (6,173)I
  181 WRITE (6,174)(SR(I,J), J=1,LN)
      N=LN-1
      DO 1060 J=1,LN
      DO 1060 I=1,LN
 1060 P(I,J)=SR(I,J)
      CALL INVERT (P,LN,DET,NNEWA,LLCODE)
      DO 1063 I=1,N
      TEMP = P(LN,LN) * P(I,I)
      IF(TEMP .LE. 0.0) GO TO 64
 1063 RES(I) = -P(I,LN) / SQRT(TEMP)
      IF(N-1) 65, 65, 67
 64   WRITE(6,641)
      GO TO 4007
   65 AA(1,1)=1.0
      GO TO 239
   67 DO 70 I=1,N
      DO 70 J=1,N
   70 AA(I,J)=SR(I,J)
      CALL INVERT(AA,N,DET,NNEWA,LLCODE)
      WRITE (6,510)
      WRITE (6,1009)
      DO 231 I=1,N
      WRITE (6,173)I
  231 WRITE (6,174)(AA(I,J),J=1,N)
  239 DO 240 I=1,N
      DO 240 J=1,N
  240 SR(I,J)=AA(I,J)*SR(I,J)/SD(I,J)
      DO 280 J=1,N
      SUM=0.0
      DO 270 I=1,N
  270 SUM=SUM+SR(I,J)*SD(I,LN)
  280 B(J)=SUM
      SSAR=0.0
      DO 290 I=1,N
  290 SSAR=SSAR+B(I)*SD(I,LN)
      SSDR=SD(LN,LN)-SSAR
      CODET=SSAR/SD(LN,LN)
      IF(CODET) 292,261,261
 292  WRITE(6,176)
      GO TO 4007
 261  COR = SQRT(CODET)
      DF=LN
      IF(OBNO-DF.EQ.1.0)SSDR=0.0
      VARI=SSDR/(OBNO-DF)
      SEE=SQRT(VARI)
      SUM=0.0
      DO 305 J=1,N
      SB(J)=SQRT(VARI*SR(J,J))
      K=NV(J)
      SUM=SUM+B(J)*AV(K)
      IF(SB(J).EQ.0.0)GO TO 303
      T(J)=B(J)/SB(J)
      GO TO 305
  303 T(J)=0.0
  305 CONTINUE
      K=NV(LN)
      BO=AV(K)-SUM
      D1=N
      D2=OBNO-DF
      D3=OBNO-1.0
      LLM=D2
      LLN=D3
      SS1=SSAR/D1
      SS2=SSDR/D2
      IF(SS2.EQ.0.0)GO TO 310
      F=SS1/SS2
      GO TO 320
  310 F=0.0
  320 NOOB=OBNO
      CALL REPORT
      IF(NRE) 4007, 4007, 348
  348 CALL RESIDU (NRE)
 4007 READ (5,10)BO,JOB,NOPR,NRE,NDP,NDEL,(ND(I),I=1,28)
      IF(JOB-JOBD)351,4008,351
 4008 IF(BO .NE. SE4) GO TO 351
      DATA SE4/6HREPDEL/
 4009 IF(MERRY)  4007,56,56
  351 MERRY=0
      IF(BO .EQ. SE4) GO TO  117
C
C
      IF(BO .EQ. SE5) GO TO 115
      IF(BO .EQ. SE1) GO TO 131
      IF(BO .NE. SE6) GO TO 119
      GO TO 352
  117 WRITE (6,1001)
  118 WRITE (6,1002)
      MERRY=-999
      GO TO 4007
  119 WRITE (6,1003)BO
      GO TO 118
  115 NVG=NRE
      NRE=NOPR
      NADD=NDP
      NSEL=NDEL
      CALL TPWD(MTAPE,NTAPE)
      GO TO 16
C
 352  IF(NTAPE-5)353,354,353
 353  REWIND NTAPE
 354  STOP
  131 NGROUP=JOB
      IND=NOPR
      NTN=NDEL/10
      NTNOOB=NRE*1000+NDP*10+NTN
      MTAPE=(NDEL-10*NTN)*10+ND(1)/10
	CALL USAGEB('BMD03R')
      KVR=ND(28)
      GO TO 132
    9 FORMAT(A6,2I2,I5,I2,53X,I2)
   10 FORMAT(A6,33I2)
   91 FORMAT(5H0SUMS)
   92 FORMAT (1H F19.5,5F20.5)
   93 FORMAT(15H0SUM OF SQUARES)
  172 FORMAT(25H0CORRELATION COEFFICIENTS)
  173 FORMAT(6H0ROW  I2)
 174  FORMAT(1H ,12G10.4)
 176  FORMAT('-THE COEFICIENT OF DETERMINATION IS NEGATIVE, BMD03R WILL   
     1SKIP TO THE NEXT REPDEL CARD,IF ANY.')
  499 FORMAT(14H0SELECTION NO.I5,1H-I2)
  500 FORMAT(14H1SELECTION NO.I5,1H-I2//)
  510 FORMAT(1H0)
 641  FORMAT('-A PARTIAL CORR. COEF. IS UNDEFINED, BMD03R WILL SKIP TO 
     1THE NEXT REPDEL CARD, IF ANY.')
 1000 FORMAT (19H0CROSS PRODUCT SUMS)
 1001 FORMAT(91H0ERRONEOUS REPDEL CARD ORDER. SELECTION NUMBER DOES NOT   
     1AGREE WITH PREVIOUS SELECTION CARD.)
 1002 FORMAT(46H PROGRAM SKIPS TO NEXT SELECTION CARD, IF ANY.)
 1003 FORMAT(73H0CONTROL CARD ORDER ERROR. PROGRAM EXPECTS A SELECT OR R
     1EPDEL CARD BUT A ,A6,18H CARD WAS PRESENT.)
 1008 FORMAT(21H SUB-SAMPLES SELECTED28I3)
 1009 FORMAT(42H0INVERSE OF CORRELATION COEFFICIENT MATRIX)
 1011 FORMAT(29H0CROSS PRODUCTS OF DEVIATIONS)
 1013 FORMAT(' A PROBLEM CARD WAS EXPECTED.')
 1014 FORMAT(' A SELECTION CARD WAS EXPECTED.')
 1015 FORMAT(' A TRANS-GENERATION CARD WAS EXPECTED.')
 1016 FORMAT(' A REPLACEMENT-DELETION CARD WAS EXPECTED.')
 1017 FORMAT(' THE CARD THAT WAS READ-IN CONTAINS AT LEAST ONE ERROR.')
 1018 FORMAT(20A4)
 1019 FORMAT(' THE LAST CARD READ-IN WAS     ',20A4)
 5004 FORMAT(1H06X,23HTRANS GENERATOR CARD(S))
 5005 FORMAT(46H0CARD    NEW     TRANS    ORIG.   ORIG. VAR(B)/45H  NO.   
     1VARIABLE   CODE    VAR(A)   OR CONSTANT)
 5006 FORMAT(A6,I3,I2,I3,F6.0)
 5007 FORMAT(2H  I2,I8,2I9,F14.4)
 5008 FORMAT(22H0SAMPLE SIZE RETAINED=I7)
 5009 FORMAT(67H0INSUFFICIENT SAMPLE,PROGRAM WILL GO TO NEXT SELECTION,E
     1TC., IF ANY)
C
      END
      SUBROUTINE DATAN
C        SUBROUTINE DATAN FOR BMD03R                   MAY  2, 1966
      COMMON  BO     , SR     , SD     , ND1    , AA     , B
      COMMON  SB     , T      , ST     , AV     , ND     , P
      COMMON  DSR    , DSD    , NV     , NNEWA  , LLCODE , LLVA
      COMMON  RES    , DATA   , IND    , NOOB   , NDP    , FRMT(180)
      COMMON  NDEL   , OBNO   , COUNT  , JOB    , NOPR   , LN
      COMMON  CODET  , COR    , SSAR   , SSDR   , VARI   , SEE
      COMMON  N      , SS1    , F      , LLM    , SS2    , LLN
      COMMON  NSEL   , NVG    , ISAMP  , LCASE  , MERRY  , NADD
      COMMON  NGROUP , NTNOOB , KVR    , NTAPE  , NT2    , NT3
      DIMENSION DATA(50),SR(40,40),SD(40,40),ND1(28),AA(40,40),B(50),
     1SB(50),T(50),ST(50),AV(180),ND(28),P(40,40),DSR(40,40),DSD(40,40),
     2NV(50),NNEWA(50),LLCODE(50),LLVA(50),RES(50)
      DOUBLE PRECISION TOTAL,SUM,SE1,BO,SIZE
      DATA TOTAL/6HTOTAL /
      DATA SUM/6H SUM  /
      DATA SE1/6HSAMSIZ/
C
      READ (5,900)BO,(ND1(I), I=1,NGROUP)
      IF(BO .EQ. SE1) GO TO 15
   10 WRITE (6,902)
      WRITE(6,1017)
C
C
C
      STOP
   15 IF(KVR.GT.0.AND.KVR.LE.10)GO TO 17
      WRITE (6,4000)
      KVR=1
   17 KVR=KVR*18
      NTOTAL=0
      DO 20 I=1,NGROUP
   20 NTOTAL=NTOTAL+ND1(I)
      SIZE=TOTAL
      IF(NTNOOB-NTOTAL)24,45,25
   24 SIZE=SUM
      NTNOOB=NTOTAL
   25 WRITE (6,903)SIZE
   45 READ (5,901) (FRMT(I), I=1,KVR)
      DO 452 I1=1,KVR
 452  AV(I1) = FRMT(I1)
      DO 55 J=1,NTNOOB
   49 READ (NTAPE,AV)(DATA(I), I=1,IND)
   55 WRITE(NT3)(DATA(I),I=1,IND)
 1000 FORMAT(20A4)
   60 END FILE NT3
      REWIND NT3
      RETURN
  900 FORMAT(A6,11I6/(6X,11I6))
 901  FORMAT(18A4)
 902  FORMAT(' A SAMPLE SIZE CARD WAS EXPECTED.')
  903 FORMAT(81H0TOTAL SAMPLE SIZE AND SUM OF SAMPLE SIZES DO NOT AGREE.
     1 THE PROGRAM ASSUMES THE A6,28H TO BE CORRECT AND PROCEEDS.)
 1017 FORMAT(' THE CARD THAT WAS READ-IN CONTAINS AT LEAST ONE ERROR.')
C
C
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
      END
      SUBROUTINE INVERT (A,N,D,L,M)
C        SUBROUTINE INVERT FOR BMD03R                  MAY  2, 1966
C     PROGRAM FOR FINDING THE INVERSE OF A NXN MATRIX
      DIMENSION A(40,40),L(50),M(50)
C     SEARCH FOR LARGEST ELEMENT
C
      D=1.0
      DO80 K=1,N
      L(K)=K
      M(K)=K
      BIGA=A(K,K)
      DO20 I=K,N
      DO20 J=K,N
      IF(ABS(BIGA)-ABS(A(I,J))) 10,20,20
   10 BIGA=A(I,J)
      L(K)=I
      M(K)=J
   20 CONTINUE
C     INTERCHANGE ROWS
      J=L(K)
      IF(L(K)-K) 35,35,25
   25 DO30 I=1,N
      HOLD=-A(K,I)
      A(K,I)=A(J,I)
   30 A(J,I)=HOLD
C     INTERCHANGE COLUMNS
   35 I=M(K)
      IF(M(K)-K) 45,45,37
   37 DO40 J=1,N
      HOLD=-A(J,K)
      A(J,K)=A(J,I)
   40 A(J,I)=HOLD
C     DIVIDE COLUMN BY MINUS PIVOT
   45 DO55 I=1,N
   46 IF(I-K)50,55,50
   50 A(I,K)=A(I,K)/(-A(K,K))
   55 CONTINUE
C     REDUCE MATRIX
      DO65 I=1,N
      DO65 J=1,N
   56 IF(I-K) 57,65,57
   57 IF(J-K) 60,65,60
   60 A(I,J)=A(I,K)*A(K,J)+A(I,J)
   65 CONTINUE
C     DIVIDE ROW BY PIVOT
      DO75 J=1,N
   68 IF(J-K)70,75,70
   70 A(K,J)=A(K,J)/A(K,K)
   75 CONTINUE
C     CONTINUED PRODUCT OF PIVOTS
      D=D*A(K,K)
C     REPLACE PIVOT BY RECIPROCAL
      A(K,K)=1.0/A(K,K)
   80 CONTINUE
C     FINAL ROW AND COLUMN INTERCHANGE
      K=N
  100 K=(K-1)
      IF(K) 150,150,103
  103 I=L(K)
      IF(I-K) 120,120,105
  105 DO110 J=1,N
      HOLD=A(J,K)
      A(J,K)=-A(J,I)
  110 A(J,I)=HOLD
  120 J=M(K)
      IF(J-K) 100,100,125
  125 DO130 I=1,N
      HOLD=A(K,I)
      A(K,I)=-A(J,I)
  130 A(J,I)=HOLD
      GO TO 100
  150 RETURN
      END
      SUBROUTINE READ
C        SUBROUTINE READ FOR BMD03R                    MAY  2, 1966
      COMMON  BO     , SR     , SD     , ND1    , AA     , B
      COMMON  SB     , T      , ST     , AV     , ND     , P
      COMMON  DSR    , DSD    , NV     , NNEWA  , LLCODE , LLVA
      COMMON  RES    , DATA   , IND    , NOOB   , NDP    , FRMT(180)
      COMMON  NDEL   , OBNO   , COUNT  , JOB    , NOPR   , LN
      COMMON  CODET  , COR    , SSAR   , SSDR   , VARI   , SEE
      COMMON  N      , SS1    , F      , LLM    , SS2    , LLN
      COMMON  NSEL   , NVG    , ISAMP  , LCASE  , MERRY  , NADD
      COMMON  NGROUP , NTNOOB , KVR    , NTAPE  , NT2    , NT3
      DIMENSION DATA(50),SR(40,40),SD(40,40),ND1(28),AA(40,40),B(50),
     1SB(50),T(50),ST(50),AV(180),ND(28),P(40,40),DSR(40,40),DSD(40,40),
     2NV(50),NNEWA(50),LLCODE(50),LLVA(50),RES(50)
      DOUBLE PRECISION BO
C
      LLM=IND+NADD
      IF((LLM-1)*(LLM-51)) 1,100,100
 1    DO 20 I=1,50
      SB(I)=0.0
      T(I)=0.0
      DO 20 J=1,50
   20 DSR(I,J)=0.0
      DO 21 I=1,NGROUP
   21 NV(I)=-ND1(I)
      DO 25 I=1,NSEL
      NQ=ND(I)
      IF(NQ)25,25,22
 22   NV(NQ)=NV(NQ)*(-1)
 25   CONTINUE
      NTAB=0
      NOOB=0
      DO 6001 I8=1,NGROUP
      I9=NV(I8)
      IF(I9) 5871, 6001, 5901
 5871 I9=I9*(-1)
      DO 5875 M9=1,I9
 5875 READ(NT3)(DATA(I),I=1,IND)
 1000 FORMAT(20A4)
      GO TO 6001
 5901 DO 6000 M9=1,I9
      READ(NT3)(DATA(I),I=1,IND)
      IF(NVG) 66, 66, 201
  201 CALL TRNGEN(DATA,IND,NVG,I9,ISAMP,LCASE,NNEWA,LLCODE,LLVA,RES,
     1MERRY,M9)
      IF(LCASE) 203, 66, 66
  203 LCASE=0
      GO TO 6000
   66 DO 67 I=1,LLM
      SB(I)=SB(I)+DATA(I)
   67 T(I)=T(I)+DATA(I)*DATA(I)
      DO 151 I=1,LLM
      DO 151 J=1,LLM
  151 DSR(I,J)=DSR(I,J)+DATA(I)*DATA(J)
      WRITE(NT2)(DATA(I),I=1,LLM)
 6000 CONTINUE
      NOOB=NOOB+I9
      NTAB=NTAB+1
      IF(NSEL-NTAB) 81, 81, 6001
 6001 CONTINUE
   81 END FILE NT2
      REWIND NT2
      REWIND NT3
      OBNO=NOOB
      ISAMP=NOOB-ISAMP
 110  RETURN
 100  IND=-LLM
      GO TO 110
C
      END
      SUBROUTINE REPORT
C        SUBROUTINE REPORT FOR BMD03R                  MAY  2, 1966
      COMMON  BO     , SR     , SD     , ND1    , AA     , B
      COMMON  SB     , T      , ST     , AV     , ND     , P
      COMMON  DSR    , DSD    , NV     , NNEWA  , LLCODE , LLVA
      COMMON  RES    , DATA   , IND    , NOOB   , NDP    , FRMT(180)
      COMMON  NDEL   , OBNO   , COUNT  , JOB    , NOPR   , LN
      COMMON  CODET  , COR    , SSAR   , SSDR   , VARI   , SEE
      COMMON  N      , SS1    , F      , LLM    , SS2    , LLN
      COMMON  NSEL   , NVG    , ISAMP  , LCASE  , MERRY  , NADD
      COMMON  NGROUP , NTNOOB , KVR    , NTAPE
      DIMENSION DATA(50),SR(40,40),SD(40,40),ND1(28),AA(40,40),B(50),
     1SB(50),T(50),ST(50),AV(180),ND(28),P(40,40),DSR(40,40),DSD(40,40),
     2NV(50),NNEWA(50),LLCODE(50),LLVA(50),RES(50)
      DOUBLE PRECISION BO
C
      WRITE (6,500)JOB, NOPR,NOOB
      WRITE (6,502)LN,NDEL
      WRITE (6,1005)NDP
      WRITE (6,503)CODET
      WRITE (6,504)COR
      WRITE (6,505)SSAR
      WRITE (6,506)SSDR
      WRITE (6,507)VARI
      WRITE (6,508)SEE
      WRITE (6,509)BO
      WRITE (6,511)
      WRITE (6,515)N,SSAR,SS1,F
      WRITE (6,516)LLM,SSDR,SS2
      WRITE (6,517)LLN,SD(LN,LN)
      DO 10 J=1,LN
 10   SR(1,J)=SD(1,J)
      DO 30 I=1,N
      II=I
      DO 20 J=II,LN
 20   P(I,J)=SR(I,J)/SR(I,I)
      K=I+1
      DO 30 J=K,LN
      SUM=0.0
      DO 25 M=1,II
 25   SUM=SUM+P(M,K)*SR(M,J)
 30   SR(K,J)=SD(K,J)-SUM
      DO 40 I=1,N
 40   AA(I,1)=SR(I,LN)*P(I,LN)
      DO 45 I=1,N
 45   AA(I,2)=AA(I,1)/SD(LN,LN)
      WRITE (6,518)
      WRITE (6,519)
      DO350 I=1,N
      K=NV(I)
 350  WRITE (6,520)NV(I),AV(K),ST(K),B(I),SB(I),T(I),RES(I), AA(I,1),AA(
     1I,2)
      K=NV(LN)
      WRITE (6,521)NV(LN),AV(K),ST(K)
      WRITE (6,522)P(N,LN)
      WRITE (6,1013)(ND(I), I=1,NDEL)
  500 FORMAT(14H1SELECTION NO.I5,1H-I2//12H SAMPLE SIZEI8)
  502 FORMAT(17H NO. OF VARIABLESI3,30H      NO. OF VARIABLES DELETEDI3,
     136H  (FOR VARIABLES DELETED, SEE BELOW))
  503 FORMAT(29H0COEFFICIENT OF DETERMINATIONF8.4)
  504 FORMAT(29H MULTIPLE CORR. COEFFICIENT  F8.4)
  505 FORMAT(44H0SUM OF SQUARES ATTRIBUTABLE TO REGRESSION  F15.5)
  506 FORMAT(44H SUM OF SQUARES OF DEVIATION FROM REGRESSIONF15.5)
  507 FORMAT(23H0VARIANCE OF ESTIMATE  F12.5)
  508 FORMAT(23H STD. ERROR OF ESTIMATEF12.5)
  509 FORMAT(20H0INTERCEPT (A VALUE)F12.5//)
  511 FORMAT(1H0,16X,37HANALYSIS OF VARIANCE FOR THE MULTIPLE/25X,18HLIN
     XEAR  REGRESSION/7X,19HSOURCE OF VARIATION,7X,4HD.F.,6X,6HSUM OF,8X
     X,4HMEAN,9X,1HF/42X,7HSQUARES,7X,7HSQUARES,5X,5HVALUE)
  515 FORMAT(30H DUE TO REGRESSION............I5,F15.5,F14.5,F11.4)
  516 FORMAT(30H DEVIATION ABOUT REGRESSION...I5,F15.5,F14.5)
  517 FORMAT(22X,8HTOTAL...I5,F15.5//)
  518 FORMAT(105H0VARIABLE     MEAN        STD.        REG.    STD.ERROR
     1    COMPUTED    PARTIAL     SUM OF SQ.  PROP. VAR.)
  519 FORMAT(101H   NO.                  DEVIATION    COEFF.   OF REG.CO
     1E.  T VALUE    CORR. COE.     ADDED       CUM.)
  520 FORMAT(I5,F15.5,F12.5,2F11.5,2F12.5,F15.5,F11.5)
  521 FORMAT(I5,F15.5,F12.5)
  522 FORMAT(32H0COMP. CHECK ON FINAL COEFF.    F11.5//)
 1005 FORMAT(30H DEPENDENT VARIABLE IS NOW NO.I3)
 1013 FORMAT(21H VARIABLES DELETED...28I3)
      RETURN
      END
      SUBROUTINE RESIDU (NRE)
C        SUBROUTINE RESIDU FOR BMD03R                  MAY  2, 1966
      COMMON  BO     , SR     , SD     , ND1    , AA     , B
      COMMON  SB     , T      , ST     , AV     , ND     , P
      COMMON  DSR    , DSD    , NV     , NNEWA  , LLCODE , LLVA
      COMMON  RES    , DATA   , IND    , NOOB   , NDP    , FRMT(180)
      COMMON  NDEL   , OBNO   , COUNT  , JOB    , NOPR   , LN
      COMMON  CODET  , COR    , SSAR   , SSDR   , VARI   , SEE
      COMMON  N      , SS1    , F      , LLM    , SS2    , LLN
      COMMON  NSEL   , NVG    , ISAMP  , LCASE  , MERRY  , NADD
      COMMON  NGROUP , NTNOOB , KVR    , NTAPE
      DIMENSION DATA(50),SR(40,40),SD(40,40),ND1(28),AA(40,40),B(50),
     1SB(50),T(50),ST(50),AV(180),ND(28),P(40,40),DSR(40,40),DSD(40,40),
     2NV(50),NNEWA(50),LLCODE(50),LLVA(50),RES(50)
      DOUBLE PRECISION BO
      SMALL=10.0**36
      BIG=-SMALL
 6    IF(NRE-1) 9,9,8
 8    WRITE (6,700)
      WRITE (6,701)
 9    DO 20 K=1,NOOB
      READ(2)(DATA(I),I=1,IND)
 1000 FORMAT(20A4)
      YEST=BO
      DO 10 I=1,N
      J=NV(I)
 10   YEST=YEST+B(I)*DATA(J)
      J=NV(LN)
      RESID=DATA(J)-YEST
      IF(RESID-SMALL) 11, 12, 12
 11   SMALL=RESID
 12   IF(RESID-BIG) 14, 14, 13
 13   BIG=RESID
   14 IF(NRE.GT.1)WRITE(6,702)K,DATA(J),YEST,RESID
 20   CONTINUE
      REWIND 2
      WRITE (6,703)
      WRITE (6,709)
 21   W=BIG-SMALL
      CRIT1=W/SEE
      WRITE (6,704)W
      WRITE (6,705)CRIT1
  700 FORMAT(1H1,19X,18HTABLE OF RESIDUALS)
 701  FORMAT(59H0OBSERVATION     Y VALUE         Y ESTIMATE        RESID
     1UAL)
 702  FORMAT(I7,2F18.5,F16.5)
 703  FORMAT(1H0)
 704  FORMAT(35H0RANGE OF RESIDUALS................F16.3)
 705  FORMAT(35H RANGE / STD. ERROR OF ESTIMATE....F16.3)
 709  FORMAT(26H0TEST OF EXTREME RESIDUALS)
 100  RETURN
      END
      SUBROUTINE TPWD(NT1,NT2)
C        SUBROUTINE TPWD FOR BMD03R                    MAY  2, 1966
C
      IF(NT1)40,10,12
 10   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   NT2=NT1
 28   RETURN
 40   WRITE (6,49)
      STOP
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      END
      SUBROUTINE TRNGEN(DATA,IND,NVG,NODATA,ISAMP,LCASE,NNEWA,LLCODE,
     1LLVA,RES,MERRY,N)
C        SUBROUTINE TRNGEN FOR BMD03R                  MAY  2, 1966
      DIMENSION DATA(50),NNEWA(50),LLCODE(50),LLVA(50),RES(50)
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
C
      SAMP=NODATA
      DO 3 J=1,NVG
  305 NEWA=NNEWA(J)
      LCODE=LLCODE(J)
      LVA=LLVA(J)
      BNEW=RES(J)
      IF(LCODE-10) 4,4,5
    5 NEWB=BNEW
    4 D=DATA(LVA)
      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, 7,17
   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, 7,27
   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, 7,33
   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)
      GO TO 3
  150 BNEW=NEWB
      IF(D-BNEW) 7,11,11
  160 IF(D-DATA(NEWB)) 7,11,11
  170 IF(D.LE.0.0)GO TO 99
      DATA(NEWA)=ALOG(D)
    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)N
      WRITE (6,1408)
      ISAMP=ISAMP+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
















































