SUBROUTINE  GLGB(N, P, PS, X, B, RHO, RHOI, RHOR, IV, LIV, LV,
     1                V, CALCRJ, UI, UR, UF)
C
C *** GENERALIZED LINEAR REGRESSION A LA NL2SOL, PLUS SIMPLE BOUNDS ***
C
C  ***  PARAMETERS  ***
C
      INTEGER N, P, PS, LIV, LV
      INTEGER IV(LIV), RHOI(*), UI(*)
      REAL B(2,P), X(P), RHOR(*), V(LV), UR(*)
      EXTERNAL CALCRJ, RHO, UF
C
C  ***  PARAMETER USAGE  ***
C
C N....... TOTAL NUMBER OF RESIDUALS.
C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED.
C PS...... NUMBER OF NON-NUISANCE PARAMETERS (THOSE INVOLVED IN S).
C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS,
C             OUTPUT = BEST VALUE FOUND).
C B....... BOUNDS TO ENFORCE... B(1,I) .LE. X(I) .LE. B(2,I).
C RHO..... SUBROUTINE FOR COMPUTING LOSS FUNCTIONS AND THEIR DERIVS.
C             SEE   RGLG FOR DETAILS ABOUT RHO.
C RHOI.... PASSED WITHOUT CHANGE TO RHO.
C RHOR.... PASSED WITHOUT CHANGE TO RHO.
C IV...... INTEGER VALUES ARRAY.
C LIV..... LENGTH OF IV, AT LEAST 82 + P.
C LV...... LENGTH OF V, AT LEAST
C              105 + P*(2*P + 20) + 4*PS
C            + N*(P + 3 + (P-PS+1)*(P-PS+2)/2).
C V....... FLOATING-POINT VALUES ARRAY.
C CALCRJ.. SUBROUTINE FOR COMPUTING RESIDUAL VECTOR AND JACOBIAN MATRIX.
C UI...... PASSED UNCHANGED TO CALCRJ.
C UR...... PASSED UNCHANGED TO CALCRJ.
C UF...... PASSED UNCHANGED TO CALCRJ.
C
C *** CALCRJ CALLING SEQUENCE...
C
C      CALL CALCRJ(N, PS, X, NF, NEED, R, RP, UI, UR, UF)
C
C PARAMETERS N, PS, X, UI, UR, AND UF ARE AS ABOVE.
C R AND RP ARE FLOATING-POINT ARRAYS DIMENSIONED R(N) AND RP(PS,N).
C NEED IS AN INTEGER ARRAY OF LENGTH 2...
C   NEED(1) = 1 MEANS CALCRJ SHOULD COMPUTE THE RESIDUAL VECTOR R,
C             AND NEED(2) IS THE VALUE NF HAD AT THE LAST X WHERE
C             CALCRJ MIGHT BE CALLED WITH NEED(1) = 2.
C   NEED(1) = 2 MEANS CALCRJ SHOULD COMPUTE THE JACOBIAN MATRIX RP,
C             WHERE RP(J,I) = DERIVATIVE OF R(I) WITH RESPECT TO X(J).
C (CALCRJ SHOULD NOT CHANGE NEED AND SHOULD CHANGE AT MOST ONE OF R
C AND RP.  IF R OR RP, AS APPROPRIATE, CANNOT BE COMPUTED, THEN CALCRJ
C SHOULD SET NF TO 0.  OTHERWISE IT SHOULD NOT CHANGE NF.)
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
C  ***  EXTERNAL SUBROUTINES  ***
C
      EXTERNAL  IVSET,   RGLGB
C
C  IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C  RGLGB... CARRIES OUT OPTIMIZATION ITERATIONS.
C
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER D1, DR1, I, IV1, NEED1(2), NEED2(2), NF, R1, RD1
C
C  ***  IV COMPONENTS  ***
C
      INTEGER D, J, NEXTV, NFCALL, NFGCAL, R, REGD, REGD0, TOOBIG, VNEED
      PARAMETER (D=27, J=70, NEXTV=47, NFCALL=6, NFGCAL=7, R=61,
     1           REGD=67, REGD0=82, TOOBIG=2, VNEED=4)
      SAVE NEED1, NEED2
      DATA NEED1(1)/1/, NEED1(2)/0/, NEED2(1)/2/, NEED2(2)/0/
C
C---------------------------------  BODY  ------------------------------
C
      IF (IV(1) .EQ. 0) CALL  IVSET(1, IV, LIV, LV, V)
      IV1 = IV(1)
      IF (IV1 .EQ. 14) GO TO 10
      IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10
      IF (IV1 .EQ. 12) IV(1) = 13
      I = (P-PS+2)*(P-PS+1)/2
      IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+1+I)
      CALL  RGLGB(B, X, V, IV, LIV, LV, N, PS, N, P, PS, V, V,
     1            RHO, RHOI,RHOR, V, X)
      IF (IV(1) .NE. 14) GO TO 999
C
C  ***  STORAGE ALLOCATION  ***
C
      IV(D) = IV(NEXTV)
      IV(R) = IV(D) + P
      IV(REGD0) = IV(R) + (P - PS + 1)*N
      IV(J) = IV(REGD0) + ((P-PS+2)*(P-PS+1)/2)*N
      IV(NEXTV) = IV(J) + N*PS
      IF (IV1 .EQ. 13) GO TO 999
C
 10   D1 = IV(D)
      DR1 = IV(J)
      R1 = IV(R)
      RD1 = IV(REGD0)
C
 20   CALL  RGLGB(B, V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS,
     1            V(R1), V(RD1), RHO, RHOI, RHOR, V, X)
      IF (IV(1)-2) 30, 50, 60
C
C  ***  NEW FUNCTION VALUE (R VALUE) NEEDED  ***
C
 30   NF = IV(NFCALL)
      NEED1(2) = IV(NFGCAL)
      CALL CALCRJ(N, PS, X, NF, NEED1, V(R1), V(DR1), UI, UR, UF)
      IF (NF .GT. 0) GO TO 40
         IV(TOOBIG) = 1
         GO TO 20
 40   IF (IV(1) .GT. 0) GO TO 20
C
C  ***  COMPUTE DR = GRADIENT OF R COMPONENTS  ***
C
 50   CALL CALCRJ(N, PS, X, IV(NFGCAL), NEED2, V(R1), V(DR1), UI, UR,UF)
      IF (IV(NFGCAL) .EQ. 0) IV(TOOBIG) = 1
      GO TO 20
C
C  ***  INDICATE WHETHER THE REGRESSION DIAGNOSTIC ARRAY WAS COMPUTED
C  ***  AND PRINT IT IF SO REQUESTED...
C
 60   IF (IV(REGD) .GT. 0) IV(REGD) = RD1
C
 999  RETURN
C
C  ***  LAST LINE OF  GLGB FOLLOWS  ***
      END
      SUBROUTINE  GLFB(N, P, PS, X, B, RHO, RHOI, RHOR, IV, LIV, LV, V,
     1                  CALCRJ, UI, UR, UF)
C
C *** GENERALIZED LINEAR REGRESSION, FINITE-DIFFERENCE JACOBIAN ***
C *** WITH SIMPLE BOUNDS ON X ***
C
C  ***  PARAMETERS  ***
C
      INTEGER N, P, PS, LIV, LV
      INTEGER IV(LIV), RHOI(*), UI(*)
      REAL B(2,P), X(P), V(LV), RHOR(*), UR(*)
      EXTERNAL CALCRJ, RHO, UF
C
C  ***  PARAMETER USAGE  ***
C
C N....... TOTAL NUMBER OF RESIDUALS.
C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED.
C PS...... NUMBER OF NON-NUISANCE PARAMETERS (THOSE INVOLVED IN S).
C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS,
C             OUTPUT = BEST VALUE FOUND).
C B....... BOUNDS TO ENFORCE... B(1,I) .LE. X(I) .LE. B(2,I).
C RHO..... SUBROUTINE FOR COMPUTING LOSS FUNCTIONS AND THEIR DERIVS.
C             SEE   RGLG FOR DETAILS ABOUT RHO.
C RHOI.... PASSED WITHOUT CHANGE TO RHO.
C RHOR.... PASSED WITHOUT CHANGE TO RHO.
C IV...... INTEGER VALUES ARRAY.
C LIV..... LENGTH OF IV, AT LEAST 82 + 5*P.
C LV...... LENGTH OF V, AT LEAST
C              105 + P*(2*P + 20) + 4*PS
C            + N*(P + 5 + (P-PS+1)*(P-PS+2)/2).
C V....... FLOATING-POINT VALUES ARRAY.
C CALCRJ.. SUBROUTINE FOR COMPUTING RESIDUAL VECTOR.
C UI...... PASSED UNCHANGED TO CALCRJ.
C UR...... PASSED UNCHANGED TO CALCRJ.
C UF...... PASSED UNCHANGED TO CALCRJ.
C
C *** CALCRJ CALLING SEQUENCE...
C
C      CALL CALCRJ(N, PS, X, NF, NEED, R, RP, UI, UR, UF)
C
C PARAMETERS N, PS, X, UI, UR, AND UF ARE AS ABOVE.
C R AND RP ARE FLOATING-POINT ARRAYS DIMENSIONED R(N) AND RP(PS,N).
C NEED MAY BE REGARDED AS AN INTEGER THAT ALWAYS HAS THE VALUE 1
C WHEN  GLFB CALLS CALCRJ.  THIS MEANS CALCRJ SHOULD COMPUTE THE
C RESIDUAL VECTOR R.  (CALCRJ SHOULD NOT CHANGE NEED OR RP.  IF R
C CANNOT BE COMPUTED, THEN CALCRJ SHOULD SET NF TO 0.  OTHERWISE IT
C SHOULD NOT CHANGE NF.  FOR COMPATIBILITY WITH    GLG, NEED IS A
C VECTOR OF LENGTH 2.)
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
C  ***  EXTERNAL SUBROUTINES  ***
C
      EXTERNAL  IVSET,   RGLGB, V7CPY
C
C  IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C  RGLGB... CARRIES OUT OPTIMIZATION ITERATIONS.
C  V7CPY.... COPIES ONE VECTOR TO ANOTHER.
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER D1, DK, DR1, I, I1, IV1, J1K, J1K0, K, NEED(2), NF,
     1        NG, RD1, R1, R21, RS1, RSN
      REAL H, H0, HLIM, NEGPT5, T, ONE, XK, XK1, ZERO
C
C  ***  IV AND V COMPONENTS  ***
C
      INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL,
     1        NGCALL, NGCOV, R, REGD0, TOOBIG, VNEED
      PARAMETER (COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70, MODE=35,
     1           NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53,
     2           R=61, REGD0=82, TOOBIG=2, VNEED=4)
      SAVE NEED
      DATA HLIM/0.1E+0/, NEGPT5/-0.5E+0/, ONE/1.E+0/, ZERO/0.E+0/
      DATA NEED(1)/1/, NEED(2)/0/
C
C---------------------------------  BODY  ------------------------------
C
      IF (IV(1) .EQ. 0) CALL  IVSET(1, IV, LIV, LV, V)
      IV(COVREQ) = -IABS(IV(COVREQ))
      IV1 = IV(1)
      IF (IV1 .EQ. 14) GO TO 10
      IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10
      IF (IV1 .EQ. 12) IV(1) = 13
      I = (P-PS+2)*(P-PS+1)/2
      IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+3+I)
      CALL  RGLGB(B, X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, RHO,
     1             RHOI, RHOR, V, X)
      IF (IV(1) .NE. 14) GO TO 999
C
C  ***  STORAGE ALLOCATION  ***
C
      IV(D) = IV(NEXTV)
      IV(R) = IV(D) + P
      IV(REGD0) = IV(R) + (P - PS + 3)*N
      IV(J) = IV(REGD0) + ((P-PS+2)*(P-PS+1)/2)*N
      IV(NEXTV) = IV(J) + N*PS
      IF (IV1 .EQ. 13) GO TO 999
C
 10   D1 = IV(D)
      DR1 = IV(J)
      R1 = IV(R)
      RD1 = IV(REGD0)
      R21 = RD1 - N
      RS1 = R21 - N
      RSN = RS1 + N - 1
C
 20   CALL  RGLGB(B, V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS,
     1            V(R1), V(RD1), RHO, RHOI, RHOR, V, X)
      IF (IV(1)-2) 30, 50, 999
C
C  ***  NEW FUNCTION VALUE (R VALUE) NEEDED  ***
C
 30   NF = IV(NFCALL)
      CALL CALCRJ(N, PS, X, NF, NEED, V(R1), V(DR1), UI, UR, UF)
      IF (NF .GT. 0) GO TO 40
         IV(TOOBIG) = 1
         GO TO 20
 40   CALL  V7CPY(N, V(RS1), V(R1))
      IF (IV(1) .GT. 0) GO TO 20
C
C  ***  COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R  ***
C
C     *** INITIALIZE D IF NECESSARY ***
C
 50   IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO)
     1        CALL  V7SCP(P, V(D1), ONE)
C
      DK = D1
      NG = IV(NGCALL) - 1
      IF (IV(1) .EQ. (-1)) IV(NGCOV) = IV(NGCOV) - 1
      J1K0 = DR1
      NF = IV(NFCALL)
      IF (NF .EQ. IV(NFGCAL)) GO TO 70
         NG = NG + 1
         CALL CALCRJ(N, PS, X, NF, NEED, V(RS1), V(DR1), UI, UR, UF)
         IF (NF .GT. 0) GO TO 70
 60          IV(TOOBIG) = 1
             IV(NGCALL) = NG
             GO TO 20
 70   DO 130 K = 1, PS
         J1K = J1K0
         J1K0 = J1K0 + 1
         IF (B(1,K) .GE. B(2,K)) GO TO 120
         XK = X(K)
         H = V(DLTFDJ) *   MAX( ABS(XK), ONE/V(DK))
         H0 = H
         DK = DK + 1
         T = NEGPT5
         XK1 = XK + H
         IF (XK - H .GE. B(1,K)) GO TO 80
            T = -T
            IF (XK1 .GT. B(2,K)) GO TO 60
 80      IF (XK1 .LE. B(2,K)) GO TO 90
            T = -T
            H = -H
            XK1 = XK + H
            IF (XK1 .LT. B(1,K)) GO TO 60
 90      X(K) = XK1
         NF = IV(NFGCAL)
         CALL CALCRJ(N, PS, X, NF, NEED, V(R21), V(DR1), UI, UR, UF)
         NG = NG + 1
         IF (NF .GT. 0) GO TO 100
              H = T * H
              XK1 = XK + H
              IF ( ABS(H/H0) .GE. HLIM) GO TO 90
                   GO TO 60
 100     X(K) = XK
         IV(NGCALL) = NG
         I1 = R21
         DO 110 I = RS1, RSN
              V(J1K) = (V(I1) - V(I)) / H
              I1 = I1 + 1
              J1K = J1K + PS
 110          CONTINUE
         GO TO 130
C        *** SUPPLY A ZERO DERIVATIVE FOR CONSTANT COMPONENTS...
 120     DO 125 I = 1, N
              V(J1K) = ZERO
              J1K = J1K + PS
 125          CONTINUE
 130     CONTINUE
      GO TO 20
C
 999  RETURN
C
C  ***  LAST LINE OF  GLFB FOLLOWS  ***
      END
      SUBROUTINE  RGLGB(B, D, DR, IV, LIV, LV, N, ND, NN, P, PS, R,
     1                  RD, RHO, RHOI, RHOR, V, X)
C
C *** ITERATION DRIVER FOR GENERALIZED (NON)LINEAR MODELS (ETC.)
C
      INTEGER LIV, LV, N, ND, NN, P, PS
      INTEGER IV(LIV), RHOI(*)
      REAL B(2,P), D(P), DR(ND,N), R(*), RD(*), RHOR(*),
     1                 V(LV), X(*)
C     DIMENSION RD(N, (P-PS)*(P-PS+1)/2 + 1)
      EXTERNAL RHO
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C B........ BOUNDS ON X.
C D........ SCALE VECTOR.
C DR....... DERIVATIVES OF R AT X.
C IV....... INTEGER VALUES ARRAY.
C LIV...... LENGTH OF IV... LIV MUST BE AT LEAST P + 82.
C LV....... LENGTH OF V...  LV  MUST BE AT LEAST 105 + P*(2*P+19).
C N........ TOTAL NUMBER OF RESIDUALS.
C ND....... LEADING DIMENSION OF DR -- MUST BE AT LEAST PS.
C NN....... LEAD DIMENSION OF R, RD.
C P........ NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED.
C PS....... NUMBER OF NON-NUISANCE PARAMETERS.
C R........ RESIDUALS (OR MEANS -- FUNCTIONS OF X) WHEN  RGLGB IS CALLED
C        WITH IV(1) = 1.
C RD....... TEMPORARY STORAGE.
C RHO...... COMPUTES INFO ABOUT OBJECTIVE FUNCTION.
C RHOI..... PASSED WITHOUT CHANGE TO RHO.
C RHOR..... PASSED WITHOUT CHANGE TO RHO.
C V........ FLOATING-POINT VALUES ARRAY.
C X........ PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS,
C             OUTPUT = BEST VALUE FOUND).
C
C *** CALLING SEQUENCE FOR RHO...
C
C  CALL RHO(NEED, F, N, NF, XN, R, RD, RHOI, RHOR, W)
C
C  PARAMETER DECLARATIONS FOR RHO...
C
C INTEGER NEED(2), N, NF, RHOI(*)
C FLOATING-POINT F, XN(*), R(*), RD(N,*), RHOR(*), W(N)
C
C    RHOI AND RHOR ARE FOR RHO TO USE AS IT SEES FIT.  THEY ARE PASSED
C TO RHO WITHOUT CHANGE.
C    F, R, RD, AND W ARE EXPLAINED BELOW WITH NEED.
C    XN IS THE VECTOR OF NUISANCE PARAMETERS (OF LENGTH P - PS).  IF
C RHO NEEDS TO KNOW THE LENGTH OF XN, THEN THIS LENGTH SHOULD BE
C COMMUNICATED THROUGH RHOI (OR THROUGH COMMON).  RHO SHOULD NOT CHANGE
C XN.
C    NEED(1) = 1 MEANS RHO SHOULD SET F TO THE SUM OF THE LOSS FUNCTION
C VALUES AT THE RESIDUALS R(I).  NF IS THE CURRENT FUNCTION INVOCATION
C COUNT (A VALUE THAT IS INCREMENTED EACH TIME A NEW PARAMETER ESTIMATE
C X IS CONSIDERED).  NEED(2) IS THE VALUE NF HAD AT THE LAST R WHERE
C RHO MIGHT BE CALLED WITH NEED(1) = 2.  IF RHO SAVES INTERMEDIATE
C RESULTS FOR USE IN CALLS WITH NEED(1) = 2, THEN IT CAN USE NF TO TELL
C WHICH INTERMEDIATE RESULTS ARE APPROPRIATE, AND IT CAN SAVE SOME OF
C THESE RESULTS IN R.
C    NEED(1) = 2 MEANS RHO SHOULD SET R(I) TO THE LOSS FUNCTION
C DERIVATIVE WITH RESPECT TO THE RESIDUALS THAT WERE PASSED TO RHO WHEN
C NF HAD THE SAME VALUE IT DOES NOW (AND NEED(1) WAS 1).  RHO SHOULD
C ALSO SET W(I) TO THE APPROXIMATION OF THE SECOND DERIVATIVE OF THE
C LOSS FUNCTION (WITH RESPECT TO THE I-TH RESIDUAL) THAT SHOULD BE USED
C IN THE GAUSS-NEWTON MODEL.  WHEN THERE ARE NUISANCE PARAMETERS (I.E.,
C WHEN PS .LT. P) RHO SHOULD ALSO SET R(I+K*N) TO THE DERIVATIVE OF THE
C LOSS FUNCTION WITH RESPECT TO THE I-TH RESIDUAL AND XN(K), AND IT
C SHOULD SET RD(I,J + K*(K+1)/2 + 1) TO THE SECOND PARTIAL DERIVATIVE
C OF THE I-TH RESIDUAL WITH RESPECT TO XN(J) AND XN(K), 0 .LE. J .LE. K
C AND 1 .LE. K .LE. P - PS, WHERE XN(0) MEANS THE I-TH RESIDUAL ITSELF.
C IN ANY EVENT, RHO SHOULD ALSO SET RD(I,1) TO THE (TRUE) SECOND
C DERIVATIVE OF THE LOSS FUNCTION WITH RESPECT TO THE I-TH RESIDUAL.
C    NF (THE FUNCTION INVOCATION COUNT WHOSE NORMAL USE IS EXPLAINED
C ABOVE) SHOULD NOT BE CHANGED UNLESS RHO CANNOT CARRY OUT THE REQUESTED
C TASK, IN WHICH CASE RHO SHOULD SET NF TO 0.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C
C+++++++++++++++++++++++++++++  DECLARATIONS  ++++++++++++++++++++++++++
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL  IVSET,  D7TPR,  D7UP5,  G7ITB, ITSUM,  L7ITV,  L7IVM,
     1         L7SRT,  L7SQR,  L7SVX,  L7SVN, L7VML, O7PRD,
     2          Q7ADR, V2AXY, V7CPY,  V7SCL,  V7SCP,  VSUM
      REAL  D7TPR,  L7SVX,  L7SVN,  VSUM
C
C  IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C  D7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS.
C  D7UP5... UPDATES SCALE VECTOR D.
C  G7ITB... PERFORMS BASIC MINIMIZATION ALGORITHM.
C  ITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X.
C  L7ITV... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR.
C  L7IVM... APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX.
C  L7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX.
C  L7SQR... COMPUTES L*(L**T) FOR LOWER TRIANG. MATRIX L.
C  L7SVX... UNDERESTIMATES LARGEST SINGULAR VALUE OF TRIANG. MATRIX.
C  L7SVN... OVERESTIMATES SMALLEST SINGULAR VALUE OF TRIANG. MATRIX.
C  L7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C  O7PRD.... ADDS OUTER PRODUCT OF VECTORS TO A MATRIX.
C  Q7ADR... ADDS ROWS TO QR FACTORIZATION.
C  V2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER.
C  V7CPY.... COPIES ONE VECTOR TO ANOTHER.
C  V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C  V7SCL... MULTIPLIES A VECTOR BY A SCALAR.
C  VSUM.... RETURNS SUM OF ELEMENTS OF A VECTOR.
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL UPDATD, ZEROG
      INTEGER G1, HN1, I, II, IV1, J, J1, JTOL1, K, LH,
     1        NEED1(2), NEED2(2),  PMPS, PS1, PSLEN, QTR1,
     2        RMAT1, STEP1, TEMP1, TEMP2, TEMP3, TEMP4, W, WI, Y1
      REAL RHMAX, RHTOL, RHO1, RHO2, T
C
      REAL ONE, ZERO
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER DINIT, DTYPE, DTINIT, D0INIT, F,
     1        F0, G, HC, IPIVOT, IVNEED, JCN, JTOL, LMAT,
     2        MODE, NEXTIV, NEXTV, NF0, NF1, NFCALL, NFGCAL,
     3        QTR, RDREQ, REGD, RESTOR, RMAT,
     4        RSPTOL, STEP, TOOBIG, VNEED
C
C  ***  IV SUBSCRIPT VALUES  ***
C
      PARAMETER (DTYPE=16, F0=13, G=28, HC=71, IPIVOT=76, IVNEED=3,
     1           JCN=66, JTOL=59, LMAT=42, MODE=35, NEXTIV=46, NEXTV=47,
     2           NFCALL=6, NF0=68, NF1=69, NFGCAL=7, QTR=77, RESTOR=9,
     3           RMAT=78, RDREQ=57, REGD=67, STEP=40, TOOBIG=2, VNEED=4)
C
C  ***  V SUBSCRIPT VALUES  ***
C
      PARAMETER (DINIT=38, DTINIT=39, D0INIT=40, F=10, RSPTOL=49)
      PARAMETER (ONE=1.E+0, ZERO=0.E+0)
      SAVE NEED1, NEED2
      DATA NEED1(1)/1/, NEED1(2)/0/, NEED2(1)/2/, NEED2(2)/0/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      LH = P * (P+1) / 2
      IF (IV(1) .EQ. 0) CALL  IVSET(1, IV, LIV, LV, V)
      PS1 = PS + 1
      IV1 = IV(1)
      IF (IV1 .GT. 2) GO TO 10
         W = IV(G) - N
         IV(RESTOR) = 0
         IF (IV(TOOBIG) .EQ. 0) GO TO (110, 120), IV1
         V(F) = V(F0)
         IF (IV1 .NE. 1) IV(1) = 2
         GO TO 40
C
C  ***  FRESH START OR RESTART -- CHECK INPUT INTEGERS  ***
C
 10   IF (ND .LT. PS) GO TO 340
      IF (PS .GT. P) GO TO 340
      IF (PS .LE. 0) GO TO 340
      IF (N .LE. 0) GO TO 340
      IF (IV1 .EQ. 14) GO TO 30
      IF (IV1 .GT. 16) GO TO 360
      IF (IV1 .LT. 12) GO TO 40
      IF (IV1 .EQ. 12) IV(1) = 13
      IF (IV(1) .NE. 13) GO TO 20
      IV(IVNEED) = IV(IVNEED) + P
      IV(VNEED) = IV(VNEED) + P*(P+13)/2 + 2*N + 4*PS
 20   CALL  G7ITB(B, D, X, IV, LIV, LV, P, PS, V, X, X)
      IF (IV(1) .NE. 14) GO TO 999
C
C  ***  STORAGE ALLOCATION  ***
C
      IV(IPIVOT) = IV(NEXTIV)
      IV(NEXTIV) = IV(IPIVOT) + P
      IV(G) = IV(NEXTV) + P + N
      IV(RMAT) = IV(G) + P + 4*PS
      IV(QTR) = IV(RMAT) + LH
      IV(JTOL) = IV(QTR) + P + N
      IV(JCN) = IV(JTOL) + 2*P
      IV(NEXTV) = IV(JCN) + P
C     *** TURN OFF COVARIANCE COMPUTATION ***
      IV(RDREQ) = 0
      IF (IV1 .EQ. 13) GO TO 999
C
 30   JTOL1 = IV(JTOL)
      IF (V(DINIT) .GE. ZERO) CALL  V7SCP(P, D, V(DINIT))
      IF (V(DTINIT) .GT. ZERO) CALL  V7SCP(P, V(JTOL1), V(DTINIT))
      I = JTOL1 + P
      IF (V(D0INIT) .GT. ZERO) CALL  V7SCP(P, V(I), V(D0INIT))
      IV(NF0) = 0
      IV(NF1) = 0
C
 40   G1 = IV(G)
      Y1 = G1 - (P + N)
      CALL  G7ITB(B, D, V(G1), IV, LIV, LV, P, PS, V, X, V(Y1))
      IF (IV(1) - 2) 50, 60, 350
C
 50   V(F) = ZERO
      IF (IV(NF1) .EQ. 0) GO TO 999
      IF (IV(RESTOR) .NE. 2) GO TO 999
      IV(NF0) = IV(NF1)
      CALL  V7CPY(N, RD, R)
      IV(REGD) = 0
      GO TO 999
C
 60   CALL  V7SCP(P, V(G1), ZERO)
      RMAT1 = IABS(IV(RMAT))
      QTR1 = IABS(IV(QTR))
      CALL  V7SCP(PS, V(QTR1), ZERO)
      IV(REGD) = 0
      CALL  V7SCP(PS, V(Y1), ZERO)
      CALL  V7SCP(LH, V(RMAT1), ZERO)
      IF (IV(RESTOR) .NE. 3) GO TO 70
         CALL  V7CPY(N, R, RD)
         IV(NF1) = IV(NF0)
 70   CALL RHO(NEED2, T, N, IV(NFGCAL), X(PS1), R, RD, RHOI, RHOR, V(W))
      IF (IV(NFGCAL) .GT. 0) GO TO 90
 80      IV(TOOBIG) = 1
         GO TO 40
 90   IF (IV(MODE) .LT. 0) GO TO 999
      DO 100 I = 1, N
 100     CALL  V2AXY(PS, V(Y1), R(I), DR(1,I), V(Y1))
      GO TO 999
C
C  ***  COMPUTE F(X)  ***
C
 110  I = IV(NFCALL)
      NEED1(2) = IV(NFGCAL)
      CALL RHO(NEED1, V(F), N, I, X(PS1), R, RD, RHOI, RHOR, V(W))
      IV(NF1) = I
      IF (I .LE. 0) GO TO 80
      GO TO 40
C
 120  G1 = IV(G)
C
C  ***  DECIDE WHETHER TO UPDATE D BELOW  ***
C
      I = IV(DTYPE)
      UPDATD = .FALSE.
      IF (I .LE. 0) GO TO 130
         IF (I .EQ. 1 .OR. IV(MODE) .LT. 0) UPDATD = .TRUE.
C
C  ***  COMPUTE RMAT AND QTR  ***
C
 130  QTR1 = IABS(IV(QTR))
      RMAT1 = IABS(IV(RMAT))
      IV(RMAT) = RMAT1
      IV(HC) = 0
      IV(NF0) = 0
      IV(NF1) = 0
      IF (IV(MODE) .LT. 0) GO TO 150
C
C  ***  ADJUST Y  ***
C
      Y1 = IV(G) - (P + N)
      WI = W
      STEP1 = IV(STEP)
      DO 140 I = 1, N
         T = V(WI) - RD(I)
         WI = WI + 1
         IF (T .NE. ZERO) CALL  V2AXY(PS, V(Y1),
     1                    T* D7TPR(PS,V(STEP1),DR(1,I)), DR(1,I), V(Y1))
 140     CONTINUE
C
C  ***  CHECK FOR NEGATIVE W COMPONENTS  ***
C
 150  J1 = W + N - 1
      DO 160 WI = W, J1
         IF (V(WI) .LT. ZERO) GO TO 230
 160     CONTINUE
C
C  ***  W IS NONNEGATIVE.  COMPUTE QR FACTORIZATION  ***
C  ***  AND, IF NECESSARY, USE SEMINORMAL EQUATIONS  ***
C
      RHMAX = ZERO
      RHTOL = V(RSPTOL)
      TEMP1 = G1 + P
      ZEROG = .TRUE.
      WI = W
      DO 190 I = 1, N
         RHO1 = R(I)
         RHO2 = V(WI)
         WI = WI + 1
         T =  SQRT(RHO2)
         IF (RHMAX .LT. RHO2) RHMAX = RHO2
         IF (RHO2 .GT. RHTOL*RHMAX) GO TO 170
C           *** SEMINORMAL EQUATIONS ***
            CALL  V2AXY(PS, V(G1), RHO1, DR(1,I), V(G1))
            RHO1 = ZERO
            ZEROG = .FALSE.
            GO TO 180
 170     RHO1 =  RHO1 / T
C        *** QR ACCUMULATION ***
 180     CALL  V7SCL(PS, V(TEMP1), T, DR(1,I))
         CALL  Q7ADR(PS, V(QTR1), V(RMAT1), V(TEMP1), RHO1)
 190     CONTINUE
C
C  ***  COMPUTE G FROM RMAT AND QTR  ***
C
      TEMP2 = TEMP1 + P
      CALL  L7VML(PS, V(TEMP1), V(RMAT1), V(QTR1))
      IF (ZEROG) GO TO 210
      IV(QTR) = -QTR1
      IF ( L7SVX(PS, V(RMAT1), V(TEMP2), V(TEMP2)) * RHTOL .GE.
     1     L7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2))) GO TO 220
         CALL  L7IVM(PS, V(TEMP2), V(RMAT1), V(G1))
C
C        *** SEMINORMAL EQUATIONS CORRECTION OF BJOERCK --
C        *** ONE CYCLE OF ITERATIVE REFINEMENT...
C
         TEMP3 = TEMP2 + PS
         TEMP4 = TEMP3 + PS
         CALL  L7ITV(PS, V(TEMP3), V(RMAT1), V(TEMP2))
         CALL  V7SCP(PS, V(TEMP4), ZERO)
         RHMAX = ZERO
         WI = W
         DO 200 I = 1, N
            RHO2 = V(WI)
            WI = WI + 1
            IF (RHMAX .LT. RHO2) RHMAX = RHO2
            RHO1 = ZERO
            IF (RHO2 .LE. RHTOL*RHMAX) RHO1 = R(I)
            T = RHO1 - RHO2* D7TPR(PS, V(TEMP3), DR(1,I))
            CALL  V2AXY(PS, V(TEMP4), T, DR(1,I), V(TEMP4))
 200        CONTINUE
         CALL  L7IVM(PS, V(TEMP3), V(RMAT1), V(TEMP4))
         CALL  V2AXY(PS, V(TEMP2), ONE, V(TEMP3), V(TEMP2))
         CALL  V2AXY(PS, V(QTR1), ONE, V(TEMP2), V(QTR1))
 210     IV(QTR) = QTR1
 220  CALL  V2AXY(PS, V(G1), ONE, V(TEMP1), V(G1))
      IF (PS .GE. P) GO TO 330
      GO TO 250
C
C  ***  INDEFINITE GN HESSIAN...  ***
C
 230  IV(RMAT) = -RMAT1
      IV(HC) = RMAT1
      CALL  O7PRD(N, LH, PS, V(RMAT1), V(W), DR, DR)
C
C  ***  COMPUTE GRADIENT  ***
C
      G1 = IV(G)
      DO 240 I = 1, N
 240     CALL  V2AXY(PS, V(G1), R(I), DR(1,I), V(G1))
      IF (PS .GE. P) GO TO 330
C
C  ***  COMPUTE GRADIENT COMPONENTS OF NUISANCE PARAMETERS ***
C
 250  K = P - PS
      J1 = 1
      G1 = G1 + PS
      DO 260 J = 1, K
         J1 = J1 + NN
         V(G1) =  VSUM(N, R(J1))
         G1 = G1 + 1
 260     CONTINUE
C
C  ***  COMPUTE HESSIAN COMPONENTS OF NUISANCE PARAMETERS  ***
C
      I = PS*PS1/2
      PSLEN = P*(P+1)/2 - I
      HN1 = RMAT1 + I
      CALL  V7SCP(PSLEN, V(HN1), ZERO)
      PMPS = P - PS
      K = HN1
      J1 = 1
      DO 290 II = 1, PMPS
         J1 = J1 + NN
         J = J1
         DO 270 I = 1, N
            CALL  V2AXY(PS, V(K), RD(J), DR(1,I), V(K))
            J = J + 1
 270        CONTINUE
         K = K + PS
         DO 280 I = 1, II
            J1 = J1 + NN
            V(K) =  VSUM(N, RD(J1))
            K = K + 1
 280        CONTINUE
 290     CONTINUE
      IF (IV(RMAT) .LE. 0) GO TO 330
      J = IV(LMAT)
      CALL  V7CPY(PSLEN, V(J), V(HN1))
      IF ( L7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2)) .LE. ZERO) GO TO 300
      CALL  L7SRT(PS1, P, V(RMAT1), V(RMAT1), I)
      IF (I .LE. 0) GO TO 310
C
C  *** HESSIAN IS NOT POSITIVE DEFINITE ***
C
 300  CALL  L7SQR(PS, V(RMAT1), V(RMAT1))
      CALL  V7CPY(PSLEN, V(HN1), V(J))
      IV(HC) = RMAT1
      IV(RMAT) = -RMAT1
      GO TO 330
C
C  *** NUISANCE PARS LEAVE HESSIAN POS. DEF.  GET REST OF QTR ***
C
 310  J = QTR1 + PS
      G1 = IV(G) + PS
      DO 320 I = PS1, P
         T =  D7TPR(I-1, V(HN1), V(QTR1))
         HN1 = HN1 + I
         V(J) = (V(G1) - T) / V(HN1-1)
         J = J + 1
         G1 = G1 + 1
 320     CONTINUE
 330  IF (UPDATD) CALL  D7UP5(D, IV, LIV, LV, P, PS, V)
      GO TO 40
C
C  ***  MISC. DETAILS  ***
C
C     ***  BAD N, ND, OR P  ***
C
 340  IV(1) = 66
      GO TO 360
C
C  ***  PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS  ***
C
 350  G1 = IV(G)
 360  CALL  ITSUM(D, V(G1), IV, LIV, LV, P, V, X)
C
 999  RETURN
C  ***  LAST LINE OF  RGLGB FOLLOWS  ***
      END
      SUBROUTINE  D7MLP(N, X, Y, Z, K)
C
C ***  SET X = DIAG(Y)**K * Z
C ***  FOR X, Z = LOWER TRIANG. MATRICES STORED COMPACTLY BY ROW
C ***  K = 1 OR -1.
C
      INTEGER N, K
      REAL X(*), Y(N), Z(*)
      INTEGER I, J, L
      REAL ONE, T
      DATA ONE/1.E+0/
C
      L = 1
      IF (K .GE. 0) GO TO 30
      DO 20 I = 1, N
         T = ONE / Y(I)
         DO 10 J = 1, I
            X(L) = T * Z(L)
            L = L + 1
 10         CONTINUE
 20      CONTINUE
      GO TO 999
C
 30   DO 50 I = 1, N
         T = Y(I)
         DO 40 J = 1, I
            X(L) = T * Z(L)
            L = L + 1
 40         CONTINUE
 50      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  D7MLP FOLLOWS  ***
      END
      SUBROUTINE  F7DHB(B, D, G, IRT, IV, LIV, LV, P, V, X)
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING
C  ***  AT V(IV(FDH)) = V(-IV(H)).  HONOR SIMPLE BOUNDS IN B.
C
C  ***  IF IV(COVREQ) .GE. 0 THEN  F7DHB USES GRADIENT DIFFERENCES,
C  ***  OTHERWISE FUNCTION DIFFERENCES.  STORAGE IN V IS AS IN  G7LIT.
C
C IRT VALUES...
C     1 = COMPUTE FUNCTION VALUE, I.E., V(F).
C     2 = COMPUTE G.
C     3 = DONE.
C
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER IRT, LIV, LV, P
      INTEGER IV(LIV)
      REAL B(2,P), D(P), G(P), V(LV), X(P)
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL OFFSID
      INTEGER GSAVE1, HES, HMI, HPI, HPM, I, K, KIND, L, M, MM1, MM1O2,
     1        NEWM1, PP1O2, STPI, STPM, STP0
      REAL DEL, DEL0, T, XM, XM1
      REAL HALF, HLIM, ONE, TWO, ZERO
C
C  ***  EXTERNAL SUBROUTINES  ***
C
      EXTERNAL  V7CPY,  V7SCP
C
C  V7CPY.... COPY ONE VECTOR TO ANOTHER.
C  V7SCP... COPY SCALAR TO ALL COMPONENTS OF A VECTOR.
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER COVREQ, DELTA, DELTA0, DLTFDC, F, FDH, FX, H, KAGQT, MODE,
     1        NFGCAL, SAVEI, SWITCH, TOOBIG, W, XMSAVE
C
      PARAMETER (HALF=0.5E+0, HLIM=0.1E+0, ONE=1.E+0, TWO=2.E+0,
     1           ZERO=0.E+0)
C
      PARAMETER (COVREQ=15, DELTA=52, DELTA0=44, DLTFDC=42, F=10,
     1           FDH=74, FX=53, H=56, KAGQT=33, MODE=35, NFGCAL=7,
     2           SAVEI=63, SWITCH=12, TOOBIG=2, W=65, XMSAVE=51)
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      IRT = 4
      KIND = IV(COVREQ)
      M = IV(MODE)
      IF (M .GT. 0) GO TO 10
         HES = IABS(IV(H))
         IV(H) = -HES
         IV(FDH) = 0
         IV(KAGQT) = -1
         V(FX) = V(F)
C        *** SUPPLY ZEROS IN CASE B(1,I) = B(2,I) FOR SOME I ***
         CALL  V7SCP(P*(P+1)/2, V(HES), ZERO)
 10   IF (M .GT. P) GO TO 999
      IF (KIND .LT. 0) GO TO 120
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND
C  ***  GRADIENT VALUES.
C
      GSAVE1 = IV(W) + P
      IF (M .GT. 0) GO TO 20
C        ***  FIRST CALL ON  F7DHB.  SET GSAVE = G, TAKE FIRST STEP  ***
         CALL  V7CPY(P, V(GSAVE1), G)
         IV(SWITCH) = IV(NFGCAL)
         GO TO 80
C
 20   DEL = V(DELTA)
      X(M) = V(XMSAVE)
      IF (IV(TOOBIG) .EQ. 0) GO TO 30
C
C     ***  HANDLE OVERSIZE V(DELTA)  ***
C
         DEL0 = V(DELTA0) *   MAX(ONE/D(M),  ABS(X(M)))
         DEL = HALF * DEL
         IF ( ABS(DEL/DEL0) .LE. HLIM) GO TO 140
C
 30   HES = -IV(H)
C
C  ***  SET  G = (G - GSAVE)/DEL  ***
C
      DEL = ONE / DEL
      DO 40 I = 1, P
         G(I) = DEL * (G(I) - V(GSAVE1))
         GSAVE1 = GSAVE1 + 1
 40      CONTINUE
C
C  ***  ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX  ***
C
      K = HES + M*(M-1)/2
      L = K + M - 2
      IF (M .EQ. 1) GO TO 60
C
C  ***  SET  H(I,M) = 0.5 * (H(I,M) + G(I))  FOR I = 1 TO M-1  ***
C
      MM1 = M - 1
      DO 50 I = 1, MM1
         IF (B(1,I) .LT. B(2,I)) V(K) = HALF * (V(K) + G(I))
         K = K + 1
 50      CONTINUE
C
C  ***  ADD  H(I,M) = G(I)  FOR I = M TO P  ***
C
 60   L = L + 1
      DO 70 I = M, P
         IF (B(1,I) .LT. B(2,I)) V(L) = G(I)
         L = L + I
 70      CONTINUE
C
 80   M = M + 1
      IV(MODE) = M
      IF (M .GT. P) GO TO 340
      IF (B(1,M) .GE. B(2,M)) GO TO 80
C
C  ***  CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE  ***
C
      DEL = V(DELTA0) *   MAX(ONE/D(M),  ABS(X(M)))
      XM = X(M)
      IF (XM .LT. ZERO) GO TO 90
         XM1 = XM + DEL
         IF (XM1 .LE. B(2,M)) GO TO 110
           XM1 = XM - DEL
           IF (XM1 .GE. B(1,M)) GO TO 100
           GO TO 280
 90    XM1 = XM - DEL
       IF (XM1 .GE. B(1,M)) GO TO 100
       XM1 = XM + DEL
       IF (XM1 .LE. B(2,M)) GO TO 110
       GO TO 280
C
 100  DEL = -DEL
 110  V(XMSAVE) = XM
      X(M) = XM1
      V(DELTA) = DEL
      IRT = 2
      GO TO 999
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY.
C
 120  STP0 = IV(W) + P - 1
      MM1 = M - 1
      MM1O2 = M*MM1/2
      HES = -IV(H)
      IF (M .GT. 0) GO TO 130
C        ***  FIRST CALL ON  F7DHB.  ***
         IV(SAVEI) = 0
         GO TO 240
C
 130  IF (IV(TOOBIG) .EQ. 0) GO TO 150
C        ***  PUNT IN THE EVENT OF AN OVERSIZE STEP  ***
 140     IV(FDH) = -2
         GO TO 350
 150  I = IV(SAVEI)
      IF (I .GT. 0) GO TO 190
C
C  ***  SAVE F(X + STP(M)*E(M)) IN H(P,M)  ***
C
      PP1O2 = P * (P-1) / 2
      HPM = HES + PP1O2 + MM1
      V(HPM) = V(F)
C
C  ***  START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H.  ***
C
      NEWM1 = 1
      GO TO 260
 160  HMI = HES + MM1O2
      IF (MM1 .EQ. 0) GO TO 180
      HPI = HES + PP1O2
      DO 170 I = 1, MM1
         T = ZERO
         IF (B(1,I) .LT. B(2,I)) T = V(FX) - (V(F) + V(HPI))
         V(HMI) = T
         HMI = HMI + 1
         HPI = HPI + 1
 170     CONTINUE
 180  V(HMI) = V(F) - TWO*V(FX)
      IF (OFFSID) V(HMI) = V(FX) - TWO*V(F)
C
C  ***  COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H.  ***
C
      I = 0
      GO TO 200
C
 190  X(I) = V(DELTA)
C
C  ***  FINISH COMPUTING H(M,I)  ***
C
      STPI = STP0 + I
      HMI = HES + MM1O2 + I - 1
      STPM = STP0 + M
      V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM))
 200  I = I + 1
      IF (I .GT. M) GO TO 230
         IF (B(1,I) .LT. B(2,I)) GO TO 210
         GO TO 200
C
 210  IV(SAVEI) = I
      STPI = STP0 + I
      V(DELTA) = X(I)
      X(I) = X(I) + V(STPI)
      IRT = 1
      IF (I .LT. M) GO TO 999
      NEWM1 = 2
      GO TO 260
 220  X(M) = V(XMSAVE) - DEL
      IF (OFFSID) X(M) = V(XMSAVE) + TWO*DEL
      GO TO 999
C
 230  IV(SAVEI) = 0
      X(M) = V(XMSAVE)
C
 240  M = M + 1
      IV(MODE) = M
      IF (M .GT. P) GO TO 330
      IF (B(1,M) .LT. B(2,M)) GO TO 250
      GO TO 240
C
C  ***  PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H.
C  ***  COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN
C  ***  F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR.
C
 250  V(XMSAVE) = X(M)
      NEWM1 = 3
 260  XM = V(XMSAVE)
      DEL = V(DLTFDC) *   MAX(ONE/D(M),  ABS(XM))
      XM1 = XM + DEL
      OFFSID = .FALSE.
      IF (XM1 .LE. B(2,M)) GO TO 270
         OFFSID = .TRUE.
         XM1 = XM - DEL
         IF (XM - TWO*DEL .GE. B(1,M)) GO TO 300
         GO TO 280
 270   IF (XM-DEL .GE. B(1,M)) GO TO 290
       OFFSID = .TRUE.
       IF (XM + TWO*DEL .LE. B(2,M)) GO TO 310
C
 280  IV(FDH) = -2
      GO TO 350
C
 290  IF (XM .GE. ZERO) GO TO 310
      XM1 = XM - DEL
 300  DEL = -DEL
 310  GO TO (160, 220, 320), NEWM1
 320  X(M) = XM1
      STPM = STP0 + M
      V(STPM) = DEL
      IRT = 1
      GO TO 999
C
C  ***  HANDLE SPECIAL CASE OF B(1,P) = B(2,P) -- CLEAR SCRATCH VALUES
C  ***  FROM LAST ROW OF FDH...
C
 330  IF (B(1,P) .LT. B(2,P)) GO TO 340
         I = HES + P*(P-1)/2
         CALL  V7SCP(P, V(I), ZERO)
C
C  ***  RESTORE V(F), ETC.  ***
C
 340  IV(FDH) = HES
 350  V(F) = V(FX)
      IRT = 3
      IF (KIND .LT. 0) GO TO 999
         IV(NFGCAL) = IV(SWITCH)
         GSAVE1 = IV(W) + P
         CALL  V7CPY(P, G, V(GSAVE1))
         GO TO 999
C
 999  RETURN
C  ***  LAST LINE OF  F7DHB FOLLOWS  ***
      END
      SUBROUTINE  G7ITB(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)
      REAL 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         G7ITB IS SIMILAR TO  G7LIT, EXCEPT FOR THE EXTRA PARAMETER B
C     --  G7ITB ENFORCES THE BOUNDS  B(1,I) .LE. X(I) .LE. B(2,I),
C     I = 1(1)P.
C         G7ITB 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  G7ITB 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      G7ITB 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,  G7ITB 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  G7ITB 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   N2GB (AND NL2SOL), EXCEPT THAT V CAN BE SHORTER
C     (SINCE THE PART OF V THAT   N2GB USES FOR STORING D, J, AND R IS
C     NOT NEEDED).  MOREOVER, COMPARED WITH   N2GB (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       N2GB (AND   N2FB), ARE NOT REFERENCED BY  G7ITB OR THE
C     SUBROUTINES IT CALLS.
C
C        WHEN  G7ITB IS FIRST CALLED, I.E., WHEN  G7ITB IS CALLED WITH
C     IV(1) = 0 OR 12, V(F), G, AND HC NEED NOT BE INITIALIZED.  TO
C     OBTAIN THESE STARTING VALUES,  G7ITB 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  G7ITB 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  G7ITB 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  G7ITB 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  G7ITB 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  G7ITB 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  G7ITB 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
      REAL E, GI, STTSST, T, T1, XI
C
C     ***  CONSTANTS  ***
C
      REAL HALF, NEGONE, ONE, ONEP2, ZERO
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      LOGICAL STOPX
      REAL  D7TPR,  RLDST,  V2NRM
      EXTERNAL  A7SST,  D7TPR,  F7DHB,  G7QSB,I7COPY, I7PNVR, I7SHFT,
     1         ITSUM,  L7MSB,  L7SQR,  L7TVM, L7VML, PARCK,  Q7RSH,
     2          RLDST,  S7DMP,  S7IPR,  S7LUP,  S7LVM, STOPX,  V2NRM,
     3         V2AXY, V7CPY,  V7IPR,  V7SCP,  V7VMP
C
C  A7SST.... ASSESSES CANDIDATE STEP.
C  D7TPR... RETURNS INNER PRODUCT OF TWO VECTORS.
C  F7DHB... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR INIT. S MATRIX).
C  G7QSB... 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  ITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X.
C  L7MSB... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL).
C  L7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L.
C  L7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C  L7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C  PARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS.
C  Q7RSH... SHIFTS A QR FACTORIZATION.
C  RLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE.
C  S7DMP... MULTIPLIES A SYM. MATRIX FORE AND AFT BY A DIAG. MATRIX.
C  S7IPR... APPLIES PERMUTATION TO (LOWER TRIANG. OF) SYM. MATRIX.
C  S7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI-
C             ANGLE OF A SYMMETRIC MATRIX.
C  S7LVM... MULTIPLIES COMPACTLY STORED SYM. MATRIX TIMES VECTOR.
C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED.
C  V2NRM... RETURNS THE 2-NORM OF A VECTOR.
C  V2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER.
C  V7CPY.... COPIES ONE VECTOR TO ANOTHER.
C  V7IPR... APPLIES A PERMUTATION TO A VECTOR.
C  V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C  V7VMP... 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
      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  ***  V SUBSCRIPT VALUES  ***
C
      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
      PARAMETER (HALF=0.5E+0, NEGONE=-1.E+0, ONE=1.E+0, ONEP2=1.2E+0,
     1           ZERO=0.E+0)
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  PARCK(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  V7SCP(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  Q7RSH(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  V7VMP(P, V(DIG1), G, D, -1)
      CALL  V7IPR(P, IV(IPI), V(DIG1))
      V(DGNORM) =  V2NRM(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  V2AXY(PP1O2, V(S1), NEGONE, V(HC1), V(H1))
         GO TO 140
 130  RMAT1 = IV(RMAT)
      LMAT1 = IV(LMAT)
      CALL  L7SQR(P, V(LMAT1), V(RMAT1))
      IPI = IV(IPIVOT)
      IPIV1 = IV(PERM) + P
      CALL I7PNVR(P, IV(IPIV1), IV(IPI))
      CALL  S7IPR(P, IV(IPIV1), V(LMAT1))
      CALL  V2AXY(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  V7SCP(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  ITSUM(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) *  V2NRM(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  V7CPY(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  L7MSB(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  G7QTS 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  L7SQR(P1, V(H1), V(RMAT1))
            GO TO 300
 290     CALL  V7CPY(P*(P+1)/2, V(H1), V(J))
         CALL  S7IPR(P, IV(IPI), V(H1))
 300     IF (IV(MODEL) .EQ. 1) GO TO 310
            LMAT1 = IV(LMAT)
            S1 = IV(S)
            CALL  V7CPY(P*(P+1)/2, V(LMAT1), V(S1))
            CALL  S7IPR(P, IV(IPI), V(LMAT1))
            CALL  V2AXY(P1LEN, V(H1), ONE, V(LMAT1), V(H1))
 310     CALL  V7CPY(P, V(TD1), D)
         CALL  V7IPR(P, IV(IPI), V(TD1))
         CALL  S7DMP(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  G7QSB(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
         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  V2AXY(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) =  RLDST(P, D, X, V(X01))
      CALL  A7SST(IV, LIV, LV, V)
      STEP1 = IV(STEP)
      LSTGST = X01 + P
      I = IV(RESTOR) + 1
      GO TO (410, 380, 390, 400), I
 380  CALL  V7CPY(P, X, V(X01))
      GO TO 410
 390   CALL  V7CPY(P, V(LSTGST), V(STEP1))
       GO TO 410
 400     CALL  V7CPY(P, V(STEP1), V(LSTGST))
         CALL  V2AXY(P, X, ONE, V(STEP1), V(X01))
         V(RELDX) =  RLDST(P, D, X, V(X01))
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  V7CPY(NVSAVE, V, V(L))
 420  CALL  V2AXY(P, V(STEP1), NEGONE, V(X01), X)
      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  S7LVM(PS, Y, V(S1), V(STEP1))
      STTSST = HALF *  D7TPR(PS, V(STEP1), Y)
      IF (IV(MODEL) .EQ. 1) STTSST = -STTSST
      IF ( ABS(E + STTSST) * V(FUZZ) .GE.  ABS(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  V7CPY(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  S7LVM(P, V(TEMP1), V(HC1), V(STEP1))
              GO TO 490
 480     RMAT1 = IV(RMAT)
         IPIV0 = IV(IPIVOT)
         CALL  V7CPY(P, V(TEMP1), V(STEP1))
         CALL  V7IPR(P, IV(IPIV0), V(TEMP1))
         CALL  L7TVM(P, V(TEMP1), V(RMAT1), V(TEMP1))
         CALL  L7VML(P, V(TEMP1), V(RMAT1), V(TEMP1))
         IPIV1 = IV(PERM) + P
         CALL I7PNVR(P, IV(IPIV1), IV(IPIV0))
         CALL  V7IPR(P, IV(IPIV1), V(TEMP1))
C
 490     IF (STPMOD .EQ. 1) GO TO 500
              S1 = IV(S)
              CALL  S7LVM(PS, V(TEMP2), V(S1), V(STEP1))
              CALL  V2AXY(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  V7CPY(P, V(G01), G)
      GO TO 690
C
C  ***  INITIALIZATIONS -- G0 = G - G0, ETC.  ***
C
 510  G01 = IV(W)
      CALL  V2AXY(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 ( V2NRM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4))  GO TO 530
              IF ( D7TPR(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  V2AXY(PS, Y, NEGONE, Y, G)
C
C  ***  DETERMINE SIZING FACTOR V(SIZE)  ***
C
C     ***  SET TEMP1 = S * STEP  ***
      S1 = IV(S)
      CALL  S7LVM(PS, V(TEMP1), V(S1), V(STEP1))
C
      T1 =  ABS( D7TPR(PS, V(STEP1), V(TEMP1)))
      T =  ABS( D7TPR(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  S7LVM(PS, V(G01), V(HC1), V(STEP1))
         GO TO 560
C
 550  RMAT1 = IV(RMAT)
      IPIV0 = IV(IPIVOT)
      CALL  V7CPY(P, V(G01), V(STEP1))
      I = G01 + PS
      IF (PS .LT. P) CALL  V7SCP(P-PS, V(I), ZERO)
      CALL  V7IPR(P, IV(IPIV0), V(G01))
      CALL  L7TVM(P, V(G01), V(RMAT1), V(G01))
      CALL  L7VML(P, V(G01), V(RMAT1), V(G01))
      IPIV1 = IV(PERM) + P
      CALL I7PNVR(P, IV(IPIV1), IV(IPIV0))
      CALL  V7IPR(P, IV(IPIV1), V(G01))
C
 560  CALL  V2AXY(PS, V(G01), ONE, Y, V(G01))
C
C  ***  UPDATE S  ***
C
      CALL  S7LUP(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  F7DHB(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  V7CPY(P*(P+1)/2, V(H1), V(HC1))
           GO TO 660
 650  RMAT1 = IV(RMAT)
      CALL  L7SQR(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  G7ITB FOLLOWS  ***
      END
      SUBROUTINE  G7QSB(B, D, DIHDI, G, IPIV, IPIV1, IPIV2, KA, L, LV,
     1                  P, P0, PC, STEP, TD, TG, V, W, X, X0)
C
C  ***  COMPUTE HEURISTIC BOUNDED NEWTON STEP  ***
C
      INTEGER KA, LV, P, P0, PC
      INTEGER IPIV(P), IPIV1(P), IPIV2(P)
      REAL B(2,P), D(P), DIHDI(1), G(P), L(1),
     1                 STEP(P,2), TD(P), TG(P), V(LV), W(P), X0(P), X(P)
C     DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2)
C
      REAL  D7TPR
      EXTERNAL  D7TPR, G7QTS,  S7BQN,  S7IPR, V7CPY,  V7IPR,
     1          V7SCP,  V7VMP
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER K, KB, KINIT, NS, P1, P10
      REAL DS0, NRED, PRED, RAD
      REAL ZERO
C
C  ***  V SUBSCRIPTS  ***
C
      INTEGER DST0, DSTNRM, GTSTEP, NREDUC, PREDUC, RADIUS
C
      PARAMETER (DST0=3, DSTNRM=2, GTSTEP=4, NREDUC=6, PREDUC=7,
     1           RADIUS=8)
      DATA ZERO/0.E+0/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      P1 = PC
      IF (KA .LT. 0) GO TO 10
         NRED = V(NREDUC)
         DS0 = V(DST0)
         GO TO 20
 10   P0 = 0
      KA = -1
C
 20   KINIT = -1
      IF (P0 .EQ. P1) KINIT = KA
      CALL  V7CPY(P, X, X0)
      PRED = ZERO
      RAD = V(RADIUS)
      KB = -1
      V(DSTNRM) = ZERO
      IF (P1 .GT. 0) GO TO 30
         NRED = ZERO
         DS0 = ZERO
         CALL  V7SCP(P, STEP, ZERO)
         GO TO 60
C
 30   CALL  V7CPY(P, TD, D)
      CALL  V7IPR(P, IPIV, TD)
      CALL  V7VMP(P, TG, G, D, -1)
      CALL  V7IPR(P, IPIV, TG)
 40   K = KINIT
      KINIT = -1
      V(RADIUS) = RAD - V(DSTNRM)
      CALL  G7QTS(TD, TG, DIHDI, K, L, P1, STEP, V, W)
      P0 = P1
      IF (KA .GE. 0) GO TO 50
         NRED = V(NREDUC)
         DS0 = V(DST0)
C
 50   KA = K
      V(RADIUS) = RAD
      P10 = P1
      CALL  S7BQN(B, D, STEP(1,2), IPIV, IPIV1, IPIV2, KB, L, LV,
     1            NS, P, P1, STEP, TD, TG, V, W, X, X0)
      IF (NS .GT. 0) CALL  S7IPR(P10, IPIV1, DIHDI)
      PRED = PRED + V(PREDUC)
      IF (NS .NE. 0) P0 = 0
      IF (KB .LE. 0) GO TO 40
C
 60   V(DST0) = DS0
      V(NREDUC) = NRED
      V(PREDUC) = PRED
      V(GTSTEP) =  D7TPR(P, G, STEP)
C
 999  RETURN
C  ***  LAST LINE OF  G7QSB FOLLOWS  ***
      END
      SUBROUTINE  H2RFA(N, A, B, X, Y, Z)
C
C  ***  APPLY 2X2 HOUSEHOLDER REFLECTION DETERMINED BY X, Y, Z TO
C  ***  N-VECTORS A, B  ***
C
      INTEGER N
      REAL A(N), B(N), X, Y, Z
      INTEGER I
      REAL T
      DO 10 I = 1, N
         T = A(I)*X + B(I)*Y
         A(I) = A(I) + T
         B(I) = B(I) + T*Z
 10      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  H2RFA FOLLOWS  ***
      END
      REAL FUNCTION  H2RFG(A, B, X, Y, Z)
C
C  ***  DETERMINE X, Y, Z SO  I + (1,Z)**T * (X,Y)  IS A 2X2
C  ***  HOUSEHOLDER REFLECTION SENDING (A,B)**T INTO (C,0)**T,
C  ***  WHERE  C = -SIGN(A)*SQRT(A**2 + B**2)  IS THE VALUE  H2RFG
C  ***  RETURNS.
C
      REAL A, B, X, Y, Z
C
      REAL A1, B1, C, T
      REAL ZERO
      DATA ZERO/0.E+0/
C
C  ***  BODY  ***
C
      IF (B .NE. ZERO) GO TO 10
         X = ZERO
         Y = ZERO
         Z = ZERO
          H2RFG = A
         GO TO 999
 10   T =  ABS(A) +  ABS(B)
      A1 = A / T
      B1 = B / T
      C =  SQRT(A1**2 + B1**2)
      IF (A1 .GT. ZERO) C = -C
      A1 = A1 - C
      Z = B1 / A1
      X = A1 / C
      Y = B1 / C
       H2RFG = T * C
 999  RETURN
C  ***  LAST LINE OF  H2RFG FOLLOWS  ***
      END
      SUBROUTINE  L7MSB(B, D, G, IERR, IPIV, IPIV1, IPIV2, KA, LMAT,
     1                  LV, P, P0, PC, QTR, RMAT, STEP, TD, TG, V,
     2                  W, WLM, X, X0)
C
C  ***  COMPUTE HEURISTIC BOUNDED NEWTON STEP  ***
C
      INTEGER IERR, KA, LV, P, P0, PC
      INTEGER IPIV(P), IPIV1(P), IPIV2(P)
      REAL B(2,P), D(P), G(P), LMAT(1), QTR(P), RMAT(1),
     1                 STEP(P,3), TD(P), TG(P), V(LV), W(P), WLM(1),
     2                 X0(P), X(P)
C     DIMENSION LMAT(P*(P+1)/2), RMAT(P*(P+1)/2), WLM(P*(P+5)/2 + 4)
C
      REAL  D7TPR
      EXTERNAL  D7MLP,  D7TPR,  L7MST,  L7TVM,  Q7RSH,  S7BQN,
     1         V2AXY, V7CPY,  V7IPR,  V7SCP,  V7VMP
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, J, K, K0, KB, KINIT, L, NS, P1, P10, P11
      REAL DS0, NRED, PRED, RAD
      REAL ONE, ZERO
C
C  ***  V SUBSCRIPTS  ***
C
      INTEGER DST0, DSTNRM, GTSTEP, NREDUC, PREDUC, RADIUS
C
      PARAMETER (DST0=3, DSTNRM=2, GTSTEP=4, NREDUC=6, PREDUC=7,
     1           RADIUS=8)
      DATA ONE/1.E+0/, ZERO/0.E+0/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      P1 = PC
      IF (KA .LT. 0) GO TO 10
         NRED = V(NREDUC)
         DS0 = V(DST0)
         GO TO 20
 10   P0 = 0
      KA = -1
C
 20   KINIT = -1
      IF (P0 .EQ. P1) KINIT = KA
      CALL  V7CPY(P, X, X0)
      CALL  V7CPY(P, TD, D)
C     *** USE STEP(1,3) AS TEMP. COPY OF QTR ***
      CALL  V7CPY(P, STEP(1,3), QTR)
      CALL  V7IPR(P, IPIV, TD)
      PRED = ZERO
      RAD = V(RADIUS)
      KB = -1
      V(DSTNRM) = ZERO
      IF (P1 .GT. 0) GO TO 30
         NRED = ZERO
         DS0 = ZERO
         CALL  V7SCP(P, STEP, ZERO)
         GO TO 90
C
 30   CALL  V7VMP(P, TG, G, D, -1)
      CALL  V7IPR(P, IPIV, TG)
      P10 = P1
 40   K = KINIT
      KINIT = -1
      V(RADIUS) = RAD - V(DSTNRM)
      CALL  V7VMP(P1, TG, TG, TD, 1)
      DO 50 I = 1, P1
 50      IPIV1(I) = I
      K0 = MAX0(0, K)
      CALL  L7MST(TD, TG, IERR, IPIV1, K, P1, STEP(1,3), RMAT, STEP,
     1            V, WLM)
      CALL  V7VMP(P1, TG, TG, TD, -1)
      P0 = P1
      IF (KA .GE. 0) GO TO 60
         NRED = V(NREDUC)
         DS0 = V(DST0)
C
 60   KA = K
      V(RADIUS) = RAD
      L = P1 + 5
      IF (K .LE. K0) CALL  D7MLP(P1, LMAT, TD, RMAT, -1)
      IF (K .GT. K0) CALL  D7MLP(P1, LMAT, TD, WLM(L), -1)
      CALL  S7BQN(B, D, STEP(1,2), IPIV, IPIV1, IPIV2, KB, LMAT,
     1            LV, NS, P, P1, STEP, TD, TG, V, W, X, X0)
      PRED = PRED + V(PREDUC)
      IF (NS .EQ. 0) GO TO 80
      P0 = 0
C
C  ***  UPDATE RMAT AND QTR  ***
C
      P11 = P1 + 1
      L = P10 + P11
      DO 70 K = P11, P10
         J = L - K
         I = IPIV2(J)
         IF (I .LT. J) CALL  Q7RSH(I, J, .TRUE., QTR, RMAT, W)
 70      CONTINUE
C
 80   IF (KB .GT. 0) GO TO 90
C
C  ***  UPDATE LOCAL COPY OF QTR  ***
C
      CALL  V7VMP(P10, W, STEP(1,2), TD, -1)
      CALL  L7TVM(P10, W, LMAT, W)
      CALL  V2AXY(P10, STEP(1,3), ONE, W, QTR)
      GO TO 40
C
 90   V(DST0) = DS0
      V(NREDUC) = NRED
      V(PREDUC) = PRED
      V(GTSTEP) =  D7TPR(P, G, STEP)
C
 999  RETURN
C  ***  LAST LINE OF  L7MSB FOLLOWS  ***
      END
      SUBROUTINE  Q7RSH(K, P, HAVQTR, QTR, R, W)
C
C  ***  PERMUTE COLUMN K OF R TO COLUMN P, MODIFY QTR ACCORDINGLY  ***
C
      LOGICAL HAVQTR
      INTEGER K, P
      REAL QTR(P), R(1), W(P)
C     DIMSNSION R(P*(P+1)/2)
C
      REAL  H2RFG
      EXTERNAL  H2RFA,  H2RFG, V7CPY
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, I1, J, JM1, JP1, J1, KM1, K1, PM1
      REAL A, B, T, WJ, X, Y, Z, ZERO
C
      DATA ZERO/0.0E+0/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      IF (K .GE. P) GO TO 999
      KM1 = K - 1
      K1 = K * KM1 / 2
      CALL  V7CPY(K, W, R(K1+1))
      WJ = W(K)
      PM1 = P - 1
      J1 = K1 + KM1
      DO 50 J = K, PM1
         JM1 = J - 1
         JP1 = J + 1
         IF (JM1 .GT. 0) CALL  V7CPY(JM1, R(K1+1), R(J1+2))
         J1 = J1 + JP1
         K1 = K1 + J
         A = R(J1)
         B = R(J1+1)
         IF (B .NE. ZERO) GO TO 10
              R(K1) = A
              X = ZERO
              Z = ZERO
              GO TO 40
 10      R(K1) =  H2RFG(A, B, X, Y, Z)
         IF (J .EQ. PM1) GO TO 30
         I1 = J1
         DO 20 I = JP1, PM1
              I1 = I1 + I
              CALL  H2RFA(1, R(I1), R(I1+1), X, Y, Z)
 20           CONTINUE
 30      IF (HAVQTR) CALL  H2RFA(1, QTR(J), QTR(JP1), X, Y, Z)
 40      T = X * WJ
         W(J) = WJ + T
         WJ = T * Z
 50      CONTINUE
      W(P) = WJ
      CALL  V7CPY(P, R(K1+1), W)
 999  RETURN
      END
      SUBROUTINE  S7BQN(B, D, DST, IPIV, IPIV1, IPIV2, KB, L, LV, NS,
     1                  P, P1, STEP, TD, TG, V, W, X, X0)
C
C  ***  COMPUTE BOUNDED MODIFIED NEWTON STEP  ***
C
      INTEGER KB, LV, NS, P, P1
      INTEGER IPIV(P), IPIV1(P), IPIV2(P)
      REAL B(2,P), D(P), DST(P), L(1),
     1                 STEP(P), TD(P), TG(P), V(LV), W(P), X(P),
     2                 X0(P)
C     DIMENSION L(P*(P+1)/2)
C
      REAL  D7TPR,  R7MDC,  V2NRM
      EXTERNAL  D7TPR, I7SHFT,  L7ITV,  L7IVM,  Q7RSH,  R7MDC,  V2NRM,
     1         V2AXY, V7CPY,  V7IPR,  V7SCP,  V7SHF
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, J, K, P0, P1M1
      REAL ALPHA, DST0, DST1, DSTMAX, DSTMIN, DX, GTS, T,
     1                 TI, T1, XI
      REAL FUDGE, HALF, MEPS2, ONE, TWO, ZERO
C
C  ***  V SUBSCRIPTS  ***
C
      INTEGER DSTNRM, GTSTEP, PHMNFC, PHMXFC, PREDUC, RADIUS, STPPAR
C
      PARAMETER (DSTNRM=2, GTSTEP=4, PHMNFC=20, PHMXFC=21, PREDUC=7,
     1           RADIUS=8, STPPAR=5)
      SAVE MEPS2
C
      DATA FUDGE/1.0001E+0/, HALF/0.5E+0/, MEPS2/0.E+0/,
     1     ONE/1.0E+0/, TWO/2.E+0/, ZERO/0.E+0/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      DSTMAX = FUDGE * (ONE + V(PHMXFC)) * V(RADIUS)
      DSTMIN = (ONE + V(PHMNFC)) * V(RADIUS)
      DST1 = ZERO
      IF (MEPS2 .LE. ZERO) MEPS2 = TWO *  R7MDC(3)
      P0 = P1
      NS = 0
      DO 10 I = 1, P
         IPIV1(I) = I
         IPIV2(I) = I
 10      CONTINUE
      DO 20 I = 1, P1
 20      W(I) = -STEP(I) * TD(I)
      ALPHA =  ABS(V(STPPAR))
      V(PREDUC) = ZERO
      GTS = -V(GTSTEP)
      IF (KB .LT. 0) CALL  V7SCP(P, DST, ZERO)
      KB = 1
C
C     ***  -W = D TIMES RESTRICTED NEWTON STEP FROM X + DST/D.
C
C     ***  FIND T SUCH THAT X - T*W IS STILL FEASIBLE.
C
 30   T = ONE
      K = 0
      DO 60 I = 1, P1
         J = IPIV(I)
         DX = W(I) / D(J)
         XI = X(J) - DX
         IF (XI .LT. B(1,J)) GO TO 40
         IF (XI .LE. B(2,J)) GO TO 60
              TI = ( X(J)  -  B(2,J) ) / DX
              K = I
              GO TO 50
 40      TI = ( X(J)  -  B(1,J) ) / DX
              K = -I
 50      IF (T .LE. TI) GO TO 60
              T = TI
 60      CONTINUE
C
      IF (P .GT. P1) CALL  V7CPY(P-P1, STEP(P1+1), DST(P1+1))
      CALL  V2AXY(P1, STEP, -T, W, DST)
      DST0 = DST1
      DST1 =  V2NRM(P, STEP)
C
C  ***  CHECK FOR OVERSIZE STEP  ***
C
      IF (DST1 .LE. DSTMAX) GO TO 80
      IF (P1 .GE. P0) GO TO 70
         IF (DST0 .LT. DSTMIN) KB = 0
         GO TO 110
C
 70   K = 0
C
C  ***  UPDATE DST, TG, AND V(PREDUC)  ***
C
 80   V(DSTNRM) = DST1
      CALL  V7CPY(P1, DST, STEP)
      T1 = ONE - T
      DO 90 I = 1, P1
 90      TG(I) = T1 * TG(I)
      IF (ALPHA .GT. ZERO) CALL  V2AXY(P1, TG, T*ALPHA, W, TG)
      V(PREDUC) = V(PREDUC) + T*((ONE - HALF*T)*GTS +
     1                        HALF*ALPHA*T* D7TPR(P1,W,W))
      IF (K .EQ. 0) GO TO 110
C
C     ***  PERMUTE L, ETC. IF NECESSARY  ***
C
      P1M1 = P1 - 1
      J = IABS(K)
      IF (J .EQ. P1) GO TO 100
         NS = NS + 1
         IPIV2(P1) = J
         CALL  Q7RSH(J, P1, .FALSE., TG, L, W)
         CALL I7SHFT(P1, J, IPIV)
         CALL I7SHFT(P1, J, IPIV1)
         CALL  V7SHF(P1, J, TG)
         CALL  V7SHF(P1, J, DST)
 100  IF (K .LT. 0) IPIV(P1) = -IPIV(P1)
      P1 = P1M1
      IF (P1 .LE. 0) GO TO 110
      CALL  L7IVM(P1, W, L, TG)
      GTS =  D7TPR(P1, W, W)
      CALL  L7ITV(P1, W, L, W)
      GO TO 30
C
C     ***  UNSCALE STEP  ***
C
 110  DO 120 I = 1, P
         J = IABS(IPIV(I))
         STEP(J) = DST(I) / D(J)
 120     CONTINUE
C
C  ***  FUDGE STEP TO ENSURE THAT IT FORCES APPROPRIATE COMPONENTS
C  ***  TO THEIR BOUNDS  ***
C
      IF (P1 .GE. P0) GO TO 150
      K = P1 + 1
      DO 140 I = K, P0
         J = IPIV(I)
         T = MEPS2
         IF (J .GT. 0) GO TO 130
            T = -T
            J = -J
            IPIV(I) = J
 130     T = T *   MAX( ABS(X(J)),  ABS(X0(J)))
         STEP(J) = STEP(J) + T
 140     CONTINUE
C
 150  CALL  V2AXY(P, X, ONE, STEP, X0)
      IF (NS .GT. 0) CALL  V7IPR(P0, IPIV1, TD)
 999  RETURN
C  ***  LAST LINE OF  S7BQN FOLLOWS  ***
      END
      SUBROUTINE  S7DMP(N, X, Y, Z, K)
C
C ***  SET X = DIAG(Z)**K * Y * DIAG(Z)**K
C ***  FOR X, Y = COMPACTLY STORED LOWER TRIANG. MATRICES
C ***  K = 1 OR -1.
C
      INTEGER N, K
      REAL X(*), Y(*), Z(N)
      INTEGER I, J, L
      REAL ONE, T
      DATA ONE/1.E+0/
C
      L = 1
      IF (K .GE. 0) GO TO 30
      DO 20 I = 1, N
         T = ONE / Z(I)
         DO 10 J = 1, I
            X(L) = T * Y(L) / Z(J)
            L = L + 1
 10         CONTINUE
 20      CONTINUE
      GO TO 999
C
 30   DO 50 I = 1, N
         T = Z(I)
         DO 40 J = 1, I
            X(L) = T * Y(L) * Z(J)
            L = L + 1
 40         CONTINUE
 50      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  S7DMP FOLLOWS  ***
      END
      SUBROUTINE  S7IPR(P, IP, H)
C
C  APPLY THE PERMUTATION DEFINED BY IP TO THE ROWS AND COLUMNS OF THE
C  P X P SYMMETRIC MATRIX WHOSE LOWER TRIANGLE IS STORED COMPACTLY IN H.
C  THUS H.OUTPUT(I,J) = H.INPUT(IP(I), IP(J)).
C
      INTEGER P
      INTEGER IP(P)
      REAL H(1)
C
      INTEGER I, J, J1, JM, K, K1, KK, KM, KMJ, L, M
      REAL T
C
C ***  BODY  ***
C
      DO 90 I = 1, P
         J = IP(I)
         IF (J .EQ. I) GO TO 90
         IP(I) = IABS(J)
         IF (J .LT. 0) GO TO 90
         K = I
 10         J1 = J
            K1 = K
            IF (J .LE. K) GO TO 20
               J1 = K
               K1 = J
 20         KMJ = K1-J1
            L = J1-1
            JM = J1*L/2
            KM = K1*(K1-1)/2
            IF (L .LE. 0) GO TO 40
               DO 30 M = 1, L
                  JM = JM+1
                  T = H(JM)
                  KM = KM+1
                  H(JM) = H(KM)
                  H(KM) = T
 30               CONTINUE
 40         KM = KM+1
            KK = KM+KMJ
            JM = JM+1
            T = H(JM)
            H(JM) = H(KK)
            H(KK) = T
            J1 = L
            L = KMJ-1
            IF (L .LE. 0) GO TO 60
               DO 50 M = 1, L
                  JM = JM+J1+M
                  T = H(JM)
                  KM = KM+1
                  H(JM) = H(KM)
                  H(KM) = T
 50               CONTINUE
 60         IF (K1 .GE. P) GO TO 80
               L = P-K1
               K1 = K1-1
               KM = KK
               DO 70 M = 1, L
                  KM = KM+K1+M
                  JM = KM-KMJ
                  T = H(JM)
                  H(JM) = H(KM)
                  H(KM) = T
 70               CONTINUE
 80         K = J
            J = IP(K)
            IP(K) = -J
            IF (J .GT. I) GO TO 10
 90      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  S7IPR FOLLOWS  ***
      END
      SUBROUTINE  V7IPR(N, IP, X)
C
C     PERMUTE X SO THAT X.OUTPUT(I) = X.INPUT(IP(I)).
C     IP IS UNCHANGED ON OUTPUT.
C
      INTEGER N
      INTEGER IP(N)
      REAL X(N)
C
      INTEGER I, J, K
      REAL T
      DO 30 I = 1, N
         J = IP(I)
         IF (J .EQ. I) GO TO 30
         IF (J .GT. 0) GO TO 10
            IP(I) = -J
            GO TO 30
 10      T = X(I)
         K = I
 20      X(K) = X(J)
         K = J
         J = IP(K)
         IP(K) = -J
         IF (J .GT. I) GO TO 20
         X(K) = T
 30      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  V7IPR FOLLOWS  ***
      END
      SUBROUTINE  V7SHF(N, K, X)
C
C  ***  SHIFT X(K),...,X(N) LEFT CIRCULARLY ONE POSITION  ***
C
      INTEGER N, K
      REAL X(N)
C
      INTEGER I, NM1
      REAL T
C
      IF (K .GE. N) GO TO 999
      NM1 = N - 1
      T = X(K)
      DO 10 I = K, NM1
 10      X(I) = X(I+1)
      X(N) = T
 999  RETURN
      END
      SUBROUTINE  V7VMP(N, X, Y, Z, K)
C
C ***  SET X(I) = Y(I) * Z(I)**K, 1 .LE. I .LE. N (FOR K = 1 OR -1)  ***
C
      INTEGER N, K
      REAL X(N), Y(N), Z(N)
      INTEGER I
C
      IF (K .GE. 0) GO TO 20
      DO 10 I = 1, N
 10      X(I) = Y(I) / Z(I)
      GO TO 999
C
 20   DO 30 I = 1, N
 30      X(I) = Y(I) * Z(I)
 999  RETURN
C  ***  LAST LINE OF  V7VMP FOLLOWS  ***
      END
      SUBROUTINE I7COPY(P, Y, X)
C
C  ***  SET Y = X, WHERE X AND Y ARE INTEGER P-VECTORS  ***
C
      INTEGER P
      INTEGER X(P), Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      Y(I) = X(I)
 999  RETURN
      END
      SUBROUTINE I7PNVR(N, X, Y)
C
C  ***  SET PERMUTATION VECTOR X TO INVERSE OF Y  ***
C
      INTEGER N
      INTEGER X(N), Y(N)
C
      INTEGER I, J
      DO 10 I = 1, N
         J = Y(I)
         X(J) = I
 10      CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF I7PNVR FOLLOWS  ***
      END
      SUBROUTINE I7SHFT(N, K, X)
C
C  ***  SHIFT X(K),...,X(N) LEFT CIRCULARLY ONE POSITION  ***
C
      INTEGER N, K
      INTEGER X(N)
C
      INTEGER I, NM1, T
C
      IF (K .GE. N) GO TO 999
      NM1 = N - 1
      T = X(K)
      DO 10 I = K, NM1
 10      X(I) = X(I+1)
      X(N) = T
 999  RETURN
      END

.