Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-05 - decus/20-0137/bmd/bmd03s.for
There is 1 other file named bmd03s.for in the archive. Click here to see a list.
C        BIOLOGICAL ASSAY - PROBIT ANALYSIS            JUNE 14, 1966
C        THIS IS A SIFTED VERSION OF BMD03S ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
C     EXTERNAL DIFF
	DOUBLE PRECISION AAA, CODE,Q000HL,Q001HL
      DIMENSION DOSE(1001), R(1001), SN(1001), ICODE(64), CONST(64)
     1, SMP(1001), P(1001), X(1001), SMY(1001)
     2,   A(3),  DC(9),   F(3)
     3,  CD(10),  DV(10),  D(10)
     4, TEMP(1000),  RX(1000)
     5,DC0(9),CD0(9)
      COMMON  CODE   , R      , SN     , ICODE  , DOSE   , CONST
      COMMON  CG     , M      , K      , IPR    , SMP    , P
      COMMON  SMY    , X      , A      , DC     , F
      ER = 0.0
	CALL USAGEB('BMD03S')
C     READ IN DATA
 998  CONTINUE
      CALL READIN
      SAVC = A(3)
      DO 2099 JK=1,K
      A(3) = SAVC
      IK = 1+8*(JK-1)
C  TRANSFORM DOSE.
      WRITE (6,9998)JK,K
      CALL  TRNGEN( DOSE,  CONST(IK),  M,  ICODE(IK),  1,  IERROR,  X)
      IF(IERROR)2098,  70,2098
 70   CONTINUE
      WRITE (6,9993)(X(I),I=1,M)
 73   CONTINUE
      IF(A(2))2000,2000,2001
C     GET LEAST SQUARES FIT
 2000 CONTINUE
      I1 = 1
      DO 2210 I = 1,M
      IF(1.0-SMP(I)) 2210,2210,2206
 2206 IF( SMP(I) -0.0) 2210,2210,2208
 2208 TEMP(I1) = SMY(I)
      RX(I1) = X(I)
      M2 = I1
      I1 = I1 +1
 2210 CONTINUE
      CALL LINREG (  RX,  TEMP,  M2,  A(2),  A(1))
 2001 CONTINUE
      CALL MNVAR (  X,  M,  XAV,  XVAR)
      WRITE (6,9988)A(1),A(2),A(3)
      IF (A(2))54,54,2201
 54   WRITE (6,46)
      GO TO 2098
 2201 CONTINUE
      DO 7002 I = 1,M
 7002 X(I) = X(I) - XAV
      A(1) = A(1) + A(2) *XAV
C  SOLVE FOR MAXIMUM LIKLIHOOD BY NEWTONS METHOD
      EPS = .00001
      NPR = IPR
C FIRST TIME GET BEST APPROX TO ALPHA, BETA
      CG1 = CG
      CG = 1.0
      DO 2376 NK = 1,2
      IF  (CG)  204, 2041, 2041
  204 NV = 3
      GO TO  2042
 2041 NV = 2
 2042 CONTINUE
      NIT = 50
      CALL NEWT
     1(DC,F,NV,NIT,TEMP,EPS,ER,NPR,A)
      CG = CG1
 2444 IF(ER)2098,2043,2044
 2044 CONTINUE
 2043 CONTINUE
      IF (CG1)2376,2377,2377
 2376 CONTINUE
 2377 CONTINUE
      CG = -1.0
      CALL DIFF
      CG = CG1
      HAMU = (5.0- A(1))/A(2)
      ALPHA = A(1)-A(2)*XAV
      ALPHA1 = A(1)
      A(1) = A(1) - A(2) *XAV
      HAMU0 = (5.0-A(1))/A(2)
      HASI = 1.0/A(2)
C  COMPUTR COVARIANCE MATRICES
      IF (CG1) 305, 207, 207
C  TWO DIMENSION CASE
 207  CONTINUE
      N1 = 2
      NN = 4
      DATA Q000HL/6H.     /
      AAA=Q000HL
      IDF = M-2
      DC(4) = DC(5)
      DC(3) = DC(2)
      DV(1) = -DC(1)/HASI**2
      DV(4) = -(HAMU*(HAMU*DC(1) + 2.0*DC(2)) + DC(5))/HASI**4 + 2.0*(HA
     1MU*F(1)+F(2))/HASI**3
      DV(2) = (HAMU*DC(1) +DC(2))/HASI**3
      DV(3) = DV(2)
      DC0(1) =-DC(1)
      DC0(2) =-DC(1)*XAV-DC(2)
      DC0(4) = -DC(4)-2.0*DC(2)*XAV -DC(1)*XAV**2
      DC0(3) = DC0(2)
      GO TO 306
C  THREE DIMENSION CASE
 305  CONTINUE
      N1 = 3
      NN = 9
      DATA Q001HL/6H, C.  /
      AAA=Q001HL
      IDF = M-3
      IF(SN(M+1)) 3051,3051,3052
 3052 CONTINUE
      IDF = IDF + 1
 3051 CONTINUE
      DV(1) = -DC(1)/HASI**2
      DV(5) = -(HAMU*(HAMU*DC(1)+2.0*DC(2))+DC(5))/HASI**4 + 2.0*(HAMU*F
     1(1)+F(2))/HASI**3
      DV(9) = -DC(9)
      DV(2) = (HAMU*DC(1)+DC(2))/HASI**3
      DV(4) = DV(2)
      DV(3) = -DC(3)/HASI
      DV(7) = DV(3)
      DV(6) = (HAMU*DC(3) +DC(6))/HASI**2
      DV(8) = DV(6)
      DC0(1) = -DC(1)
      DC0(2) = - DC(1)*XAV-DC(4)
      DC0(3) = -DC(3)
      DC0(5) = -DC(1)*XAV**2 -2.0*DC(4)*XAV -DC(5)
      DC0(6) = -DC(3)*XAV-DC(6)
      DC0(9) = -DC(9)
      DC0(4) = DC0(2)
      DC0(7) = DC0(3)
      DC0(8) = DC0(6)
 306  CONTINUE
      DO 3061 I=1,NN
 3061 DC(I) = -DC(I)
      CALL INVERT( DC, CD,  N1,   DET,TEMP(20) )
      CALL INVERT( DV,   D,  N1,  DET,TEMP(20))
      CALL INVERT (DC0,CD0,N1,DET,TEMP(20))
C  COMPUTE CHI SQUARE
      M1 = M
      IF (CG1)308,3081,3081
 308  CONTINUE
      M1 = M + 1
 3081 CONTINUE
      CHI = 0.0
      DO 307 I = 1,M
      RI = R(I)
      PP = A(3) + P(I)*(1.0-A(3))
      CHI = CHI + (RI-SN(I)*PP)**2 / (PP-PP**2) /SN(I)
 307  CONTINUE
C  PRINT RESULTS
      WRITE (6,9885)
      WRITE (6,9980)ALPHA1, F(1), ALPHA, A(2), F(2),A(2), A(3), F(3), A(
     13), HAMU0, HASI
      WRITE (6,9981)AAA ,AAA
      DO 3071  I= 1,N1
      I1 = 1 + N1*(I-1)
      I3 = N1*I
      GO TO (3071,3170,3171),N1
 3170 WRITE (6,9883)(CD(IJ),IJ=I1,I3),(CD0(IJ),IJ=I1,I3)
      GO TO 3071
 3171 WRITE (6,9884)(CD(IJ),IJ=I1,I3),(CD0(IJ),IJ=I1,I3)
 3071 CONTINUE
      WRITE (6,9982)AAA
      DO 3072  I= 1,N1
      I1 = 1 + N1*(I-1)
      I3 = N1*I
 3072 WRITE (6,9983)( D(IJ),IJ=I1,I3)
      WRITE (6,9984)CHI,IDF
 2098 CONTINUE
      WRITE (6,9989)JK,K
      A(2) = 0.0
 2099 CONTINUE
      WRITE (6,38)CODE
      GO TO 998
C NORMAL OUTPUT FORMAT
 9883 FORMAT(3X, 2F12.4, 35X, 1H*, 2F12.4)
 9884 FORMAT(3X,3F12.4,23X, 1H*, 3F12.4 )
 9885 FORMAT (67H0THE FINAL VALUES OF THE PARAMETERS AND THE PARTIAL DER
     1IVATIVES ARE   /
     2 5X,40HFOR PROBIT LINE (ALPHA1+BETA*(X-XBAR))  17X,1H*,5X,
     3 35HFOR PROBIT LINE (ALPHA+BETA*(X))       )
 9980 FORMAT(
     15X, 7HALPHA1= F8.3, 10X,6HDL/DA=3X,E10.4 ,13X,1H*,5X,7HALPHA =F8.3
     2/5X,7HBETA  =F8.3, 10X,6HDL/DB= 3X,E10.4, 13X, 1H*, 5X,7HBETA  =F8
     3.3/  5X, 1HC 5X,1H=, F8.3, 10X, 6HDL/DC= 3X, E10.4,13X,1H*, 5X,
     41HC,5X, 1H= F8.3   /
     5  120H0THE ESTIMATES MU-HAT AND SIGMA-HAT OF THE MEAN AND STANDARD
     6 DEVIATION OF THE TRANSFORMED DOSE TOLERANCE ARE AS FOLLOWS.     /
     7 5X, 10H   MU-HAT= 3X, F9.4 /
     8  5X,  10HSIGMA-HAT= 3X, F9.4    )
 9981 FORMAT(39H0THE COVARIANCE MATRICES ARE AS FOLLOWS   /
     1,  5X, 16HFOR ALPHA1, BETA  A6,35X, 1H*,5X,15HFOR ALPHA, BETA  A6)
 9982 FORMAT ( 5X,  13HFOR MU, SIGMA       A6     )
 9983 FORMAT (3X, 9F12.4)
 9984 FORMAT (12H0CHI-SQUARE=3X, F12.4, 6H, WITHI5, 20H DEGREES OF FREED
     1OM.  )
 9985 FORMAT (50H THE MEAN AND VARIANCE OF THE TRANSFORMED DOSE ARE /
     1  7H MEAN=  F8.4,  7H  VAR=   F8.4)
 9988 FORMAT (54H0THE INITIAL VALUES AT START OF ITERATION ARE, ALPHA=
     1  F9.3,  3X,5HBETA= F9.3,  3X,  2HC= F9.3)
 9993 FORMAT (3X,9F12.2 )
 9998 FORMAT (///
     1       45H0TRANSFORMED DOSE FOR TRANSFORMATION CARD NO  I3,3H OF
     2 ,I3)
 9989 FORMAT (46H END FOR OUTPUT USING TRANSFORMATION CARD NO  I3, 3H OF
     1,  I3  //)
   38 FORMAT(1H05X,26HEND OF OUTPUT FOR PROBLEM  A6 )
C  ERROR OUTPUT FORMAT
 46   FORMAT (48H0THE DATA GIVEN IS UNWORKABLE, BETA IS NEGATIVE.    )
      END
C        SUBROUTINE DERIV2 FOR BMD03S                  JUNE 14, 1966
C DERIVATIVES WITH RESPECT TO A AND B.
      SUBROUTINE DERIV2 ( X,  YY,  Z,  P,  R,  SN,  M,  C,  F,  D)
      DIMENSION X(1),Z(1),YY(1),P(1),SN(1)
     1,  F(1),  D(1),  R(1)
      CQ = 1.0-C
      DO 10 I=1,4
      D(I) = 0.0
 10   CONTINUE
      F(1) = 0.0
      F(2) = 0.0
 12   CONTINUE
      PCC = 0.0
      DO 100I=1,M
      PP = C+P(I)*CQ
      Q = 1.0- PP
      PQ = PP*Q
      XLP = (R(I) - PP*SN(I)) /PQ
      XLPP = -R(I)/PP**2 -(SN(I)-R(I))/Q**2
      PA = CQ*Z(I)
      PB = PA*X(I)
      PC = 1.0-P(I)
      PAA = -PA*YY(I)
      PAB = -PB*YY(I)
      PBB = PAB*X(I)
      PAC = -Z(I)
      PBC = PAC*X(I)
      F(1) = F(1) + XLP*PA
      F(2) = F(2) + XLP*PB
      D(1) = D(1) + XLPP*PA*PA + XLP*PAA
      D(2) = D(2) +XLPP*PA*PB + XLP*PAB
      D(4) = D(4) + XLPP*PB*PB + XLP*PBB
  100   CONTINUE
      D(3) = D(2)
       RETURN
       END
C        SUBROUTINE DERIV3 FOR BMD03S                  JUNE 14, 1966
C DERIVATIVES WITH RESPECT TO A, B, AND C.
      SUBROUTINE DERIV3 ( X,  YY,  Z,  P,  R,  SN,  M,  C,  F,  D)
      DIMENSION X(1),Z(1),YY(1),P(1),SN(1)
     1,  F(1),  D(1),  R(1)
      CQ = 1.0-C
      DO 10 I=1,3
      F(I) = 0.0
      D(I) = 0.0
 10   CONTINUE
      DO 12 I=4,9
      D(I) = 0.0
 12   CONTINUE
      PCC = 0.0
      M1 = M
      IF (C) 14,13,14
 14   M1 = M+1
 13   CONTINUE
      DO 100I=1,M1
      PP = C+P(I)*CQ
      Q = 1.0- PP
      PQ = PP*Q
      XLP = (R(I) - PP*SN(I)) /PQ
      XLPP = -R(I)/PP**2 -(SN(I)-R(I))/Q**2
      PA = CQ*Z(I)
      PB = PA*X(I)
      PC = 1.0-P(I)
      PAA = -PA*YY(I)
      PAB = -PB*YY(I)
      PBB = PAB*X(I)
      PAC = -Z(I)
      PBC = PAC*X(I)
      F(1) = F(1) + XLP*PA
      F(2) = F(2) + XLP*PB
      F(3) = F(3) + XLP*PC
      D(1) = D(1) + XLPP*PA*PA + XLP*PAA
      D(2) = D(2) +XLPP*PA*PB + XLP*PAB
      D(5) = D(5) + XLPP*PB*PB + XLP*PBB
      D(3) = D(3) + XLPP*PA*PC + XLP*PAC
      D(6) = D(6) + XLPP*PB*PC + XLP*PBC
      D(9) = D(9) + XLPP*PC*PC + XLP*PCC
  100   CONTINUE
      D(4) = D(2)
      D(8) = D(6)
      D(7) = D(3)
 900  RETURN
       END
C        SUBROUTINE DIFF FOR BMD03S                    JUNE 14, 1966
C  DIFFERENTIAL COEFICIENTS AND FUNCTION AUXILIARY SUBROUTINE FOR NEWT
      SUBROUTINE DIFF
	DOUBLE PRECISION CODE
      DIMENSION DOSE(1001), R(1001), SN(1001), ICODE(64), CONST(64)
     1, SMP(1001), P(1001), Z(1001), X(1001), YY(1001),SMY(1001),Y(1001)
     2,   A( 3),   DC(9),   F(3)
      COMMON  CODE   , R      , SN     , ICODE  , DOSE   , CONST
      COMMON  CG     , M      , K      , IPR    , SMP    , P
      COMMON  SMY    , X      , A      , DC     , F
C     PREDICT PROBIT YY AND PROBABILITY P.
         DO 200 I=1,M
      Y(I) = A(1) + A(2)*X(I)
       YY(I)=Y(I) -5.0
      CALL PHI2(YY(I),  P(I),  Z(I))
 200   CONTINUE
      M1 = M+1
      P(M1) = 0.0
      YY(M1)=0.0
      X(M1) = 0.0
      Z(M1)=0.0
 2049 IF (CG)300,2050,2050
C     C CONSTANT
 2050 CONTINUE
      CALL       DERIV2 ( X,  YY,  Z,  P,  R,  SN,  M,  A(3),  F,  DC)
      RETURN
C     OPTIMIZE C
 300  CONTINUE
      CALL       DERIV3 ( X,  YY,  Z,  P,  R,  SN,  M,  A(3),  F,  DC)
 900  RETURN
      END
C        SUBROUTINE INVERT FOR BMD03S                  JUNE 14, 1966
      SUBROUTINE INVERT (  G,  R,  N,  DET,  ER)
      DIMENSION G(1), R(1)
      ER = 0.0
 1    K = N
      GO TO (10,20,30),K
 10   R(1)= 1.0/G(1)
      RETURN
   20 CONTINUE
      DET = G(1)*G(4) - G(3)*G(2)
      IF (DET)80,99,80
 99   ER = 1.0
 199  RETURN
 80   CONTINUE
      R(1) = G(4)/DET
      R(2) = -G(2)/DET
      R(3) = -G(3)/DET
      R(4) = G(1)/DET
      RETURN
   30 CONTINUE
      R(1) =  (G(5)*G(9) - G(6)*G(8))
      R(4) = -(G(4)*G(9) - G(6)*G(7))
      R(7) =  (G(4)*G(8) - G(5)*G(7))
      R(2) = -(G(2)*G(9) - G(3)*G(8))
      R(5) =  (G(1)*G(9) - G(3)*G(7))
      R(8) = -(G(1)*G(8) - G(2)*G(7))
      R(3) =  (G(2)*G(6) - G(3)*G(5))
      R(6) = -(G(1)*G(6) - G(3)*G(4))
      R(9) =  (G(1)*G(5) - G(2)*G(4))
      DET = R(1)*G(1) + R(4)*G(2) + R(7)*G(3)
 35   IF(DET)84,99,84
 84   CONTINUE
      DO  32  I = 1, 9
   32 R(I) = R(I)/DET
 88   RETURN
 9999 FORMAT (1H  E12.4   )
      END
C        SUBROUTINE LINREG FOR BMD03S                  JUNE 14, 1966
C  REGRESS Y AS A FUNCTION OF X.  Y=AX + B.
      SUBROUTINE LINREG (  X,  Y,  M,  A,  B)
      DIMENSION  X(1),   Y(1)
      XM =M
      SX = 0.0
      SY = 0.0
      SXX= 0.0
      SXY= 0.0
      DO 10 I=1,M
      SXY = SXY +X(I) * Y(I)
      SX = SX+ X(I)
      SY=  SY+Y(I)
      SXX = SXX+ X(I)**2
 10   CONTINUE
      A = (XM*SXY -SX*SY)/(XM*SXX-SX**2)
      B = (SY -A*SX)/XM
 1000 RETURN
      END
C        SUBROUTINE MNVAR FOR BMD03S                   JUNE 14, 1966
C        COMPUTES MEAN AND VARIANCE
      SUBROUTINE MNVAR (  X,  N,  XB,  XV)
      DIMENSION X(1)
      XN1 = N-1
      XN = N
      XV=0.0
      XB=0.0
      DO 10 I=1,N
      XB = XB+X(I)
 10   CONTINUE
      XB = XB/XN
      DO 20 I=1,N
      XV = XV + (X(I) - XB)**2
 20   CONTINUE
      XV = XV/XN1
      RETURN
      END
C        SUBROUTINE NEWT FOR BMD03S                    JUNE 14, 1966
      SUBROUTINE NEWT
     1(D,F,N,NI,TEMP,EPS,ER,NPR,A)
      DIMENSION D(1),F(1),TEMP(1)
     1, A(1)
      ER=0.0
      NM = N+1
      NF = N
      NN = N**2
      IPR = NI+1
      IF (NPR)8,8,5
 5    IPR = 1
 8    CONTINUE
      DO 200 K =1,NI
      CALL DIFF
      CALL INVERT( D, TEMP(NM), NF,DET, TEMP(20) )
 9    DEL = 0.0
      IF(DET)10,99,10
 10   CONTINUE
      DO  35  I = 1, N
      TEMP(I) = 0.0
      IJ = I
      DO  30  J = 1, N
      IJ = IJ+N
      TEMP(I) = TEMP(I) - F(J)*TEMP(IJ)
   30 CONTINUE
      DEL = DEL +  TEMP(I)**2
      A(I) = A(I) + TEMP(I)
 35   CONTINUE
      IF(K-IPR)18,15,18
 15   CONTINUE
      IPR = IPR + NPR
      WRITE (6,800)K
      WRITE (6,805)(F(J), J=1,N)
      WRITE (6,801)
      WRITE (6,805)(A(J), J=1,N)
      WRITE (6,802)
      DO 16J=1,N
      JL = N*(J-1)+1
      JU = JL + N-1
      WRITE (6,805)(D(IJ),IJ=JL,JU)
 16   CONTINUE
 18   CONTINUE
      IF  (DEL - EPS)  80, 80, 190
  190 CONTINUE
  200 CONTINUE
      WRITE (6,806)NI
      ER = 1.0
      GO TO 80
 99   CONTINUE
      WRITE (6,807)
      ER = -1.0
      NI = K
   80 RETURN
 800  FORMAT (16H0ITERATION COUNT I5, 20H. VALUE OF FUNCTION.  )
 801  FORMAT ( 27H APPROXIMATION TO SOLUTION.)
 805  FORMAT ( 3H   , 7E16.6)
 806  FORMAT ( 43H0CRITERION FOR CONVERGENCE IS NOT MET AFTER I5,
     1  11H ITERATIONS  )
 802  FORMAT (34H DIFFERENTIAL COEFFICIENT MATRIX.  )
 807  FORMAT (78H0MATRIX OF DIFFERENTIAL COEFICIENTS IS SINGULAR, PROBLE
     1M CANNOT BE CONTINUED.     )
      END
C        FUNCTION PHINV FOR BMD03S                     JUNE 14, 1966
C        USES HASTINGS APPROXIMATIONS
      FUNCTION PHINV (P)
      IF(P-.5)10,100,12
 100  PHINV = 0.0
      RETURN
 12   CONTINUE
      Q = 1.0-P
      T = -1.0
      GO TO 13
 10   Q = P
      T = 1.0
 13   CONTINUE
      IF(Q)121,121,122
 121  CONTINUE
      PHINV = -T*999.0
      RETURN
 122  CONTINUE
      E=SQRT(ALOG(1.0/Q**2))
      X = -E+ ((.010328*E+.802853)*E +2.515517) /(((.001308*E+.189269)*E
     1 + 1.432788)*E+1.0)
      PHINV=T*X
      RETURN
      END
C        SUBROUTINE PHI2 FOR BMD03S                    JUNE 14, 1966
C        COMPUTES PHI(X) AND PHI-PRIME(X)
      SUBROUTINE PHI2( X, PHI, PP)
      SQ2 = .141421356E+01
      Y = ABS(X)/SQ2
      YY = X*X/2.0
      PP = EXP(-YY)*.39894228
      P1 = PP*2.0*SQ2
      ET = 1.0/ (1.0+.3275911*Y)
      PH = 1.0-((((.94064607*ET-.128782245E+01)*ET+1.25969513)*ET-.25212
     18668)*ET+.225836846)*ET*P1
      IF(X) 1, 4, 3
 1    PH = -PH
 3    PHI = .5 + PH/2.0
      RETURN
 4    PHI = 0.0
      RETURN
      END
C        SUBROUTINE READIN FOR BMD03S                  JUNE 14, 1966
      SUBROUTINE READIN
      DOUBLE PRECISION PPP,AAA,TRG,AA,PRBLM,TODE,CODE,FMT
     1,Q003HL,Q004HL,Q005HL,Q006HL,Q007HL,Q008HL
      DIMENSION DOSE(1001), R(1001), SN(1001), ICODE(64), CONST(64)
     1, SMP(1001), P(1001), X(1001), SMY(1001)
     2,   A(3),  DC(9),   F(3)
     3,  AA(3),  FMT(120)
      COMMON  CODE   , R      , SN     , ICODE  , DOSE   , CONST
      COMMON  CG     , M      , K      , IBG    , SMP    , P
      COMMON  SMY    , X      , A      , DC     , F
C
 11   FORMAT(53H1BMD03S BIOLOGICAL ASSAY, PROBIT ANALYSIS - REVISED  ,2X
     1,'FEBRUARY 11, 1969'/
     242H HEALTH SCIENCES COMPUTING FACILITY, UCLA.//)
C
      DATA PPP,AAA,TRG/6HPROBLM,6HFINISH,6HSPECTG/
      DATA Q003HL/6HOPTIMI/
      DATA Q004HL/6HZE C. /
      DATA Q005HL/6H      /
      DATA Q006HL/6HC HELD/
      DATA Q007HL/6H CONST/
      DATA Q008HL/6HANT.  /
 1    NTAP=5
 2    ER = 0.0
      READ (5,10)PRBLM, CODE, IC, CG, M, R(1), SN(1), A(1), A(2), IBG, K
     1, IPR , MTAP, IVF
      IF(PRBLM.NE.PPP) GO TO 990
 49   CONTINUE
      IF(M*(M-1001))495,993,993
  495 M1=M+1
      R(M1) = R(1)
      SN(M1) = AMAX1(SN(1),0.0)
 50   CONTINUE
      K=MIN0(K,8)
      IF(K)  991,991,51
 51   DO 501 I=1,K
      I2 = 1+8*(I-1)
      I3 = I2+7
      READ(5,19)TODE,K1,(ICODE(I1),CONST(I1),I1=I2,I3)
      IF(TODE.NE.TRG)GO TO 994
 501  CONTINUE
      IF(-K)502,887,887
  502 IF(IVF.GT.0.AND.IVF.LE.10)GO TO 52
      IVF=1
      WRITE (6,4000)
   52 IVF=IVF*12
      READ (5,47)(FMT(I),I=1, IVF )
      CALL TPWD(MTAP,NTAP)
 53   READ (NTAP,FMT)(DOSE(I),I=1,M)
      READ (NTAP,FMT)(R(I), I=1,M)
      READ (NTAP,FMT)(SN(I), I =1,M)
C     GET SAMPLE P ( = R/S )
C     GET INVERSE NORMAL OF SAMPLE P.
      DO 60 I=1,M
      IF(SN(I)) 54,54,540
 540  CONTINUE
      SMP(I)=R(I)/SN(I)
      IF( SMP(I)*(SMP(I)-1.0))58,58,54
 58   CONTINUE
   60 CONTINUE
C  CHECK OPTION FOR C
      IF(ABS(CG)-1.0) 261,28,28
 261  IF(CG) 29,262,262
 262  IF(IC)263,30,263
 263  IF(CG) 29,28,29
C  WE GUESS C
 28   IF (SN(M1)) 281,281,282
 281  CONTINUE
      SN(M1) = 0.0
      R(M1) = 0.0
      SMALN = SN(1)
      SMALP = SMP(1)
      GO TO 2821
 282  CONTINUE
      SMALP = R(M1) /SN(M1)
      SMALN = SN(M1)
 2821 CONTINUE
      DO 301 I=1,M
      IF(SMP(I) - SMALP) 300,301,301
 300  SMALN = SN(I)
      SMALP = SMP(I)
 301  CONTINUE
      A(3) = SMALP
      IF(A(3)) 54,302,303
 302  A(3) = .5/SMALN
 303  CONTINUE
      GO TO 291
C  USER GUESS C
 29   A(3) = ABS(CG)
 291  CONTINUE
      AA(1)=Q003HL
      AA(2)=Q004HL
      AA(3)=Q005HL
      CG = -1.0
      GO TO 32
C  HOLD C CONSTANT
   30 AA(1)=Q006HL
      AA(2)=Q007HL
      AA(3)=Q008HL
      A(3) = CG
 32   CONTINUE
      IF (SN(M1))611,611,610
 611  SMP(M1) = 0.0
      SMY(M1) = 5.0+ PHINV(0.0)
      GO TO 612
 610  CONTINUE
      SMP(M1) = R(M1)/SN(M1)
      SMY(M1) = 5.0+ PHINV(SMP(M1))
 612  CONTINUE
      QC = 1.0- A(3)
      DO 6110 I = 1,M
      SMP(I) = ( SMP(I) - A(3))/QC
      SMY(I) = 5.0+ PHINV(SMP(I))
 6110 CONTINUE
C  PRINT INPUT DATA
 210  CONTINUE
      WRITE (6,11)
      WRITE (6,9990)CODE
      WRITE (6,9991)SN(M1), R(M1), AA
      IF (IPR)25,21,25
 21   CONTINUE
      WRITE (6,9992)
      WRITE (6,9993)(DOSE(I) ,I=1,M)
      WRITE (6,9994)
      WRITE (6,9993)(R(I),I=1,M)
      WRITE (6,9996)
      WRITE (6,9993)(SN(I), I=1,M)
      IF(ER)41,25,41
 41   WRITE (6,9985)
      WRITE (6,38)CODE
      GO TO 1
 25   CONTINUE
      WRITE (6,9981)
      WRITE (6,9993)(SMY(I),I= 1,M)
      WRITE (6,9982)
      WRITE (6,9993)(SMP(I),I= 1,M)
      RETURN
 990  CONTINUE
      IF(PRBLM.EQ.AAA)GO TO 999
      GO TO 992
C  DATA ERROR EXIT
 54   CONTINUE
 995  ER = 1.0
      GO TO 210
C LAST PROBLEM COMPLETED.
 999  CONTINUE
      WRITE (6,9995)
  887 IF(NTAP.EQ.5)GO TO 889
      REWIND NTAP
  889 STOP
C PROBLEM CARD ERROR
 992  CONTINUE
      WRITE (6,11)
      WRITE (6,9997)
      GO TO 887
 991  CONTINUE
      WRITE (6,11)
      WRITE (6,9998)CODE,K
      GO TO 887
  993 WRITE (6,11)
      WRITE (6,9999)
      GO TO 887
  994 WRITE (6,11)
      WRITE (6,9989)KODE
      K=-1001
      GO TO 501
 10   FORMAT( A6, A6, I1, F5.5, I4, 4F6.0, 18X, I2,2I1,2I2)
 19   FORMAT (  A6,  I1,  8( I2,F6.0))
 47   FORMAT (12A6)
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
 9989 FORMAT(48H0CONTROL CARD ERROR. SPECTG CARD EXPECTED BUT A A6,35H C
     1ARD WAS FOUND. EXECUTION DELETED.)
 9990 FORMAT (8H PROBLEM 1X,A6,  1H. )
 9991 FORMAT (35H CONTROL GROUP DATA. SAMPLE SIZE N= F6.0, 13H, RESPONSE
     1 R= ,  F6.0,  1H.
     2  3X,  3A6)
 9992 FORMAT (20H UNTRANSFORMED DOSE. )
 9993 FORMAT (3X,9F12.2  )
 9994 FORMAT (10H RESPONSE.  )
   38 FORMAT(1H05X,26HEND OF OUTPUT FOR PROBLEM  A6 )
 9981 FORMAT (14H SAMPLE PROBIT  )
 9982 FORMAT (19H SAMPLE PROBABILITY )
 9985 FORMAT (33H0DATA INCORRECT,PROBLEM DELETED.  )
 9995 FORMAT (52H0NO NEW PROBLEMS, PROGRAM TERMINATED BY FINISH CARD.  )
 9996 FORMAT (13H SAMPLE SIZE. )
 9997 FORMAT (66H0NON PROBLEM CARD ENCOUNTERED CHECK DATA DECK. EXECUTIO
     1N DELETED.  )
 9998 FORMAT(75H THE NUMBER OF SPECIAL TRANSGENERATION CARDS IS ILLEGAL.
     1 PROBLEM CARD CODE=  3X,  A6,  3X,  2HT=,  3X,  I3,  1H.  )
 9999 FORMAT(54H0ILLEGAL NUMBER OF DOSES SPECIFIED. EXECUTION DELETED.)
      END
C        SUBROUTINE TPWD FOR BMD03S                    JUNE 14, 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
C        SUBROUTINE TRNGEN FOR BMD03S                  JUNE 14, 1966
      SUBROUTINE TRNGEN (DOSE,CONST,M,ICODE,K,IERROR,X)
      DIMENSION DOSE(1000),CONST(10),ICODE(10),X(1000)
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
      IERROR = 0
      DO 2 J2 = 1,M
 2    X(J2) = DOSE(J2)
      IF (ICODE(1)) 1000,1000,5
 5    CONTINUE
      DO 210 JK = 1,8
      JUMP = ICODE(JK)
      FM = CONST(JK)
      IF(JUMP*(JUMP-11)) 200,1000,1000
 200    DO 8 J=1,M
      D = X(J)
        GO TO (10,20,30,40,50,60,70,80,90,110), JUMP
C  1  SQRT(DOSE)
  10    IF(D)99,11,12
 11   TD = 0.0
        GO TO 8
 12   TD = SQRT(D)
        GO TO 8
C  2  SQRT(DOSE) + SQRT(DOSE + 1)
 20        IF( D) 99,21,22
 21   TD = 1.0
        GO TO 8
 22   TD = SQRT(D) + SQRT(D+1.0)
        GO TO 8
C  3  LOG10( DOSE)
   30  IF(D) 99,99,31
 31   TD =ALOG10(D)
        GO TO 8
C  4  EXP(DOSE).
 40   TD = EXP(D)
        GO TO 8
C  5  ARCSINE(DOSE)
   50 IF(-D)52,11,99
  52   IF(D-1.0)   53,54,99
 54   TD = .157079632E+01
        GO TO 8
 53   CONTINUE
      TD = ASN(A)
        GO TO 8
C  6  ARCSIN(SQRT(DOSE/(M+1))) + ARCSIN(SQRT(A+10(M+1)))
  60   SAMP = M
       A=D/(SAMP+1.0)
       B=A+1.0/(SAMP+1.0)
      IF (A)99,61,61
 61   TD = ASN(SQRT(A))
 62   IF (B) 99,65,65
 65   TD = TD + ASN(SQRT(B))
        GO TO 8
C  7   1/ DOSE
  70    IF(D) 71,99,71
 71   TD = 1.0/D
        GO TO 8
C  8  DOSE + CONST
 80   TD = D + FM
        GO TO 8
C  9  DOSE X CONST
 90   TD = D*FM
        GO TO 8
C  10  DOSETO POWER OF CONST
 110  TD = D**FM
 8    X(J) = TD
 210   CONTINUE
 1000   RETURN
  99   IERROR =-999
      WRITE (6,9980)JUMP,DOSE(J),D
       GO TO 1000
 9980 FORMAT(43H0TRANSFORMATION CANNOT BE PERFORMED. CODE =   I4,
     1 6H DOSE=  F12.4,  5H ARG=  F12.4)
      END