Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - decus/20-0137/bmd/bmdx63.for
There is 1 other file named bmdx63.for in the archive. Click here to see a list.
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