Trailing-Edge
-
PDP-10 Archives
-
-
There are no other files named in the archive.
CID BMDX63
DOUBLE PRECISION PP,P1,P2,NA
DOUBLE PRECISION XX,X
DIMENSION NA(9)
COMMON XX(5000),X(100),FMT(324),NI,MM,ND,NH,MT,NO,
1KA,KB,NT,NMP,QA,QC,QD
DATA P1,P2/6HFINISH,6HPROBLM/
CALL USAGEB('BMDX63')
20 READ (5,1) PP,NI,ND,NH,NO,MT,NF,NA
1 FORMAT(A6,3I2,I4,2I2,8A6,A4)
IF(PP.EQ.P1)GO TO 22
23 IF(PP.EQ.P2)GO TO 21
24 WRITE (6,25)
CALL EXIT
25 FORMAT(30H0ILLEGAL PROBLM OR FINISH CARD )
21 NT=NI+ND
WRITE (6,82) NA,NI,ND,NO,NH
82 FORMAT('1BMDX63 - MULTIVARIATE GENERAL LINEAR HYPOTHESIS - REVISED
1 ','APRIL 14, 1969'/
X41H HEALTH SCIENCES COMPUTING FACILITY, UCLA//1H0,8X,9A6//
129H NO. OF INDEPENDENT VARIABLES I5/
229H NO. OF DEPENDENT VARIABLES I5/
329H NO. OF OBSERVATIONS I5/
429H NO. OF HYPOTHESES I5///)
IF(NT*NT+ND*ND- 5000)3,2,2
4 FORMAT(1H1,9A6//22H TOO MANY VARIABLES )
2 WRITE (6,4) NA
22 CALL EXIT
3 LLLL=0
IF(MT)40,40,41
40 MT=5
GO TO 42
41 IF(MT-5)43,42,43
43 LLLL=1
REWIND MT
42 NF=18*MAX0(NF,1)
WRITE(6,50)
50 FORMAT(///38H0D-TYPE VARIABLE FORMAT CARDS FOR DATA)
READ (5,8) (FMT(I),I=1,NF)
8 FORMAT(18A4)
WRITE(6,9)(FMT(I),I=1,NF)
9 FORMAT(1X,18A4)
L=NI+1
M=NI*NT+1
N=M+NI
K=NT*NT+1
CALL DOIT1(XX,XX(M),XX(L),XX(N),XX(K),NT,ND)
IF(LLLL)20,20,46
46 REWIND MT
GO TO 20
END
SUBROUTINE DOIT1(XX,XY,YX,YY,P,NT1,NT2)
DOUBLE PRECISION XX(NT1,NT1),XY(NT1,NT1)
DOUBLE PRECISION YX(NT1,NT1),YY(NT1,NT1),P(NT2,NT2)
DOUBLE PRECISION XXX,X
1,D,TT
DOUBLE PRECISION H,UIQ,HYP
COMMON XXX(5000),X(100),FMT(324),NI,MM,ND,NH,MT,NO,
1KA,KB,NT,NMP,QA,QC,QD
DATA HYP/6HHYPOTH/
DO 85 I=1,NT
NNN=I
DO 85 J=1,NNN
85 XX(I,J)=0.0
DO 6 N=1,NO
READ (MT,FMT) (X(I),I=1,NT)
DO 6 I=1,NT
DO 6 J=1,I
6 XX(I,J)=XX(I,J)+X(I)*X(J)
DO 7 I=1,NT
DO 7 J=1,I
7 XX(J,I)=XX(I,J)
WRITE (6,5000)
5000 FORMAT(34H-CROSS PRODUCT MATRIX (X,Y)'(X,Y)//)
DO 400 I=1,NT
400 WRITE (6,99) I,(XX(I,J),J=1,NT)
CALL MATINV(XX,NI,D,NT,TT)
WRITE (6,5005) TT
5005 FORMAT(1H-,10X,2H-1/6X,5H(X'X)5X,10H TOLERANCE F12.8//)
DO 405 I=1,NI
405 WRITE (6,99) I,(XX(I,J),J=1,NI)
50 CALL MPY(XX,XY,XY,NI,NI,ND,NT,NT,NT,1)
CALL MPY(YX,XY,P,ND,NI,ND,NT,NT,ND,1)
DO 10 I=1,ND
DO 10 J=1,ND
YY(I,J)=YY(I,J)-P(I,J)
10 P(I,J)=YY(I,J)
WRITE (6,5002)
5002 FORMAT(1H-32X,2H-1/38H REGRESSION COEFFICIENTS B=(X'X) X'Y//)
DO 402 I=1,NI
402 WRITE (6,99) I,(XY(I,J),J=1,ND)
99 FORMAT(I4,10F12.5/(4X,10F12.5))
WRITE (6,5003)
5003 FORMAT(43H-RESIDUAL SUM OF PRODUCT MATRIX E=Y'Y-Y'XB//)
DO 71 I=1,ND
71 WRITE (6,99) I,(YY(J,I),J=1,ND)
NMP=NO-NI
READ (5,66) (FMT(I),I=1,18)
66 FORMAT(18A4)
WRITE(6,6666)
6666 FORMAT(39H0D-TYPE VARIABLE FORMAT CARD FOR HYPOTH)
WRITE(6,77)(FMT(I),I=1,18)
77 FORMAT(1X,18A4)
IF(NH)80,80,82
80 CALL EXIT
82 DO 2 IQ=1,NH
READ (5,1) H,UIQ,MM,MNM,NMN
1 FORMAT(2A6,3I2)
QA=MM
QC=MNM
QD=NMN
IF(MM.NE.0) GO TO 300
MM=NI
300 IF(MNM.NE.0) GO TO 301
MNM=ND
301 IF(H.NE.HYP)GO TO 1111
302 IF(MM.GT.NI)GO TO 2222
303 IF(MNM.GT.ND)GO TO 3333
KB=MAX0(ND,NI)
KA=NI
L=NT*NT+1
M=KA*MAX0(ND,MM)+L
N=KB*MM+M
NNNN=N+ND*MAX0(MNM,MM)
IF(ND*MNM+NNNN-5000)4,3,3
1111 WRITE(6,1119)UIQ
1119 FORMAT(//40H0HYPOTH IS PUNCHED WRONG ON HYPOTH CARD ,A6)
STOP
2222 WRITE(6,2229)UIQ
2229 FORMAT(//46H COL.13,14 ARE SPECIFIED WRONG ON HYPOTH CARD ,A6)
STOP
3333 WRITE(6,3339)UIQ
3339 FORMAT(//46H COL.15,16 ARE SPECIFIED WRONG ON HYPOTH CARD ,A6)
STOP
8 FORMAT(1X,I4,2X,A6,5X,30HTHIS HYPOTHESIS IS TOO LARGE )
3 WRITE (6,8) IQ,UIQ
DO 200 I=1,MM
200 READ (5,FMT) (X(J),J=1,NI)
DO 201 I=1,MNM
201 READ (5,FMT) (X(J),J=1,ND)
DO 202 I=1,MM
202 READ (5,FMT) (X(J),J=1,MNM)
GO TO 2
4 WRITE (6,6000) IQ,UIQ
6000 FORMAT(1H080(1H*)/29H-OUTPUT FOR HYPOTHESIS NUMBERI4,4X,A6)
CALL DOIT2(XX,XY,YX,YY,XXX(N),XXX(L),XXX(M),XXX(NNNN),NT,ND,KA,KB,
1MNM)
2 CONTINUE
RETURN
END
SUBROUTINE DOIT2(XX,XY,YX,YY,C,A,B,CEC,NT1,ND1,KA1,KB1,MMM)
DOUBLE PRECISION XXX,X,DE,XX,XY,YX,YY
1,C,A,B,DEH,TT,CEC
COMMON XXX(5000),X(100),FMT(324),NI,MM,ND,NH,MT,NO,
1KA,KB,NT,NMP,QA,QC,QD
DIMENSION XX(NT1,NT1),XY(NT1,NT1), YX(NT1,NT1),YY(NT1,NT1),C(ND1,N
1D1),A(KA1,KA1),B(KB1,KB1),CEC(ND1,ND1)
DATA Z1,Z2,Z3/1HA,1HC,1HD/
WRITE (6,201) Z1
201 FORMAT(1H-A2,6HMATRIX//)
DO 3 I=1,MM
IF(QA)500,500,501
500 DO 502 J=1,NI
502 A(J,I)=0.
A(I,I)=1.
GO TO 3
501 READ (5,FMT) (A(J,I),J=1,NI)
3 WRITE (6,200) I,(A(J,I),J=1,NI)
200 FORMAT(I4,10F12.5/(4X,10F12.5))
CALL MPY(XY,A,B,ND,NI,MM,NT,KA,KB,2)
WRITE (6,201) Z2
DO 20 I=1,MMM
IF(QC)508,508,509
508 DO 510 J=1,ND
510 C(J,I)=0.
C(I,I)=1.
GO TO 20
509 READ (5,FMT) (C(J,I),J=1,ND)
20 WRITE (6,200) I,(C(J,I),J=1,ND)
CALL MPY(YY,C,CEC,ND,ND,MMM,NT,ND,ND,2)
CALL MPY(C,CEC,CEC,MMM,ND,MMM,ND,ND,ND,2)
CALL MPY(C,B,B,MMM,ND,MM,ND,KB,KB,2)
WRITE (6,201) Z3
DO 21 J=1,MM
IF(QD)564,564,566
564 DO 562 I=1,MMM
562 C(I,J)=0.
GO TO 21
566 READ (5,FMT) (C(I,J),I=1,MMM)
21 WRITE (6,200) J,(C(I,J),I=1,MMM)
NDD=ND
ND=MMM
DO 1 I=1,ND
DO 1 J=1,MM
1 C(I,J)=B(I,J)-C(I,J)
WRITE (6,6001)
6001 FORMAT(17H-CONTRAST ABC'-D//)
DO 801 I=1,MM
801 WRITE (6,200) I,(C(J,I),J=1,ND)
CALL MPY(XX,A,B,NI,NI,MM,NT,KA,KB,1)
CALL MPY(A,B,B,MM,NI,MM,KA,KB,KB,2)
WRITE (6,6000)
6000 FORMAT(1H-10X,2H-1/5X,10HA(X'X) A'//)
DO 800 I=1,MM
800 WRITE (6,200) I,(B(I,J),J=1,MM)
CALL MATINV(B,MM,D,KB,TT)
80 CALL MPY(B,C,A,MM,MM,ND,KB,NDD,KA,3)
DO 10 I=1,ND
DO 11 J=1,MM
11 X(J)=C(I,J)
DO 12 J=1,ND
12 C(I,J)=0.
DO 10 K=1,ND
DO 10 J=1,MM
10 C(I,K)=C(I,K)+X(J)*A(J,K)
WRITE (6,6009)
6009 FORMAT(1H-4X,4HCEC'//)
DO 301 I=1,ND
301 WRITE (6,200) I,(CEC(I,J),J=1,ND)
WRITE (6,6007)
6007 FORMAT(1H-52X,7H-1 -1/68H HYPOTHESIS SUM OF PRODUCT MATRIX H=(A
1BC'-D)'(A(X'X) A') (ABC'-D)//)
DO 810 I=1,ND
810 WRITE (6,200) I,(C(I,J),J=1,ND)
DO 2 I=1,ND
DO 2 J=1,ND
2 C(I,J)=C(I,J)+CEC(I,J)
CALL MATINV(C,ND,DEH,NDD,TT)
CALL MATINV(CEC,ND,DE,NDD,TT)
V=DE/DEH
WRITE (6,310) DE,DEH,V,ND,MM,NMP
310 FORMAT(20H-DETERMINANT(CEC') F20.6/20H0DETERMINANT(CEC'+H)
1F20.6/20H0U-STATISTIC F20.6,3X,18HDEGREES OF FREEDOM3I4)
ZS=ND
ZR=MM
ZN=NMP
ZD=ZR*ZR+ZS*ZS-5.
ZT=1.
IF(ZD.EQ.0.) GO TO 50
ZT=SQRT((ZR*ZR*ZS*ZS-4.)/ZD)
50 ZH=(ZN-.5*(ZS-ZR+1.))*ZT-ZR*ZS*.5+1.
IDF1=ZR*ZS
ZY=V**(1./ZT)
ZF=(1.-ZY)/ZY*ZH/(ZR*ZS)
WRITE (6,55) ZF,IDF1,ZH
55 FORMAT(20H0F-STATISTIC F20.6,3X,18HDEGREES OF FREEDOM
XI4,F7.0)
ND=NDD
RETURN
END
SUBROUTINE MATINV(A,N,D,ND,T)
DOUBLE PRECISION A,U,C,T,D,X,Y,G
DIMENSION A(ND,ND),IND(100),U(100),C(100)
DO 1 I=1,N
C(I)=A(I,I)
1 IND(I)=0
K=1
T=1.
D=1.
Y=1.
DO 7 L=1,N
IF(Y.LT.T)T=Y
IF(K)9,9,10
10 IND(K)=1
DO 2 I=1,K
U(I)=A(K,I)
2 A(K,I)=0.
X=U(K)
DO 3 I=K,N
U(I)=A(I,K)
3 A(I,K)=0.
U(K)=1.
D=D*X
Y=0.
DO 7 I=1,N
DO 4 J=1,I
4 A(I,J)=A(I,J)-U(I)*U(J)/X
IF(IND(I))5,5,7
5 G=A(I,I)/C(I)
IF(Y-G)6,7,7
6 Y=G
K=I
7 CONTINUE
DO 11 I=1,N
DO 11 J=1,I
A(I,J)=-A(I,J)
11 A(J,I)=A(I,J)
RETURN
9 D=0.
RETURN
END
SUBROUTINE MPY(A,B,C,II,JJ,KK,L,M,N,NNN)
DOUBLE PRECISION A,B,C,U
DIMENSION A(L,L),B(M,M),C(N,N),U(100)
DO 1 K=1,KK
GO TO (2,2,7),NNN
2 DO 4 J=1,JJ
4 U(J)=B(J,K)
3 GO TO (5,6,7),NNN
5 DO 8 I=1,II
C(I,K)=0.
DO 8 J=1,JJ
8 C(I,K)=C(I,K)+A(I,J)*U(J)
GO TO 1
6 DO 9 I=1,II
C(I,K)=0.
DO 9 J=1,JJ
9 C(I,K)=C(I,K)+A(J,I)*U(J)
GO TO 1
7 DO 10 I=1,II
C(I,K)=0.
DO 10 J=1,JJ
10 C(I,K)=C(I,K)+A(I,J)*B(K,J)
1 CONTINUE
RETURN
END