PROGRAM PMAIN
C *** MAIN PROGRAM FOR RUNNING PREG EXAMPLES USING DGLG ***
      INTEGER LIV, LV, MMAX, NMAX, NW, NR0, PMAX
      PARAMETER (LIV=200, LV=8000, NW=6, MMAX = 18, NMAX=200, NR0=8,
     1           PMAX=20)
      CHARACTER*72 FNAME
      CHARACTER*6 ALGNAM(4)
      INTEGER ALG, I, IV(LIV), J, J0, J1, K, KDIAG, M, MDL(6), MODEL,
     1        N, NIN, NR, NRUN, P, P0, PS, RHOI(NMAX+6), UI(7)
      DOUBLE PRECISION A((MMAX+6)*NMAX), B(2,PMAX),
     1                 RHOR((17+PMAX)*NMAX+4), T, T1, V(LV), X(PMAX+3),
     1                 X0(PMAX+3), YN(2,7*NMAX+3)
      EQUIVALENCE (RHOI(1), MDL(1)), (RHOR(1), YN(1,1))
      CHARACTER*96 DESC, FMT
      CHARACTER*8 WNAME(4)
      DOUBLE PRECISION DR7MDC
      EXTERNAL BRJ, CHKDER, DEVIAN, DGLF, DGLFB, DGLG, DGLGB, DIVSET,
     1         DR7MDC, DV7CPY, DV7SCP, LOUCHK, POIX0, RHPOIL, RPOIL0
      DOUBLE PRECISION ONE
      INTEGER BS, BSSTR, F, FLO, FLOSTR, LOO, NB, NFIX, RDREQ, XNOTI
      PARAMETER (BS=85, BSSTR=86, F=10, FLO=88, FLOSTR=89, LOO=84,
     1           NB=87, NFIX=83, RDREQ=57, XNOTI=90)
      DATA ALG/1/, KDIAG/0/, NIN/5/
      DATA ALGNAM(1)/'DGLG'/,  ALGNAM(2)/'DGLF'/
      DATA ALGNAM(3)/'DGLGB'/, ALGNAM(4)/'DGLFB'/
      DATA ONE/1.D+0/
      DATA WNAME(1)/'  RHO"  '/, WNAME(2)/'  IRLS  '/,
     1     WNAME(3)/' SCORE  '/, WNAME(4)/'DEVIANCE'/
C
C *** BODY ***
C
      CALL DIVSET(1, IV, LIV, LV, V)
      IV(FLO) = 16*NMAX + 5
      IV(XNOTI) = IV(FLO) + NMAX
      IV(BS) = 7
      IV(BSSTR) = 1
      IV(FLOSTR) = 1
      IV(LOO) = 1
      IV(NB) = 5
      IV(NFIX) = 0
      CALL DV7SCP(NMAX, RHOR(IV(FLO)), ONE)
      CALL DV7SCP(NMAX, RHOR(IV(XNOTI)), -2.D+0)
      DO 10 I = IV(BS), IV(BS) + NMAX - 1
 10      RHOI(I) = 1
      T = DR7MDC(6)
      DO 20 I = 1, PMAX
         B(1,I) = -T
         B(2,I) = T
 20      CONTINUE
      NRUN = 0
      MDL(6) = 1
 30   READ(NIN,*,END=210) K
      WRITE(NW,*) '*', K
      GO TO (40, 50, 60, 70, 80, 90, 100, 110, 170, 180, 220,
     1       230, 240, 250, 260, 270, 300, 310, 320, 340,
     2       350, 360, 370, 380, 390, 430, 440, 450),  K
      WRITE(NW,*) '/// Invalid command', K
 40   WRITE(NW,*) '1 = LIST MENU'
      WRITE(NW,*) '2 = READ IV'
      WRITE(NW,*) '3 = READ V'
      WRITE(NW,*)
     1 '4 = READ ALG: 1 = DGLG, 2 = DGLF, 3 = DGLGB, 4 = DGLFB'
      WRITE(NW,*) '5 = READ ALL OF X0'
      WRITE(NW,*) '6 = COPY X TO X0'
      WRITE(NW,*) '7 = START'
      WRITE(NW,*) '8 = CONTINUE'
      WRITE(NW,*) '9 = READ COMMANDS FROM SPECIFIED FILE'
      WRITE(NW,*) '10 = READ PROBLEM'
      WRITE(NW,*) '11 = READ RHO'
      WRITE(NW,*) '12 = READ MODEL'
      WRITE(NW,*) '13 = CHECK RHO DERIVATIVES'
      WRITE(NW,*) '14 = READ P'
      WRITE(NW,*) '15 = READ X0 COMPONENTWISE'
      WRITE(NW,*) '16 = read new Y'
      WRITE(NW,*)
     1 '17 = negate RHO (negative ==> use weights; see KW = 19)'
      WRITE(NW,*) '18 = read KDIAG: 1 = from X*, 2 = from X0, 3 = both'
      WRITE(NW,*)
     1 '19 = read KW: 1 = RHO", 2 = IRLS, 3 = score, 4 = deviance'
      WRITE(NW,*) '20 = READ B (format i, b(1,i), b(2,i))'
      WRITE(NW,*) '21,22 = Read,Show RHOI (componentwise)'
      WRITE(NW,*) '23,24 = Read,Show RHOR        "'
      WRITE(NW,*) '25 = Show range of RHOR components'
      WRITE(NW,*) '26,27 = Show IV, V components'
      WRITE(NW,*) '28 = Read and echo comment'
      GO TO 30
 50   READ(NIN,*,END=210) I, J
      IF (I .LE. 0) GO TO 30
      IV(I) = J
      GO TO 50
 60   READ(NIN,*,END=210) I, T
      IF (I .LE. 0) GO TO 30
      V(I) = T
      GO TO 60
 70   READ(NIN,*,END=210) ALG
      GO TO 30
 80   READ(NIN,*,END=210) (X0(I), I = 1, P0)
      GO TO 30
 90   CALL DV7CPY(P0+3, X0, X)
      GO TO 30
 100  CALL DV7CPY(P0+3, X, X0)
      IV(1) = 12
 110  UI(6) = M
      NRUN = NRUN + 1
      IF (IV(1) .EQ. 0 .OR. IV(1) .EQ. 12) THEN
         WRITE(NW,'(/'' Run'',I5,'':  calling '',A,'' with PS ='',I5)')
     1      NRUN, ALGNAM(ALG), PS
       ELSE
         WRITE(NW,'(/'' Run'',I5,'':  continuing '',A,'', PS ='',I5)')
     1      NRUN, ALGNAM(ALG), PS
         END IF
      IF (KDIAG .GT. 0) IV(RDREQ) = 2
      GO TO (120,130,140,150), ALG
 120  CALL DGLG(N, P, PS, X, RHPOIL, RHOI, YN,
     1       IV, LIV, LV, V, BRJ, UI, A, BRJ)
      GO TO 160
 130  CALL DGLF(N, P, PS, X, RHPOIL, RHOI, YN,
     1       IV, LIV, LV, V, BRJ, UI, A, BRJ)
      GO TO 160
 140  CALL DGLGB(N, P, PS, X, B, RHPOIL, RHOI, YN,
     1       IV, LIV, LV, V, BRJ, UI, A, BRJ)
      GO TO 30
 150  CALL DGLFB(N, P, PS, X, B, RHPOIL, RHOI, YN,
     1       IV, LIV, LV, V, BRJ, UI, A, BRJ)
      GO TO 30
 160  IF (IV(1) .LT. 8) THEN
         CALL DEVIAN(V(F), MDL(1), N, NW, X(PS+1), YN)
         IF (ALG .EQ. 1) CALL LOUCHK(KDIAG,   DGLG, X0, N, P, PS, X,
     1       RHPOIL, MDL, YN, IV, LIV, LV, V, BRJ, UI, A, BRJ)
         IF (ALG .EQ. 2) CALL LOUCHK(KDIAG,   DGLF, X0, N, P, PS, X,
     1       RHPOIL, MDL, YN, IV, LIV, LV, V, BRJ, UI, A, BRJ)
         END IF
      GO TO 30
 170  IF (NIN .LE. 1) THEN
         WRITE(NW,*) '*** TOO MANY FILES OPEN'
         GO TO 30
         END IF
      READ(NIN,'(A)',END=200) FNAME
      NIN = NIN - 1
      OPEN(NIN,FILE=FNAME,STATUS='OLD',ERR=410)
      REWIND NIN
      GO TO 30
 180  READ(NIN,'(A)',END=200) FNAME
      IF (FNAME .EQ. '-') THEN
         NR = NIN
      ELSE
         OPEN(NR0,FILE=FNAME,STATUS='OLD',ERR=410)
         REWIND NR0
         NR = NR0
         END IF
      READ(NR, '(A)', END=200) DESC
      WRITE(NW,*) DESC
      READ(NR, '(9I4)', END=200) N, P, MODEL, M, MDL(1), I, J, PS
      P0 = P
      IF (PS .EQ. 0) PS = P
      IF (MODEL .LE. 2) M = PS
      IF (MIN(MDL(1),M,N,PS,P-PS+1,MODEL+1) .LE. 0 .OR. P .GT. PMAX
     1          .OR. M .GT. MMAX) THEN
         WRITE(NW,*) 'INVALID PROBLEM DIMENSIONS: M, N, P, MODEL  =',
     1                  M, N, P, MODEL
         STOP
         END IF
      MDL(2) = P
      MDL(3) = PS
      UI(1) = M
      UI(2) = MODEL
      UI(3) = 2
      UI(4) = 0
      UI(5) = 0
      UI(7) = PS
      CALL DV7SCP(3, X0(P+1), ONE)
      IF (MODEL .GT. 2) THEN
        READ(NR, *, END=200) (X0(I), I = 1, P)
       ELSE IF (PS .LT. P) THEN
        READ(NR, *, END=200) (X0(I), I = PS+1, P)
        END IF
      READ(NR, '(A)', END=200) FMT
      J1 = 0
      DO 190 I = 1, N
         J0 = J1 + 1
         J1 = J1 + M
         READ(NR, FMT, END=200) YN(1,I), YN(2,I), (A(J), J = J0, J1)
C        FROME*S DOCUMENTATION CLAIMS Y(I) IS YBAR(I), BUT HIS PROGRAM
C        ASSUMES IT IS THE TOTAL COUNT AND TURNS Y(I) INTO YBAR(I)
C        BY THE EQUIVALENT OF THE FOLLOWING STATEMENT...
C        YN(1,I) = YN(1,I) / YN(2,I)
 190     CONTINUE
      IF (MODEL .LE. 2) THEN
          CALL POIX0(A, IV, PS, LIV, LV, MODEL, N, PS, V, X0, YN)
          END IF
      GO TO 30
 200  WRITE(NW,*) '*** PREMATURE END OF FILE'
      IF (NR .NE. NIN) GO TO 30
 210  IF (NIN .GE. 5) STOP
      NIN = NIN + 1
      GO TO 30
 220  READ(NIN,*,END=210) I
      IF (I .LE. 0) I = MDL(1)
      WRITE(NW,*) 'Changing RHO from ', MDL(1), ' to ', I
      MDL(1) = I
      GO TO 30
 230  READ(NIN,*,END=210) I
      IF (I .EQ. 0) I = MODEL
      WRITE(NW,*) 'Changing MODEL from ', MODEL, ' to ', I
      MODEL = I
      UI(2) = MODEL
      GO TO 30
 240  CALL CHKDER(MDL, N, P-PS, X0(PS+1), V(200), RHPOIL, RPOIL0, YN)
      GO TO 30
 250  READ(NIN,*,END=210) I
      IF (I .GT. P0 .OR. I .LT. P0-3) THEN
         WRITE(NW,*) 'INVALID P = ', I, ' -- P REMAINS ', P
       ELSE
         P = I
         MDL(2) = I
         END IF
      GO TO 30
 260  READ(NIN,*,END=210) I, T
      IF (I .LE. 0) GO TO 30
      X0(I) = T
      GO TO 260
 270  DO 280 I = 1, N
 280     READ(NIN, FMT, END=290) YN(1,I), YN(2,I)
      GO TO 30
 290  WRITE(NW,*) 'Premature end of file!'
      GO TO 210
 300  I = 1
      IF (MDL(6) .EQ. 1) I = 2
      GO TO 330
 310  READ(NIN,*,END=210) KDIAG
      GO TO 30
 320  READ(NIN,*,END=210) I
      I =  MIN(4, MAX0(I,1))
 330  WRITE(NW,*) 'KW changed from ', MDL(6), ' = ', WNAME(MDL(6)),
     1 ' to ', I, ' = ', WNAME(I)
      MDL(6) = I
      GO TO 30
 340  READ(NIN,*,END=210) I, T, T1
      IF (I .LE. 0) GO TO 30
      B(1,I) = T
      B(2,I) = T1
      GO TO 340
 350  READ(NIN,*,END=210) I, J
      IF (I .LE. 0) GO TO 30
      RHOI(I) = J
      GO TO 350
 360  READ(NIN,*,END=210) I
      IF (I .LE. 0) GO TO 30
      WRITE(*,*) 'RHOI(',I,') = ', RHOI(I)
      GO TO 360
 370  READ(NIN,*,END=210) I, T
      IF (I .LE. 0) GO TO 30
      RHOR(I) = T
      GO TO 370
 380  READ(NIN,*,END=210) I
      IF (I .LE. 0) GO TO 30
      WRITE(*,*) 'RHOR(',I,') = ', RHOR(I)
      GO TO 380
 390  READ(NIN,*,END=210) I, J
      IF (I .LE. 0) GO TO 30
      WRITE(*,*) (RHOR(K), K = I, J)
      GO TO 390
 410  WRITE(*,420) FNAME
 420  FORMAT(' Can''t open ',A)
      GO TO 30
 430  READ(NIN,*,END=210) I
      IF (I .LE. 0) GO TO 30
      WRITE(*,*) 'IV(',I,') = ', IV(I)
      GO TO 430
 440  READ(NIN,*,END=210) I
      IF (I .LE. 0) GO TO 30
      WRITE(*,*) 'V(',I,') = ', V(I)
      GO TO 440
 450  READ(NIN,'(A)',END=200) FNAME
      WRITE(NW,*) FNAME
      GO TO 30
      END
      SUBROUTINE BRJ(N, P, X, NF, NEED, R, RP, UI, A, UF)
      INTEGER N, P, NF, NEED(2), UI(7)
      DOUBLE PRECISION X(P), R(N), RP(P,N), A(*)
      EXTERNAL UF
      EXTERNAL BRJ1
      INTEGER M
C
C *** BODY ***
C
      M = UI(6)
      CALL BRJ1(M, N, UI(7), X, NF, NEED, R, RP, UI, A, A(M*N+1), UF)
 999  RETURN
      END
      SUBROUTINE BRJ1(M, N, P, X, NF, NEED, R, RP, UI, A, UR, UF)
      INTEGER M, N, P, NF, NEED(2), UI(5)
      DOUBLE PRECISION X(P), R(N), RP(P,N), A(M,N), UR(N,6)
      EXTERNAL UF
      EXTERNAL DD7TPR, DR7MDC
      DOUBLE PRECISION DD7TPR, DR7MDC
C
C *** LOCAL VARIABLES ***
C
      INTEGER I, J, J2, J4, MODEL
      DOUBLE PRECISION ALPHA, BETA1, BETA2, DI, E, EMX, PHI, T, T1,
     1                 THETA, TI, X1, X1INV, X2, X3, X3M1, X4
      DOUBLE PRECISION EXPMAX, EXPMIN, ONE, TWO, ZERO
      DATA EXPMAX/0.D+0/, EXPMIN/0.D+0/, ONE/1.D+0/, TWO/2.D+0/,
     1     ZERO/0.D+0/
C
C *** BODY ***
C
      MODEL = IABS(UI(2))
      IF (MODEL .LE. 0) GO TO 520
      IF (MODEL .GT. 11) GO TO 520
      IF (EXPMAX .GT. ZERO) GO TO 10
         EXPMAX = TWO * DLOG(DR7MDC(5))
         EXPMIN = TWO * DLOG(DR7MDC(2))
 10   IF (NEED(1) .EQ. 2) GO TO 260
      J = 3 - UI(3)
      IF (UI(3+J) .EQ. NEED(2)) J = UI(3)
      UI(3) = J
      UI(3+J) = NF
      J2 = J + 2
      J4 = J + 4
      GO TO (20, 40, 60, 60, 80, 100, 120, 170, 190, 210, 230), MODEL
C
C *** LINEAR MODEL ***
C
 20   DO 30 I = 1, N
 30      R(I) = DD7TPR(P, X, A(1,I))
      GO TO 999
C
C *** EXPONENTIAL OF LINEAR ***
C
 40   DO 50 I = 1, N
         T = DD7TPR(P, X, A(1,I))
         IF (T .GE. EXPMAX) GO TO 520
         E = ZERO
         IF (T .GT. EXPMIN) E = DEXP(T)
         R(I) = E
         UR(I,J) = E
 50      CONTINUE
      GO TO 999
C
C *** NONLINEAR POISSON EXAMPLE FROM FROME*S PREG MANUAL ***
C
 60   X1 = X(1)
      X2 = X(2)
      X3 = X(3)
      DO 70 I = 1, N
         E = DEXP(-X2*A(2,I))
         UR(I,J2) = E
         T = (ONE - E) ** X3
         UR(I,J4) = T
         T = X1*A(1,I) * (ONE - T)
         IF (T .LE. ZERO) GO TO 520
         UR(I,J) = T
         IF (MODEL .EQ. 3) T = DLOG(T)
         R(I) = T
 70      CONTINUE
      GO TO 999
C
C *** CAESIUM DOSE EFFECT MODEL ***
C
 80   X1 = X(1)
      X2 = X(2)
      X3 = X(3)
      DO 90 I = 1, N
         DI = A(1,I)
         TI = A(2,I)
         IF (X3 .EQ. ZERO) GO TO 520
         IF (TI .EQ. ZERO) GO TO 520
         T = -TI / X3
         IF (T .GE. EXPMAX) GO TO 520
         E = ZERO
         IF (T .GT. EXPMIN) E = DEXP(T)
         UR(I,J) = E
         T = X3 / TI
         T = DI * (X2 + TWO*T*DI*(ONE - T*(ONE - E)))
         UR(I,J2) = T
         R(I) = X1 * T
 90      CONTINUE
      GO TO 999
C
C *** LUNG CANCER MODEL ***
C
 100  X1 = X(1)
      X2 = X(2)
      X3 = X(3)
      X4 = X(4)
      EMX = EXPMAX - 10.D+0
      DO 110 I = 1, N
         T1 = X1 * A(1,I)
         T = X2 + X3*A(2,I) + T1
         IF (T .GE. EMX) GO TO 520
         E = ZERO
         IF (T .GT. EXPMIN) E = DEXP(T)
         T = X4 + T1
         IF (T .GE. EMX) GO TO 520
         T1 = ZERO
         IF (T .GT. EXPMIN) T1 = DEXP(T)
         T = E + T1
         R(I) = T
         UR(I,J) = E
         UR(I,J2) = T1
         UR(I,J4) = T
 110     CONTINUE
      GO TO 999
C
C *** LOGISTIC OF LINEAR ***
C
 120  DO 160 I = 1, N
         T = DD7TPR(P, A(1,I), X)
         IF (T .LE. EXPMIN) GO TO 130
         IF (T .GE. EXPMAX) GO TO 140
         E = DEXP(T)
         T1 = ONE / (ONE + E)
         T = E * T1
         T1 = T * T1
         GO TO 150
 130     T = ZERO
         T1 = ZERO
         GO TO 150
 140     T = ONE
         T1 = ZERO
 150     R(I) = T
         UR(I,J) = T1
 160     CONTINUE
      GO TO 999
C
C *** LOG OF LINEAR ***
C
 170  DO 180 I = 1, N
         T = DD7TPR(P, X, A(1,I))
         IF (T .LE. ZERO) GO TO 520
         R(I) = DLOG(T)
         UR(I,J) = T
 180     CONTINUE
      GO TO 999
C
C *** EXAMPLE ON P. 204 OF MCCULLAGH AND NELDER ***
C
 190  ALPHA = X(1)
      BETA1 = X(2)
      BETA2 = X(3)
      PHI = X(4)
      DO 200 I = 1, N
         X2 = A(2,I)
         R(I) = ALPHA + BETA1*DLOG(A(1,I)) + BETA2*X2/(PHI + X2)
 200     CONTINUE
      GO TO 999
C
C *** EXAMPLE ON P. 205 OF MCCULLAGH AND NELDER ***
C
 210  ALPHA = X(1)
      BETA1 = X(2)
      BETA2 = X(3)
      PHI = X(4)
      THETA = X(5)
      DO 220 I = 1, N
         X2 = A(2,I)
         T = A(1,I) - THETA
         IF (T .LE. ZERO) GO TO 520
         R(I) = ALPHA + BETA1*DLOG(T) + BETA2*X2/(PHI + X2)
 220     CONTINUE
      GO TO 999
C
C *** EXAMPLE P. 202 OF MCCULLAGH AND NELDER ***
C
 230  DO 250 I = 1, N
         T = X(1)
         DO 240 J = 1, 3
            T1 = A(J,I) + X(2*J+1)
            IF (T1 .LE. ZERO) GO TO 520
 240        T = T + X(2*J)/T1
         R(I) = T
 250     CONTINUE
      GO TO 999
C
C *** JACOBIAN EVALUATIONS...
C
 260  J = UI(3)
      IF (NF .EQ. UI(J+3)) GO TO 270
      J = 3 - J
      IF (NF .EQ. UI(J+3)) GO TO 270
      WRITE(6,*) 'HELP! UNAVAILABLE INTERMEDIATE INFO!'
      GO TO 520
 270  J2 = J + 2
      J4 = J + 4
      GO TO (280, 290, 310, 340, 370, 390, 410, 430, 450, 470, 490),
     1           MODEL
C
C *** LINEAR MODEL ***
C
C
 280  CALL DV7CPY(N*P, RP, A)
      GO TO 999
C
C *** EXPONENTIAL OF LINEAR MODEL ***
C
 290  DO 300 I = 1, N
 300     CALL DV7SCL(P, RP(1,I), UR(I,J), A(1,I))
      GO TO 999
C
C *** LOG OF NONLINEAR POISSON EXAMPLE FROM FROME*S PREG MANUAL ***
C
 310  X1 = X(1)
      X2 = X(2)
      X3 = X(3)
      X3M1 = X3 - ONE
      X1INV = ONE / X1
      DO 330 I = 1, N
         RP(1,I) = X1INV
         E = UR(I,J2)
         T1 = ONE - E
         T = -A(1,I) * X1 / UR(I,J)
         RP(2,I) = T * X3 * A(2,I) * E * T1**X3M1
         IF (T1 .LE. ZERO) GO TO 320
         RP(3,I) = T * UR(I,J4) * DLOG(T1)
         GO TO 330
 320     RP(3,I) = ZERO
 330     CONTINUE
      GO TO 999
C
C *** NONLINEAR POISSON EXAMPLE FROM FROME*S PREG MANUAL ***
C
 340  X1 = X(1)
      X2 = X(2)
      X3 = X(3)
      X3M1 = X3 - ONE
      X1INV = ONE / X1
      DO 360 I = 1, N
         RP(1,I) = A(1,I) * (ONE - UR(I,J4))
         E = UR(I,J2)
         T1 = ONE - E
         T = -A(1,I) * X1
         RP(2,I) = T * X3 * A(2,I) * E * T1**X3M1
         IF (T1 .LE. ZERO) GO TO 350
         RP(3,I) = T * UR(I,J4) * DLOG(T1)
         GO TO 360
 350     RP(3,I) = ZERO
 360     CONTINUE
      GO TO 999
C
C *** CAESIUM DOSE EFFECT MODEL ***
C
 370  X1 = X(1)
      X3 = X(3)
      DO 380 I = 1, N
         RP(1,I) = UR(I,J2)
         DI = A(1,I)
         TI = A(2,I)
         RP(2,I) = X1 * DI
         E = UR(I,J)
         T = TWO * X3 / TI
         RP(3,I) = TWO * X1 * (DI/TI) * DI * (ONE - T + E*(T + ONE))
 380     CONTINUE
      GO TO 999
C
C *** LUNG CANCER MODEL ***
C
 390  DO 400 I = 1, N
         RP(1,I) = UR(I,J4) * A(1,I)
         T = UR(I,J)
         RP(2,I) = T
         RP(3,I) = T * A(2,I)
         RP(4,I) = UR(I,J2)
 400     CONTINUE
      GO TO 999
C
C *** LOGISTIC OF LINEAR ***
C
 410  DO 420 I = 1, N
 420     CALL DV7SCL(P, RP(1,I), UR(I,J), A(1,I))
      GO TO 999
C
C *** LOG OF LINEAR ***
C
 430  DO 440 I = 1, N
 440     CALL DV7SCL(P, RP(1,I), ONE/UR(I,J), A(1,I))
      GO TO 999
C
C *** EXAMPLE ON P. 204 OF MCCULLAGH AND NELDER ***
C
 450  ALPHA = X(1)
      BETA1 = X(2)
      BETA2 = X(3)
      PHI = X(4)
      DO 460 I = 1, N
         X2 = A(2,I)
C        R(1,I) = ALPHA + BETA1*DLOG(A(1,I)) + BETA2*X2/(PHI + X2)
         RP(1,I) = ONE
         RP(2,I) = DLOG(A(1,I))
         RP(3,I) = X2/(PHI + X2)
         RP(4,I) = -BETA2*X2/(PHI + X2)**2
         RP(1,I) = ONE
 460     CONTINUE
      GO TO 999
C
C
C *** EXAMPLE ON P. 205 OF MCCULLAGH AND NELDER ***
C
 470  ALPHA = X(1)
      BETA1 = X(2)
      BETA2 = X(3)
      PHI = X(4)
      THETA = X(5)
      DO 480 I = 1, N
         X2 = A(2,I)
C        R(I) = ALPHA + BETA1*DLOG(A(1,I) - THETA) + BETA2*X2/(PHI + X2)
         RP(1,I) = ONE
         RP(2,I) = DLOG(A(1,I) - THETA)
         RP(3,I) = X2/(PHI + X2)
         RP(4,I) = -BETA2*X2/(PHI + X2)**2
         RP(5,I) = -BETA1/(A(1,I) - THETA)
 480     CONTINUE
      GO TO 999
C
C *** EXAMPLE P. 202 OF MCCULLAGH AND NELDER ***
C
 490  DO 510 I = 1, N
C        DO 453 J = 1, 3
C453        RI = RI + X(2*J)/(A(J,I) + X(2*J+1))
         RP(1,I) = ONE
         DO 500 J = 1, 3
            T = ONE / (A(J,I) + X(2*J+1))
            RP(2*J,I) = T
            RP(2*J+1,I) = -X(2*J)*T*T
 500        CONTINUE
 510     CONTINUE
      GO TO 999
 520  NF = 0
 999  RETURN
      END
      SUBROUTINE CHKDER(MDL, N, NPT, PT, R, RHO, RHO0, YN)
      INTEGER MDL(1), N, NPT
C     DOUBLE PRECISION PT(NPT) -- BUT NPT MAY BE 0
      DOUBLE PRECISION PT(1), R(N,20), YN(2,N)
      EXTERNAL RHO, RHO0
      EXTERNAL DV2NRM
      DOUBLE PRECISION DV2NRM
      INTEGER I, J
      DOUBLE PRECISION F, H, T
      REAL FOO(10), FAC
      DATA FOO/.1, -.1, .2, -.2, .4, -.4, .6, -.6, .8, -.9/, H/.001D0/
C
C *** BODY ***
C
      J = 1
      FAC = 1.0
      DO 10 I = 1, N
         T = FAC * FOO(J)
         R(I,1) = T
         R(I,10) = T + H
         J = J + 1
         IF (J .LE. 10) GO TO 10
                J = 1
                FAC = 10. * FAC
 10      CONTINUE
      CALL RHO0(MDL, N, PT, R, R(1,4), YN)
      CALL RHO0(MDL, N, PT, R(1,10), R(1,13), YN)
      DO 20 I = 1, N
         T = R(I,10) - R(I,1)
         IF (T .NE. 0.D0) T = 1.D0 / T
         R(I,20) = T
 20      CONTINUE
      CALL DV2AXY(N, R(1,13), -1.D0, R(1,4), R(1,13))
      CALL DV7VMP(N, R(1,13), R(1,13), R(1,20), 1)
      J = 1
      CALL RHO(0, F, N, J, PT, R, R(1,4), MDL, YN)
      CALL RHO(1, F, N, J, PT, R, R(1,4), MDL, YN)
      CALL DV2AXY(N, R(1,19), -1.D0, R(1,13), R)
      T = DV2NRM(N,R(1,19))/(DV2NRM(N,R(1,13)) + DV2NRM(N,R))
      WRITE(6,*) '1ST DERIV RELATIVE DIFFERENCE =', T
      IF (T .GT. .01) THEN
        WRITE(6,*) 'I   FD(I)   AN(I)'
        WRITE(6,'(I5,2G13.4)') (I, R(I,13), R(I,1), I = 1, N)
        END IF
      CALL RHO(0, F, N, J, PT, R(1,10), R(1,13), MDL, YN)
      CALL RHO(1, F, N, J, PT, R(1,10), R(1,13), MDL, YN)
      CALL DV2AXY(N, R(1,19), -1.D0, R, R(1,10))
      CALL DV7VMP(N, R(1,19), R(1,19), R(1,20), 1)
      CALL DV2AXY(N, R(1,13), -1.D0, R(1,19), R(1,4))
      T = DV2NRM(N,R(1,13))/(DV2NRM(N,R(1,4)) + DV2NRM(N,R(1,19)))
      WRITE(6,*) '2ND DERIV RELATIVE DIFFERENCE =', T
      IF (T .GT. .01) THEN
        WRITE(6,*) 'I   FD(I)   AN(I)'
        WRITE(6,'(I5,2G13.4)') (I, R(I,19), R(I,4), I = 1, N)
        END IF
 999  RETURN
      END
      SUBROUTINE RPOIL0(MDL, N, PT, R, RHO, YN)
      INTEGER N, MDL(1)
      DOUBLE PRECISION PT(1), R(N), RHO(N), YN(2,N)
      EXTERNAL LPN, DR7MDC
      DOUBLE PRECISION LPN, DR7MDC
      INTEGER I, MODEL
      DOUBLE PRECISION E, RI, T, YI
      DOUBLE PRECISION DEXP, DLOG
      DOUBLE PRECISION EXPMAX, EXPMIN, HALF, ONE, TWO, ZERO
      DATA EXPMAX/0.D+0/, EXPMIN/0.D+0/,
     1     HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C
C *** BODY ***
C
      MODEL = MDL(1)
      I = MODEL + 2
      IF (I .LE. 0 .OR. I .GT. 11) THEN
        WRITE(6,*) 'HELP! RPOIL0 HAS MODEL =', MODEL
        STOP
        END IF
      IF (EXPMAX .GT. ZERO) GO TO 10
         EXPMAX = TWO * DLOG(DR7MDC(5))
         EXPMIN = TWO * DLOG(DR7MDC(2))
 10   GO TO (20, 20, 40, 60, 80, 80, 100, 120, 140, 160, 180), I
C
C *** POISSON RHO (AND CONVENTIONAL IRLS) ***
C
 20   DO 30 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) THEN
                RI = ONE
                R(I) = ONE
                END IF
         RHO(I) = YN(2,I)*RI - YN(1,I)*DLOG(RI)
 30      CONTINUE
      GO TO 999
C
C *** LOG LINEAR ***
C
 40   DO 50 I = 1, N
         E = ZERO
         RI = R(I)
         IF (RI .GT. EXPMAX) THEN
                RI = HALF * EXPMAX
                R(I) = RI
                END IF
         IF (RI .GT. EXPMIN) E = EXP(RI)
         RHO(I) = YN(2,I)*E - YN(1,I)*RI
 50      CONTINUE
      GO TO 999
C
C *** SQUARE-ROOT LINEAR POISSON ***
C
 60   DO 70 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) THEN
                RI = ONE
                R(I) = RI
                END IF
         RHO(I) = YN(2,I)*RI**2 - TWO*YN(1,1)*DLOG(RI)
 70      CONTINUE
      GO TO 999
C
C *** BINOMIAL RHO (AND CONVENTIONAL IRLS) ***
C
 80   DO 90 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO .OR. RI .GE. ONE) THEN
                RI = HALF
                R(I) = RI
                END IF
         RHO(I) = -YN(1,I)*DLOG(RI) - (YN(2,I) - YN(1,I))*DLOG(ONE-RI)
 90      CONTINUE
      GO TO 999
C
C *** BINOMIAL LOGISTIC RHO ***
C
 100  DO 110 I = 1, N
         RI = R(I)
         IF (RI .GT. EXPMAX) THEN
                RI = HALF * EXPMAX
                R(I) = RI
                END IF
         E = ZERO
         IF (RI .GT. EXPMIN) E = DEXP(RI)
         RHO(I) = YN(2,I)*DLOG(ONE + E) - YN(1,I)*RI
 110     CONTINUE
      GO TO 999
C
C *** PROBIT ***
C
 120  DO 130 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         RHO(I) = -YI*LPN(RI) - (YN(2,I)-YI)*LPN(-RI)
 130     CONTINUE
      GO TO 999
C
C *** WEIBULL ***
C
 140  DO 150 I = 1, N
         RI = R(I)
         IF (RI .GT. EXPMAX) THEN
                RI = HALF * EXPMAX
                R(I) = RI
                END IF
         E = ZERO
         IF (RI .GT. EXPMIN) E = DEXP(RI)
         T = ZERO
         IF (-E .GT. EXPMIN) T = DEXP(-E)
         RHO(I) = (YN(2,I) - YN(1,I))*E - YN(1,I)*DLOG(ONE - T)
 150     CONTINUE
      GO TO 999
C
C  *** GAMMA ERRORS ***
C
 160  DO 170 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) THEN
                WRITE(6,*) 'HELP! CHKDER HAS R(',I,') =', RI,' < 0'
                STOP
                END IF
         RHO(I) = YN(2,I) * (YN(1,I)*RI - DLOG(RI))
 170     CONTINUE
      GO TO 999
C
C  ***  PREGIBON ERRORS ***
C
C      *** IN THIS CASE, YN(1,I) = Y(I), YN(2,I) = LOG(Y(I))
C      *** AND YN(I,J), J = N+1(1)2*N, I = 1 OR 2 = SCRATCH
C
 180  DO 190 I = 1, N
         IF (R(I) .LT. ZERO) R(I) = -R(I)
 190     CONTINUE
      CALL PRGRH1(N, PT, R, RHO, MDL, YN)
C
 999  RETURN
      END
      SUBROUTINE DEVIAN(F, MODEL0, N, NW, PT, YN)
      INTEGER MODEL0, N, NW
      DOUBLE PRECISION F, PT(2), YN(2,N)
      DOUBLE PRECISION DATAN, DLOG
      INTEGER I, MODEL
      DOUBLE PRECISION CI, D, S, T, T1, YI
      DOUBLE PRECISION EIGHT, HALF, ONE, TWO, ZERO
      DATA EIGHT/8.D+0/, HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/,
     1     ZERO/0.D+0/
C
C *** BODY ***
C
      D = F
      MODEL = IABS(MODEL0)
      IF (MODEL .LT. 5) GO TO 20
      IF (MODEL .GT. 9) GO TO (40, 60, 999, 80) MODEL - 9
C
C *** BINOMIAL DEVIANCE ***
C
      DO 10 I = 1, N
         YI = YN(1,I)
         CI = YN(2,I)
         T = YI / CI
         IF (T .GT. ZERO) D = D + YI*DLOG(T)
         IF (T .LT. ONE) D = D + (CI-YI)*DLOG(ONE-T)
 10      CONTINUE
      GO TO 100
C
C *** POISSON DEVIANCE ***
C
 20   DO 30 I = 1, N
         YI = YN(1,I)
         IF (YI .GT. ZERO) D = D + YI*(DLOG(YI/YN(2,I)) - ONE)
 30      CONTINUE
      GO TO 100
C
C *** GAMMA DEVIANCE ***
C
 40   DO 50 I = 1, N
         YI = YN(1,I)
         IF (YI .LE. ZERO) GO TO 999
         D = D - YN(2,I)*(ONE + DLOG(YI))
 50      CONTINUE
      GO TO 100
C
C  *** PREGIBON DEVIANCE, REPLICATE WEIGHTS ***
C
 60   T = PT(2)
      T1 = DLOG(EIGHT*DATAN(ONE)*PT(1))
      S = ZERO
      DO 70 I = 1, N
 70      S = S + YN(2,I) * (T*DLOG(DBLE(YN(1,I))) + T1)
      D = PT(1) * (D - HALF*S)
      GO TO 100
C
C  *** PREGIBON DEVIANCE, VARIANCE WEIGHTS ***
C
 80   S = ZERO
      T = ZERO
      DO 90 I = 1, N
         S = S + DLOG(DBLE(YN(1,I)))
         T = T + DLOG(DBLE(YN(2,I)))
 90      CONTINUE
      D = PT(1) * (D -
     1     HALF*(PT(2)*S - T + N*DLOG(EIGHT*DATAN(ONE)*PT(1))))
C
 100  WRITE(NW,*) 'DEVIANCE = ', TWO*D
 999  RETURN
      END
      DOUBLE PRECISION FUNCTION DZERO(F,A,B,T)
C *** THE PORT ROUTINE, MODIFIED TO STOP RATHER THAN CALLING SETERR ***
C *** AND TO CALL DR7MDC RATHER THAN D1MACH ***
C
C  FINDS THE REAL ROOT OF THE FUNCTION F LYING BETWEEN A AND B
C  TO WITHIN A TOLERANCE OF
C
C         6*D1MACH(3) *  ABS(DZERO) + 2 * T
C
C  F(A) AND F(B) MUST HAVE OPPOSITE SIGNS
C
C  THIS IS BRENTS ALGORITHM
C
C  A, STORED IN SA, IS THE PREVIOUS BEST APPROXIMATION (I.E. THE OLD B)
C  B, STORED IN SB, IS THE CURRENT BEST APPROXIMATION
C  C IS THE MOST RECENTLY COMPUTED POINT SATISFYING F(B)*F(C) .LT. 0
C  D CONTAINS THE CORRECTION TO THE APPROXIMATION
C  E CONTAINS THE PREVIOUS VALUE OF D
C  M CONTAINS THE BISECTION QUANTITY (C-B)/2
C
      DOUBLE PRECISION F,A,B,T,TT,SA,SB,C,D,E,FA,FB,FC,TOL,M,P,Q,R,S
      EXTERNAL F
      DOUBLE PRECISION DR7MDC
C
      TT = T
      IF (T .LE. 0.0D0) TT = 10.D0*DR7MDC(1)
C
      SA = A
      SB = B
      FA = F(SA)
      FB = F(SB)
      IF (FA .NE. 0.0D0) GO TO 5
      DZERO = SA
      RETURN
  5   IF (FB .EQ. 0.0D0) GO TO 140
        IF (DSIGN(FA,FB) .EQ. FA) THEN
                WRITE(*,*) 'DZERO: F(A) = ', FA, '; F(B) = ', FB
                STOP
                END IF
C
 10   C  = SA
      FC = FA
      E  = SB-SA
      D  = E
C
C  INTERCHANGE B AND C IF  ABS F(C) .LT.  ABS F(B)
C
 20   IF ( ABS(FC).GE. ABS(FB)) GO TO 30
      SA = SB
      SB = C
      C  = SA
      FA = FB
      FB = FC
      FC = FA
C
 30   TOL = 2.0D0*DR7MDC(3)* ABS(SB)+TT
      M = 0.5D0*(C-SB)
C
C  SUCCESS INDICATED BY M REDUCES TO UNDER TOLERANCE OR
C  BY F(B) = 0
C
      IF (( ABS(M).LE.TOL).OR.(FB.EQ.0.0D0)) GO TO 140
C
C  A BISECTION IS FORCED IF E, THE NEXT-TO-LAST CORRECTION
C  WAS LESS THAN THE TOLERANCE OR IF THE PREVIOUS B GAVE
C  A SMALLER F(B).  OTHERWISE GO TO 40.
C
      IF (( ABS(E).GE.TOL).AND.( ABS(FA).GE. ABS(FB))) GO TO 40
      E = M
      D = E
      GO TO 100
 40   S = FB/FA
C
C  QUADRATIC INTERPOLATION CAN ONLY BE DONE IF A (IN SA)
C  AND C ARE DIFFERENT POINTS.
C  OTHERWISE DO THE FOLLOWING LINEAR INTERPOLATION
C
      IF (SA.NE.C) GO TO 50
      P = 2.0D0*M*S
      Q = 1.0D0-S
      GO TO 60
C
C  INVERSE QUADRATIC INTERPOLATION
C
 50   Q = FA/FC
      R = FB/FC
      P = S*(2.0D0*M*Q*(Q-R)-(SB-SA)*(R-1.0D0))
      Q = (Q-1.0D0)*(R-1.0D0)*(S-1.0D0)
 60   IF (P.LE.0.0D0) GO TO 70
      Q = -Q
      GO TO 80
 70   P = -P
C
C  UPDATE THE QUANTITIES USING THE NEWLY COMPUTED
C  INTERPOLATE UNLESS IT WOULD EITHER FORCE THE
C  NEW POINT TOO FAR TO ONE SIDE OF THE INTERVAL
C  OR WOULD REPRESENT A CORRECTION GREATER THAN
C  HALF THE PREVIOUS CORRECTION.
C
C  IN THESE LAST TWO CASES - DO THE BISECTION
C  BELOW (FROM STATEMENT 90 TO 100)
C
 80   S = E
      E = D
      IF ((2.0D0*P.GE.3.0D0*M*Q- ABS(TOL*Q)).OR.
     1    (P.GE. ABS(0.5D0*S*Q))) GO TO 90
      D = P/Q
      GO TO 100
 90   E = M
      D = E
C
C  SET A TO THE PREVIOUS B
C
 100  SA = SB
      FA = FB
C
C  IF THE CORRECTION TO BE MADE IS SMALLER THAN
C  THE TOLERANCE, JUST TAKE A  DELTA STEP  (DELTA=TOLERANCE)
C         B = B + DELTA * SIGN(M)
C
      IF ( ABS(D).LE.TOL) GO TO 110
      SB = SB+D
      GO TO 130
C
 110  IF (M.LE.0.0D0) GO TO 120
      SB = SB+TOL
      GO TO 130
C
 120  SB = SB-TOL
 130  FB = F(SB)
C
C  IF F(B) AND F(C) HAVE THE SAME SIGN ONLY
C  LINEAR INTERPOLATION (NOT INVERSE QUADRATIC)
C  CAN BE DONE
C
      IF ((FB.GT.0.0D0).AND.(FC.GT.0.0D0)) GO TO 10
      IF ((FB.LE.0.0D0).AND.(FC.LE.0.0D0)) GO TO 10
      GO TO 20
C
C***SUCCESS***
 140  DZERO = SB
      RETURN
      END
        DOUBLE PRECISION FUNCTION INVCN(X, ERRFLG)
        DOUBLE PRECISION X
        INTEGER ERRFLG
        COMMON /INVCMN/ XC, TOL, NCALL
        DOUBLE PRECISION XC, TOL
        INTEGER NCALL

        DOUBLE PRECISION CNERR, DZERO, PNORMS, DR7MDC
        EXTERNAL CNERR, PNORMS, DR7MDC

        DOUBLE PRECISION A, B
        DOUBLE PRECISION HALF, ONE, ZERO
        LOGICAL FIRST
        DOUBLE PRECISION HUGE
        PARAMETER (HALF = 0.5D+0, ONE = 1.D+0, ZERO = 0.D+0)
        SAVE FIRST, HUGE
        DATA FIRST/.TRUE./, HUGE/0.D+0/

        IF (FIRST) THEN
                TOL = 10.D+0 * DR7MDC(1)
                HUGE = 0.1D+0 * DR7MDC(6)
                FIRST = .FALSE.
                END IF

        NCALL = 0
        ERRFLG = 0
        IF (X .LE. ZERO) THEN
C               IF (X .EQ. ZERO) THEN
C                       INVCN = -HUGE
C                       GO TO 999
C                       END IF
                ERRFLG = 1
                INVCN = ZERO
                GO TO 999
                END IF
        IF (X .GE. ONE) THEN
C               IF (X .EQ. ONE) THEN
C                       INVCN = HUGE
C                       GO TO 999
C                       END IF
                ERRFLG = 1
                INVCN = ZERO
                GO TO 999
                END IF
        IF (X .GE. HALF) THEN
                A = ZERO
                B = ONE
 10             IF (PNORMS(B) .LT. X) THEN
                        B = B + ONE
                        GO TO 10
                        END IF
        ELSE
                B = ZERO
                A = -ONE
 20             IF (PNORMS(A) .GT. X) THEN
                        A = A - ONE
                        GO TO 20
                        END IF
                END IF
        XC = X
        INVCN = DZERO(CNERR,A,B,TOL)
 999    RETURN
        END

        DOUBLE PRECISION FUNCTION CNERR(X)
        DOUBLE PRECISION X

        COMMON /INVCMN/ XC, TOL, NCALL
        DOUBLE PRECISION XC, TOL
        INTEGER NCALL

        DOUBLE PRECISION PNORMS
        EXTERNAL PNORMS
        NCALL = NCALL + 1
        CNERR = XC - PNORMS(X)
        END
      SUBROUTINE LOUCHK(KDIAG, DGLG, X0, N, P, PS, X, RHPOIL, MDL, YN,
     1                  IV, LIV, LV, V, BRJ, UI, A, BRJ1)
      EXTERNAL DGLG, RHPOIL, BRJ, BRJ1
      INTEGER KDIAG, N, P, PS, LIV, LV
      INTEGER IV(LIV), MDL(2), UI(*)
      DOUBLE PRECISION X0(P), X(P), V(LV), A(*), YN(N)
C
C *** DUMMY REPLACEMENT FOR C ROUTINE (USED FOR DEBUGGING) ***
C
      END
      DOUBLE PRECISION FUNCTION PNORMS(X)
      DOUBLE PRECISION X

      EXTERNAL MECDF
      DOUBLE PRECISION D(1), PROB, RHO(1)
      INTEGER IER

      D(1) = X
      CALL MECDF(1, D, RHO, PROB, IER)
      PNORMS = 1.D+0 - PROB
      END
      SUBROUTINE POISX0(A, C, LA, LC, MODEL, N, P, QTR, X, YN)
      INTEGER LA, LC, MODEL, N, P
      DOUBLE PRECISION A(LA,N), C(LC), QTR(P), X(P), YN(2,N)
      EXTERNAL DL7ITV, DL7SVX, DL7SVN,DQ7ADR, DR7MDC, DV7SCL, DV7SCP
      DOUBLE PRECISION DL7SVX, DL7SVN, DR7MDC
      INTEGER I
      DOUBLE PRECISION SX, W, WRT, WY, YN1
      DOUBLE PRECISION HALF, ONE, ZERO
      DATA HALF/0.5D+0/, ONE/1.D+0/, ZERO/0.D+0/
C
C *** BODY ***
C
      CALL DV7SCP(LC, C, ZERO)
      CALL DV7SCP(P, QTR, ZERO)
      DO 30 I = 1, N
         W = YN(2,I)
         IF (W .LE. ZERO) GO TO 40
         WRT =  SQRT(W)
         YN1 = YN(1,I) / YN(2,I)
         IF (MODEL .EQ. 2) GO TO 10
            WY = WRT * YN1
            GO TO 20
 10      WY = WRT * DLOG(  MAX(YN1, HALF/W))
 20      CALL DV7SCL(P, X, WRT, A(1,I))
         CALL DQ7ADR(P, QTR, C, X, WY)
 30      CONTINUE
      SX = DL7SVX(P, C, X, X)
      IF (SX .LE. ZERO) GO TO 40
      IF (DL7SVN(P, C, X, X)/SX .LE. DR7MDC(3)) GO TO 40
      CALL DL7ITV(P, X, C, QTR)
      GO TO 999
 40   W = ONE
      IF (MODEL .EQ. 2) W = ZERO
      CALL DV7SCP(P, X, W)
C
 999  RETURN
      END
      SUBROUTINE POIX0(A, IV, LA, LIV, LV, MODEL, N, P, V, X, YN)
C
C *** COMPUTE INITIAL X OF E. L. FROME ***
C
      INTEGER LA, LIV, LV, MODEL, N, P
      INTEGER IV(LIV)
      DOUBLE PRECISION X(P), A(LA,N), V(LV), YN(2,N)
C
      EXTERNAL DIVSET, POISX0, DV7SCP
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER C1, PP1O2, QTR1, TEMP1
      DOUBLE PRECISION ONE, ZERO
C
C  ***  IV COMPONENTS  ***
C
      INTEGER LMAT
      PARAMETER (LMAT=42)
      DATA ONE/1.D+0/, ZERO/0.D+0/
C
C---------------------------------  BODY  ------------------------------
C
      IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V)
C
      C1 = IV(LMAT)
      PP1O2 = P * (P + 1) / 2
      QTR1 = C1 + PP1O2
      TEMP1 = QTR1 + P
      IF (TEMP1 .GT. LV) GO TO 10
      CALL POISX0(A, V(C1), LA, P*(P+1)/2, MODEL, N, P, V(QTR1), X, YN)
      GO TO 999
C
 10   IF (MODEL .GT. 1) GO TO 20
      CALL DV7SCP(P, X, ONE)
      GO TO 999
 20   CALL DV7SCP(P, X, ZERO)
C
 999   RETURN
       END
      SUBROUTINE PREGRH(DERIV, F, N, NF, PT, R, RD, RHOI, YLOG, YN, ZN)
C
C  ***  RHO FOR PREGIBON ERROR MODELS WITH REPLICATE WEIGHTS ***
C  ***  SEE PREGRV FOR THE RIGHT WEIGHTING FOR THE INSURANCE EXAMPLE ***
C
      INTEGER DERIV, N, NF, RHOI(*)
      DOUBLE PRECISION F, PT(3), R(*), RD(*), YLOG(*), YN(2,N), ZN(3,N)
      EXTERNAL DR7MDC
      DOUBLE PRECISION DR7MDC
C
C *** LOCAL VARIABLES ***
C
      INTEGER I, K, KMP, KMPS, KMT, KPP, KPPS, KPSPS, KPT, KTT, KTPS
      DOUBLE PRECISION F1, MU, PHI, PHII2, PHII3, PHIINV, PSI, PSPHII,
     1                 RI, RL, RP0, RPP0, RT1, RT1L, RT2, RT2L, RTOL, T,
     2                 T1, T1INV, T1INV2, T2, T2INV, T2INV2, THETA, TT,
     3                 WI, WOVPHI, YI, YL, YT1, YT1L, YT2, YT2L
C
      DOUBLE PRECISION BIG, BIGH, TWOPI
      DOUBLE PRECISION BTOL, EIGHT, HALF, ONE, THREE, TWO, ZERO
      DATA BIG/0.D+0/, BIGH/0.D+0/, TWOPI/0.D+0/
      DATA BTOL/1.01D+0/, EIGHT/8.D+0/, HALF/0.5D+0/, ONE/1.D+0/,
     1     THREE/3.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C
C *** BODY ***
C
      IF (NF .GT. 1) GO TO 20
      IF (DERIV .GT. 0) GO TO 20
      DO 10 I = 1, N
 10      YLOG(I) = DLOG(YN(1,I))
 20   PHI = PT(1)
      PSI = PT(3)
      IF (PHI .LE. ZERO) GO TO 240
      THETA = PT(2)
      IF (TWOPI .GT. ZERO) GO TO 30
         TWOPI = EIGHT * DATAN(ONE)
         BIGH = DR7MDC(5)
         BIG = DR7MDC(6)
 30   T2 = TWO - THETA
      T1 = ONE - THETA
      IF (DERIV .GT. 0) GO TO 120
      RTOL = BIG
      IF (T2 .LT. BTOL) GO TO 40
         RTOL = BIGH**(ONE/T2)
         RTOL = RTOL*RTOL
 40   T = DLOG(TWOPI * PHI)
      F = ZERO
      DO 50 I = 1, N
 50      F = F + YN(2,I)*(T + THETA*YLOG(I))
      F1 = ZERO
      IF (THETA .EQ. ONE) GO TO 70
      IF (THETA .EQ. TWO) GO TO 90
      T1INV = ONE / T1
      T2INV = ONE / T2
      DO 60 I = 1, N
         RI = R(I)
         IF (RI .GE. RTOL) GO TO 240
         IF (RI .LE. ZERO) GO TO 240
         YI = YN(1,I)
         RT1 = RI**(T1*PSI)
         ZN(2,I) = RT1
         YT1 = YI**T1
         ZN(3,I) = YT1
         T = T2INV*(RI**(T2*PSI) - YI*YT1) + YI*T1INV*(YT1 - RT1)
         F1 = F1 + T*YN(2,I)
         ZN(1,I) = T
 60      CONTINUE
      GO TO 110
C
C *** THETA == 1 ***
C
 70   DO 80 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 240
         MU = RI**PSI
         YI = YN(1,I)
         T = MU - YI - YI*DLOG(MU/YI)
         F1 = F1 + T*YN(2,I)
         ZN(1,I) = T
         ZN(2,I) = ONE
 80      CONTINUE
      GO TO 110
C
C *** THETA == 2 ***
C
 90   DO 100 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 240
         T1 = RI**(-PSI)
         YI = YN(1,I) * T1
         T = YI - DLOG(YI) - ONE
         F1 = F1 + T*YN(2,I)
         ZN(1,I) = T
         ZN(2,I) = T1
 100     CONTINUE
 110  F = HALF*F + F1/PHI
      GO TO 999
C
C  ***  GRADIENT COMPUTATIONS  ***
C
 120  PHIINV = ONE / PHI
      PHII2 = PHIINV * PHIINV
      RP0 = HALF * PHIINV
      RPP0 = -PHIINV * RP0
      PHII3 = TWO * PHIINV * PHII2
      KMP = N
      KPP = N + N
      T1 = ONE - THETA
      T2 = TWO - THETA
      IF (RHOI(2) .LE. RHOI(3)+2) GO TO 140
C
C  *** PSI DERIVATIVES ***
C
      K = KPP + N
      KMPS = 6*N
      KPPS = KMPS + N
      KTPS = KPPS + N
      KPSPS = KTPS + N
      DO 130 I = 1, N
         WI = YN(2,I)
         RI = R(I)
         MU = RI**PSI
         RL = DLOG(RI)
         RT1 = WI * ZN(2,I)
         RT2 = RT1 * MU
         YI = YN(1,I)
         T = (RL/PHI) * (RT2 - YI*RT1)
         K = K + 1
         R(K) = T
         KMPS = KMPS + 1
         TT = RL * (T2*RT2 - YI*T1*RT1)
         RD(KMPS) = (RT2 - YI*RT1 + PSI*TT) / (RI*PHI)
         KPPS = KPPS + 1
         RD(KPPS) = -T / PHI
         KTPS = KTPS + 1
         RD(KTPS) = -PSI * RL * T
         KPSPS = KPSPS + 1
         RD(KPSPS) = TT * RL / PHI
 130     CONTINUE
C
 140  IF (RHOI(2) .LE. RHOI(3)) GO TO 220
      IF (RHOI(2) .EQ. RHOI(3)+1) GO TO 200
C
C  *** THETA DERIVATIVES ***
C
      K = KPP
      KMT = K + N
      KPT = KMT + N
      KTT = KPT + N
      IF (THETA .EQ. ONE) GO TO 160
      IF (THETA .EQ. TWO) GO TO 180
      T1INV = ONE / T1
      T1INV2 = T1INV + T1INV
      T2INV = ONE / T2
      T2INV2 = T2INV + T2INV
      DO 150 I = 1, N
         WI = YN(2,I)
         WOVPHI = WI * PHIINV
         RI = R(I)
         MU = RI**PSI
         RT1 = ZN(2,I)
         RT2 = RT1 * MU
         RL = DLOG(MU)
         RT1L = RT1 * RL
         RT2L = RT2 * RL
         YI = YN(1,I)
         YT1 = ZN(3,I)
         YT2 = YT1 * YI
         YL = YLOG(I)
         YT1L = YT1 * YL
         YT2L = YT2 * YL
         T = PHIINV * (YI * T1INV * (RL*RT1 - YL*YT1 +
     1                          T1INV*(YT1 - RT1))
     2                  + T2INV * (YL*YT2 - RL*RT2 +
     3                          T2INV*(RT2 - YT2)))
         K = K + 1
         R(K) = WI * (HALF*YL + T)
         KMT = KMT + 1
         RD(KMT) = PSI * WOVPHI * RL * (YI*RT1 - RT2) / RI
         KPT = KPT + 1
         RD(KPT) = -WOVPHI * T
         KTT = KTT + 1
         RD(KTT) = WOVPHI*(T1INV*YI*(YT1L*YL - RT1L*RL +
     1                       T1INV2*(RT1L - YT1L +
     2                        T1INV*(YT1 - RT1))) +
     3                        T2INV*(RT2L*RL - YT2L*YL +
     4                       T2INV2*(YT2L - RT2L +
     5                        T2INV*(RT2 - YT2))))
 150     CONTINUE
      GO TO 200
C
C *** THETA DERIVATIVES AT THETA == 1 ***
C
 160  DO 170 I = 1, N
         WI = YN(2,I)
         WOVPHI = WI * PHIINV
         YI = YN(1,I)
         YL = YLOG(I)
         RI = R(I)
         MU = RI**PSI
         RL = DLOG(MU)
         K = K + 1
         T = HALF*YI*(RL*RL - YL*YL) + YI*YL - MU*RL + MU - YI
         R(K) =  WI*(HALF*YL + T)
         KMT = KMT + 1
         RD(KMT) = PSI * WOVPHI * RL * (YI - MU) / RI
         KPT = KPT + 1
         RD(KPT) = -WOVPHI * T
         KTT = KTT + 1
         T = RL * RL
         RD(KTT) = WOVPHI * ( MU * (TWO - TWO*RL + T)
     1          -YI*(TWO - T*RL/THREE + YL*(YL - TWO + YL*YL/THREE)))
 170     CONTINUE
      GO TO 200
C
C *** THETA DERIVATIVES AT THETA == 2 ***
C
 180  DO 190 I = 1, N
         WI = YN(2,I)
         WOVPHI = WI * PHIINV
         YI = YN(1,I)
         YL = YLOG(I)
         RI = R(I)
         MU = RI**PSI
         RL = DLOG(MU)
         K = K + 1
         T = HALF*(YL*YL - RL*RL) + YL + ONE - (YI + YI*RL)/MU
         R(K) =  WI*(HALF*YL + T)
         KMT = KMT + 1
         RD(KMT) = PSI * WOVPHI * RL * (YI/MU - ONE) / RI
         KPT = KPT + 1
         RD(KPT) = -WOVPHI * T
         KTT = KTT + 1
         T = RL * RL
         RD(KTT) = WOVPHI * ((YL/MU)*(T + TWO*RL + TWO) - TWO
     1                  - YL*(TWO + YL*(ONE + YL/THREE)) + T*RL/THREE)
 190     CONTINUE
C
C *** PHI AND MU DERIVATIVES ***
C
 200  K = N
      THETA = ONE - PSI*T1
      T1 = PSI*T2 - ONE
      PSPHII = PSI * PHIINV
      PHIINV = -PHIINV
      DO 210 I = 1, N
         WI = YN(2,I)
         WOVPHI = WI * PSPHII
         RI = R(I)
         MU = RI**PSI
         YI = YN(1,I)
         RT1 = ZN(2,I)/RI
         T2 = WOVPHI * RT1 * (MU - YI)
         R(I) = T2
         RD(I) = WOVPHI * RT1 * (T1*MU + YI*THETA) / RI
         T = ZN(1,I)
         K = K + 1
         R(K) = WI * (RP0 - PHII2*T)
         KMP = KMP + 1
         RD(KMP) = PHIINV * T2
         KPP = KPP + 1
         RD(KPP) = WI * (RPP0 + PHII3*T)
 210     CONTINUE
      GO TO 999
C
C *** JUST MU DERIVATIVES ***
C
 220  PHIINV = PHIINV * PSI
      THETA = ONE - PSI*T1
      T1 = PSI*T2 - ONE
      DO 230 I = 1, N
         WOVPHI = YN(2,I) * PHIINV
         RI = R(I)
         MU = RI**PSI
         YI = YN(1,I)
         RT1 = ZN(2,I)/RI
         R(I) = WOVPHI * RT1 * (MU - YI)
         RD(I) = WOVPHI * RT1 * (T1*MU + YI*THETA) / RI
 230     CONTINUE
      GO TO 999
C
 240  NF = 0
C
 999  RETURN
      END
      SUBROUTINE PREGRV(DERIV, F, N, NF, PT, R, RD, RHOI, YLOG, YN, ZN)
C
C  ***  RHO FOR PREGIBON ERROR MODELS WITH VARIANCE WEIGHTS ***
C
      INTEGER DERIV, N, NF, RHOI(*)
      DOUBLE PRECISION F, PT(3), R(*), RD(*), YLOG(N+2),YN(2,N),ZN(3,N)
      EXTERNAL DR7MDC
      DOUBLE PRECISION DR7MDC
C
C *** LOCAL VARIABLES ***
C
      INTEGER I, K, KMP, KMPS, KMT, KPP, KPPS, KPSPS, KPT, KTT, KTPS
      DOUBLE PRECISION F1, MU, PHI, PHII2, PHII3, PHIINV, PSI, PSPHII,
     1                 RI, RL, RP0, RPP0, RT1, RT1L, RT2, RT2L, RTOL, T,
     2                 T1, T1INV, T1INV2, T2, T2INV, T2INV2, THETA, TT,
     3                 WI, WOVPHI, YI, YL, YT1, YT1L, YT2, YT2L
C
      DOUBLE PRECISION BIG, BIGH, TWOPI
      DOUBLE PRECISION BTOL, EIGHT, HALF, ONE, THREE, TWO, ZERO
      DATA BIG/0.D+0/, BIGH/0.D+0/, TWOPI/0.D+0/
      DATA BTOL/1.01D+0/, EIGHT/8.D+0/, HALF/0.5D+0/, ONE/1.D+0/,
     1     THREE/3.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C
C *** BODY ***
C
      PHI = PT(1)
      IF (PHI .LE. ZERO) GO TO 230
      IF (TWOPI .GT. ZERO) GO TO 10
         TWOPI = EIGHT * DATAN(ONE)
         BIGH = DR7MDC(5)
         BIG = DR7MDC(6)
 10   IF (NF .GT. 1) GO TO 30
      IF (DERIV .GT. 0) GO TO 30
      T1 = ZERO
      T2 = ZERO
      DO 20 I = 1, N
         T = DLOG(YN(1,I))
         YLOG(I) = T
         T1 = T1 + T
         T2 = T2 + DLOG(YN(2,I))
 20      CONTINUE
      YLOG(N+1) = T1
      YLOG(N+2) = -T2
 30   PSI = PT(3)
      THETA = PT(2)
      T2 = TWO - THETA
      T1 = ONE - THETA
      IF (DERIV .GT. 0) GO TO 110
      RTOL = BIG
      IF (T2 .LT. BTOL) GO TO 40
         RTOL = BIGH**(ONE/T2)
         RTOL = RTOL*RTOL
 40   F = N*DLOG(TWOPI*PHI) + YLOG(N+2) + THETA*YLOG(N+1)
      F1 = ZERO
      IF (THETA .EQ. ONE) GO TO 60
      IF (THETA .EQ. TWO) GO TO 80
      T1INV = ONE / T1
      T2INV = ONE / T2
      DO 50 I = 1, N
         RI = R(I)
         IF (RI .GE. RTOL) GO TO 230
         IF (RI .LE. ZERO) GO TO 230
         YI = YN(1,I)
         RT1 = RI**(T1*PSI)
         ZN(2,I) = RT1
         YT1 = YI**T1
         ZN(3,I) = YT1
         T = T2INV*(RI**(T2*PSI) - YI*YT1) + YI*T1INV*(YT1 - RT1)
         F1 = F1 + T*YN(2,I)
         ZN(1,I) = T
 50      CONTINUE
      GO TO 100
C
C *** THETA == 1 ***
C
 60   DO 70 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 230
         MU = RI**PSI
         YI = YN(1,I)
         T = MU - YI - YI*DLOG(MU/YI)
         F1 = F1 + T*YN(2,I)
         ZN(1,I) = T
         ZN(2,I) = ONE
 70      CONTINUE
      GO TO 100
C
C *** THETA == 2 ***
C
 80   DO 90 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 230
         T1 = RI**(-PSI)
         YI = YN(1,I) * T1
         T = YI - DLOG(YI) - ONE
         F1 = F1 + T*YN(2,I)
         ZN(1,I) = T
         ZN(2,I) = T1
 90      CONTINUE
 100  F = HALF*F + F1/PHI
      GO TO 999
C
C  ***  GRADIENT COMPUTATIONS  ***
C
 110  PHIINV = ONE / PHI
      PHII2 = PHIINV * PHIINV
      RP0 = HALF * PHIINV
      RPP0 = -PHIINV * RP0
      PHII3 = TWO * PHIINV * PHII2
      KMP = N
      KPP = N + N
      T1 = ONE - THETA
      T2 = TWO - THETA
      IF (RHOI(2) .LE. RHOI(3)+2) GO TO 130
C
C  *** PSI DERIVATIVES ***
C
      K = KPP + N
      KMPS = 6*N
      KPPS = KMPS + N
      KTPS = KPPS + N
      KPSPS = KTPS + N
      DO 120 I = 1, N
         WI = YN(2,I)
         RI = R(I)
         MU = RI**PSI
         RL = DLOG(RI)
         RT1 = WI * ZN(2,I)
         RT2 = RT1 * MU
         YI = YN(1,I)
         T = (RL/PHI) * (RT2 - YI*RT1)
         K = K + 1
         R(K) = T
         KMPS = KMPS + 1
         TT = RL * (T2*RT2 - YI*T1*RT1)
         RD(KMPS) = (RT2 - YI*RT1 + PSI*TT) / (RI*PHI)
         KPPS = KPPS + 1
         RD(KPPS) = -T / PHI
         KTPS = KTPS + 1
         RD(KTPS) = -PSI * RL * T
         KPSPS = KPSPS + 1
         RD(KPSPS) = TT * RL / PHI
 120     CONTINUE
C
 130  IF (RHOI(2) .LE. RHOI(3)) GO TO 210
      IF (RHOI(2) .EQ. RHOI(3)+1) GO TO 190
C
C  *** THETA DERIVATIVES ***
C
      K = KPP
      KMT = K + N
      KPT = KMT + N
      KTT = KPT + N
      IF (THETA .EQ. ONE) GO TO 150
      IF (THETA .EQ. TWO) GO TO 170
      T1INV = ONE / T1
      T1INV2 = T1INV + T1INV
      T2INV = ONE / T2
      T2INV2 = T2INV + T2INV
      DO 140 I = 1, N
         WI = YN(2,I)
         WOVPHI = WI * PHIINV
         RI = R(I)
         MU = RI**PSI
         RT1 = ZN(2,I)
         RT2 = RT1 * MU
         RL = DLOG(MU)
         RT1L = RT1 * RL
         RT2L = RT2 * RL
         YI = YN(1,I)
         YT1 = ZN(3,I)
         YT2 = YT1 * YI
         YL = YLOG(I)
         YT1L = YT1 * YL
         YT2L = YT2 * YL
         T = PHIINV * (YI * T1INV * (RL*RT1 - YL*YT1 +
     1                          T1INV*(YT1 - RT1))
     2                  + T2INV * (YL*YT2 - RL*RT2 +
     3                          T2INV*(RT2 - YT2)))
         K = K + 1
         R(K) = HALF*YL + WI*T
         KMT = KMT + 1
         RD(KMT) = PSI * WOVPHI * RL * (YI*RT1 - RT2) / RI
         KPT = KPT + 1
         RD(KPT) = -WOVPHI * T
         KTT = KTT + 1
         RD(KTT) = WOVPHI*(T1INV*YI*(YT1L*YL - RT1L*RL +
     1                       T1INV2*(RT1L - YT1L +
     2                        T1INV*(YT1 - RT1))) +
     3                        T2INV*(RT2L*RL - YT2L*YL +
     4                       T2INV2*(YT2L - RT2L +
     5                        T2INV*(RT2 - YT2))))
 140     CONTINUE
      GO TO 190
C
C *** THETA DERIVATIVES AT THETA == 1 ***
C
 150  DO 160 I = 1, N
         WI = YN(2,I)
         WOVPHI = WI * PHIINV
         YI = YN(1,I)
         YL = YLOG(I)
         RI = R(I)
         MU = RI**PSI
         RL = DLOG(MU)
         K = K + 1
         T = HALF*YI*(RL*RL - YL*YL) + YI*YL - MU*RL + MU - YI
         R(K) = HALF*YL + WI*T
         KMT = KMT + 1
         RD(KMT) = PSI * WOVPHI * RL * (YI - MU) / RI
         KPT = KPT + 1
         RD(KPT) = -WOVPHI * T
         KTT = KTT + 1
         T = RL * RL
         RD(KTT) = WOVPHI * ( MU * (TWO - TWO*RL + T)
     1          -YI*(TWO - T*RL/THREE + YL*(YL - TWO + YL*YL/THREE)))
 160     CONTINUE
      GO TO 190
C
C *** THETA DERIVATIVES AT THETA == 2 ***
C
 170  DO 180 I = 1, N
         WI = YN(2,I)
         WOVPHI = WI * PHIINV
         YI = YN(1,I)
         YL = YLOG(I)
         RI = R(I)
         MU = RI**PSI
         RL = DLOG(MU)
         K = K + 1
         T = HALF*(YL*YL - RL*RL) + YL + ONE - (YI + YI*RL)/MU
         R(K) = HALF*YL + WI*T
         KMT = KMT + 1
         RD(KMT) = PSI * WOVPHI * RL * (YI/MU - ONE) / RI
         KPT = KPT + 1
         RD(KPT) = -WOVPHI * T
         KTT = KTT + 1
         T = RL * RL
         RD(KTT) = WOVPHI * ((YL/MU)*(T + TWO*RL + TWO) - TWO
     1                  - YL*(TWO + YL*(ONE + YL/THREE)) + T*RL/THREE)
 180     CONTINUE
C
C *** PHI AND MU DERIVATIVES ***
C
 190  K = N
      THETA = ONE - PSI*T1
      T1 = PSI*T2 - ONE
      PSPHII = PSI * PHIINV
      PHIINV = -PHIINV
      DO 200 I = 1, N
         WI = YN(2,I)
         WOVPHI = WI * PSPHII
         RI = R(I)
         MU = RI**PSI
         YI = YN(1,I)
         RT1 = ZN(2,I)/RI
         T2 = WOVPHI * RT1 * (MU - YI)
         R(I) = T2
         RD(I) = WOVPHI * RT1 * (T1*MU + YI*THETA) / RI
         T = ZN(1,I)
         K = K + 1
         R(K) = RP0 - WI*PHII2*T
         KMP = KMP + 1
         RD(KMP) = PHIINV * T2
         KPP = KPP + 1
         RD(KPP) = RPP0 + WI*PHII3*T
 200     CONTINUE
      GO TO 999
C
C *** JUST MU DERIVATIVES ***
C
 210  PHIINV = PHIINV * PSI
      THETA = ONE - PSI*T1
      T1 = PSI*T2 - ONE
      DO 220 I = 1, N
         WOVPHI = YN(2,I) * PHIINV
         RI = R(I)
         MU = RI**PSI
         YI = YN(1,I)
         RT1 = ZN(2,I)/RI
         R(I) = WOVPHI * RT1 * (MU - YI)
         RD(I) = WOVPHI * RT1 * (T1*MU + YI*THETA) / RI
 220     CONTINUE
      GO TO 999
C
 230  NF = 0
C
 999  RETURN
      END
      SUBROUTINE PRGRH1(N, PT, R, RHO, RHOI, YN)
C
C  ***  RHO FOR PREGIBON ERROR MODELS ***
C
      INTEGER N, RHOI(3)
      DOUBLE PRECISION PT(2), R(*), RHO(N), YN(2,N)
C *** LOCAL VARIABLES ***
C
      INTEGER I
      DOUBLE PRECISION HTHETA, PHI, RI, RT1, T, T1, T1INV, T2, T2INV,
     1                  THETA, YI, YT1
C
      DOUBLE PRECISION HALF, ONE, TWO
      DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/
C
C *** BODY ***
C
      PHI = PT(1)
      THETA = PT(2)
      HTHETA = HALF * THETA
      DO 10 I = 1, N
 10      RHO(I) = HTHETA*DLOG(PHI*YN(1,I))
      IF (THETA .EQ. ONE) GO TO 30
      IF (THETA .EQ. TWO) GO TO 50
      T1 = ONE - THETA
      T1INV = ONE / T1 / PHI
      T2 = TWO - THETA
      T2INV = ONE / T2 / PHI
      DO 20 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         RT1 = RI**T1
         YT1 = YI**T1
         RHO(I) = RHO(I) + T2INV*(RI*RT1 - YI*YT1) + YI*T1INV*(YT1- RT1)
 20      CONTINUE
      GO TO 999
 30   DO 40 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         T = RI - YI - YI*DLOG(RI/YI)
         RHO(I) = RHO(I) + T / PHI
 40      CONTINUE
      GO TO 999
 50   DO 60 I = 1, N
         YI = YN(1,I) / R(I)
         T = YI - DLOG(YI) - ONE
         RHO(I) = RHO(I) + T / PHI
 60      CONTINUE
 999  RETURN
      END
      SUBROUTINE RHPOIL(NEED, F, N, NF, PT, R, RD, RHOI, YN, W)
      COMMON /FUDGE/ NFUDGE
      INTEGER NFUDGE
      INTEGER NEED(2), N, NF, RHOI(6)
      DOUBLE PRECISION F, PT(3), R(*), RD(*), W(N), YN(2,N)
C PT = PHI AND THETA (WHEN PS == P, I.E. RHOI(2) == RHOI(3))
C
      DOUBLE PRECISION INVCN, LPN, PNORMS, DR7MDC
      EXTERNAL INVCN, LPN, PNORMS, DR7MDC
      INTEGER ERRFLG, I, IM, WCOMP
      DOUBLE PRECISION CI, E, PHI, PHIRI, PHIMRI, PSI, PSI1, PSI2,
     1                 RI, T, T1, T2, THETA, YI
      DOUBLE PRECISION DATAN, DEXP, DLOG,  SQRT
      DOUBLE PRECISION CNN, EIGHT, EXPMAX, EXPMIN, FOUR, HALF, ONE, TWO,
     1                 TWOPI, ZERO
      DATA CNN/0.D+0/, EXPMAX/0.D+0/, EIGHT/8.D+0/, EXPMIN/0.D+0/,
     1     FOUR/4.0D+0/, HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/,
     2     TWOPI/0.D+0/, ZERO/0.D+0/
C
C *** BODY ***
C
      IM = RHOI(1)
      WCOMP = RHOI(6)
      IF (IM .LE. 0) GO TO 800
      IF (IM .GT. 13) GO TO 800
      IF (EXPMAX .GT. ZERO) GO TO 10
         EXPMAX = TWO * DLOG(DR7MDC(5))
         EXPMIN = TWO * DLOG(DR7MDC(2))
         TWOPI = EIGHT * DATAN(ONE)
 10   IF (NEED(1) .EQ. 2) GO TO 240
      F = ZERO
      GO TO (20,20,40,60,80,80,100,120,140,160,180,220,180), IM
C
C *** POISSON RHO (AND CONVENTIONAL IRLS) ***
C
 20   DO 30 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         F = F + YN(2,I)*RI - YN(1,I)*DLOG(RI)
 30      CONTINUE
      GO TO 999
C
C *** LOG LINEAR POISSON ***
C
 40   DO 50 I = 1, N
         E = ZERO
         RI = R(I)
         IF (RI .GT. EXPMAX) GO TO 800
         IF (RI .GT. EXPMIN) E = EXP(RI)
         F = F + YN(2,I)*E - YN(1,I)*RI
         R(I) = E
 50      CONTINUE
      GO TO 999
C
C *** SQUARE-ROOT LINEAR POISSON ***
C
 60   DO 70 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         F = F + YN(2,I)*RI**2 - TWO*YN(1,1)*DLOG(RI)
 70      CONTINUE
      GO TO 999
C
C *** BINOMIAL RHO (AND CONVENTIONAL IRLS) ***
C
 80   DO 90 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         IF (RI .GE. ONE) GO TO 800
         F = F - YN(1,I)*DLOG(RI) - (YN(2,I) - YN(1,I))*DLOG(ONE-RI)
 90      CONTINUE
      GO TO 999
C
C *** BINOMIAL LOGISTIC RHO ***
C
 100  DO 110 I = 1, N
         RI = R(I)
         IF (RI .GE. EXPMAX) GO TO 800
         E = ZERO
         IF (RI .GT. EXPMIN) E = DEXP(RI)
         F = F + YN(2,I)*DLOG(ONE + E) - YN(1,I)*RI
         R(I) = E
 110     CONTINUE
      GO TO 999
C
C *** PROBIT ***
C
 120  DO 130 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         F = F - YI*LPN(RI) - (YN(2,I)-YI)*LPN(-RI)
 130     CONTINUE
        IF (NFUDGE .GT. 0) WRITE(*,*) 'NFUDGE =', NFUDGE
        NFUDGE = 0
      GO TO 999
C
C *** WEIBULL ***
C
 140  DO 150 I = 1, N
         RI = R(I)
         IF (RI .GE. EXPMAX) GO TO 800
         E = ZERO
         IF (RI .GT. EXPMIN) E = DEXP(RI)
         R(I) = E
         T = ZERO
         IF (-E .GT. EXPMIN) T = DEXP(-E)
         F = F + (YN(2,I) - YN(1,I))*E - YN(1,I)*DLOG(ONE - T)
 150     CONTINUE
      GO TO 999
C
C  *** GAMMA ERRORS ***
C
 160  DO 170 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         F = F + YN(1,I)*RI - YN(2,I)*DLOG(RI)
 170     CONTINUE
      GO TO 999
C
C  ***  PREGIBON ERRORS ***
C
C      *** IN THIS CASE, YN(1,I) = Y(I), YN(2,I) = LOG(Y(I))
C      *** AND YN(I,J), J = N+1(1)2*N, I = 1 OR 2 = SCRATCH
C
 180  IF (NF .GT. 1) GO TO 190
      RHOI(4) = 0
      RHOI(5) = 0
 190  I = N + N + 3
C     *** THE YLOG ARRAY PASSED TO PREGRV MUST BE AT LEAST N+2 LONG
      IF (NEED(2) .NE. RHOI(4)) GO TO 200
         I = I + 3*N
         RHOI(5) = NF
         GO TO 210
 200  RHOI(4) = NF
 210  IF (IM .EQ. 11) THEN
        CALL PREGRH(0, F, N, NF, PT, R, RD, RHOI, YN(1,N+1), YN,
     1            YN(1,I))
      ELSE
        CALL PREGRV(0, F, N, NF, PT, R, RD, RHOI, YN(1,N+1), YN,
     1            YN(1,I))
        END IF
      GO TO 999
C
C *** LEAST-SQUARES ***
C
 220  DO 230 I = 1, N
        E = R(I) - YN(1,I)
        F = F + E*E
 230    CONTINUE
      F = HALF * F
      GO TO 999
C
 240  GO TO (250,270,310,350,400,420,460,500,570,620,660,780,660), IM
C
C *** IRLS POISSON DERIVATIVES ***
C
 250  DO 260 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         R(I) = YN(2,I) - YN(1,I) / RI
         RD(I) = YN(2,I) / RI
 260     CONTINUE
      GO TO 820
C
C *** POISSON DERIVATIVES ***
C
 270  DO 300 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         YI = YN(1,I)
         CI = YN(2,I)
         E = YI / RI
         R(I) = CI - E
         RD(I) = E / RI
         GO TO (300, 280, 280, 290), WCOMP
 280     W(I) = CI / RI
         GO TO 300
 290     IF (YI .LE. ZERO) THEN
             W(I) = HALF * CI / RI
         ELSE
            T1 = CI*RI + YI*(DLOG(E/CI) - ONE)
            IF (T1 .NE. ZERO) THEN
               T = R(I)
               W(I) = T*T / (T1+T1)
            ELSE
               W(I) = RD(I)
               END IF
            END IF
 300     CONTINUE
      GO TO 810
C
C *** LOG LINEAR POISSON ***
C
 310  DO 340 I = 1, N
         YI = YN(1,I)
         CI = YN(2,I)
         RI = CI*R(I)
         R(I) = RI - YI
         RD(I) = RI
         GO TO (340,340,320,330), WCOMP
 320     T = RI/YI
         IF (T .EQ. ONE) THEN
            W(I) = YI
         ELSE
            W(I) = YI * ((T - ONE) / DLOG(T))
            ENDIF
         GO TO 340
 330     T1 = RI + YI*(DLOG(YI/RI) - ONE)
         IF (T1 .NE. ZERO) THEN
            T = RI - YI
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 340     CONTINUE
      IF (WCOMP .LE. 2) GO TO 820
      GO TO 999
C
C *** SQUARE-ROOT LINEAR POISSON ***
C
 350  DO 390 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         YI = YN(1,I)
         CI = YN(2,I)
         E = YI / RI
         R(I) = TWO * (CI*RI - E)
         RD(I) = TWO * (CI + E/RI)
         GO TO (390, 360, 370, 380), WCOMP
 360     W(I) = FOUR * CI
         GO TO 390
 370     T1 = RI -  SQRT(YI/CI)
         IF (T1 .NE. ZERO) THEN
            T = CI*RI - YI/RI
            W(I) = (T+T) / T1
         ELSE
            W(I) = RD(I)
            END IF
         GO TO 390
 380     T1 = CI*RI*RI - YI + YI*DLOG(YI/(CI*RI*RI))
         IF (T1 .NE. ZERO) THEN
            T = CI*RI - YI/RI
            T = T / T1
            W(I) = T + T
         ELSE
            W(I) = RD(I)
            END IF
 390     CONTINUE
      GO TO 810
C
C *** IRLS BINOMIAL ***
C
 400  DO 410 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         IF (RI .GE. ONE) GO TO 800
         YI = YN(1,I)
         CI = YN(2,I)
         T = ONE / (ONE - RI)
         R(I) = (CI - YI) * T  -  YI / RI
         RD(I) = T * CI / RI
 410     CONTINUE
      GO TO 820
C
C *** BINOMIAL ***
C
 420  DO 450 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         IF (RI .GE. ONE) GO TO 800
         YI = YN(1,I)
         T = ONE / (ONE - RI)
         CI = (YN(2,I) - YI) * T
         YI = YI / RI
         R(I) = CI - YI
         RD(I) = T*CI + YI/RI
         GO TO (450,430,430,440), WCOMP
 430     W(I) = T*YN(2,I) / RI
         GO TO 450
 440     YI = YN(1,I)
         CI = YN(2,I)
         T2 = YI / CI
         T1 = (YI - CI)*DLOG((ONE - RI)/(ONE - T2)) + YI*DLOG(T2/RI)
         IF (T1 .NE. ZERO) THEN
            T = (CI*RI - YI)/(RI * (ONE - RI))
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 450     CONTINUE
      GO TO 810
C
C *** BINOMIAL LOGISTIC ***
C
 460  DO 490 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         CI = YN(2,I)
         T = ONE / (ONE + RI)
         T1 = T * RI * CI
         R(I) = T1 - YI
         RD(I) = T * T1
         GO TO (490,490,470,480), WCOMP
 470     T1 = (ONE + RI)*DLOG(RI*(CI-YI)/YI)
         IF (T1 .NE. ZERO) THEN
            W(I) = ((CI - YI)*RI - YI) / T1
         ELSE
            W(I) = RD(I)
            END IF
         GO TO 490
 480     T1 = CI*DLOG((ONE+RI)*(ONE - YI/CI)) + YI*DLOG(YI/(RI*(CI-YI)))
         IF (T1 .NE. ZERO) THEN
            T = ((CI - YI)*RI - YI) / (ONE + RI)
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 490     CONTINUE
      IF (WCOMP .LE. 2) GO TO 820
      GO TO 999
C
C *** PROBIT ***
C
 500  IF (CNN .LE. ZERO) CNN = ONE /  SQRT(TWOPI)
      DO 560 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         CI = YN(2,I) - YI
         E = ZERO
         T = -HALF * RI**2
         IF (T .GT. EXPMIN) E = CNN * DEXP(T)
         PHIRI = PNORMS(RI)
         IF (WCOMP .EQ. 2)
     1          W(I) = YN(2,I) * (E / PHIRI) * (E / (ONE - PHIRI))
         IF (PHIRI .LE. ZERO) GO TO 510
            PHIRI = ONE / PHIRI
            T1 = E*PHIRI*YI
            T2 = T1*(RI + PHIRI*E)
            T1 = -T1
            GO TO 520
 510     T1 = YI * (RI + ONE/RI)
         T2 = YI * (ONE - ONE/RI**2)
 520     PHIMRI = PNORMS(-RI)
         IF (PHIMRI .LE. ZERO) GO TO 530
            PHIMRI = ONE / PHIMRI
            T = E*CI*PHIMRI
            R(I) = T + T1
            RD(I) = T*(PHIMRI*E - RI) + T2
            GO TO (560,560,540,550), WCOMP
 530     R(I) = CI*(RI + ONE/RI) + T1
         RD(I) = CI*(ONE - ONE/RI**2) + T2
         GO TO (560,560,540,550), WCOMP
 540     T = RI - INVCN(YI/YN(2,I), ERRFLG)
         IF (ERRFLG .NE. 0) THEN
            WRITE(*,*) 'ERROR FROM INVCN: I, YI, YN(1,I), YN(2,I) ='
     1                  , I, YI, YN(1,I), YN(2,I)
             GO TO 800
             END IF
         IF (T .NE. ZERO) THEN
             W(I) = R(I) / T
         ELSE
             W(I) = RD(I)
             END IF
         GO TO 560
 550     T2 = CI
         CI = YN(2,I)
         T1 = T2*(DLOG(T2/CI) - LPN(-RI))
         IF (YI .GT. ZERO) T1 = T1 + YI*(DLOG(YI/CI) - LPN(RI))
         IF (T1 .NE. ZERO) THEN
             T = R(I)
             W(I) = T*T / (T1+T1)
         ELSE
             W(I) = RD(I)
             END IF
 560     CONTINUE
      GO TO 810
C
C *** WEIBULL ***
C
 570  DO 610 I = 1, N
         RI = R(I)
         E = ZERO
         IF (-RI .GT. EXPMIN) E = DEXP(-RI)
         T = RI / (ONE - E)
         CI = YN(2,I)*RI
         YI = YN(1,I)*T
         R(I) = CI - YI
         RD(I) = CI - YI*(ONE - E*T)
         GO TO (570,580,590,600), WCOMP
 580     W(I) = E*CI*RI / (ONE - E)
         GO TO 610
 590     T1 = DLOG(-RI / DLOG(ONE - YN(1,I)/YN(2,I)))
         IF (T1 .NE. ZERO) THEN
            W(I) = (CI - YI) / T1
         ELSE
            W(I) = RD(I)
            END IF
         GO TO 610
 600     YI = YN(1,I)
         CI = YN(2,I)
         T2 = YI / CI
         CI = CI - YI
         T1 = CI*(RI + DLOG(ONE - T2)) + YI*(DLOG(T2/(ONE - E)))
         IF (T1 .NE. ZERO) THEN
            T = CI - YI
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 610     CONTINUE
      GO TO 810
C
C  *** GAMMA ERRORS ***
C
 620  DO 650 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
C        F = F + YN(1,I)*RI - YN(2,I)*DLOG(RI)
         T = YN(2,I)/RI
         T1 = ONE
         R(I) = YN(1,I) - T
         RD(I) = T/RI
         GO TO (650,650,630,640), WCOMP
 630     W(I) = YN(1,I) / RI
         GO TO 650
 640     T2 = YN(1,I) * RI / YN(2,I)
         T1 = T2 - ONE
         T = T1*RD(I)*T1
         IF (T .GT. ZERO) THEN
            T2 = T1 - DLOG(T2)
            T = T / (T2+T2)
            END IF
         W(I) = T
 650     CONTINUE
      IF (WCOMP .LE. 2) GO TO 820
      GO TO 999
C
C ***  PREGIBON ERRORS ***
C
 660  IF (WCOMP .GE. 2) CALL DV7CPY(N, W, R)
      I = N + N + 3
      IF (RHOI(4) .EQ. NF) GO TO 670
         I = I + 3*N
         IF (RHOI(5) .EQ. NF) GO TO 670
         WRITE(6,*) 'HELP! NF =', NF, ' BUT RHOI =', RHOI
         GO TO 800
 670  IF (IM .EQ. 11) THEN
         CALL PREGRH(1, F, N, NF, PT, R, RD, RHOI, YN(1,N+1), YN,
     1            YN(1,I))
      ELSE
         CALL PREGRV(1, F, N, NF, PT, R, RD, RHOI, YN(1,N+1), YN,
     1            YN(1,I))
         END IF
      IF (NF .EQ. 0) GO TO 999
      GO TO (820,680,700,720), WCOMP
 680  PSI = PT(3)
      T = (TWO - PT(2))*PSI - TWO
      T1 = PSI*PSI
      DO 690 I = 1, N
 690        W(I) = YN(2,I) * T1 * W(I)**T
      GO TO 999
 700  T = ONE / PT(3)
      DO 710 I = 1, N
         T1 = W(I) - ONE
         IF (T1 .NE. ZERO) THEN
            YI = YN(1,I)
            W(I) = R(I) / (W(I) - YI**T)
         ELSE
            W(I) = RD(I)
            END IF
 710     CONTINUE
      GO TO 999
 720  PHI = PT(1)
      THETA = PT(2)
      PSI = PT(3)
      IF (THETA .EQ. ONE) GO TO 740
      IF (THETA .EQ. TWO) GO TO 760
      T1 = ONE - THETA
      T2 = TWO - THETA
      PSI1 = PSI * T1
      PSI2 = PSI * T2
      DO 730 I = 1, N
         RI = W(I)
         YI = YN(1,I)
         T = YI**T2
         E = YN(2,I)/PHI * ((T - YI*RI**PSI1)/T1 - (T - RI**PSI2)/T2)
         IF (E .NE. ZERO) THEN
            T = R(I)
            W(I) = T*T / (E+E)
         ELSE
            W(I) = RD(I)
            END IF
 730     CONTINUE
      GO TO 999
 740  DO 750 I = 1, N
         RI = W(I)
         YI = YN(1,I)
         T1 = YN(2,I)/PHI * (RI**PSI - YI + YI*(DLOG(YI)-PSI*DLOG(RI)))
         IF (T1 .NE. ZERO) THEN
            T = R(I)
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 750     CONTINUE
      GO TO 999
 760  DO 770 I = 1, N
         RI = W(I)
         YI = YN(1,I)
         T1 = YI*RI**(-PSI) - ONE + PSI*DLOG(RI) - DLOG(YI)
         IF (T1 .NE. ZERO) THEN
            T1 = T1 * YN(2,I) / PHI
            T = R(I)
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 770     CONTINUE
      GO TO 999
C
C *** LEAST SQUARES ***
C
 780  DO 790 I = 1, N
         R(I) = R(I) - YN(1,I)
         RD(I) = ONE
 790     CONTINUE
      GO TO 820
C
 800  NF = 0
      GO TO 999
C
 810  IF (WCOMP .GT. 1) GO TO 999
 820  CALL DV7CPY(N, W, RD)
C
 999  RETURN
      END
      DOUBLE PRECISION FUNCTION LPN(X)
      COMMON /FUDGE/ NFUDGE
      INTEGER NFUDGE
      DOUBLE PRECISION X
      EXTERNAL PNORMS
      DOUBLE PRECISION PNORMS
      DOUBLE PRECISION T
      DOUBLE PRECISION DLOG
      DOUBLE PRECISION HALF, ZERO
      DATA HALF/0.5D+0/, ZERO/0.D+0/
C
C *** BODY ***
C
      T = PNORMS(X)
      IF (T .GT. ZERO) THEN
         LPN = DLOG(T)
      ELSE
         NFUDGE = NFUDGE + 1
         LPN = -HALF*X**2 - DLOG(-X)
         END IF
 999  RETURN
      END

.