PDP-10 Archive: decus/20-0137/bmd/bmd06v.for from decuslib20-05

Trailing-Edge - PDP-10 Archives - decuslib20-05 - decus/20-0137/bmd/bmd06v.for
There is 1 other file named bmd06v.for in the archive. Click here to see a list.
C     GENERAL LINEAR HYPOTHESIS WITH CONTRASTS         AUGUST  8, 1966
C     THIS IS A SYSTEM/360 FORTRAN IV (LEVEL H) VERSION OF BMD06V
C     ORIGINALLY WRITTEN IN FORTRAN II.  THE PROGRAM HAS BEEN SIFTED
C     AND SLIGHTLY MODIFIED TO MAKE IT OPERABLE ON THE 7094.
C     IT HAS BEEN FURTHER MODIFIED TO ALLOW ITS EXECUTION ON THE 360.
C
      DIMENSION DATA(100),FMT(10),B(61),C(61,61),CI(61,61),NV(61),
     1MV(61),NN(61),P(6),A(61,61),XM(61),GUNK(180),Q(6)
      DIMENSION DUMFMT(120)
      DOUBLE PRECISION P,RPR,GOBLE0,GOBLE1,GOBLE2,RRD,FMT,
     1 RRE,URGLY,COEFFS
      COMMON  DATA   , GUNK   , B      , C      , CI     , NV
      COMMON  MV     , NN     , Q      , A      , PR     , RD
      COMMON  NP     , NO     , NT     , NEW    , NH     , NC
      COMMON  ND     , ON     , RE     , N1     , N2     , N3
      COMMON  N4     , NNT    , NPP    , MMT    , MD     , OON
      COMMON  NTAPE  , XM     , NCOV   , NT1    , NT2
C
  902 FORMAT('1BMD06V - GENERAL LINEAR HYPOTHESIS WITH CONTRASTS - ',
     1'REVISED SEPTEMBER 18, 1968'/
     241H HEALTH SCIENCES COMPUTING FACILITY, UCLA)
C
      DATA  P            /6H66F1.0,6H33F2.0,6H22F3.0,6H16F4.0,6H13F5.0,
     16H11F6.0/
      NTAPE=5
	CALL USAGEB('BMD06V')
      NT1=2
      NT2=3
      REWIND NT1
      REWIND NT2
   10 READ (5,900) RPR,RRD,NP,NO,NT,NEW,NH,NC,MD,MTAPE,NCOV,ND
      DATA GOBLE0/6HPROBLM/
      IF(RPR.EQ. GOBLE0) GO TO 30
      DATA GOBLE1/6HFINISH/
 15   IF(RPR.EQ. GOBLE1) GO TO 1000
      WRITE(6,10000)
10000 FORMAT(87H0PROBLM OR FINISH IS PUNCHED WRONG ON THE PROPER CARD OR
     1 EITHER CARD IS OUT OF SEQUENCE)
      STOP
6000  WRITE(6,7000)
7000  FORMAT(55H0COVARS IS PUNCHED WRONG OR COVARS CARD OUT OF SEQUENCE)
      STOP
8000  WRITE(6,9000)
9000  FORMAT(55H0COEFFS IS PUNCHED WRONG OR COEFFS CARD OUT OF SEQUEBCE)
      STOP
 20   WRITE (6,901)
      STOP
C
C
C
22    WRITE(6,1022)
 1022 FORMAT(55H0VIOLATING THE LIMITATION FOR TOTAL NUMBER OF VARIABLES)
      STOP
 66   WRITE(6,666)
666   FORMAT(41H0WRONG SET UP IN COL.13-72 ON CNTRST CARD)
      GO TO 170
 2000 WRITE(6,3000)
 3000 FORMAT(39H0ILLEGAL NUMBER OF CASES ON PRMBLM CARD)
      STOP
 4000 WRITE(6,5000)
 5000 FORMAT(45H0ILLEGALLY SPECIFIED NUMBER OF CNTRST CARD(S))
      STOP
 30   IF((MTAPE-NT1)*(MTAPE-NT2))8,20,8
 8    CALL TPWD(MTAPE,NTAPE)
 34   NNT=NP+NEW
      MMT=NNT-1
      ON=NO
      IF(NNT.LE.0.OR.NNT.GT.60)GO TO 22
C
      WRITE(6,902)
      IF(NO.LE.0) GO TO 2000
      IF(NH.LE.0) GO TO 4000
      WRITE(6,903)RRD,NNT,NO
      IF(NC) 40, 40, 35
 35   NPP=NP+1
      GO TO 45
 40   NPP=NP
 45   DO 47 I=1,NNT
 47   XM(I)=0.0
      CALL READ
      IF(NEW) 50, 50, 1000
 50   IF(NCOV) 53, 53, 52
 52   READ(5,924)RRD,(NN(J), J=1,NCOV)
      DATA URGLY/6HCOVARS/
      DATA COEFFS/6HCOEFFS/
      IF(RRD.NE.URGLY) GO TO 6000
 53   CALL ADJUST
      DO 170 IJ=1,NH
      READ(5,904)RRD,NT,ND,MD,(NV(I),I=1,MMT)
      DATA GOBLE2/6HCNTRST/
      IF(RRD.EQ. GOBLE2) GO TO 60
 55   WRITE (6,905) IJ
      GO TO 170
 60   WRITE (6,906)IJ,(NV(I),I=1,MMT)
      NP=0
      DO 70 I=1,MMT
      IF(NV(I)) 70, 70, 65
 65   NP=NP+1
      MV(NP)=I
 70   CONTINUE
      IF(NP) 66, 66, 75
 75   DO 80 J=1,NP
      N1=MV(J)
      DATA(J)=C(N1,NNT)
      DO 80 I=1,NP
      N2=MV(I)
 80   CI(I,J)=C(N2,N1)
      WRITE (6,907)
      DO 85 I=1,NP
      WRITE (6,908)I
 85   WRITE (6,909)(CI(I,J),J=1,NP),DATA(I)
      CALL INVERT (CI,NP,RD,NV,NN)
      WRITE (6,910)
      DO 90 I=1,NP
      WRITE (6,908)I
 90   WRITE (6,909)(CI(I,J),J=1,NP)
      DO 95 J=1,NP
      B(J)=0.0
      DO 95 I=1,NP
 95   B(J)=B(J)+CI(I,J)*DATA(I)
      IF(MD) 98, 98, 96
 96   WRITE (6,921)
 98   RD=0.0
      NDF=0
      DO 108 I=1,NO
      READ(NTAPE)(DATA(J),J=1,NNT),PR
      RE=0.0
      DO 100 J=1,NP
      N1=MV(J)
 100  RE=RE+DATA(N1)*B(J)
      IF(RE) 101, 108, 101
 101  OON=DATA(NNT)-RE
      NDF=NDF+1
      IF(MD) 105, 105, 102
 102  WRITE (6,922)I,DATA(NNT),RE,OON
 105  RD=RD+(OON**2)*PR
 108  CONTINUE
      REWIND NTAPE
      NDF=NDF-NP
      RE=NDF
      RD=RD/RE
      RE=SQRT(RD)
      DO 110 I=1,NP
 110  DATA(I)=SQRT(RD*CI(I,I))
      WRITE (6,911)
      WRITE (6,909)(B(I),I=1,NP)
      WRITE (6,912)
      WRITE (6,909)(DATA(I),I=1,NP)
      WRITE (6,913)RD
      WRITE (6,914)RE
      WRITE (6,923)NDF
      DATA FMT(1)/'(     '/
      DATA FMT(2)/'A6,   '/
      FMT(3)=P(ND)
      DATA FMT(4)/'/     '/
      FMT(5) = FMT(1)
      DATA FMT(6)/'6X,   '/
      FMT(7)=P(ND)
      DATA FMT(8)/',     '/
      DATA  FMT(9)/')     '/
      FMT(10) = FMT(9)
      DO 115 I=1,NT
 115  READ(5,FMT)RRE,(A(I,J),J=1,NP)
      IF(RRE.NE.COEFFS) GO TO 8000
      DO 120 I=1,NT
      DATA(I)=0.0
      DO 120 J=1,NP
 120  DATA(I)=DATA(I)+A(I,J)*B(J)
      DO 130 I=1,NT
      DUMFMT(I)=0.0
      DO 125 K=1,NP
      DO 125 J=1,NP
  125 DUMFMT(I)=DUMFMT(I)+A(I,K)*A(I,J)*CI(J,K)
  130 DUMFMT(I)=SQRT(DUMFMT(I)*RD)
      DO 132 I=1,NT
 132  NN(I)=I
      WRITE (6,915)
      N1=1
      N2=0
 135  IF(NT-7) 140, 140, 145
 140  N2=N2+NT
      GO TO 150
 145  N2=N2+7
 150  WRITE (6,920)(NN(J),J=N1,N2)
      DO 160 I=1,NP
 160  WRITE (6,916)MV(I),(A(J,I),J=N1,N2)
      WRITE (6,917)(DATA(J),J=N1,N2)
      WRITE (6,918) (DUMFMT(J),J=N1,N2)
      WRITE (6,919)
      NT=NT-7
      IF(NT) 170, 170, 165
 165  N1=N1+7
      GO TO 135
 170  CONTINUE
      GO TO 10
 900  FORMAT(A6,A2,I2,I4,7I2,42X,I2)
 901  FORMAT(33H1ERROR ON ASSIGNMENT OF TAPE UNIT)
 903  FORMAT(17H0PROBLEM NUMBER  A2//17H NO. OF VARIABLESI4/15H SAMPLE S
     1IZE   I6//)
 904  FORMAT(A6,3I2,60I1)
905   FORMAT(43H0CNTRST IS PUNCHED WRONG ON CNTRST CARD NO.,I4,31H OR TH
     1E CARD IS OUT OF SEQUENCE)
 906  FORMAT(1H0//9H0CONTRAST/9H CARD NO.I3,5X,60I1)
 907  FORMAT(1H0//23H SUMS OF CROSS PRODUCTS)
 908  FORMAT(4H0ROWI3)
 909  FORMAT(8F15.5)
 910  FORMAT(1H0//31H INVERTED CROSS PRODUCTS MATRIX)
 911  FORMAT(1H0//24H REGRESSION COEFFICIENTS)
 912  FORMAT(47H0STANDARD DEVIATIONS OF REGRESSION COEFFICIENTS)
 913  FORMAT(1H0//21H VARIANCE OF ESTIMATEF15.5)
 914  FORMAT(23H STD. ERROR OF ESTIMATEF13.5//)
 915  FORMAT(24H0CONTRASTS AND ESTIMATES)
 916  FORMAT(1H I5,3X,7F15.4)
 917  FORMAT(10H0ESTIMATE 7F15.5)
 918  FORMAT(9H0STANDARD/10H ERROR OF 7F15.5)
 919  FORMAT(9H ESTIMATE//)
 920  FORMAT(10H0PARAMETER7X,I4,6(11X,I4))
 921  FORMAT(1H0/1H019X,18HTABLE OF RESIDUALS/12H OBSERVATION5X,7HY VALU
     1E9X,10HY ESTIMATE8X,8HRESIDUAL)
 922  FORMAT(I7,2F18.5,F16.5)
 923  FORMAT(19H DEGREES OF FREEDOM5X,I6///)
 924  FORMAT(A6,33I2,/(6X,33I2))
 9876 FORMAT(20A4)
 1000 IF(NTAPE-5)1001,1002,1001
 1001 REWIND NTAPE
 1002 STOP
      END
C            SUBROUTINE ADJUST FOR BMD06V              MAY 10, 1966
      SUBROUTINE ADJUST
      DIMENSION DATA(100),FMT(180),B(61),C(61,61),CI(61,61),NV(61),
     1MV(61),NN(61),P(6),A(61,61),XM(61)
      COMMON  DATA   , FMT    , B      , C      , CI     , NV
      COMMON  MV     , NN     , P      , A      , PR     , RD
      COMMON  NP     , NO     , NT     , NEW    , NH     , NC
      COMMON  ND     , ON     , RE     , N1     , N2     , N3
      COMMON  N4     , NNT    , NPP    , MMT    , MD     , OON
      COMMON  NTAPE  , XM     , NCOV   , NT1    , NT2
C
      MMMT=MMT-NCOV
      DO 3 I=1,MMT
    3 MV(I)=I
      IF(NCOV) 40, 40, 10
 10   WRITE (6,901)
      I=0
      DO 5 J=1,MMT
      DO 1 K=1,NCOV
      IF(NN(K)-J)1,5,1
    1 CONTINUE
      I=I+1
      MV(I)=J
    5 CONTINUE
      DO 20 I=1,NCOV
      K=NN(I)
      DO 15 L=1,MMMT
      J=MV(L)
   15 DATA(L)=C(J,K)/C(J,J)
      WRITE (6,902)K
 20   WRITE (6,900)(DATA(J), J=1,MMMT)
 30   WRITE (6,903)
      GO TO 45
 40   WRITE (6,904)
   45 DO 50 L=1,MMMT
      J=MV(L)
   50 DATA(L)=C(J,NNT)/C(J,J)
      WRITE (6,900)(DATA(J), J=1,MMMT)
      IF(NCOV) 1000, 1000, 60
 60   DO 70 J=1,NNT
      DO 70 I=1,NNT
 70   C(I,J)=0.0
      DO 80 I=1,NCOV
      K=NN(I)
 80   XM(K)=XM(K)/ON
      DO 100 I=1,NO
      READ(NT1)(DATA(J),J=1,NNT),PR
      DO 85 J=1,NCOV
      K=NN(J)
 85   DATA(K)=DATA(K)-XM(K)
      DO 90 J=1,NNT
      DO 90 L=J,NNT
 90   C(J,L)=C(J,L)+DATA(J)*DATA(L)*PR
  100 WRITE(NT2)(DATA(J),J=1,NNT),PR
      REWIND NT1
      END FILE NT2
      REWIND NT2
      NTAPE = NT2
      N1=NNT-1
      DO 110 I=1,N1
      N2=I+1
      DO 110 J=N2,NNT
 110  C(J,I)=C(I,J)
 900  FORMAT(8F15.5)
 901  FORMAT(1H0//36H ORIGINAL MEANS FOR DESIGN VARIABLES)
 902  FORMAT(13H0VARIABLE NO.I3,13H  (COVARIATE))
 903  FORMAT(19H0DEPENDENT VARIABLE)
 904  FORMAT(1H0//11H CELL MEANS)
 9876 FORMAT(20A4)
 1000 RETURN
      END
C            SUBROUTINE INVERT FOR BMD06V              MAY 10, 1966
      SUBROUTINE INVERT (A,N,D,L,M)
      DIMENSION A(61,61),L(61),M(61)
C
      D=1.0
      DO80 K=1,N
      L(K)=K
      M(K)=K
      BIGA=A(K,K)
      DO20 I=K,N
      DO20 J=K,N
      IF(ABS(BIGA)-ABS(A(I,J))) 10,20,20
   10 BIGA=A(I,J)
      L(K)=I
      M(K)=J
   20 CONTINUE
      J=L(K)
      IF(L(K)-K) 35,35,25
   25 DO30 I=1,N
      HOLD=-A(K,I)
      A(K,I)=A(J,I)
   30 A(J,I)=HOLD
   35 I=M(K)
      IF(M(K)-K) 45,45,37
   37 DO40 J=1,N
      HOLD=-A(J,K)
      A(J,K)=A(J,I)
   40 A(J,I)=HOLD
   45 DO55 I=1,N
   46 IF(I-K)50,55,50
   50 A(I,K)=A(I,K)/(-A(K,K))
   55 CONTINUE
      DO65 I=1,N
      DO65 J=1,N
   56 IF(I-K) 57,65,57
   57 IF(J-K) 60,65,60
   60 A(I,J)=A(I,K)*A(K,J)+A(I,J)
   65 CONTINUE
      DO75 J=1,N
   68 IF(J-K)70,75,70
   70 A(K,J)=A(K,J)/A(K,K)
   75 CONTINUE
      D=D*A(K,K)
      A(K,K)=1.0/A(K,K)
   80 CONTINUE
      K=N
  100 K=(K-1)
      IF(K) 150,150,103
  103 I=L(K)
      IF(I-K) 120,120,105
  105 DO110 J=1,N
      HOLD=A(J,K)
      A(J,K)=-A(J,I)
  110 A(J,I)=HOLD
  120 J=M(K)
      IF(J-K) 100,100,125
  125 DO130 I=1,N
      HOLD=A(K,I)
      A(K,I)=-A(J,I)
  130 A(J,I)=HOLD
      GO TO 100
  150 RETURN
      END
C          SUBROUTINE READ FOR BMD06V                  MAY 10, 1966
      SUBROUTINE READ
      DIMENSION DATA(100),FMT(180),B(61),C(61,61),CI(61,61),NV(61),
     1MV(61),NN(61),P(6),A(61,61),XM(61)
      DOUBLE PRECISION GOBLE3,RRD
      COMMON  DATA   , FMT    , B      , C      , CI     , NV
      COMMON  MV     , NN     , P      , A      , PR     , RD
      COMMON  NP     , NO     , NT     , NEW    , NH     , NC
      COMMON  ND     , ON     , RE     , N1     , N2     , N3
      COMMON  N4     , NNT    , NPP    , MMT    , MD     , OON
      COMMON  NTAPE  , XM     , NCOV   , NT1    , NT2
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
C
      NEW=0
      IF(NT) 20, 20, 8
 8    OON=ON+1.0
      WRITE (6,900)
      WRITE (6,901)
      DATA GOBLE3/6HTRNGEN/
      DO 18 I=1,NT
      READ(5,905)RRD,(A(I,J), J=1,4)
      IF(RRD.EQ. GOBLE3) GO TO 12
 10   WRITE (6,906)
   11 STOP
 12   N1=A(I,1)
      N2=A(I,2)
      N3=A(I,3)
      N4=A(I,4)
      IF(N2-8) 13, 14, 14
 13   WRITE (6,902)I,N1,N2,N3
      GO TO 18
 14   IF(N2-10) 15, 16, 16
 15   WRITE (6,903)I,N1,N2,N3,A(I,4)
      GO TO 18
 16   WRITE (6,902)I,N1,N2,N3,N4
 18   CONTINUE
   20 IF(ND.GT.0.AND.ND.LE.10)GO TO 25
      ND=1
      WRITE(6,4000)
   25 ND=ND*18
      WRITE(6,5000)
5000  FORMAT(24H0VARIABLE FORMAT CARD(S))
      READ (5,912)(FMT(I), I=1,ND)
      WRITE(6,6000)(FMT(I),I=1,ND)
6000  FORMAT(1X,18A4)
      DO 70 J=1,NNT
      DO 70 I=1,NNT
 70   C(I,J)=0.0
      ITAPE=5
      IF(NTAPE.GT.0)ITAPE=NTAPE
      DO 280 I=1,NO
 75   READ (ITAPE,FMT)(DATA(J), J=1,NPP)
 77   IF(NC) 81, 81, 82
 81   PR=1.0
      GO TO 84
 82   PR=DATA(NPP)
      IF(PR) 84, 81, 84
 84   IF(NT) 260, 260, 90
 90   DO 255 J=1,NT
      N1=A(J,1)
      N2=A(J,2)
      N3=A(J,3)
      D3=DATA(N3)
      N4=A(J,4)
      GO TO (110,110,110,140,150,160,170,180,190,200,210,220,230,240),N2
  110 IF(D3)113,111,115
 111  IF(N2-2) 114, 112, 114
  112 D1=1.0
      GO TO 250
 113  WRITE (6,904)N2,N3,I,D3
      NEW=NEW+1
  114 D1=0.0
      GO TO 250
 115  GO TO (118,120,130), N2
  118 D1=SQRT(D3)
      GO TO 250
  120 D1=SQRT(D3)+SQRT(D3+1.0)
      GO TO 250
  130 D1=ALOG10(D3)
      GO TO 250
  140 D1=EXP(D3)
      GO TO 250
  150 IF(-D3)154,154,113
  154 RD=SQRT(D3)
      IF(RD-1.0) 155, 156, 113
  155 D1=ASN(RD)
      GO TO 250
  156 D1=1.5708
      GO TO 250
  160 IF(-D3)161,161,113
  161 RD=SQRT(D3/OON)
      RE=SQRT((D3+1.0)/OON)
      IF(RD-1.0) 162, 163, 113
 162  RD=ASN(RD)
      GO TO 164
 163  RD=1.5708
 164  IF(RE-1.0) 165, 166, 113
 165  RE=ASN(RE)
      GO TO 167
 166  RE=1.5708
  167 D1=RD+RE
      GO TO 250
  170 IF(D3)175,113,175
  175 D1=1.0/D3
      GO TO 250
  180 D1=D3+A(J,4)
      GO TO 250
  190 D1=D3*A(J,4)
      GO TO 250
  200 D1=D3**N4
      GO TO 250
  210 D1=D3+DATA(N4)
      GO TO 250
  220 D1=D3-DATA(N4)
      GO TO 250
  230 D1=D3*DATA(N4)
      GO TO 250
 240  IF(DATA(N4)) 242, 241, 242
 241  WRITE (6,904)N2,N4,I,DATA(N4)
      NEW=NEW+1
      GO TO 255
  242 D1=D3/DATA(N4)
  250 DATA(N1)=D1
  255 CONTINUE
 260  DO 270 J=1,NNT
      XM(J)=XM(J)+DATA(J)
      DO 270 L=J,NNT
 270  C(J,L)=C(J,L)+DATA(J)*DATA(L)*PR
  280 WRITE(NT1)(DATA(J),J=1,NNT),PR
      END FILE NT1
      REWIND NT1
      IF(NEW) 285, 285, 295
 285  N1=NNT-1
      DO 290 I=1,N1
      N2=I+1
      DO 290 J=N2,NNT
 290  C(J,I)=C(I,J)
 295  IF(MD) 1000, 1000, 296
 296  WRITE (6,908)
      IF(NC) 297, 297, 298
 297  N3=NNT
      GO TO 299
 298  N3=NNT+1
 299  DO 350 I=1,NO
      READ(NT1)(DATA(J),J=1,NNT),PR
      IF(NC) 305, 305, 300
 300  DATA(N3)=PR
 305  N1=1
      N2=0
 310  IF((N3-N2)-8) 315, 315, 320
 315  N2=N3
      GO TO 325
 320  N2=N2+8
 325  IF(N1-1) 330, 330, 335
 330  WRITE (6,909)I,(DATA(J), J=N1,N2)
      GO TO 340
 335  WRITE (6,910)(DATA(J), J=N1,N2)
 340  IF(N3-N2) 350, 350, 345
 345  N1=N1+8
      GO TO 310
 350  CONTINUE
      REWIND NT1
 900  FORMAT(15H0TRANSFORMATION4X,12HVARIABLE NO.6X,7HTYPE OF9X,8HVARIAB
     1LE6X,12H2ND VARIABLE)
 901  FORMAT(12H    CARD NO.9X,8HASSIGNED5X,14HTRANSFORMATION4X,11HTRANS
     1FORMED4X,12HOR  CONSTANT)
 902  FORMAT(1H 4X,I4,13X,I4,13X,I3,13X,I3,10X,I6)
 903  FORMAT(1H 4X,I4,13X,I4,13X,I3,13X,I3,10X,F12.5)
 904  FORMAT(1H0//27H TRANSFORMATION OF THE TYPEI3,24H  CANNOT BE PERFOR
     1MED ON/9H VARIABLEI3,6H  CASEI5,15H.  THE VALUE ISF12.5)
 905  FORMAT(A6,F3.0,F2.0,F3.0,F6.0)
 906  FORMAT(29H0ERROR ON TRANSFORMATION CARD)
 907  FORMAT(12I6)
 908  FORMAT(1H0//1H 10X,34HINPUT  DATA  AS  READ  FROM  CARDS/9H CASE N
     1O.)
 909  FORMAT(I6,8F13.5)
 910  FORMAT(6X,8F13.5)
  912 FORMAT (18A4)
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
 9876 FORMAT(20A4)
 1000 NTAPE=NT1
      RETURN
      END
C        SUBROUTINE TPWD FOR BMD06V                    MAY 10, 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 18
 17   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