C             CORRELATION WITH ITEM DELETION         MARCH  1, 1966
C        THIS IS A SIFTED VERSION OF BMD03D ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
C               HEALTH SCIENCES COMPUTING FACILITY
C                      UCLA MEDICAL SCHOOL
      DIMENSION NOC(90),NVAR(90),CODE(90,10),FMT(180),
     1SUMX(2000),SUMY(2000),SUMX2(2000),SUMY2(2000),SUMXY(2000),
     2ITEM(2000),X(90),LL(90),JUNK(90),          NEW(99),LUMP(99),NA(99)
     3,BN(99),SUM(90),SUM2(90),COUNT(90)
      COMMON  FMT
      COMMON  SUMX   , SUMY   , SUMXY  , ITEM   , NVNEW  , JUNK
      COMMON  NV     , STAR
C
      EQUIVALENCE (FMT,LL)
      DOUBLE PRECISION TODE,PROB,TRAN,A123,B123,C123,D123,E123,F123
C
  308 FORMAT('1BMD03D - CORRELATION WITH ITEM DELETION - REVISED ',
     1'JUNE 26, 1969'
     2/41H HEALTH SCIENCES COMPUTING FACILITY, UCLA //
     313H PROBLEM CODE21(1H.),A6/20H NUMBER OF VARIABLES14(1H.),I6/16H N
     4UMBER OF CASES18(1H.),I6/34H NUMBER OF TRANSGENERATION CARDS..I6/
     526H NUMBER OF VARIABLES ADDED 8(1H.),I6/34H NUMBER OF VARIABLE FOR
     6MAT CARD(S)I6//33H TRANSGENERATION (IF ANY) OCCURS A6,14H ITEM DEL
     7ETION///)
C
      DATA A123,B123,C123,D123,E123,F123,STRR/'PROBLM','FINISH','DELETE'
     X,'TRNGEN','AFTER ','BEFORE','$'/
      DATA BLANK,EMPTY,TEN/'    ','    ',10.0/
      DATA RE / 'NO' /
      STAR=STRR
      NTAPE=5
	CALL USAGEB('BMD03D')
   18 READ(5,10) TODE,PROB,NV,N,MISC,NTG,NADD,MMQ,REW,
     * TRAN,MTAPE,KVR
      IF(TODE.EQ.B123)GO TO 301
      IF(TODE.EQ.A123)GO TO 20
      PRINT 1000, TODE
  302 WRITE (6,304)
  301 IF(NTAPE.EQ.5)GO TO 310
      REWIND NTAPE
  310 STOP
 1001 PRINT 1002
      GO TO 302
 1003 PRINT 1004
      GO TO 302
 1007 PRINT 1005, TODE
      GO TO 302
 1006 PRINT 1008, I
      GO TO 302
 1009 PRINT 1010
      GO TO 302
 1011 PRINT 1012
      GO TO 302
 1013 PRINT 1014
      GO TO 302
 1015 PRINT 1016, TODE
      GO TO 302
   20 IF(REW.NE.RE.OR.MTAPE.NE.NTAPE)CALL TPWD(MTAPE,NTAPE)
   27 IF(NV*(91-NV))1001,1001,307
  307 JUMP=3
      ITRAN=1
      IF(TRAN.EQ.E123)GO TO 3077
      TRAN=F123
      ITRAN=2
 3077 WRITE (6,308)PROB,NV,N,NTG,NADD,KVR, TRAN
      IF(MISC)1003,65,500
  500 IF(MISC-90)25,25,1003
   25 DO 60 I=1,MISC
      READ (5,12)TODE,NVAR(I),NOC(I),(CODE(I,J),J=1,10)
      IF(TODE.NE.C123)GO TO 1007
  306 IF(NVAR(I))1006,30,35
   30 JUMP=2
   35 KK=NOC(I)
      IF(KK)1009,60,36
   36 IF(KK-10)37,37,1009
 37   DO 50 J=1,KK
      IF(CODE(I,J))50,40,50
   40 D=CODE(I,J)
      TEST=SIGNE(TEN,D)
      IF(TEST)45,50,50
   45 CODE(I,J)=BLANK
   50 CONTINUE
   60 CONTINUE
      GO TO 100
   65 JUMP=1
 100  NVNEW=NV+NADD
      IF(NVNEW-90)101,101,1011
 101  DO 400 I=1,NVNEW
  400 JUNK(I)=0
C
C     READ IN TRANSGENERATION CARDS
C
      IF(NTG*(NTG-100))275,277,1013
  275 WRITE (6,104)
      WRITE (6,108)
      NN=NTG
      DO 411 I=1,NTG
      READ (5,106)TODE,NEW(I),LUMP(I),NA(I),BN(I)
      IF(TODE.NE.D123)GO TO 1015
 28   WRITE (6,107)I,NEW(I),LUMP(I),NA(I),BN(I)
      IF(LUMP(I).GT.0.AND.LUMP(I).LE.10)GO TO 411
 413  WRITE (6,903)I
      NN=-NTG
  411 CONTINUE
      IF(NN)301,412,412
  277 ITRAN=3
  412 IF(KVR.GT.0.AND.KVR.LE.10)GO TO 105
      WRITE(6,4000)
      KVR=1
  105 KVR=KVR*18
      READ (5,16)(FMT(I),I=1,KVR)
      PRINT 1017, (FMT(I), I=1,KVR)
      NN=(NVNEW*NVNEW-NVNEW)/2
      DO 120 I=1,NN
      SUMX(I)=0.0
      SUMX2(I)=0.0
      SUMY(I)=0.0
      SUMY2(I)=0.0
      SUMXY(I)=0.0
  120 ITEM(I)=0
      DO 119 I=1,NVNEW
      SUM(I)=0.0
      SUM2(I)=0.0
  119 COUNT(I)=0.0
C
C        THE CODING USING THE MOD FUNCTION IS DONE TO ALLOW THE TOTAL
C        NUMBER OF CASES TO BE GREATER THAN 2**15 -1 (32767).
C
      N1=N
 1195 N2=MOD(N1,32767)
      IF(N2.EQ.0)N2=32767
      DO 215 III=1,N2
      READ (NTAPE,FMT)(X(I),I=1,NV)
      GO TO(122,121,122),ITRAN
  121 CALL TRANS(X,NEW,LUMP,NA,BN,N,III,NTG,ITRAN)
 122  GO TO (195,165,125),JUMP
  125 DO 160 K=1,MISC
      KK=NVAR(K)
      NO=NOC(K)
      IF(NO.LE.0)GO TO 160
 126  DO 155 J=1,NO
      IF(CODE(K,J).EQ.BLANK)GO TO 135
      IF(CODE(K,J).EQ.X(KK))GO TO 150
      GO TO 155
  135 IF(X(KK))155,140,155
  140 TEST=SIGNE(TEN,X(KK))
      IF(TEST)150,155,155
  150 X(KK)=BLANK
  155 CONTINUE
  160 CONTINUE
      GO TO 195
  165 NO=NOC(1)
      IF(NO.LE.0)GO TO 195
 166  DO 191 J=1,NVNEW
      DO 190 K=1,NO
      IF(CODE(1,K).EQ.BLANK)GO TO 175
      IF(CODE(1,K).EQ.X(J))GO TO 185
      GO TO 190
  175 IF(X(J))190,180,190
  180 TEST=SIGNE(TEN,X(J))
      IF(TEST)185,190,190
  185 X(J)=BLANK
      GO TO 191
  190 CONTINUE
  191 CONTINUE
  195 GO TO(1955,1958,1958),ITRAN
 1955 CALL TRANS(X,NEW,LUMP,NA,BN,N,III,NTG,ITRAN)
 1958 DO 197 JK=1,NVNEW
      IF(X(JK).EQ.BLANK)GO TO 197
  196 JUNK(JK)=JUNK(JK)+1
  197 CONTINUE
      K=0
      NV1=NVNEW-1
      DO 213 I=1,NV1
      IF(X(I).EQ.BLANK)GO TO 212
  200 SUM(I)=SUM(I)+X(I)
      SUM2(I)=SUM2(I)+(X(I)*X(I))
      COUNT(I)=COUNT(I)+1.0
      IX=I+1
      DO 210 J=IX,NVNEW
      K=K+1
      IF(X(J).EQ.BLANK)GO TO 210
  205 SUMX(K)=SUMX(K)+X(I)
      SUMX2(K)=SUMX2(K)+X(I)*X(I)
      SUMY(K)=SUMY(K)+X(J)
      SUMY2(K)=SUMY2(K)+X(J)*X(J)
      SUMXY(K)=SUMXY(K)+X(I)*X(J)
      ITEM(K)=ITEM(K)+1
  210 CONTINUE
      GO TO 213
  212 K=K+NVNEW-I
  213 CONTINUE
 2105 IF(X(NVNEW).EQ.BLANK)GO TO 215
  211 SUM(NVNEW)=SUM(NVNEW)+X(NVNEW)
      SUM2(NVNEW)=SUM2(NVNEW)+(X(NVNEW)*X(NVNEW))
      COUNT(NVNEW)=COUNT(NVNEW)+1.0
  215 CONTINUE
      N1=N1-N2
      IF(N1.GT.0)GO TO 1195
C     DETERMINE MEANS AND STANDARD DEVIATIONS
      WRITE (6,904)
      DO 2151 I=1,NVNEW
      Y=0.0
      IF(COUNT(I).NE.0.0)Y=SUM(I)/COUNT(I)
      K=COUNT(I)
      B=0.0
      IF(COUNT(I).NE.1.0)B=SQRT((SUM2(I)-(Y*SUM(I)))/(COUNT(I)-1.0))
 2151 WRITE (6,905)I,Y,B,K
      DO 220 K=1,NN
      FNTOT=ITEM(K)
      TOP=FNTOT*SUMXY(K)-SUMX(K)*SUMY(K)
      BO=(FNTOT*SUMX2(K)-SUMX(K)**2)*(FNTOT*SUMY2(K)-SUMY(K)**2)
      IF(BO)217,217,218
 217  SUMXY(K)=STAR
      BO=0.0
      SUMX(K)=0.0
      GO TO 220
 218  SUMXY(K)=EMPTY
      BO=SQRT(BO)
      SUMX(K)=TOP/BO
 220  CONTINUE
      DO 230 I=1,NVNEW
  230 LL(I)=I
      CALL PRINT(1,NVNEW,MMQ,1)
      WRITE (6,900)
      WRITE (6,901)
      MAX=10
      NF=1
      IF(NVNEW-MAX)235,235,240
  235 NL=NVNEW
      CALL PRINT(NF,NL,0,0)
      GO TO 18
  240 NL=MAX
      CALL PRINT(NF,NL,0,0)
      NO=NVNEW
  243 NO=NO-MAX
      NF=NF+MAX
      WRITE (6,900)
      WRITE (6,902)
      IF(NO-MAX)250,250,245
  245 NL=NL+MAX
      CALL PRINT(NF,NL,0,0)
      GO TO 243
  250 NL=NL+NO
      CALL PRINT(NF,NL,0,0)
  270 GO TO 18
C
   10 FORMAT(2A6,I2,I5,2I2,I3,I2,33X,A2,A5,2I2)
   12 FORMAT(A6,2I2,10F6.0)
   16 FORMAT(18A4)
  104 FORMAT(23H TRANS-GENERATION CARDS//)
  106 FORMAT(A6,I3,I2,I3,F6.0)
  107 FORMAT(I4,I5,I6,I5,4X,F9.4)
  108 FORMAT(30H CARD NEW  TRANS OLD    CONST./23H  NO. VAR.  CODE VAR(A
     1)//)
  304 FORMAT(45H0CONTROL CARDS INCORRECTLY ORDERED OR PUNCHED)
  900 FORMAT(1H047X,18HCORRELATION MATRIX/43X29H(SAMPLE SIZES IN PARENTH
     1ESES))
  901 FORMAT(1H08X,12HVARIABLE NO.)
  902 FORMAT(1H08X,24HVARIABLE NO. (CONTINUED))
 903  FORMAT(1H026X19HTRANSGENERATOR CARDI4,42H HAS ILLEGAL CODE. PROGRA
     1M CANNOT PROCEED.)
  904 FORMAT(45X,30HMEANS AND STANDARD DEVIATIONS //5X,8HVARIABLE,6X,4HM
     1EAN,2X,8HSTANDARD3X9HNUMBER OF/25X,9HDEVIATION,4X,5HITEMS,//)
  905 FORMAT(6X,I4,4X,2F10.4,3X,I6)
 1000 FORMAT(' PROGRAM EXPECTED PROBLM OR FINISH CARD INSTEAD READ THE
     1 FOLLOWING'/1X,A6)
 1002 FORMAT(' NUMBER OF VARIABLES CANNOT EXCEED 90')
 1004 FORMAT(' NUMBER OF DELETE CARDS INCORRECTLY SPECIFIED')
 1005 FORMAT(' PROGRAM EXPECTED DELETE CARD INSTEAD READ THE FOLLOWING'/
     11X,A6)
 1008 FORMAT(' COLUMNS 7 AND 8 OF DELETE CARD',I4,' CONTAIN A NEGATIVE N
     1UMBER')
 1010 FORMAT(' NUMBER OF DELETION CODES INCORRECTLY SPECIFIED')
 1012 FORMAT(' NUMBER OF VARIABLES AFTER TRANSGENERATION CANNOT EXCEED 9
     10')
 1014 FORMAT(' NUMBER OF TRANSGENERATION CARDS INCORRECTLY SPECIFIED')
 1016 FORMAT(' PROGRAM EXPECTED TRNGEN CARD INSTEAD READ THE FOLLOWING'/
     11X,A6)
 1017 FORMAT(' VARIABLE FORMAT CARD(S)'/1X,18A4)
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
      END
      FUNCTION SIGNE(A,B)
      EXTERNAL SIGN
      SIGNE=SIGN(A,B)
      RETURN
      END
C             FUNCTION KEY FOR BMD03D                MARCH  1, 1965
      FUNCTION KEY (I,J,N)
      IF(I-J)5,5,7
    5 II=I
      JJ=J
      GO TO 10
    7 II=J
      JJ=I
   10 KEY=0
      IX=II-1
      IF(IX.LE.0)GO TO 20
   15 NN=N
      DO 17 IJ=1,IX
      NN=NN-1
   17 KEY=KEY+NN
   20 KEY=KEY+JJ-II
      RETURN
      END
C             SUBROUTINE PRINT FOR BMD03D            MARCH  1, 1966
      SUBROUTINE PRINT(NF,NL,N1,N2)
      DIMENSION SUMX(2000),SUMY(2000), SUMXY(2000),ITEM(2000),LL(90),
     1ICASE(90),JUNK(90),FMT(180),Q(90),FM(18)
      COMMON  FMT
      COMMON  SUMX   , SUMY   , SUMXY  , ITEM   , NVNEW  , JUNK
      COMMON  NV     , STAR
      EQUIVALENCE (FMT,LL)
      DATA EMPTY/'    '/
C
      N3=1
      IF(N2)1,1,2
 1    WRITE (6,902)(LL(I),I= NF,NL)
      GO TO 70
 2    IF(N1)3,60,4
 3    N3=N3+1
 4    N3=N3+1
      READ (5, 9999)FM
 9999 FORMAT(18A4)
 70   ASSIGN 27 TO ISKIP
      ASSIGN 60 TO JSKIP
      DO 40 I=1,NVNEW
      K=0
      DO 30 J=NF,NL
      K=K+1
      IF(I-J)25,20,25
   20 SUMY(K)=1.0
      Q(K)=EMPTY
      ICASE(K)=JUNK(I)
      GO TO 30
   25 KK=KEY(I,J,NVNEW)
      SUMY(K)=SUMX(KK)
      ICASE(K)=ITEM(KK)
      Q(K)=SUMXY(KK)
      GO TO ISKIP,(27,30)
   27 IF(SUMXY(KK).NE.STAR)GO TO 30
 28   ASSIGN 30 TO ISKIP
      ASSIGN 50 TO JSKIP
   30 CONTINUE
      GO TO (5,6,7),N3
 5    WRITE (6,900)I,(SUMY(IJ),Q(IJ),IJ=1,K)
      WRITE (6,901)(ICASE(IJ),IJ=1,K)
      GO TO 40
    6 WRITE (N1,FM)(SUMY(IJ),IJ=1,K)
      GO TO 40
    7 PUNCH FM ,(SUMY(IJ),IJ=1,K)
   40 CONTINUE
      N3=1
      GO TO JSKIP,(50,60)
 50   WRITE (6,903)
 60   RETURN
  900 FORMAT(1H0I3,3H  *F10.5,9(A1,F9.5),A1)
  901 FORMAT(1H 6X,10(3X,1H(I5,1H)))
  902 FORMAT(1H03X,10I10)
 903  FORMAT(1H0,14X,90HA $ INDICATES THE COEFFICIENT IS NOT COMPUTED DU
     1E TO A ZERO DIVISOR, A ZERO IS INSERTED.  )
      END
C        SUBROUTINE TPWD FOR BMD03D                  MARCH  1, 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
   15 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 TRANS FOR BMD03D            MARCH  1, 1966
      SUBROUTINE TRANS(DATA,NEW,JUMP,NA,BN,NCASE,II,NTG,ITRANS)
      DIMENSION DATA(90),NEW(99),JUMP(99),NA(99),BN(99)
      DATA BLANK/'    '/
      DATA PI,PI2/3.14159265,1.57079633/
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
C
      MAX=91
      SAMPX=NCASE
      DO 1010 I=1,NTG
      JOKE=JUMP(I)
      NW=NEW(I)
      IF(NW-MAX)80,97,97
 80   NAA=NA(I)
      IF(NAA-MAX)100,97,97
 100  BBN=BN(I)
      D1=DATA(NAA)
      GO TO(85,33),ITRANS
   85 IF(D1.EQ.BLANK)GO TO 96
   33 GO TO (1,2,3,4,5,6,7,8,9,10),JOKE
    1 IF(D1)99,11,12
   11 D2=0.0
      GO TO 1000
   12 D2=SQRT(D1)
      GO TO 1000
    2 IF(D1)99,13,14
   13 D2=1.0
      GO TO 1000
   14 D2=SQRT(D1)+SQRT(D1+1.0)
      GO TO 1000
 3    IF(D1.LE.0.0) GO TO 99
   15 D2=ALOG10(D1)
      GO TO 1000
    4 D2=EXP(D1)
      GO TO 1000
    5 IF(-D1)17,11,99
   17 IF(D1-1.0)18,19,99
   19 D2=PI2
      GO TO 1000
   18 D2=ASN(SQRT(D1))
      GO TO 1000
    6 A=D1/(SAMPX+1.0)
      B=A+1.0/(SAMPX+1.0)
      IF(A)99,27,24
   27 D2=ASN(SQRT(B))
      GO TO 1000
   24 IF(B)99,28,29
 28   IF(A-1.0)285,287,99
  285 D2=ASN(SQRT(A))*2.0
      GO TO 1000
  287 D2=PI
      GO TO 1000
   29 A=SQRT(A)
      B=SQRT(B)
      IF(A.GT.1.0.OR.B.GT.1.0)GO TO 99
  292 D2=ASN(A)+ASN(B)
      GO TO 1000
    7 IF(D1.EQ.0.0)GO TO 99
   31 D2=1.0/D1
      GO TO 1000
    8 D2=D1+BBN
      GO TO 1000
    9 D2=D1*BBN
      GO TO 1000
   10 IF(D1)99,11,32
   32 D2=D1**BBN
      GO TO 1000
   96 D2=BLANK
      GO TO 1000
 97   WRITE (6,4001)I
      GO TO 991
 99   WRITE (6,4000)NAA,JUMP(I),II
 991  WRITE (6,4002)
      GO TO 1010
 1000 DATA(NW)=D2
 1010 CONTINUE
  999 RETURN
 4000 FORMAT(17H0DATA OF VARIABLEI4,30H VIOLATES RESTRICTION FOR CODEI3,
     119H. THIS WAS FOR CASEI6)
 4001 FORMAT(20H0TRANSGENERATOR CARDI3,44H HAS VARIABLE NUMBERS TOO LARG
     1E FOR PROGRAM.)
 4002 FORMAT(32H THE DATA WILL REMAIN UNCHANGED.)
      END

























































