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