SUBROUTINE DG7ITB(B, D, G, IV, LIV, LV, P, PS, V, X, Y)
C
C  ***  CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR   ***
C  ***  REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE)     ***
C  ***  HAVING SIMPLE BOUNDS ON THE PARAMETERS BEING ESTIMATED.   ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, P, PS
      INTEGER IV(LIV)
      DOUBLE PRECISION B(2,P), D(P), G(P), V(LV), X(P), Y(P)
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C B.... VECTOR OF LOWER AND UPPER BOUNDS ON X.
C D.... SCALE VECTOR.
C IV... INTEGER VALUE ARRAY.
C LIV.. LENGTH OF IV.  MUST BE AT LEAST 80.
C LH... LENGTH OF H = P*(P+1)/2.
C LV... LENGTH OF V.  MUST BE AT LEAST P*(3*P + 19)/2 + 7.
C G.... GRADIENT AT X (WHEN IV(1) = 2).
C HC... GAUSS-NEWTON HESSIAN AT X (WHEN IV(1) = 2).
C P.... NUMBER OF PARAMETERS (COMPONENTS IN X).
C PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S.
C V.... FLOATING-POINT VALUE ARRAY.
C X.... PARAMETER VECTOR.
C Y.... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE).
C
C  ***  DISCUSSION  ***
C
C        DG7ITB IS SIMILAR TO DG7LIT, EXCEPT FOR THE EXTRA PARAMETER B
C     -- DG7ITB ENFORCES THE BOUNDS  B(1,I) .LE. X(I) .LE. B(2,I),
C     I = 1(1)P.
C        DG7ITB PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF
C     REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES
C     IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED
C     FIRST-ORDER TERM AND A SECOND-ORDER TERM.  THE CALLER SUPPLIES
C     THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED
C     COMPACTLY BY ROWS), AND DG7ITB BUILDS AN APPROXIMATION, S, TO THE
C     SECOND-ORDER TERM.  THE CALLER ALSO PROVIDES THE FUNCTION VALUE,
C     GRADIENT, AND PART OF THE YIELD VECTOR USED IN UPDATING S.
C     DG7ITB DECIDES DYNAMICALLY WHETHER OR NOT TO USE S WHEN CHOOSING
C     THE NEXT STEP TO TRY...  THE HESSIAN APPROXIMATION USED IS EITHER
C     HC ALONE (GAUSS-NEWTON MODEL) OR HC + S (AUGMENTED MODEL).
C     IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT
C     CONSTANT.  THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO
C     1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM, OR THE CALLER SETS
C     IV(INITS) TO 3 OR 4 AND THE FINITE-DIFFERENCE INITIAL S THEN
C     COMPUTED HAS NONZERO VALUES IN THESE ROWS.
C
C        IF IV(INITS) IS 3 OR 4, THEN THE INITIAL S IS COMPUTED BY
C     FINITE DIFFERENCES.  3 MEANS USE FUNCTION DIFFERENCES, 4 MEANS
C     USE GRADIENT DIFFERENCES.  FINITE DIFFERENCING IS DONE THE SAME
C     WAY AS IN COMPUTING A COVARIANCE MATRIX (WITH IV(COVREQ) = -1, -2,
C     1, OR 2).
C
C        FOR UPDATING S, DG7ITB ASSUMES THAT THE GRADIENT HAS THE FORM
C     OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE
C     GRADIENT WITH RESPECT TO X.  THE TRUE SECOND-ORDER TERM THEN IS
C     THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)).  IF X = X0 + STEP,
C     THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF
C     RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))).  THE CALLER MUST SUPPLY
C     PART OF THIS IN Y, NAMELY THE SUM OVER I OF
C     RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING DG7ITB WITH IV(1) = 2 AND
C     IV(MODE) = 0 (WHERE MODE = 38).  G THEN CONTANS THE OTHER PART,
C     SO THAT THE DESIRED YIELD VECTOR IS G - Y.  IF PS .LT. P, THEN
C     THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF
C     GRAD(R(I,X)), STEP, AND Y.
C
C        PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING
C     ONES TO  DN2GB (AND NL2SOL), EXCEPT THAT V CAN BE SHORTER
C     (SINCE THE PART OF V THAT  DN2GB USES FOR STORING D, J, AND R IS
C     NOT NEEDED).  MOREOVER, COMPARED WITH  DN2GB (AND NL2SOL), IV(1)
C     MAY HAVE THE TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE
C     EXPLAINED BELOW, AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL).
C     THE VALUES IV(D), IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM
C      DN2GB (AND  DN2FB), ARE NOT REFERENCED BY DG7ITB OR THE
C     SUBROUTINES IT CALLS.
C
C        WHEN DG7ITB IS FIRST CALLED, I.E., WHEN DG7ITB IS CALLED WITH
C     IV(1) = 0 OR 12, V(F), G, AND HC NEED NOT BE INITIALIZED.  TO
C     OBTAIN THESE STARTING VALUES, DG7ITB RETURNS FIRST WITH IV(1) = 1,
C     THEN WITH IV(1) = 2, WITH IV(MODE) = -1 IN BOTH CASES.  ON
C     SUBSEQUENT RETURNS WITH IV(1) = 2, IV(MODE) = 0 IMPLIES THAT
C     Y MUST ALSO BE SUPPLIED.  (NOTE THAT Y IS USED FOR SCRATCH -- ITS
C     INPUT CONTENTS ARE LOST.  BY CONTRAST, HC IS NEVER CHANGED.)
C     ONCE CONVERGENCE HAS BEEN OBTAINED, IV(RDREQ) AND IV(COVREQ) MAY
C     IMPLY THAT A FINITE-DIFFERENCE HESSIAN SHOULD BE COMPUTED FOR USE
C     IN COMPUTING A COVARIANCE MATRIX.  IN THIS CASE DG7ITB WILL MAKE
C     A NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE.
C     WHEN IV(MODE) IS POSITIVE, Y SHOULD NOT BE CHANGED.
C
C IV(1) = 1 MEANS THE CALLER SHOULD SET V(F) (I.E., V(10)) TO F(X), THE
C             FUNCTION VALUE AT X, AND CALL DG7ITB AGAIN, HAVING CHANGED
C             NONE OF THE OTHER PARAMETERS.  AN EXCEPTION OCCURS IF F(X)
C             CANNOT BE EVALUATED (E.G. IF OVERFLOW WOULD OCCUR), WHICH
C             MAY HAPPEN BECAUSE OF AN OVERSIZED STEP.  IN THIS CASE
C             THE CALLER SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL
C             CAUSE DG7ITB TO IGNORE V(F) AND TRY A SMALLER STEP.  NOTE
C             THAT THE CURRENT FUNCTION EVALUATION COUNT IS AVAILABLE
C             IN IV(NFCALL) = IV(6).  THIS MAY BE USED TO IDENTIFY
C             WHICH COPY OF SAVED INFORMATION SHOULD BE USED IN COM-
C             PUTING G, HC, AND Y THE NEXT TIME DG7ITB RETURNS WITH
C             IV(1) = 2.  SEE MLPIT FOR AN EXAMPLE OF THIS.
C IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT OF F AT
C             X.  THE CALLER SHOULD ALSO SET HC TO THE GAUSS-NEWTON
C             HESSIAN AT X.  IF IV(MODE) = 0, THEN THE CALLER SHOULD
C             ALSO COMPUTE THE PART OF THE YIELD VECTOR DESCRIBED ABOVE.
C             THE CALLER SHOULD THEN CALL DG7ITB AGAIN (WITH IV(1) = 2).
C             THE CALLER MAY ALSO CHANGE D AT THIS TIME, BUT SHOULD NOT
C             CHANGE X.  NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE
C             VALUE THAT IV(NFCALL) HAD DURING THE RETURN WITH
C             IV(1) = 1 IN WHICH X HAD THE SAME VALUE AS IT NOW HAS.
C             IV(NFGCAL) IS EITHER IV(NFCALL) OR IV(NFCALL) - 1.  MLPIT
C             IS AN EXAMPLE WHERE THIS INFORMATION IS USED.  IF G OR HC
C             CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET
C             IV(NFGCAL) TO 0, IN WHICH CASE DG7ITB WILL RETURN WITH
C             IV(1) = 15.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C
C        (SEE NL2SOL FOR REFERENCES.)
C
C+++++++++++++++++++++++++++  DECLARATIONS  ++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL HAVQTR, HAVRM
      INTEGER DUMMY, DIG1, G01, H1, HC1, I, I1, IPI, IPIV0, IPIV1,
     1        IPIV2, IPN, J, K, L, LMAT1, LSTGST, P1, P1LEN, PP1, PP1O2,
     2        QTR1, RMAT1, RSTRST, STEP1, STPMOD, S1, TD1, TEMP1, TEMP2,
     3        TG1, W1, WLM1, X01
      DOUBLE PRECISION E, GI, STTSST, T, T1, XI
C
C     ***  CONSTANTS  ***
C
      DOUBLE PRECISION HALF, NEGONE, ONE, ONEP2, ZERO
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      LOGICAL STOPX
      DOUBLE PRECISION DD7TPR, DRLDST, DV2NRM
      EXTERNAL DA7SST, DD7TPR, DF7DHB, DG7QSB,I7COPY, I7PNVR, I7SHFT,
     1        DITSUM, DL7MSB, DL7SQR, DL7TVM,DL7VML,DPARCK, DQ7RSH,
     2         DRLDST, DS7DMP, DS7IPR, DS7LUP, DS7LVM, STOPX, DV2NRM,
     3        DV2AXY,DV7CPY, DV7IPR, DV7SCP, DV7VMP
C
C DA7SST.... ASSESSES CANDIDATE STEP.
C DD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS.
C DF7DHB... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR INIT. S MATRIX).
C DG7QSB... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL).
C I7COPY.... COPIES ONE INTEGER VECTOR TO ANOTHER.
C I7PNVR... INVERTS PERMUTATION ARRAY.
C I7SHFT... SHIFTS AN INTEGER VECTOR.
C DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X.
C DL7MSB... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL).
C DL7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L.
C DL7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C DPARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS.
C DQ7RSH... SHIFTS A QR FACTORIZATION.
C DRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE.
C DS7DMP... MULTIPLIES A SYM. MATRIX FORE AND AFT BY A DIAG. MATRIX.
C DS7IPR... APPLIES PERMUTATION TO (LOWER TRIANG. OF) SYM. MATRIX.
C DS7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI-
C             ANGLE OF A SYMMETRIC MATRIX.
C DS7LVM... MULTIPLIES COMPACTLY STORED SYM. MATRIX TIMES VECTOR.
C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED.
C DV2NRM... RETURNS THE 2-NORM OF A VECTOR.
C DV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER.
C DV7CPY.... COPIES ONE VECTOR TO ANOTHER.
C DV7IPR... APPLIES A PERMUTATION TO A VECTOR.
C DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C DV7VMP... MULTIPLIES (DIVIDES) VECTORS COMPONENTWISE.
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER CNVCOD, COSMIN, COVMAT, COVREQ, DGNORM, DIG,
     1        DSTNRM, F, FDH, FDIF, FUZZ, F0, GTSTEP, H, HC, IERR,
     2        INCFAC, INITS, IPIVOT, IRC, IVNEED, KAGQT, KALM, LMAT,
     3        LMAX0, LMAXS, MODE, MODEL, MXFCAL, MXITER, NEXTIV, NEXTV,
     4        NFCALL, NFGCAL, NFCOV, NGCOV, NGCALL, NITER, NVSAVE, P0,
     5        PC, PERM, PHMXFC, PREDUC, QTR, RADFAC, RADINC, RADIUS,
     6        RAD0, RDREQ, REGD, RELDX, RESTOR, RMAT, S, SIZE, STEP,
     7        STGLIM, STPPAR, SUSED, SWITCH, TOOBIG, TUNER4, TUNER5,
     8        VNEED, VSAVE, W, WSCALE, XIRC, X0
C
C  ***  IV SUBSCRIPT VALUES  ***
C
C  ***  (NOTE THAT P0 AND PC ARE STORED IN IV(G0) AND IV(STLSTG) RESP.)
C
C/6
C     DATA CNVCOD/55/, COVMAT/26/, COVREQ/15/, DIG/37/, FDH/74/, H/56/,
C    1     HC/71/, IERR/75/, INITS/25/, IPIVOT/76/, IRC/29/, IVNEED/3/,
C    2     KAGQT/33/, KALM/34/, LMAT/42/, MODE/35/, MODEL/5/,
C    3     MXFCAL/17/, MXITER/18/, NEXTIV/46/, NEXTV/47/, NFCALL/6/,
C    4     NFGCAL/7/, NFCOV/52/, NGCOV/53/, NGCALL/30/, NITER/31/,
C    5     P0/48/, PC/41/, PERM/58/, QTR/77/, RADINC/8/, RDREQ/57/,
C    6     REGD/67/, RESTOR/9/, RMAT/78/, S/62/, STEP/40/, STGLIM/11/,
C    7     SUSED/64/, SWITCH/12/, TOOBIG/2/, VNEED/4/, VSAVE/60/, W/65/,
C    8     XIRC/13/, X0/43/
C/7
      PARAMETER (CNVCOD=55, COVMAT=26, COVREQ=15, DIG=37, FDH=74, H=56,
     1           HC=71, IERR=75, INITS=25, IPIVOT=76, IRC=29, IVNEED=3,
     2           KAGQT=33, KALM=34, LMAT=42, MODE=35, MODEL=5,
     3           MXFCAL=17, MXITER=18, NEXTIV=46, NEXTV=47, NFCALL=6,
     4           NFGCAL=7, NFCOV=52, NGCOV=53, NGCALL=30, NITER=31,
     5           P0=48, PC=41, PERM=58, QTR=77, RADINC=8, RDREQ=57,
     6           REGD=67, RESTOR=9, RMAT=78, S=62, STEP=40, STGLIM=11,
     7           SUSED=64, SWITCH=12, TOOBIG=2, VNEED=4, VSAVE=60, W=65,
     8           XIRC=13, X0=43)
C/
C
C  ***  V SUBSCRIPT VALUES  ***
C
C/6
C     DATA COSMIN/47/, DGNORM/1/, DSTNRM/2/, F/10/, FDIF/11/, FUZZ/45/,
C    1     F0/13/, GTSTEP/4/, INCFAC/23/, LMAX0/35/, LMAXS/36/,
C    2     NVSAVE/9/, PHMXFC/21/, PREDUC/7/, RADFAC/16/, RADIUS/8/,
C    3     RAD0/9/, RELDX/17/, SIZE/55/, STPPAR/5/, TUNER4/29/,
C    4     TUNER5/30/, WSCALE/56/
C/7
      PARAMETER (COSMIN=47, DGNORM=1, DSTNRM=2, F=10, FDIF=11, FUZZ=45,
     1           F0=13, GTSTEP=4, INCFAC=23, LMAX0=35, LMAXS=36,
     2           NVSAVE=9, PHMXFC=21, PREDUC=7, RADFAC=16, RADIUS=8,
     3           RAD0=9, RELDX=17, SIZE=55, STPPAR=5, TUNER4=29,
     4           TUNER5=30, WSCALE=56)
C/
C
C
C/6
C     DATA HALF/0.5D+0/, NEGONE/-1.D+0/, ONE/1.D+0/, ONEP2/1.2D+0/,
C    1     ZERO/0.D+0/
C/7
      PARAMETER (HALF=0.5D+0, NEGONE=-1.D+0, ONE=1.D+0, ONEP2=1.2D+0,
     1           ZERO=0.D+0)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      I = IV(1)
      IF (I .EQ. 1) GO TO 50
      IF (I .EQ. 2) GO TO 60
C
      IF (I .LT. 12) GO TO 10
      IF (I .GT. 13) GO TO 10
         IV(VNEED) = IV(VNEED) + P*(3*P + 25)/2 + 7
         IV(IVNEED) = IV(IVNEED) + 4*P
 10   CALL DPARCK(1, D, IV, LIV, LV, P, V)
      I = IV(1) - 2
      IF (I .GT. 12) GO TO 999
      GO TO (360, 360, 360, 360, 360, 360, 240, 190, 240, 20, 20, 30), I
C
C  ***  STORAGE ALLOCATION  ***
C
 20   PP1O2 = P * (P + 1) / 2
      IV(S) = IV(LMAT) + PP1O2
      IV(X0) = IV(S) + PP1O2
      IV(STEP) = IV(X0) + 2*P
      IV(DIG) = IV(STEP) + 3*P
      IV(W) = IV(DIG) + 2*P
      IV(H) = IV(W) + 4*P + 7
      IV(NEXTV) = IV(H) + PP1O2
      IV(IPIVOT) = IV(PERM) + 3*P
      IV(NEXTIV) = IV(IPIVOT) + P
      IF (IV(1) .NE. 13) GO TO 30
         IV(1) = 14
         GO TO 999
C
C  ***  INITIALIZATION  ***
C
 30   IV(NITER) = 0
      IV(NFCALL) = 1
      IV(NGCALL) = 1
      IV(NFGCAL) = 1
      IV(MODE) = -1
      IV(STGLIM) = 2
      IV(TOOBIG) = 0
      IV(CNVCOD) = 0
      IV(COVMAT) = 0
      IV(NFCOV) = 0
      IV(NGCOV) = 0
      IV(RADINC) = 0
      IV(PC) = P
      V(RAD0) = ZERO
      V(STPPAR) = ZERO
      V(RADIUS) = V(LMAX0) / (ONE + V(PHMXFC))
C
C  ***  CHECK CONSISTENCY OF B AND INITIALIZE IP ARRAY  ***
C
      IPI = IV(IPIVOT)
      DO 40 I = 1, P
         IV(IPI) = I
         IPI = IPI + 1
         IF (B(1,I) .GT. B(2,I)) GO TO 680
 40      CONTINUE
C
C  ***  SET INITIAL MODEL AND S MATRIX  ***
C
      IV(MODEL) = 1
      IV(1) = 1
      IF (IV(S) .LT. 0) GO TO 710
      IF (IV(INITS) .GT. 1) IV(MODEL) = 2
      S1 = IV(S)
      IF (IV(INITS) .EQ. 0 .OR. IV(INITS) .GT. 2)
     1   CALL DV7SCP(P*(P+1)/2, V(S1), ZERO)
      GO TO 710
C
C  ***  NEW FUNCTION VALUE  ***
C
 50   IF (IV(MODE) .EQ. 0) GO TO 360
      IF (IV(MODE) .GT. 0) GO TO 590
C
      IF (IV(TOOBIG) .EQ. 0) GO TO 690
         IV(1) = 63
         GO TO 999
C
C  ***  MAKE SURE GRADIENT COULD BE COMPUTED  ***
C
 60   IF (IV(TOOBIG) .EQ. 0) GO TO 70
         IV(1) = 65
         GO TO 999
C
C  ***  NEW GRADIENT  ***
C
 70   IV(KALM) = -1
      IV(KAGQT) = -1
      IV(FDH) = 0
      IF (IV(MODE) .GT. 0) GO TO 590
      IF (IV(HC) .LE. 0 .AND. IV(RMAT) .LE. 0) GO TO 670
C
C  ***  CHOOSE INITIAL PERMUTATION  ***
C
      IPI = IV(IPIVOT)
      IPN = IPI + P - 1
      IPIV2 = IV(PERM) - 1
      K = IV(PC)
      P1 = P
      PP1 = P + 1
      RMAT1 = IV(RMAT)
      HAVRM = RMAT1 .GT. 0
      QTR1 = IV(QTR)
      HAVQTR = QTR1 .GT. 0
C     *** MAKE SURE V(QTR1) IS LEGAL (EVEN WHEN NOT REFERENCED) ***
      W1 = IV(W)
      IF (.NOT. HAVQTR) QTR1 = W1 + P
C
      DO 100 I = 1, P
         I1 = IV(IPN)
         IPN = IPN - 1
         IF (B(1,I1) .GE. B(2,I1)) GO TO 80
         XI = X(I1)
         GI = G(I1)
         IF (XI .LE. B(1,I1) .AND. GI .GT. ZERO) GO TO 80
         IF (XI .GE. B(2,I1) .AND. GI .LT. ZERO) GO TO 80
C           *** DISALLOW CONVERGENCE IF X(I1) HAS JUST BEEN FREED ***
            J = IPIV2 + I1
            IF (IV(J) .GT. K) IV(CNVCOD) = 0
            GO TO 100
 80      IF (I1 .GE. P1) GO TO 90
            I1 = PP1 - I
            CALL I7SHFT(P1, I1, IV(IPI))
            IF (HAVRM)
     1          CALL DQ7RSH(I1, P1, HAVQTR, V(QTR1), V(RMAT1), V(W1))
 90      P1 = P1 - 1
 100     CONTINUE
      IV(PC) = P1
C
C  ***  COMPUTE V(DGNORM) (AN OUTPUT VALUE IF WE STOP NOW)  ***
C
      V(DGNORM) = ZERO
      IF (P1 .LE. 0) GO TO 110
      DIG1 = IV(DIG)
      CALL DV7VMP(P, V(DIG1), G, D, -1)
      CALL DV7IPR(P, IV(IPI), V(DIG1))
      V(DGNORM) = DV2NRM(P1, V(DIG1))
 110  IF (IV(CNVCOD) .NE. 0) GO TO 580
      IF (IV(MODE) .EQ. 0) GO TO 510
      IV(MODE) = 0
      V(F0) = V(F)
      IF (IV(INITS) .LE. 2) GO TO 170
C
C  ***  ARRANGE FOR FINITE-DIFFERENCE INITIAL S  ***
C
      IV(XIRC) = IV(COVREQ)
      IV(COVREQ) = -1
      IF (IV(INITS) .GT. 3) IV(COVREQ) = 1
      IV(CNVCOD) = 70
      GO TO 600
C
C  ***  COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S  ***
C
 120  H1 = IV(FDH)
      IF (H1 .LE. 0) GO TO 660
      IV(CNVCOD) = 0
      IV(MODE) = 0
      IV(NFCOV) = 0
      IV(NGCOV) = 0
      IV(COVREQ) = IV(XIRC)
      S1 = IV(S)
      PP1O2 = PS * (PS + 1) / 2
      HC1 = IV(HC)
      IF (HC1 .LE. 0) GO TO 130
         CALL DV2AXY(PP1O2, V(S1), NEGONE, V(HC1), V(H1))
         GO TO 140
 130  RMAT1 = IV(RMAT)
      LMAT1 = IV(LMAT)
      CALL DL7SQR(P, V(LMAT1), V(RMAT1))
      IPI = IV(IPIVOT)
      IPIV1 = IV(PERM) + P
      CALL I7PNVR(P, IV(IPIV1), IV(IPI))
      CALL DS7IPR(P, IV(IPIV1), V(LMAT1))
      CALL DV2AXY(PP1O2, V(S1), NEGONE, V(LMAT1), V(H1))
C
C     *** ZERO PORTION OF S CORRESPONDING TO FIXED X COMPONENTS ***
C
 140  DO 160 I = 1, P
         IF (B(1,I) .LT. B(2,I)) GO TO 160
         K = S1 + I*(I-1)/2
         CALL DV7SCP(I, V(K), ZERO)
         IF (I .GE. P) GO TO 170
         K = K + 2*I - 1
         I1 = I + 1
         DO 150 J = I1, P
            V(K) = ZERO
            K = K + J
 150        CONTINUE
 160     CONTINUE
C
 170  IV(1) = 2
C
C
C-----------------------------  MAIN LOOP  -----------------------------
C
C
C  ***  PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT  ***
C
 180  CALL DITSUM(D, G, IV, LIV, LV, P, V, X)
 190  K = IV(NITER)
      IF (K .LT. IV(MXITER)) GO TO 200
         IV(1) = 10
         GO TO 999
 200  IV(NITER) = K + 1
C
C  ***  UPDATE RADIUS  ***
C
      IF (K .EQ. 0) GO TO 220
      STEP1 = IV(STEP)
      DO 210 I = 1, P
         V(STEP1) = D(I) * V(STEP1)
         STEP1 = STEP1 + 1
 210     CONTINUE
      STEP1 = IV(STEP)
      T = V(RADFAC) * DV2NRM(P, V(STEP1))
      IF (V(RADFAC) .LT. ONE .OR. T .GT. V(RADIUS)) V(RADIUS) = T
C
C  ***  INITIALIZE FOR START OF NEXT ITERATION  ***
C
 220  X01 = IV(X0)
      V(F0) = V(F)
      IV(IRC) = 4
      IV(H) = -IABS(IV(H))
      IV(SUSED) = IV(MODEL)
C
C     ***  COPY X TO X0  ***
C
      CALL DV7CPY(P, V(X01), X)
C
C  ***  CHECK STOPX AND FUNCTION EVALUATION LIMIT  ***
C
 230  IF (.NOT. STOPX(DUMMY)) GO TO 250
         IV(1) = 11
         GO TO 260
C
C     ***  COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX.
C
 240  IF (V(F) .GE. V(F0)) GO TO 250
         V(RADFAC) = ONE
         K = IV(NITER)
         GO TO 200
C
 250  IF (IV(NFCALL) .LT. IV(MXFCAL) + IV(NFCOV)) GO TO 270
         IV(1) = 9
 260     IF (V(F) .GE. V(F0)) GO TO 999
C
C        ***  IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH
C        ***  IMPROVED V(F), EVALUATE THE GRADIENT AT X.
C
              IV(CNVCOD) = IV(1)
              GO TO 500
C
C. . . . . . . . . . . . .  COMPUTE CANDIDATE STEP  . . . . . . . . . .
C
 270  STEP1 = IV(STEP)
      TG1 = IV(DIG)
      TD1 = TG1 + P
      X01 = IV(X0)
      W1 = IV(W)
      H1 = IV(H)
      P1 = IV(PC)
      IPI = IV(PERM)
      IPIV1 = IPI + P
      IPIV2 = IPIV1 + P
      IPIV0 = IV(IPIVOT)
      IF (IV(MODEL) .EQ. 2) GO TO 280
C
C        ***  COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE...
C
         RMAT1 = IV(RMAT)
         IF (RMAT1 .LE. 0) GO TO 280
         QTR1 = IV(QTR)
         IF (QTR1 .LE. 0) GO TO 280
         LMAT1 = IV(LMAT)
         WLM1 = W1 + P
         CALL DL7MSB(B, D, G, IV(IERR), IV(IPIV0), IV(IPIV1),
     1               IV(IPIV2), IV(KALM), V(LMAT1), LV, P, IV(P0),
     2               IV(PC), V(QTR1), V(RMAT1), V(STEP1), V(TD1),
     3               V(TG1), V, V(W1), V(WLM1), X, V(X01))
C        *** H IS STORED IN THE END OF W AND HAS JUST BEEN OVERWRITTEN,
C        *** SO WE MARK IT INVALID...
         IV(H) = -IABS(H1)
C        *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO
C        *** MARK INVALID THE INFORMATION DG7QTS MAY HAVE STORED IN V...
         IV(KAGQT) = -1
         GO TO 330
C
 280  IF (H1 .GT. 0) GO TO 320
C
C     ***  SET H TO  D**-1 * (HC + T1*S) * D**-1.  ***
C
         P1LEN = P1*(P1+1)/2
         H1 = -H1
         IV(H) = H1
         IV(FDH) = 0
         IF (P1 .LE. 0) GO TO 320
C        *** MAKE TEMPORARY PERMUTATION ARRAY ***
         CALL I7COPY(P, IV(IPI), IV(IPIV0))
         J = IV(HC)
         IF (J .GT. 0) GO TO 290
            J = H1
            RMAT1 = IV(RMAT)
            CALL DL7SQR(P1, V(H1), V(RMAT1))
            GO TO 300
 290     CALL DV7CPY(P*(P+1)/2, V(H1), V(J))
         CALL DS7IPR(P, IV(IPI), V(H1))
 300     IF (IV(MODEL) .EQ. 1) GO TO 310
            LMAT1 = IV(LMAT)
            S1 = IV(S)
            CALL DV7CPY(P*(P+1)/2, V(LMAT1), V(S1))
            CALL DS7IPR(P, IV(IPI), V(LMAT1))
            CALL DV2AXY(P1LEN, V(H1), ONE, V(LMAT1), V(H1))
 310     CALL DV7CPY(P, V(TD1), D)
         CALL DV7IPR(P, IV(IPI), V(TD1))
         CALL DS7DMP(P1, V(H1), V(H1), V(TD1), -1)
         IV(KAGQT) = -1
C
C  ***  COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP  ***
C
 320  LMAT1 = IV(LMAT)
      CALL DG7QSB(B, D, V(H1), G, IV(IPI), IV(IPIV1), IV(IPIV2),
     1            IV(KAGQT), V(LMAT1), LV, P, IV(P0), P1, V(STEP1),
     2            V(TD1), V(TG1), V, V(W1), X, V(X01))
      IF (IV(KALM) .GT. 0) IV(KALM) = 0
C
 330  IF (IV(IRC) .NE. 6) GO TO 340
         IF (IV(RESTOR) .NE. 2) GO TO 360
         RSTRST = 2
         GO TO 370
C
C  ***  CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE  ***
C
 340  IV(TOOBIG) = 0
      IF (V(DSTNRM) .LE. ZERO) GO TO 360
      IF (IV(IRC) .NE. 5) GO TO 350
      IF (V(RADFAC) .LE. ONE) GO TO 350
      IF (V(PREDUC) .GT. ONEP2 * V(FDIF)) GO TO 350
         STEP1 = IV(STEP)
         X01 = IV(X0)
         CALL DV2AXY(P, V(STEP1), NEGONE, V(X01), X)
         IF (IV(RESTOR) .NE. 2) GO TO 360
         RSTRST = 0
         GO TO 370
C
C  ***  COMPUTE F(X0 + STEP)  ***
C
 350  X01 = IV(X0)
      STEP1 = IV(STEP)
      CALL DV2AXY(P, X, ONE, V(STEP1), V(X01))
      IV(NFCALL) = IV(NFCALL) + 1
      IV(1) = 1
      GO TO 710
C
C. . . . . . . . . . . . .  ASSESS CANDIDATE STEP  . . . . . . . . . . .
C
 360  RSTRST = 3
 370  X01 = IV(X0)
      V(RELDX) = DRLDST(P, D, X, V(X01))
      CALL DA7SST(IV, LIV, LV, V)
      STEP1 = IV(STEP)
      LSTGST = X01 + P
      I = IV(RESTOR) + 1
      GO TO (410, 380, 390, 400), I
 380  CALL DV7CPY(P, X, V(X01))
      GO TO 410
 390   CALL DV7CPY(P, V(LSTGST), V(STEP1))
       GO TO 410
 400     CALL DV7CPY(P, V(STEP1), V(LSTGST))
         CALL DV2AXY(P, X, ONE, V(STEP1), V(X01))
         V(RELDX) = DRLDST(P, D, X, V(X01))
         IV(RESTOR) = RSTRST
C
C  ***  IF NECESSARY, SWITCH MODELS  ***
C
 410  IF (IV(SWITCH) .EQ. 0) GO TO 420
         IV(H) = -IABS(IV(H))
         IV(SUSED) = IV(SUSED) + 2
         L = IV(VSAVE)
         CALL DV7CPY(NVSAVE, V, V(L))
 420  L = IV(IRC) - 4
      STPMOD = IV(MODEL)
      IF (L .GT. 0) GO TO (440,450,460,460,460,460,460,460,570,510), L
C
C  ***  DECIDE WHETHER TO CHANGE MODELS  ***
C
      E = V(PREDUC) - V(FDIF)
      S1 = IV(S)
      CALL DS7LVM(PS, Y, V(S1), V(STEP1))
      STTSST = HALF * DD7TPR(PS, V(STEP1), Y)
      IF (IV(MODEL) .EQ. 1) STTSST = -STTSST
      IF (DABS(E + STTSST) * V(FUZZ) .GE. DABS(E)) GO TO 430
C
C     ***  SWITCH MODELS  ***
C
         IV(MODEL) = 3 - IV(MODEL)
         IF (-2 .LT. L) GO TO 470
              IV(H) = -IABS(IV(H))
              IV(SUSED) = IV(SUSED) + 2
              L = IV(VSAVE)
              CALL DV7CPY(NVSAVE, V(L), V)
              GO TO 230
C
 430  IF (-3 .LT. L) GO TO 470
C
C     ***  RECOMPUTE STEP WITH DIFFERENT RADIUS  ***
C
 440  V(RADIUS) = V(RADFAC) * V(DSTNRM)
      GO TO 230
C
C  ***  COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST
C
 450  V(RADIUS) = V(LMAXS)
      GO TO 270
C
C  ***  CONVERGENCE OR FALSE CONVERGENCE  ***
C
 460  IV(CNVCOD) = L
      IF (V(F) .GE. V(F0)) GO TO 580
         IF (IV(XIRC) .EQ. 14) GO TO 580
              IV(XIRC) = 14
C
C. . . . . . . . . . . .  PROCESS ACCEPTABLE STEP  . . . . . . . . . . .
C
 470  IV(COVMAT) = 0
      IV(REGD) = 0
C
C  ***  SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS  ***
C
      IF (IV(IRC) .NE. 3) GO TO 500
         STEP1 = IV(STEP)
         TEMP1 = STEP1 + P
         TEMP2 = IV(X0)
C
C     ***  SET  TEMP1 = HESSIAN * STEP  FOR USE IN GRADIENT TESTS  ***
C
         HC1 = IV(HC)
         IF (HC1 .LE. 0) GO TO 480
              CALL DS7LVM(P, V(TEMP1), V(HC1), V(STEP1))
              GO TO 490
 480     RMAT1 = IV(RMAT)
         IPIV0 = IV(IPIVOT)
         CALL DV7CPY(P, V(TEMP1), V(STEP1))
         CALL DV7IPR(P, IV(IPIV0), V(TEMP1))
         CALL DL7TVM(P, V(TEMP1), V(RMAT1), V(TEMP1))
         CALL DL7VML(P, V(TEMP1), V(RMAT1), V(TEMP1))
         IPIV1 = IV(PERM) + P
         CALL I7PNVR(P, IV(IPIV1), IV(IPIV0))
         CALL DV7IPR(P, IV(IPIV1), V(TEMP1))
C
 490     IF (STPMOD .EQ. 1) GO TO 500
              S1 = IV(S)
              CALL DS7LVM(PS, V(TEMP2), V(S1), V(STEP1))
              CALL DV2AXY(PS, V(TEMP1), ONE, V(TEMP2), V(TEMP1))
C
C  ***  SAVE OLD GRADIENT AND COMPUTE NEW ONE  ***
C
 500  IV(NGCALL) = IV(NGCALL) + 1
      G01 = IV(W)
      CALL DV7CPY(P, V(G01), G)
      GO TO 690
C
C  ***  INITIALIZATIONS -- G0 = G - G0, ETC.  ***
C
 510  G01 = IV(W)
      CALL DV2AXY(P, V(G01), NEGONE, V(G01), G)
      STEP1 = IV(STEP)
      TEMP1 = STEP1 + P
      TEMP2 = IV(X0)
      IF (IV(IRC) .NE. 3) GO TO 540
C
C  ***  SET V(RADFAC) BY GRADIENT TESTS  ***
C
C     ***  SET  TEMP1 = D**-1 * (HESSIAN * STEP  +  (G(X0) - G(X)))  ***
C
         K = TEMP1
         L = G01
         DO 520 I = 1, P
              V(K) = (V(K) - V(L)) / D(I)
              K = K + 1
              L = L + 1
 520          CONTINUE
C
C        ***  DO GRADIENT TESTS  ***
C
         IF (DV2NRM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4))  GO TO 530
              IF (DD7TPR(P, G, V(STEP1))
     1                  .GE. V(GTSTEP) * V(TUNER5))  GO TO 540
 530               V(RADFAC) = V(INCFAC)
C
C  ***  COMPUTE Y VECTOR NEEDED FOR UPDATING S  ***
C
 540  CALL DV2AXY(PS, Y, NEGONE, Y, G)
C
C  ***  DETERMINE SIZING FACTOR V(SIZE)  ***
C
C     ***  SET TEMP1 = S * STEP  ***
      S1 = IV(S)
      CALL DS7LVM(PS, V(TEMP1), V(S1), V(STEP1))
C
      T1 = DABS(DD7TPR(PS, V(STEP1), V(TEMP1)))
      T = DABS(DD7TPR(PS, V(STEP1), Y))
      V(SIZE) = ONE
      IF (T .LT. T1) V(SIZE) = T / T1
C
C  ***  SET G0 TO WCHMTD CHOICE OF FLETCHER AND AL-BAALI  ***
C
      HC1 = IV(HC)
      IF (HC1 .LE. 0) GO TO 550
         CALL DS7LVM(PS, V(G01), V(HC1), V(STEP1))
         GO TO 560
C
 550  RMAT1 = IV(RMAT)
      IPIV0 = IV(IPIVOT)
      CALL DV7CPY(P, V(G01), V(STEP1))
      I = G01 + PS
      IF (PS .LT. P) CALL DV7SCP(P-PS, V(I), ZERO)
      CALL DV7IPR(P, IV(IPIV0), V(G01))
      CALL DL7TVM(P, V(G01), V(RMAT1), V(G01))
      CALL DL7VML(P, V(G01), V(RMAT1), V(G01))
      IPIV1 = IV(PERM) + P
      CALL I7PNVR(P, IV(IPIV1), IV(IPIV0))
      CALL DV7IPR(P, IV(IPIV1), V(G01))
C
 560  CALL DV2AXY(PS, V(G01), ONE, Y, V(G01))
C
C  ***  UPDATE S  ***
C
      CALL DS7LUP(V(S1), V(COSMIN), PS, V(SIZE), V(STEP1), V(TEMP1),
     1            V(TEMP2), V(G01), V(WSCALE), Y)
      IV(1) = 2
      GO TO 180
C
C. . . . . . . . . . . . . .  MISC. DETAILS  . . . . . . . . . . . . . .
C
C  ***  BAD PARAMETERS TO ASSESS  ***
C
 570  IV(1) = 64
      GO TO 999
C
C
C  ***  CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE  ***
C
 580  IF (IV(RDREQ) .EQ. 0) GO TO 660
      IF (IV(FDH) .NE. 0) GO TO 660
      IF (IV(CNVCOD) .GE. 7) GO TO 660
      IF (IV(REGD) .GT. 0) GO TO 660
      IF (IV(COVMAT) .GT. 0) GO TO 660
      IF (IABS(IV(COVREQ)) .GE. 3) GO TO 640
      IF (IV(RESTOR) .EQ. 0) IV(RESTOR) = 2
      GO TO 600
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE  ***
C
 590  IV(RESTOR) = 0
 600  CALL DF7DHB(B, D, G, I, IV, LIV, LV, P, V, X)
      GO TO (610, 620, 630), I
 610  IV(NFCOV) = IV(NFCOV) + 1
      IV(NFCALL) = IV(NFCALL) + 1
      IV(1) = 1
      GO TO 710
C
 620  IV(NGCOV) = IV(NGCOV) + 1
      IV(NGCALL) = IV(NGCALL) + 1
      IV(NFGCAL) = IV(NFCALL) + IV(NGCOV)
      GO TO 690
C
 630  IF (IV(CNVCOD) .EQ. 70) GO TO 120
      GO TO 660
C
 640  H1 = IABS(IV(H))
      IV(FDH) = H1
      IV(H) = -H1
      HC1 = IV(HC)
      IF (HC1 .LE. 0) GO TO 650
           CALL DV7CPY(P*(P+1)/2, V(H1), V(HC1))
           GO TO 660
 650  RMAT1 = IV(RMAT)
      CALL DL7SQR(P, V(H1), V(RMAT1))
C
 660  IV(MODE) = 0
      IV(1) = IV(CNVCOD)
      IV(CNVCOD) = 0
      GO TO 999
C
C  ***  SPECIAL RETURN FOR MISSING HESSIAN INFORMATION -- BOTH
C  ***  IV(HC) .LE. 0 AND IV(RMAT) .LE. 0
C
 670  IV(1) = 1400
      GO TO 999
C
C  ***  INCONSISTENT B  ***
C
 680  IV(1) = 82
      GO TO 999
C
C  *** SAVE, THEN INITIALIZE IPIVOT ARRAY BEFORE COMPUTING G ***
C
 690  IV(1) = 2
      J = IV(IPIVOT)
      IPI = IV(PERM)
      CALL I7PNVR(P, IV(IPI), IV(J))
      DO 700 I = 1, P
         IV(J) = I
         J = J + 1
 700     CONTINUE
C
C  ***  PROJECT X INTO FEASIBLE REGION (PRIOR TO COMPUTING F OR G)  ***
C
 710  DO 720 I = 1, P
         IF (X(I) .LT. B(1,I)) X(I) = B(1,I)
         IF (X(I) .GT. B(2,I)) X(I) = B(2,I)
 720     CONTINUE
      IV(TOOBIG) = 0
C
 999  RETURN
C
C  ***  LAST LINE OF DG7ITB FOLLOWS  ***
      END

.