Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - decus/20-0137/bmd/bmdx84.for
There is 1 other file named bmdx84.for in the archive. Click here to see a list.
C        ASYMETRICAL CORRELATION WITH MISSING DATA OCTOBER 21, 1965
      DOUBLE PRECISION P,FIN,SUBPRO,PROBLM,VARSEM
      DIMENSION F(162),DATE(4),A(12000),Q(400)
      DATA MXST/12000/
C
      DATA PER,ONO,FIN,YES,SUBPRO,PROBLM,VARSEM/1H.,2HNO,6HFINISH,3HYES,
     X6HSUBPRO,6HPROBLM,6HVARSEL/
      DIMENSION PC(2)
      LOGICAL NBL
      NPR=0
	CALL USAGEB('BMDX84')
      NERR=0
      GO TO 100
  203 NERR=NERR+1
 202  NC=NC1
      NV=NV1
      NVA=NVA1
      IT=IT1
      NF=NF1
      ON=ON1
      GO TO 105
 100  READ (5,1) P,PC,NC,NV,NVA,IT,NF,ON
    1 FORMAT(A6,2A3,5I6,4X,A2)
 105  IF(P.EQ.FIN) STOP
      IF(P.NE.PROBLM)GO TO 399
C
C
C
C
      NERR=0
  405 NPR=NPR+1
      IF(IT.EQ.1.OR.IT.EQ.6)WRITE(6, 502)
      IF(IT.EQ.0) IT=5
      IF(IT.EQ.5) ON=ONO
      IF(ON.NE.ONO) ON=YES
      IF(ON.NE.ONO) REWIND IT
      IF(NF.GE.1.AND.NF.LE.9)GO TO 500
      NF=1
      WRITE(6, 501)
  500 NF=18*NF
      IF(IT.GE.0) READ (5,10) (F(I),I=1,NF)
10    FORMAT(18A4)
C
C          WRITE(6, HEADING
       WRITE (6,2) PC,NC,NV,NVA,IT,ON
 2    FORMAT ('1BMDX84 - ASYMMETRICAL CORRELATION WITH MISSING DATA',
     * ' - REVISED MAY 10, 1968' /
     X41H0HEALTH SCIENCES COMPUTING FACILITY, UCLA//
     X31H0PROBLEM CODE                            2A3/
     X31H0NUMBER OF CASES                         I6/
     X31H0NUMBER OF VARIABLES READ IN             I6/
     X31H0NUMBER OF VARIABLES ADDED               I6/
     X31H0INPUT TAPE NUMBER                       I6/
     X35H0REWIND INPUT TAPE                       A4)
      WRITE (6,22) (F(I),I=1,NF)
22    FORMAT(31H0INPUT FORMAT                     18A4/(31X,18A4))
C
C          SET ARGUMENTS AND CALL PASS1
      NV1=NV+NVA
      L1=1+NV1
      L11=L1+MAX0(NV1,NV)
      L2=L11+NV1
      L21=L2+NV1
      CALL PASS1(A,A(L11),A(L1),A(L2),A(L21),NV,NV1,NC,IT,F,NPR)
      NV=NV1
      IPP=0
C
C          READ SUBPROBLEM CARD
 205  READ (5,1) P,PC,NC1,NV1,NVA1,IT1,NF1,ON1
      IF(P.NE.SUBPRO)GO TO 203
      IPP=IPP+1
      N1=NC1
      N2=NV1
      SSS=PC(2)
      IF(SSS.NE.YES) SSS=ONO
C          WRITE(6, SUBPROBLEM INFORMATION
      WRITE (6,49) IPP,SSS
 49   FORMAT(31H1SUBPROBLEM NUMBER                I6/
     X35H0BLANKS TREATED AS MISSING               A4)
C
C          SET DIAGONAL LIMITS
      KH=1
      IF(NBL(N1).OR.NBL(N2)) GO TO 52
      N2=NV
      NRZ=0
 52   MM1=66
      M1=0
 99   M0=M1+1
      M1=M1+MM1
C
C          READ IN VARIABLE SELECTION CARDS
      READ (5,76) P,(Q(M),M=M0,M1)
   76 FORMAT(A6,66A1)
      IF(P.NE.VARSEM)GO TO 402
      NRZ=NRZ+1
      DO 98 M=M0,M1
 98   IF(Q(M).EQ.PER) GO TO (50,89),KH
      GO TO 99
 50   CALL VARSEL(A(L2),NXX,NV,Q,NRZ)
      L3=L2+NXX
      KH=2
      K2=L3
      J1=1
      WRITE (6,60) (Q(I),I=1,M1)
 60   FORMAT(31H0VARIABLE SELECTION               100A1/(31X,100A1))
      GO TO 52
 89   WRITE (6,61) (Q(I),I=1,M1)
 61   FORMAT(/(31X,100A1))
      WRITE (6,685) N1,N2
 685  FORMAT(31H0LOWER DIAGONAL LIMIT             I6/
     X31H0UPPER DIAGONAL LIMIT               I6)
 88   NP=0
      CALL VARSEL(A(L3),NYY,NV,Q,NRZ)
      II0=MIN0(MAX0(1,J1+N1),NXX)
      DO 13 J=J1,NYY
      II1=MIN0(MAX0(1,J+N1),NXX)
      II2=MIN0(MAX0(1,J+N2),NXX)
      NP1=NP+II2-II1+1
      NX=II2-II0+1
      NY=J-J1+1
      NST=3*NV+NXX+2*NX+3*NY+NP1
      NST1=NST+5*NP1+NX+NY
      IF(SSS.NE.YES) GO TO 41
      IF(NST1.GT.MXST) GO TO 30
      GO TO 13
 41   IF(NST.GT.MXST) GO TO 30
 13   NP=NP1
      J=NYY+1
 30   NY=J-J1
      II2=MIN0(MAX0(1,J-1+N2),NXX)
      NX=II2-II0+1
      L22=L2+II0-1
      M1=N1+J1-II0
      M2=N2+J1-II0
      J1=J
      K1=L3
      DO 18 I=1,NY
      A(K1)=A(K2)
      K1=K1+1
 18   K2=K2+1
      L4=L3+NY
      L5=L4+NX
      L6=L5+NY
      L7=L6+NX
      L8=L7+NY
      L9=L8+NX
      L10=L9+NY
      L13=L10+NP
      L14=L13+NP
      L15=L14+NP
      L16=L15+NP
      L17=L16+NP
      REWIND 1
      IF(SSS.NE.YES) GO TO 44
      CALL CORR(A,A(L1),A(L22),A(L3),A(L4),A(L5),A(L6),A(L7),A(L8),A(L9)
     X,A(L10),A(L13),A(L14),A(L15),A(L16),A(L17),NX,NY,M1,M2,NV,NC,NP,SS
     XS)
      GO TO 46
 44   CALL COR(A(L11),A(L1),A(L22),A(L3),A(L4),A(L5),A(L6),A(L7),A(L8),N
     XX,NY,M1,M2,NV,NC,NP,SSS)
 46   REWIND 1
      IF(J1.LE.NYY) GO TO 88
      GO TO 205
  399 IF(NERR)400,400,403
  400 WRITE(6, 4000)P
 4000 FORMAT(' PROGRAM EXPECTED PROBLM CARD INSTEAD READ THE FOLLOWING'/
     11X,A6)
  501 FORMAT(' NUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIFIED,ASS
     1UMED TO BE 1')
  502 FORMAT(' INPUT TAPE CANNOT EQUAL 1 OR 6')
      STOP
  402 WRITE(6, 4002)NRZ,P
 4002 FORMAT(' VARSEL CARD',1X,I4,' READS AS FOLLOWS',1X,A6)
      STOP
  403 WRITE(6, 4003)P
 4003 FORMAT(' PROGRAM EXPECTED SUBPRO OR PROBLM CARD INSTEAD READ THE
     1FOLLOWING '/1X,A6)
      STOP
      END
C        SUBROUTINE CON FOR BMDX84                 OCTOBER 21, 1965
      SUBROUTINE CON(N,B)
      DIMENSION A(10),B(6)
      DATA A,BL,ALP/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H ,1H(/
      M=N
      DO 1 L=1,6
 1    B(L)=BL
      L=6
      DO 2 I=1,5
      J=MOD(M,10)
      B(L)=A(J+1)
      M=M/10
      L=L-1
 2    IF(M.EQ.0) GO TO 7
 7    B(L)=ALP
      RETURN
      END
C        SUBROUTINE COR FOR BMDX84                 OCTOBER 21, 1965
      SUBROUTINE COR(EAN,H,IX,IY,X,Y,SSX,SSY,SXY,NX,NY,N1,N2,NV,NC,NP,SS
     XS)
      DIMENSION EAN(2),H(2),IX(2),IY(2),X(2),Y(2),SSX(2),SSY(2),SXY(2)
      DATA BL/1H /
      DO 1 I=1,NX
 1    SSX(I)=0.
      DO 2 J=1,NY
 2    SSY(J)=0.
      DO 3 K=1,NP
 3    SXY(K)=0.
      DO 6 L=1,NC
       READ (1) (H(I),I=1,NV)
  400  FORMAT ( 20A4 )
      DO 5 I=1,NX
      I1=IX(I)
      IF(H(I1).EQ.BL) H(I1)=0.
      X(I)=H(I1)-EAN(I1)
 5    SSX(I)=SSX(I)+X(I)*X(I)
      K=1
      DO 6 J=1,NY
      II1=MIN0(MAX0(1,J+N1),NX)
      II2=MIN0(MAX0(1,J+N2),NX)
      J1=IY(J)
      IF(H(J1).EQ.BL) H(J1)=0.
      Y(J)=H(J1)-EAN(J1)
      SSY(J)=SSY(J)+Y(J)*Y(J)
      DO 6 I=II1,II2
      SXY(K)=SXY(K)+X(I)*Y(J)
 6    K=K+1
      K=1
      DO 7 J=1,NY
      II1=MIN0(MAX0(1,J+N1),NX)
      II2=MIN0(MAX0(1,J+N2),NX)
      DO 7 I=II1,II2
      A=-0.
      IF(SSX(I)*SSY(J).LE.0.0) GO TO 8
      A=SXY(K)/SQRT(SSX(I)*SSY(J))
 8    SXY(K)=A
 7    K=K+1
      CALL OUTPUT(SXY,MXY,IX,IY,NX,NY,N1,N2,NP,NC,SSS)
      RETURN
      END
C        SUBROUTINE CORR FOR BMDX84                OCTOBER 21, 1965
      SUBROUTINE CORR(EAN,H,IX,IY,X,Y,SSX,SSY,MX,MY,SXMY,SYMX,SSXMY,SSYM
     XX,SXY,MXY,NX,NY,N1,N2,NV,NC,NP,SSS)
      DIMENSION IX(2),IY(2),X(2),Y(2),SSX(2),SSY(2),MX(2),MY(2),SXMY(2),
     XSYMX(2),SSXMY(2),SSYMX(2),SXY(2),MXY(2),EAN(2),H(2)
      DATA AMISS/1H /
      DO 1 I=1,NX
      SSX(I)=0.
 1    MX(I)=0
      DO 2 I=1,NY
      SSY(I)=0.
 2    MY(I)=0
      DO 3 I=1,NP
      SXY(I)=0.
      MXY(I)=0
      SXMY(I)=0.
      SYMX(I)=0.
      SSXMY(I)=0.
 3    SSYMX(I)=0.
      DO 8 L=1,NC
       READ (1) (H(I),I=1,NV)
  400  FORMAT ( 20A4 )
      DO 6 I=1,NX
      I1=IX(I)
      X(I)=H(I1)
      IF(X(I).EQ.AMISS) GO TO 7
      X(I)=X(I)-EAN(I1)
      SSX(I)=SSX(I)+X(I)*X(I)
      GO TO 6
 7    MX(I)=MX(I)+1
 6    CONTINUE
      K=1
      DO 8 J=1,NY
      II1=MIN0(MAX0(1,J+N1),NX)
      II2=MIN0(MAX0(1,J+N2),NX)
      J1=IY(J)
      Y(J)=H(J1)
      IF(Y(J).EQ.AMISS) GO TO 9
      Y(J)=Y(J)-EAN(J1)
      SSY(J)=SSY(J)+Y(J)*Y(J)
      DO 10 I=II1,II2
      IF(X(I).EQ.AMISS) GO TO 11
      SXY(K)=SXY(K)+X(I)*Y(J)
      GO TO 10
 11   SYMX(K)=SYMX(K)+Y(J)
      SSYMX(K)=SSYMX(K)+Y(J)*Y(J)
 10   K=K+1
      GO TO 8
 9    MY(J)=MY(J)+1
      DO 12 I=II1,II2
      IF(X(I).EQ.AMISS) GO TO 13
      SXMY(K)=SXMY(K)+X(I)
      SSXMY(K)=SSXMY(K)+X(I)*X(I)
      GO TO 12
 13   MXY(K)=MXY(K)+1
 12   K=K+1
 8    CONTINUE
      K=1
      DO 16 J=1,NY
      II1=MIN0(MAX0(1,J+N1),NX)
      II2=MIN0(MAX0(1,J+N2),NX)
      DO 15 I=II1,II2
      MXY(K)=NC-MX(I)-MY(J)+MXY(K)
      XY=MXY(K)
      A=SSX(I)-SSXMY(K)
      B=SSY(J)-SSYMX(K)
      C=SXY(K)
      IF(XY.LE.0.0)GO TO 14
      A=A-SXMY(K)*SXMY(K)/XY
      B=B-SYMX(K)*SYMX(K)/XY
      C=C-SXMY(K)*SYMX(K)/XY
   14 SXY(K)=-0.0
      IF(A*B.LE.0.) GO TO 15
      SXY(K)=C/SQRT(A*B)
 15   K=K+1
 16   CONTINUE
      CALL OUTPUT(SXY,MXY,IX,IY,NX,NY,N1,N2,NP,NC,SSS)
      RETURN
      END
C        FUNCTION NBL FOR BMDX84                   OCTOBER 21, 1965
      LOGICAL FUNCTION NBL(X)
      EXTERNAL SIGN
      NBL=.TRUE.
      IF(SIGN(1.0,X).NE.1.0.AND.X.EQ.0.0)NBL=.FALSE.
      RETURN
      END
C        SUBROUTINE OUTPUT FOR BMDX84              OCTOBER 21, 1965
      SUBROUTINE OUTPUT(R,MC,IX,IY,NX,NY,N1,N2,NP,NC,SSS)
      DATA YES/3HYES/
      DIMENSION R(2),MC(2),IX(2),IY(2),II(2),JJ(2)
      DIMENSION RR(10),I1(10),J1(10),A1(6,10)
      JJ(1)=0
      II(1)=0
      II2=0
      K1=0
      L1=0
 6    L1=MIN0(L1+200,NP)
      WRITE (6,3)
 3    FORMAT(1H1,4X,10(12H    I   CORR)/5X,10(12H    J  COUNT))
 5    K0=K1+1
      K1=MIN0(K1+10,L1)
      I=0
      DO 8 K=K0,K1
      I=I+1
      II(1)=II(1)+1
      IF(II(1).LE.II2) GO TO 9
      JJ(1)=JJ(1)+1
      JJJ=JJ(1)
      II(1)=MIN0(MAX0(1,JJJ+N1),NX)
      II2=MIN0(MAX0(1,JJJ+N2),NX)
9     I1(I)=IX(II(1))
      J1(I)=IY(JJ(1))
      RR(I)=R(K)
      MCK=NC
      IF(SSS.EQ.YES) MCK=MC(K)
 8    CALL CON(MCK,A1(1,I))
      WRITE (6,4) (I1(K),RR(K),K=1,I)
 4    FORMAT(1H0,4X,10(I5,F7.3))
      WRITE (6,7) (J1(K),(A1(J,K),J=1,6),K=1,I)
 7    FORMAT(5X,10(I5,6A1,1H)))
      IF(K1.LT.L1) GO TO 5
      IF(L1.LT.NP) GO TO 6
      RETURN
      END
C        SUBROUTINE PASS1 FOR BMDX84               OCTOBER 21, 1965
      SUBROUTINE PASS1(XM,XM1,X,SD,SD1,NV,NV1,NC,IT,F1,NPR)
C
C          READ IN DATA MATRIX, INITIATE TRANSGENERATIONS AND SELECTIONS
C          WRITE RESULTS ON INTERMEDIATE TAPE
C          COMPUTE AND WRITE(6, MEANS AND STANDARD DEVIATIONS
C
      DIMENSION XM(2),XM1(2),X(2),F1(162),SD(2),SD1(2)
      LOGICAL NBL
      DATA BIN,BL/3HBIN,1H /
      BBCD=1.
      IF(IT.LT.0) BBCD=BIN
      IT=IABS(IT)
      REWIND 1
      DO 11 I=1,NV1
      SD(I)=0.
      XM1(I)=0.
 11   XM(I)=0.
      MY=0
      DO 12 J=1,NC
      SELECT=1.
      IF(BBCD.EQ.BIN) READ (IT) (X(I),I=1,NV)
      IF(BBCD.NE.BIN) READ (IT,F1) (X(I),I=1,NV)
      CALL TRANS(X,NV1,J,NPR,SELECT)
      IF(SELECT.EQ.0.) GO TO 12
      MY=MY+1
      DO 14 I=1,NV1
      IF(NBL(X(I))) GO TO 16
      XM1(I)=XM1(I)+1.
      X(I)=BL
      GO TO 14
 16   XM(I)=XM(I)+X(I)
      SD(I)=SD(I)+X(I)*X(I)
 14   CONTINUE
       WRITE (1) (X(I),I=1,NV1)
  410  FORMAT ( 20A4 )
 12   CONTINUE
      NC=MY
      H=NC
      DO 15 I=1,NV1
      S=XM(I)
      X(I)=H-XM1(I)
      XM(I)=0.0
      IF(X(I).GT.0.0)XM(I)=S/X(I)
      XM1(I)=S/H
      SD1(I)=SQRT((SD(I)-S*XM1(I))/(H-1.))
 15   SD(I)=SQRT((SD(I)-S*XM(I))/(X(I)-1.))
      WRITE (6,1) (I,XM1(I),SD1(I),H,XM(I),SD(I),X(I),I=1,NV1)
 1    FORMAT(76H1                BLANKS TREATED AS ZEROS           BLANK
     XS CONSIDERED MISSING/
     X77H0VARIABLE       MEAN       STD. DEV.  COUNT       MEAN       ST
     XD. DEV.  COUNT/(I6,3X,2F13.5,F8.0,2F13.5,F8.0))
      ENDFILE 1
      REWIND 1
      RETURN
      END
C        SUBROUTINE VARSEL FOR BMDX84              OCTOBER 21, 1965
      SUBROUTINE VARSEL(LL,NX,NV,Q,NRZ)
C
C          SELECT VARIABLES SPECIFIED FOR SUBPROBLEM
C
      DIMENSION Q(2), LL(2), R(18)
      DATA R/1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H0,1H,,1H-,1H.,1H ,1H(
     X,1H),1H(,1H)/
      L=0
      I=0
 1    NN=0
 19   M=0
      M1=1
 6    L=L+1
      DO 2 J=1,18
 2    IF(Q(L).EQ.R(J)) GO TO 3
      L=L-(NRZ-1)*66+6
      WRITE(6, 101)L
  101 FORMAT(' CHARACTER IN COLUMN',I2,' IS ILLEGAL')
 30   WRITE (6,4) NRZ
    4 FORMAT(21H1ERROR ON VARSEL CARDI4)
   60 STOP
 3    IF(J.GT.10) GO TO 5
      J=MOD(J,10)
      M=M*10+J
      GO TO 6
 9    M1=M
      M=0
      GO TO 6
 5    J=J-10
      GO TO (11,11,11,6,11,9,11,9), J
   40 WRITE(6, 100)
  100 FORMAT(' INDEX OF VARIABLE SELECTED EXCEEDS TOTAL NUMBER OF VARIAB
     1LES')
      GO TO 30
   50 WRITE(6,102),NRZ
  102 FORMAT(' INDICES ARE REVERSED ON VARSEL CARD',I4)
      STOP
 11   IF(M.GT.NV) GO TO 40
      I=I+1
      LL(I)=M
      IF(NN.EQ.0) GO TO 10
      I=I-1
      L1=LL(I)+M1
      IF(M.LT.L1) GO TO 50
      DO 13 K=L1,M,M1
      I=I+1
 13   LL(I)=K
 10   NN=1
      GO TO (1,19,20,6,19,9,19,9), J
 20   NX=I
      RETURN
      END
C        SUBROUTINE TRANS FOR BMDX84               OCTOBER 21, 1965
      SUBROUTINE TRANS(X,NV,NC,NPR,SELECT)
      DIMENSION X(NV)
      RETURN
      END