SUBROUTINE    GLG(N, P, PS, X, RHO, RHOI, RHOR, IV, LIV, LV, V,
     1                  CALCRJ, UI, UR, UF)
C
C *** GENERALIZED LINEAR REGRESSION A LA NL2SOL ***
C
C  ***  PARAMETERS  ***
C
      INTEGER N, P, PS, LIV, LV
      INTEGER IV(LIV), RHOI(*), UI(*)
      REAL X(*), 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 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 90 + P.
C LV...... LENGTH OF V, AT LEAST
C              105 + P*(2*P + 17) + 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,   RGLG
C
C  IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C  RGLG ... 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   RGLG(X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, RHO, RHOI,
     1            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   RGLG(V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS, V(R1),
     1            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    GLG FOLLOWS  ***
      END
      SUBROUTINE    GLF(N, P, PS, X, RHO, RHOI, RHOR, IV, LIV, LV, V,
     1                  CALCRJ, UI, UR, UF)
C
C *** GENERALIZED LINEAR REGRESSION, FINITE-DIFFERENCE JACOBIAN ***
C
C  ***  PARAMETERS  ***
C
      INTEGER N, P, PS, LIV, LV
      INTEGER IV(LIV), RHOI(*), UI(*)
      REAL X(*), 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 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 90 + P.
C LV...... LENGTH OF V, AT LEAST
C              105 + P*(2*P + 17) + 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    GLF 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,   RGLG, V7CPY
C
C  IVSET... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C  RGLG... 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, RN, RS1
      REAL H, H0, HLIM, NEGPT5, ONE, XK, ZERO
C
C  ***  IV AND V COMPONENTS  ***
C
      INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL,
     1        NGCALL, NGCOV, R, RDREQ, REGD, 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, RDREQ=57, REGD=67, 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))
      IF (IV(COVREQ) .EQ. 0 .AND. IV(RDREQ) .GT. 0) IV(COVREQ) = -1
      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   RGLG(X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, RHO, RHOI,
     1            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
      RN = RS1 + N - 1
C
 20   CALL   RGLG(V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS, V(R1),
     1            V(RD1), RHO, RHOI, RHOR, V, X)
      IF (IV(1)-2) 30, 50, 120
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 110 K = 1, PS
         XK = X(K)
         H = V(DLTFDJ) *   MAX( ABS(XK), ONE/V(DK))
         H0 = H
         DK = DK + 1
 80      X(K) = XK + H
         NG = NG + 1
         NF = -NG
         CALL CALCRJ(N, PS, X, NF, NEED, V(R21), V(DR1), UI, UR, UF)
         IF (NF .LT. 0) GO TO 90
              H = NEGPT5 * H
              IF ( ABS(H/H0) .GE. HLIM) GO TO 80
                   GO TO 60
 90      X(K) = XK
         IV(NGCALL) = NG
         I1 = R21
         J1K = J1K0
         J1K0 = J1K0 + 1
         DO 100 I = RS1, RN
              V(J1K) = (V(I1) - V(I)) / H
              I1 = I1 + 1
              J1K = J1K + PS
 100          CONTINUE
 110     CONTINUE
      GO TO 20
C
 120  IF (IV(REGD) .GT. 0) IV(REGD) = RD1
C
 999  RETURN
C
C  ***  LAST LINE OF    GLF FOLLOWS  ***
      END
      SUBROUTINE   RGLG(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 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 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 + 90.
C LV...... LENGTH OF V...  LV  MUST BE AT LEAST
C              105 + P*(2*P+16) + 2*N + 4*PS.
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   RGLG IS CALLED
C          WITH IV(1) = 1.
C RD...... RD(I) = HALF * (G(I)**T * H(I)**-1 * G(I)) ON OUTPUT WHEN
C          IV(RDREQ) IS 2, 3, 4, OR 5.    RGLG SETS IV(REGD) = 1 IF RD
C          IS SUCCESSFULLY COMPUTED, TO 0 IF NO ATTEMPT WAS MADE
C          TO COMPUTE IT, AND TO -1 IF H (THE FINITE-DIFFERENCE HESSIAN)
C          WAS INDEFINITE.  BEFORE CONVERGENCE, RD IS ALSO USED AS
C          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.  IF IV(RDREQ) IS AT LEAST 4 (MOD 6), I.E., IF
C MORE THAN THE SIMPLEST REGRESSION DIAGNOSTIC INFORMATION IS TO BE
C COMPUTED, THEN SOME COMPONENTS OF RHOI AND RHOR MUST CONVEY SOME EXTRA
C DETAILS, AS DESCRIBED BELOW.
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
C  ***  REGRESSION DIAGNOSTICS  ***
C
C IV(RDREQ) INDICATES WHETHER A COVARIANCE MATRIX AND REGRESSION
C DIAGNOSTIC VECTOR ARE TO BE COMPUTED.  IV(RDREQ) (MOD 6) HAS THE FORM
C IV(RDREQ) = CVR +2*RDR, WHERE CVR = 0 OR 1 AND RDR = 0, 1, OR 2,
C SO THAT
C
C      CVR = MOD(IV(RDREQ), 2)
C      RDR = MOD(IV(RDREQ)/2, 3).
C
C    CVR = 0 FOR NO COVARIANCE MATRIX
C        = 1 IF A COVARIANCE MATRIX ESTIMATE IS DESIRED
C
C    RDR = 0 FOR NO LEAVE-ONE-OUT DIAGNOSTIC INFORMATION.
C        = 1 TO HAVE ONE-STEP ESTIMATES OF F(X(I)) - F(X*) STORED IN RD,
C            WHERE X(I) MINIMIZES F (THE OBJECTIVE FUNCTION) WITH
C            COMPONENT I OF R REMOVED AND X* MINIMIZES THE FULL F.
C        = 2 FOR MORE DETAILED ONE-STEP LEAVE-ONE-OUT INFORMATION, AS
C            DICTATED BY THE IV COMPONENTS DESCRIBED BELOW.
C
C FOR RDR = 2, THE FOLLOWING COMPONENTS OF IV ARE RELEVANT...
C
C  NFIX = IV(83) = NUMBER OF TRAILING NUISANCE PARAMETERS TO TREAT
C          AS FIXED WHEN COMPUTING DIAGNOSTIC VECTORS (0 .LE. NFIX .LE.
C          P - PS, SO X(I) IS KEPT FIXED FOR P - NFIX .LT. I .LE. P).
C
C   LOO = IV(84) TELLS WHAT TO LEAVE OUT...
C       = 1 MEANS LEAVE OUT EACH COMPONENT OF R SEPARATELY, AND
C       = 2 MEANS LEAVE OUT CONTIGUOUS BLOCKS OF R COMPONENTS.
C           FOR LOO = 2, IV(85) IS THE STARTING SUBSCRIPT IN RHOI
C           OF AN ARRAY BS OF BLOCK SIZES, IV(86) IS THE STRIDE FOR BS,
C           AND IV(87) = NB IS THE NUMBER OF BLOCKS, SO THAT
C           BS(I) = RHOI(IV(85) + (I-1)*IV(86)), 1 .LE. I .LE. NB.
C           NOTE THAT IF ALL BLOCKS ARE THE SAME SIZE, THEN IT SUFFICES
C           TO SET RHOI(IV(85)) = BLOCKSIZE AND IV(86) = 0.
C           NOTE THAT LOO = 1 IS EQUIVALENT TO LOO = 2 WITH
C           RHOI(IV(85)) = 1, IV(86) = 0, IV(87) = N.
C       = 3,4 ARE SIMILAR TO LOO = 1,2, RESPECTIVELY, BUT LEAVING A
C           FRACTION OUT.  IN THIS CASE, IV(88) IS THE STARTING
C           SUBSCRIPT IN RHOR OF AN ARRAY FLO OF FRACTIONS TO LEAVE OUT,
C           AND IV(89) IS THE STRIDE FOR FLO...
C           FLO(I) = RHOR(IV(88) + (I-1)*IV(89)), 1 .LE. I .LE. NB.
C
C XNOTI = IV(90) TELLS WHAT DIAGNOSTIC INFORMATION TO STORE...
C       = 0  MEANS JUST STORE ONE-STEP ESTIMATES OF F(X(I)) - F(X*) IN
C            RD(I), 1 .LE. I .LE. NB.
C       .GT. 0 MEANS ALSO STORE ONE-STEP ESTIMATES OF X(I) ESTIMATES
C            IN RHOR, STARTING AT RHOR(XNOTI)...
C              X(I)(J) = RHOR((I-1)*(P-NFIX) + J + XNOTI-1),
C              1 .LE. I .LE. NB, 1 .LE. J .LE. P - NFIX.
C
C    SOMETIMES ONE-STEP ESTIMATES OF X(I) DO NOT EXIST, BECAUSE THE
C APPROXIMATE UPDATED HESSIAN MATRIX IS INDEFINITE.  IN SUCH CASES,
C THE CORRESPONDING RD COMPONENT IS SET TO -1, AND, IF XNOTI IS
C POSITIVE, THE SOLUTION X IS RETURNED AS X(I).  WHEN ONE-STEP ESTIMATES
C OF X(I) DO EXIST, THE CORRESPONDING COMPONENT OF RD IS POSITIVE.
C
C SUMMARY OF RHOI COMPONENTS (FOR RDR = MOD(IV(RDREQ)/2, 3) = 2)...
C
C IV(83) = NFIX
C IV(84) = LOO
C IV(85) = START IN RHOI OF BS
C IV(86) = STRIDE FOR BS
C IV(87) = NB
C IV(88) = START IN RHOR OF FLO
C IV(89) = STRIDE FOR FLO
C IV(90) = XNOTI (START IN RHOR OF X(I)).
C
C
C  ***  COVARIANCE MATRIX ESTIMATE  ***
C
C IF IV(RDREQ) INDICATES THAT A COVARIANCE MATRIX IS TO BE COMPUTED,
C THEN IV(COVREQ) = IV(15) DETERMINES THE FORM OF THE COMPUTED
C COVARIANCE MATRIX ESTIMATE AND, SIMULTANEOUSLY, THE FORM OF
C APPROXIMATE HESSIAN MATRIX USED IN COMPUTING REGRESSION DIAGNOSTIC
C INFORMATION.  IN ALL CASES, SOME APPROXIMATE FINAL HESSIAN MATRIX
C IS OBTAINED, AND ITS INVERSE IS THE COVARIANCE MATRIX ESTIMATE
C (WHICH MAY HAVE TO BE SCALED APPROPRIATELY -- THAT IS UP TO YOU).
C IF IV(COVREQ) IS AT MOST 2 IN ABSOLUTE VALUE, THEN THE FINAL
C HESSIAN APPROXIMATION IS COMPUTED BY FINITE DIFFERENCES -- GRADIENT
C DIFFERENCES IF IV(COVREQ) IS NONNEGATIVE, FUNCTION DIFFERENCES
C OTHERWISE.  IF (IV(COVREQ)) IS AT LEAST 3 IN ABSOLUTE VALUE, THEN THE
C CURRENT GAUSS-NEWTON HESSIAN APPROXIMATION IS TAKEN AS THE FINAL
C HESSIAN APPROXIMATION.  FOR SOME PROBLEMS THIS SAVES TIME AND YIELDS
C THE SAME OR NEARLY THE SAME HESSIAN APPROXIMATION AS DO FINITE
C DIFFERENCES.  FOR OTHER PROBLEMS, THE TWO KINDS OF HESSIAN
C APPROXIMATIONS MAY GIVE DECIDEDLY DIFFERENT REGRESSION DIAGNOSTICS AND
C COVARIANCE MATRIX ESTIMATES.
C
C SET IV(RDREQ) TO 2, 4 OR 6 TO REQUEST A CHOLESKY FACTOR FOR THE FINAL
C APPROXIMATE HESSIAN MATRIX.  IF IV(FDH) IS POSITIVE WHEN THIS ROUTINE
C RETURNS, THEN THE CHOLESKY FACTOR IS STORED IN V, STARTING AT
C V(IV(FDH)), IN THE SEQUENCE L(1,1), L(2,1), L(2,2), L(3,1), ... .
C
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C
C+++++++++++++++++++++++++++++  DECLARATIONS  ++++++++++++++++++++++++++
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL  D7UP5, IVSET,  G2LRD,  N3RDP,  D7TPR,  Q7ADR,  VSUM,
     1         G7LIT, ITSUM,  L7NVR,  L7ITV,  L7IVM, L7SRT,  L7SQR,
     2          L7SVX,  L7SVN,  L7TSQ, L7VML, O7PRD, V2AXY, V7CPY,
     3          V7SCL,  V7SCP
      REAL  D7TPR,  L7SVX,  L7SVN, VSUM
C
C  D7UP5...  UPDATES SCALE VECTOR D.
C  IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C  G2LRD.... COMPUTES REGRESSION DIAGNOSTIC.
C  N3RDP... PRINTS REGRESSION DIAGNOSTIC.
C  D7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS.
C  Q7ADR.... ADDS ROWS TO QR FACTORIZATION.
C  VSUM..... RETURNS SUM OF ELEMENTS OF A VECTOR.
C  G7LIT.... PERFORMS BASIC MINIMIZATION ALGORITHM.
C  ITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X.
C  L7NVR... INVERTS COMPACTLY STORED TRIANGULAR MATRIX.
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  L7TSQ... COMPUTES (L**T)*L FOR LOWER TRIANG. MATRIX L.
C  L7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C  O7PRD.... ADDS OUTER PRODUCT OF VECTORS TO A MATRIX.
C  V2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER.
C  V7CPY.... COPIES ONE VECTOR TO ANOTHER.
C  V7SCL... MULTIPLIES A VECTOR BY A SCALAR.
C  V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL JUSTG, UPDATD, ZEROG
      INTEGER G1, HN1, I, II, IV1, J, J1, JTOL1, K, L, 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 CNVCOD, COVMAT, DINIT, DTYPE, DTINIT, D0INIT, F,
     1        F0, FDH, G, H, HC, IPIVOT, IVNEED, JCN, JTOL, LMAT,
     2        MODE, NEXTIV, NEXTV, NF0, NF1, NFCALL, NFCOV, NFGCAL,
     3        NGCALL, NGCOV, PERM, QTR, RDREQ, REGD, RESTOR,
     4        RMAT, RSPTOL, STEP, TOOBIG, VNEED, XNOTI, Y
C
C  ***  IV SUBSCRIPT VALUES  ***
C
      PARAMETER (CNVCOD=55, COVMAT=26, DTYPE=16, F0=13, FDH=74, G=28,
     1           H=56, HC=71, IPIVOT=76, IVNEED=3, JCN=66, JTOL=59,
     2           LMAT=42, MODE=35, NEXTIV=46, NEXTV=47, NFCALL=6,
     3           NFCOV=52, NF0=68, NF1=69, NFGCAL=7, NGCALL=30,
     4           NGCOV=53, PERM=58, QTR=77, RESTOR=9, RMAT=78, RDREQ=57,
     5           REGD=67, STEP=40, TOOBIG=2, VNEED=4, XNOTI=90, Y=48)
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(Y) + P
         IV(RESTOR) = 0
         I = IV1 + 2
         IF (IV(TOOBIG) .EQ. 0) GO TO (120, 110, 110, 130), I
         V(F) = V(F0)
         IF (I .NE. 3) IV(1) = 2
         GO TO 40
C
C  ***  FRESH START OR RESTART -- CHECK INPUT INTEGERS  ***
C
 10   IF (ND .LT. PS) GO TO 360
      IF (PS .GT. P) GO TO 360
      IF (PS .LE. 0) GO TO 360
      IF (N .LE. 0) GO TO 360
      IF (IV1 .EQ. 14) GO TO 30
      IF (IV1 .GT. 16) GO TO 420
      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
C     *** ADJUST IV(PERM) TO MAKE ROOM FOR IV INPUT COMPONENTS
C     *** NEEDED WHEN IV(RDREQ) IS 4 OR 5...
      I = XNOTI + 1
      IF (IV(PERM) .LT. I) IV(PERM) = I
C
 20   CALL  G7LIT(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(Y) = IV(NEXTV)
      IV(G) = IV(Y) + 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
      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 = IV(Y)
      CALL  G7LIT(D, V(G1), IV, LIV, LV, P, PS, V, X, V(Y1))
      IF (IV(1) - 2) 50, 60, 380
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   IF (IV(MODE) .GT. 0) GO TO 370
      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
C  ***  COMPUTE GRADIENT INFORMATION FOR FINITE-DIFFERENCE HESSIAN  ***
C
 120  IV(1) = 2
      JUSTG = .TRUE.
      I = IV(NFCALL)
      CALL RHO(NEED1, T, N, I, X(PS1), R, RD, RHOI, RHOR, V(W))
      IF (I .LE. 0) GO TO 80
      CALL RHO(NEED2, T, N, I, X(PS1), R, RD, RHOI, RHOR, V(W))
      IF (I .LE. 0) GO TO 80
      GO TO 250
C
C  ***  PREPARE TO COMPUTE GRADIENT INFORMATION WHILE ITERATING  ***
C
 130  JUSTG = .FALSE.
      G1 = IV(G)
C
C  ***  DECIDE WHETHER TO UPDATE D BELOW  ***
C
      I = IV(DTYPE)
      UPDATD = .FALSE.
      IF (I .LE. 0) GO TO 140
         IF (I .EQ. 1 .OR. IV(MODE) .LT. 0) UPDATD = .TRUE.
C
C  ***  COMPUTE RMAT AND QTR  ***
C
 140  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 160
C
C  ***  ADJUST Y  ***
C
      Y1 = IV(Y)
      WI = W
      STEP1 = IV(STEP)
      DO 150 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))
 150     CONTINUE
C
C  ***  CHECK FOR NEGATIVE W COMPONENTS  ***
C
 160  J1 = W + N - 1
      DO 170 WI = W, J1
         IF (V(WI) .LT. ZERO) GO TO 240
 170     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 200 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 180
C           *** SEMINORMAL EQUATIONS ***
            CALL  V2AXY(PS, V(G1), RHO1, DR(1,I), V(G1))
            RHO1 = ZERO
            ZEROG = .FALSE.
            GO TO 190
 180     RHO1 =  RHO1 / T
C        *** QR ACCUMULATION ***
 190     CALL  V7SCL(PS, V(TEMP1), T, DR(1,I))
         CALL  Q7ADR(PS, V(QTR1), V(RMAT1), V(TEMP1), RHO1)
 200     CONTINUE
C
C  ***  COMPUTE G FROM RMAT AND QTR  ***
C
      TEMP2 = TEMP1 + PS
      CALL  L7VML(PS, V(TEMP1), V(RMAT1), V(QTR1))
      IF (ZEROG) GO TO 220
      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 230
         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 210 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))
 210        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))
 220     IV(QTR) = QTR1
 230  CALL  V2AXY(PS, V(G1), ONE, V(TEMP1), V(G1))
      IF (PS .GE. P) GO TO 350
      GO TO 270
C
C  ***  INDEFINITE GN HESSIAN...  ***
C
 240  IV(RMAT) = -RMAT1
      IV(HC) = RMAT1
      CALL  O7PRD(N, LH, PS, V(RMAT1), V(W), DR, DR)
C
C  ***  COMPUTE GRADIENT  ***
C
 250  G1 = IV(G)
      CALL  V7SCP(P, V(G1), ZERO)
      DO 260 I = 1, N
 260     CALL  V2AXY(PS, V(G1), R(I), DR(1,I), V(G1))
      IF (PS .GE. P) GO TO 350
C
C  ***  COMPUTE GRADIENT COMPONENTS OF NUISANCE PARAMETERS ***
C
 270  K = P - PS
      J1 = 1
      G1 = G1 + PS
      DO 280 J = 1, K
         J1 = J1 + NN
         V(G1) = VSUM(N, R(J1))
         G1 = G1 + 1
 280     CONTINUE
      IF (JUSTG) GO TO 390
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 310 II = 1, PMPS
         J1 = J1 + NN
         J = J1
         DO 290 I = 1, N
            CALL  V2AXY(PS, V(K), RD(J), DR(1,I), V(K))
            J = J + 1
 290        CONTINUE
         K = K + PS
         DO 300 I = 1, II
            J1 = J1 + NN
            V(K) = VSUM(N, RD(J1))
            K = K + 1
 300        CONTINUE
 310     CONTINUE
      IF (IV(RMAT) .LE. 0) GO TO 350
      J = IV(LMAT)
      CALL  V7CPY(PSLEN, V(J), V(HN1))
      IF ( L7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2)) .LE. ZERO) GO TO 320
      CALL  L7SRT(PS1, P, V(RMAT1), V(RMAT1), I)
      IF (I .LE. 0) GO TO 330
C
C  *** HESSIAN IS NOT POSITIVE DEFINITE ***
C
 320  CALL  L7SQR(PS, V(RMAT1), V(RMAT1))
      CALL  V7CPY(PSLEN, V(HN1), V(J))
      IV(HC) = RMAT1
      IV(RMAT) = -RMAT1
      GO TO 350
C
C  *** NUISANCE PARS LEAVE HESSIAN POS. DEF.  GET REST OF QTR ***
C
 330  J = QTR1 + PS
      G1 = IV(G) + PS
      DO 340 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
 340     CONTINUE
 350  IF (JUSTG) GO TO 390
      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
 360  IV(1) = 66
      GO TO 420
C
C  ***  COVARIANCE OR INITIAL S COMPUTATION  ***
C
 370  IV(NFCOV) = IV(NFCOV) + 1
      IV(NFCALL) = IV(NFCALL) + 1
      IV(NFGCAL) = IV(NFCALL)
      IV(1) = -1
      GO TO 999
C
C  ***  CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE  ***
C
 380  IF (IV(COVMAT) .NE. 0) GO TO 410
      IF (IV(REGD) .NE. 0) GO TO 410
C
C     ***  SEE IF CHOLESKY FACTOR OF HESSIAN IS AVAILABLE  ***
C
      K = IV(FDH)
      IF (K .LE. 0) GO TO 400
      IF (IV(RDREQ) .LE. 0) GO TO 410
C
C     ***  COMPUTE REGRESSION DIAGNOSTICS AND DEFAULT COVARIANCE IF
C          DESIRED  ***
C
      IV(MODE) = P + 1
      IV(NGCALL) = IV(NGCALL) + 1
      IV(NGCOV) = IV(NGCOV) + 1
      IV(CNVCOD) = IV(1)
      IV(NFCOV) = IV(NFCOV) + 1
      IV(NFCALL) = IV(NFCALL) + 1
      IV(NFGCAL) = IV(NFCALL)
      IV(1) = -1
      GO TO 999
C
 390  IF (IV(MODE) .LE. P) GO TO 40
C     *** SAVE RD IN W FOR POSSIBLE USE IN OTHER DIAGNOSTICS ***
      CALL  V7CPY(N, V(W), RD)
C     *** OVERWRITE RD WITH REGRESSION DIAGNOSTICS ***
      L = IV(LMAT)
      I = IV(JCN)
      STEP1 = IV(STEP)
      CALL  G2LRD(DR, IV, V(L), LH, LIV, LV, ND, N, P, PS, R, RD,
     1            RHOI, RHOR, V, V(STEP1), X, V(I))
      IV(1) = IV(CNVCOD)
      IV(CNVCOD) = 0
      IF (MOD(IV(RDREQ),2) .EQ. 0) GO TO 410
C
C        *** FINISH COVARIANCE COMPUTATION ***
C
         I = IABS(IV(H))
         IV(FDH) = 0
         CALL  L7NVR(P, V(I), V(L))
         CALL  L7TSQ(P, V(I), V(I))
         IV(COVMAT) = I
         GO TO 410
C
C  ***  COME HERE FOR INDEFINITE FINITE-DIFFERENCE HESSIAN  ***
C
 400  IV(COVMAT) = K
      IV(REGD) = K
C
C  ***  PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS  ***
C
 410  G1 = IV(G)
 420  CALL  ITSUM(D, V(G1), IV, LIV, LV, P, V, X)
      IF (IV(1) .LE. 6 .AND. IV(RDREQ) .GT. 0)
     1     CALL  N3RDP(IV, LIV, LV, N, P, RD, RHOI, RHOR, V)
C
 999  RETURN
C  ***  LAST LINE OF   RGLG FOLLOWS  ***
      END
      SUBROUTINE  F7HES(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)).
C
C  ***  IF IV(COVREQ) .GE. 0 THEN  F7HES 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 D(P), G(P), V(LV), X(P)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER GSAVE1, HES, HMI, HPI, HPM, I, K, KIND, L, M, MM1, MM1O2,
     1        PP1O2, STPI, STPM, STP0
      REAL DEL, HALF, NEGPT5, ONE, TWO, ZERO
C
C  ***  EXTERNAL SUBROUTINES  ***
C
      EXTERNAL  V7CPY
C
C  V7CPY.... COPY ONE VECTOR TO ANOTHER.
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, NEGPT5=-0.5E+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
         IV(H) = -IABS(IV(H))
         IV(FDH) = 0
         IV(KAGQT) = -1
         V(FX) = V(F)
 10   IF (M .GT. P) GO TO 999
      IF (KIND .LT. 0) GO TO 110
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  F7HES.  SET GSAVE = G, TAKE FIRST STEP  ***
         CALL  V7CPY(P, V(GSAVE1), G)
         IV(SWITCH) = IV(NFGCAL)
         GO TO 90
C
 20   DEL = V(DELTA)
      X(M) = V(XMSAVE)
      IF (IV(TOOBIG) .EQ. 0) GO TO 40
C
C     ***  HANDLE OVERSIZE V(DELTA)  ***
C
         IF (DEL*X(M) .GT. ZERO) GO TO 30
C             ***  WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT  ***
              IV(FDH) = -2
              GO TO 220
C
C        ***  TRY SHRINKING V(DELTA)  ***
 30      DEL = NEGPT5 * DEL
         GO TO 100
C
 40   HES = -IV(H)
C
C  ***  SET  G = (G - GSAVE)/DEL  ***
C
      DO 50 I = 1, P
         G(I) = (G(I) - V(GSAVE1)) / DEL
         GSAVE1 = GSAVE1 + 1
 50      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 70
C
C  ***  SET  H(I,M) = 0.5 * (H(I,M) + G(I))  FOR I = 1 TO M-1  ***
C
      MM1 = M - 1
      DO 60 I = 1, MM1
         V(K) = HALF * (V(K) + G(I))
         K = K + 1
 60      CONTINUE
C
C  ***  ADD  H(I,M) = G(I)  FOR I = M TO P  ***
C
 70   L = L + 1
      DO 80 I = M, P
         V(L) = G(I)
         L = L + I
 80      CONTINUE
C
 90   M = M + 1
      IV(MODE) = M
      IF (M .GT. P) GO TO 210
C
C  ***  CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE  ***
C
      DEL = V(DELTA0) *   MAX(ONE/D(M),  ABS(X(M)))
      IF (X(M) .LT. ZERO) DEL = -DEL
      V(XMSAVE) = X(M)
 100  X(M) = X(M) + DEL
      V(DELTA) = DEL
      IRT = 2
      GO TO 999
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY.
C
 110  STP0 = IV(W) + P - 1
      MM1 = M - 1
      MM1O2 = M*MM1/2
      IF (M .GT. 0) GO TO 120
C        ***  FIRST CALL ON  F7HES.  ***
         IV(SAVEI) = 0
         GO TO 200
C
 120  I = IV(SAVEI)
      HES = -IV(H)
      IF (I .GT. 0) GO TO 180
      IF (IV(TOOBIG) .EQ. 0) GO TO 140
C
C     ***  HANDLE OVERSIZE STEP  ***
C
         STPM = STP0 + M
         DEL = V(STPM)
         IF (DEL*X(XMSAVE) .GT. ZERO) GO TO 130
C             ***  WE ALREADY TRIED SHRINKING THE STEP, SO QUIT  ***
              IV(FDH) = -2
              GO TO 220
C
C        ***  TRY SHRINKING THE STEP  ***
 130     DEL = NEGPT5 * DEL
         X(M) = X(XMSAVE) + DEL
         V(STPM) = DEL
         IRT = 1
         GO TO 999
C
C  ***  SAVE F(X + STP(M)*E(M)) IN H(P,M)  ***
C
 140  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
      HMI = HES + MM1O2
      IF (MM1 .EQ. 0) GO TO 160
      HPI = HES + PP1O2
      DO 150 I = 1, MM1
         V(HMI) = V(FX) - (V(F) + V(HPI))
         HMI = HMI + 1
         HPI = HPI + 1
 150     CONTINUE
 160  V(HMI) = V(F) - TWO*V(FX)
C
C  ***  COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H.  ***
C
      I = 1
C
 170  IV(SAVEI) = I
      STPI = STP0 + I
      V(DELTA) = X(I)
      X(I) = X(I) + V(STPI)
      IF (I .EQ. M) X(I) = V(XMSAVE) - V(STPI)
      IRT = 1
      GO TO 999
C
 180  X(I) = V(DELTA)
      IF (IV(TOOBIG) .EQ. 0) GO TO 190
C        ***  PUNT IN THE EVENT OF AN OVERSIZE STEP  ***
         IV(FDH) = -2
         GO TO 220
C
C  ***  FINISH COMPUTING H(M,I)  ***
C
 190  STPI = STP0 + I
      HMI = HES + MM1O2 + I - 1
      STPM = STP0 + M
      V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM))
      I = I + 1
      IF (I .LE. M) GO TO 170
      IV(SAVEI) = 0
      X(M) = V(XMSAVE)
C
 200  M = M + 1
      IV(MODE) = M
      IF (M .GT. P) GO TO 210
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
      DEL = V(DLTFDC) *   MAX(ONE/D(M),  ABS(X(M)))
      IF (X(M) .LT. ZERO) DEL = -DEL
      V(XMSAVE) = X(M)
      X(M) = X(M) + DEL
      STPM = STP0 + M
      V(STPM) = DEL
      IRT = 1
      GO TO 999
C
C  ***  RESTORE V(F), ETC.  ***
C
 210  IV(FDH) = HES
 220  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  F7HES FOLLOWS  ***
      END
      SUBROUTINE  G2LRD(DR, IV, L, LH, LIV, LV, ND, N, P, PS, R, RD,
     1                  RHOI, RHOR, V, W, X, Z)
C
C  ***  COMPUTE REGRESSION DIAGNOSTIC FOR   RGLG  ***
C
C  ***  PARAMETERS  ***
C
      INTEGER LH, LIV, LV, ND, N, P, PS
      INTEGER IV(LIV), RHOI(*)
      REAL DR(ND,P), L(LH), R(N), RD(N), RHOR(*), V(LV),
     1                 W(P), X(P), Z(P)
C
C  ***  CODED BY DAVID M. GAY (SPRING 1986, SUMMER 1991)  ***
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      EXTERNAL  D7TPR,  L7ITV,  L7IVM, L7SRT,  L7SQR,  S7LVM,
     1         V2AXY, V7CPY,  V7SCP
      REAL  D7TPR
C
C  D7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS.
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  S7LVM... MULTIPLIES COMPACTLY STORED SYM. MATRIX TIMES VECTOR.
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
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL USEFLO
      INTEGER BS1, BSINC, FLO1, FLOINC, H1, HPS1, I,
     1        J, J1, K, KI, KI1, KID, L1, LE, LL, LOO1, N1,
     2        PMPS, PP1O2, PS1, PX, RDR, XNI, ZAP1, ZAPLEN
      REAL FRAC, HI, RI, S, T, T1
C
C  ***  CONSTANTS  ***
C
      REAL HALF, NEGONE, ONE, ZERO
C
C
C  ***  IV SUBSCRIPTS  ***
C
      INTEGER BS, BSSTR, COVREQ, FDH, FLO, FLOSTR, LOO, NB, NFIX,
     1        RDREQ, REGD, XNOTI
      PARAMETER (BS=85, BSSTR=86, COVREQ=15, FDH=74, FLO=88, FLOSTR=89,
     1           LOO=84, NB=87, NFIX=83, RDREQ=57, REGD=67, XNOTI=90)
      PARAMETER (HALF=0.5E+0, NEGONE=-1.E+0, ONE=1.E+0, ZERO=0.E+0)
C
C++++++++++++++++++++++++++++++++  BODY  +++++++++++++++++++++++++++++++
C
      I = IV(RDREQ)
      RDR = MOD(I/2, 3)
      IF (RDR .EQ. 0) GO TO 999
      H1 = IV(FDH)
      USEFLO = .FALSE.
      PX = P
      N1 = N
      FRAC = ONE
      XNI = 0
      IF (RDR .EQ. 1) GO TO 120
      LOO1 = IV(LOO)
      IF (LOO1 .LE. 0 .OR. LOO1 .GT. 6) THEN
         IV(REGD) = -1
         GO TO 999
         ENDIF
      IF (LOO1 .GT. 3) THEN
         USEFLO = .TRUE.
         FLO1 = IV(FLO)
         FLOINC = IV(FLOSTR)
         LOO1 = LOO1 - 3
         ENDIF
      XNI = IV(XNOTI)
      PX = P - IV(NFIX)
      IF (PX .LT. PS .OR. PX .GT. P) THEN
         IV(REGD) = -2
         GO TO 999
         ENDIF
      IF (LOO1 .EQ. 1) GO TO 120
      N1 = IV(NB)
      IF (N1 .LE. 0 .OR. N1 .GT. N) THEN
         IV(REGD) = -3
         GO TO 999
         ENDIF
      BS1 = IV(BS)
      BSINC = IV(BSSTR)
      IF (H1 .LE. 0) GO TO 190
      IF (IABS(IV(COVREQ)) .GE. 3) CALL  L7SQR(P, V(H1), L)
      PP1O2 = PX*(PX+1)/2
      PS1 = PS + 1
      ZAP1 = PS*(PS1)/2 + 1
      LE = 0
      DO 100 I = 1, N1
         IF (USEFLO) THEN
            FRAC = RHOR(FLO1)
            FLO1 = FLO1 + FLOINC
            ENDIF
         L1 = LE + 1
         IF (L1 .GT. N) GO TO 110
         LE = LE + RHOI(BS1)
         IF (LE .GT. N) LE = N
         BS1 = BS1 + BSINC
         CALL  V7CPY(PP1O2, L, V(H1))
         IF (PS .GE. PX) GO TO 50
            K = ZAP1
            KI = L1
            DO 40 J = PS1, P
               KI = KI + N
               KI1 = KI
               DO 10 LL = L1, LE
                  CALL  V2AXY(PS, L(K), -FRAC*RD(KI1), DR(1,LL), L(K))
                  KI1 = KI1 + 1
 10               CONTINUE
               K = K + PS
               DO 30 J1 = PS1, J
                  KI = KI + N
                  KI1 = KI
                  T = ZERO
                  DO 20 LL = L1, LE
                     T = T + RD(KI1)
                     KI1 = KI1 + 1
 20                  CONTINUE
                  L(K) = L(K) - FRAC*T
                  K = K + 1
 30               CONTINUE
 40            CONTINUE
 50      DO 70 LL = L1, LE
            T = -FRAC*RD(LL)
            K = 1
            DO 60 J = 1, PS
               CALL  V2AXY(J, L(K), T*DR(J,LL), DR(1,LL), L(K))
               K = K + J
 60            CONTINUE
 70         CONTINUE
         CALL  L7SRT(1, PX, L, L, J)
         IF (J .EQ. 0) THEN
            CALL  V7SCP(PX, W, ZERO)
            DO 90 LL = L1, LE
               CALL  V2AXY(PS, W, R(LL), DR(1,LL), W)
               IF (PS1 .GT. PX) GO TO 90
               K = L1
               DO 80 J = PS1, P
                  K = K + N
                  W(J) = W(J) + R(K)
 80               CONTINUE
 90            CONTINUE
            CALL  L7IVM(PX, W, L, W)
            CALL  L7ITV(PX, W, L, W)
            CALL  S7LVM(PX, Z, V(H1), W)
            RD(I) = HALF * FRAC *  D7TPR(PX, W, Z)
            IF (XNI .GT. 0) THEN
               CALL  V2AXY(PX, RHOR(XNI), FRAC, W, X)
               XNI = XNI + PX
               ENDIF
         ELSE
            RD(I) = NEGONE
            IF (XNI .GT. 0) THEN
               CALL  V7CPY(PX, RHOR(XNI), X)
               XNI = XNI + PX
               ENDIF
            ENDIF
 100     CONTINUE
 110  IV(REGD) = 1
C     *** RESTORE L ***
      CALL  L7SRT(1, P, L, V(H1), J)
      GO TO 999
C
 120  IF (H1 .LE. 0) GO TO 190
      IF (IABS(IV(COVREQ)) .GE. 3) CALL  L7SQR(P, V(H1), L)
      IF (PS .GE. PX) GO TO 170
      PS1 = PS + 1
      PMPS = PX - PS
      ZAP1 = PS*(PS1)/2
      ZAPLEN = PX*(PX+1)/2 - ZAP1
      HPS1 = H1 + ZAP1
      ZAP1 = ZAP1 + 1
      DO 160 I = 1, N
         IF (USEFLO) THEN
            FRAC = RHOR(FLO1)
            FLO1 = FLO1 + FLOINC
            ENDIF
         CALL  V7CPY(ZAPLEN, L(ZAP1), V(HPS1))
         CALL  V7SCP(PS, W, ZERO)
         K = ZAP1
         KI = I
         KID = KI
         DO 140 J = PS1, PX
            KI = KI + N
            CALL  V2AXY(PS, L(K), -FRAC*RD(KI), DR(1,I), L(K))
            K = K + PS
            KID = KID + N
            W(J) = FRAC*R(KID)
            DO 130 J1 = PS1, J
               KI = KI + N
               L(K) =  L(K) - FRAC*RD(KI)
               K = K + 1
 130           CONTINUE
 140        CONTINUE
         CALL  L7SRT(PS1, PX, L, L, J)
         IF (J .NE. 0) GO TO 150
         CALL  V7CPY(PS, Z, DR(1,I))
         CALL  V7SCP(PMPS, Z(PS1), ZERO)
         CALL  L7IVM(PX, Z, L, Z)
         HI =  D7TPR(PX, Z, Z)
         CALL  L7IVM(PX, W, L, W)
         RI = FRAC*R(I)
C        *** FIRST PS ELEMENTS OF W VANISH ***
         T =  D7TPR(PMPS, W(PS1), Z(PS1))
         S = FRAC*RD(I)
         T1 = ONE - S*HI
         IF (T1 .LE. ZERO) GO TO 150
         CALL  V2AXY(PX, W, (RI + S*T)/T1, Z, W)
         CALL  L7ITV(PX, W, L, W)
         CALL  S7LVM(PX, Z, V(H1), W)
         RD(I) = HALF *  D7TPR(PX, W, Z)
         IF (XNI .GT. 0) THEN
            CALL  V2AXY(PX, RHOR(XNI), ONE, W, X)
            XNI = XNI + PX
            ENDIF
         GO TO 160
 150     RD(I) = NEGONE
         IF (XNI .GT. 0) THEN
            CALL  V7CPY(PX, RHOR(XNI), X)
            XNI = XNI + PX
            ENDIF
 160     CONTINUE
C
C     *** RESTORE L ***
C
      CALL  V7CPY(ZAPLEN, L(ZAP1), V(HPS1))
      CALL  L7SRT(PS1, PX, L, L, J)
      GO TO 200
C
 170  DO 180 I = 1, N
         IF (USEFLO) THEN
            FRAC = RHOR(FLO1)
            FLO1 = FLO1 + FLOINC
            ENDIF
         CALL  L7IVM(PX, Z, L, DR(1,I))
         S =  D7TPR(PX, Z, Z)
         T = ONE - FRAC*RD(I) * S
         IF (T .LE. ZERO) THEN
            RD(I) = NEGONE
            IF (XNI .GT. 0) THEN
               CALL  V7CPY(PX, RHOR(XNI), X)
               XNI = XNI + PX
               ENDIF
         ELSE
            RD(I) = HALF * FRAC * (R(I)/T)**2 * S
            IF (XNI .GT. 0) THEN
               CALL  L7ITV(PX, Z, L, Z)
               CALL  V2AXY(PX, RHOR(XNI), FRAC*R(I)/T, Z, X)
               XNI = XNI + PX
               ENDIF
            ENDIF
 180     CONTINUE
      GO TO 200
C
 190  CALL  V7SCP(N1, RD, NEGONE)
 200  IV(REGD) = 1
C
 999  RETURN
C  ***  LAST LINE OF  G2LRD FOLLOWS  ***
      END
      SUBROUTINE  G7LIT(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
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LIV, LV, P, PS
      INTEGER IV(LIV)
      REAL D(P), G(P), V(LV), X(P), Y(P)
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C D.... SCALE VECTOR.
C IV... INTEGER VALUE ARRAY.
C LIV.. LENGTH OF IV.  MUST BE AT LEAST 82.
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 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        G7LIT 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 IN V, STARTING AT IV(HC)), AND  G7LIT BUILDS AN
C     APPROXIMATION, S, TO THE SECOND-ORDER TERM.  THE CALLER ALSO
C     PROVIDES THE FUNCTION VALUE, GRADIENT, AND PART OF THE YIELD
C     VECTOR USED IN UPDATING S.  G7LIT DECIDES DYNAMICALLY WHETHER OR
C     NOT TO USE S WHEN CHOOSING THE NEXT STEP TO TRY...  THE HESSIAN
C     APPROXIMATION USED IS EITHER HC ALONE (GAUSS-NEWTON MODEL) OR
C     HC + S (AUGMENTED MODEL).
C
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, G7LIT 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  G7LIT 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 NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER
C     (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS
C     NOT NEEDED).  MOREOVER, COMPARED WITH NL2SOL, IV(1) MAY HAVE THE
C     TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW,
C     AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL).  THE VALUES IV(D),
C     IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM NL2SOL (AND
C     NL2SNO), ARE NOT REFERENCED BY  G7LIT OR THE SUBROUTINES IT CALLS.
C
C        WHEN  G7LIT IS FIRST CALLED, I.E., WHEN  G7LIT IS CALLED WITH
C     IV(1) = 0 OR 12, V(F), G, AND HC NEED NOT BE INITIALIZED.  TO
C     OBTAIN THESE STARTING VALUES, G7LIT 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  G7LIT WILL MAKE A
C     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  G7LIT 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  G7LIT 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  G7LIT 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  G7LIT 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(TOOBIG) TO 1, IN WHICH CASE  G7LIT WILL RETURN WITH
C             IV(1) = 15.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED IN PART BY D.O.E. GRANT EX-76-A-01-2295 TO MIT/CCREMS.
C
C        (SEE NL2SOL FOR REFERENCES.)
C
C+++++++++++++++++++++++++++  DECLARATIONS  ++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER DUMMY, DIG1, G01, H1, HC1, I, IPIV1, J, K, L, LMAT1,
     1        LSTGST, PP1O2, QTR1, RMAT1, RSTRST, STEP1, STPMOD, S1,
     2        TEMP1, TEMP2, W1, X01
      REAL E, STTSST, T, T1
C
C     ***  CONSTANTS  ***
C
      REAL HALF, NEGONE, ONE, ONEP2, ZERO
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      LOGICAL STOPX
      REAL  D7TPR,  L7SVX,  L7SVN,  RLDST,  R7MDC,  V2NRM
      EXTERNAL  A7SST,  D7TPR, F7HES, G7QTS, ITSUM,  L7MST, L7SRT,
     1          L7SQR,  L7SVX,  L7SVN,  L7TVM, L7VML, PARCK,  RLDST,
     2          R7MDC,  S7LUP,  S7LVM, STOPX, V2AXY, V7CPY,  V7SCP,
     3          V2NRM
C
C  A7SST.... ASSESSES CANDIDATE STEP.
C  D7TPR... RETURNS INNER PRODUCT OF TWO VECTORS.
C  F7HES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE).
C  G7QTS.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL).
C  ITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X.
C  L7MST... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL).
C  L7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX.
C  L7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L.
C  L7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C  L7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX.
C  L7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX.
C  L7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
C  PARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS.
C  RLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE.
C  R7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS.
C  S7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI-
C             ANGLE OF A SYMMETRIC MATRIX.
C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED.
C  V2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER.
C  V7CPY.... COPIES ONE VECTOR TO ANOTHER.
C  V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR.
C  V2NRM... RETURNS THE 2-NORM OF A VECTOR.
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER CNVCOD, COSMIN, COVMAT, COVREQ, DGNORM, DIG, DSTNRM, F,
     1        FDH, FDIF, FUZZ, F0, GTSTEP, H, HC, IERR, INCFAC, INITS,
     2        IPIVOT, IRC, KAGQT, KALM, LMAT, LMAX0, LMAXS, MODE, MODEL,
     3        MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL, NFCOV, NGCOV,
     4        NGCALL, NITER, NVSAVE, PHMXFC, PREDUC, QTR, RADFAC,
     5        RADINC, RADIUS, RAD0, RCOND, RDREQ, REGD, RELDX, RESTOR,
     6        RMAT, S, SIZE, STEP, STGLIM, STLSTG, STPPAR, SUSED,
     7        SWITCH, TOOBIG, TUNER4, TUNER5, VNEED, VSAVE, W, WSCALE,
     8        XIRC, X0
C
C  ***  IV SUBSCRIPT VALUES  ***
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, KAGQT=33,
     2           KALM=34, LMAT=42, MODE=35, MODEL=5, MXFCAL=17,
     3           MXITER=18, NEXTV=47, NFCALL=6, NFGCAL=7, NFCOV=52,
     4           NGCOV=53, NGCALL=30, NITER=31, QTR=77, RADINC=8,
     5           RDREQ=57, REGD=67, RESTOR=9, RMAT=78, S=62, STEP=40,
     6           STGLIM=11, STLSTG=41, SUSED=64, SWITCH=12, TOOBIG=2,
     7           VNEED=4, VSAVE=60, W=65, 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, RCOND=53, RELDX=17, SIZE=55, STPPAR=5,
     4           TUNER4=29, 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 40
      IF (I .EQ. 2) GO TO 50
C
      IF (I .EQ. 12 .OR. I .EQ. 13)
     1     IV(VNEED) = IV(VNEED) + P*(3*P + 19)/2 + 7
      CALL  PARCK(1, D, IV, LIV, LV, P, V)
      I = IV(1) - 2
      IF (I .GT. 12) GO TO 999
      GO TO (290, 290, 290, 290, 290, 290, 170, 120, 170, 10, 10, 20), I
C
C  ***  STORAGE ALLOCATION  ***
C
 10   PP1O2 = P * (P + 1) / 2
      IV(S) = IV(LMAT) + PP1O2
      IV(X0) = IV(S) + PP1O2
      IV(STEP) = IV(X0) + P
      IV(STLSTG) = IV(STEP) + P
      IV(DIG) = IV(STLSTG) + P
      IV(W) = IV(DIG) + P
      IV(H) = IV(W) + 4*P + 7
      IV(NEXTV) = IV(H) + PP1O2
      IF (IV(1) .NE. 13) GO TO 20
         IV(1) = 14
         GO TO 999
C
C  ***  INITIALIZATION  ***
C
 20   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(RESTOR) = 0
      IV(FDH) = 0
      V(RAD0) = ZERO
      V(STPPAR) = ZERO
      V(RADIUS) = V(LMAX0) / (ONE + V(PHMXFC))
C
C  ***  SET INITIAL MODEL AND S MATRIX  ***
C
      IV(MODEL) = 1
      IF (IV(S) .LT. 0) GO TO 999
      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)
      IV(1) = 1
      J = IV(IPIVOT)
      IF (J .LE. 0) GO TO 999
      DO 30 I = 1, P
         IV(J) = I
         J = J + 1
 30      CONTINUE
      GO TO 999
C
C  ***  NEW FUNCTION VALUE  ***
C
 40   IF (IV(MODE) .EQ. 0) GO TO 290
      IF (IV(MODE) .GT. 0) GO TO 520
C
      IV(1) = 2
      IF (IV(TOOBIG) .EQ. 0) GO TO 999
         IV(1) = 63
         GO TO 999
C
C  ***  NEW GRADIENT  ***
C
 50   IV(KALM) = -1
      IV(KAGQT) = -1
      IV(FDH) = 0
      IF (IV(MODE) .GT. 0) GO TO 520
C
C  ***  MAKE SURE GRADIENT COULD BE COMPUTED  ***
C
      IF (IV(TOOBIG) .EQ. 0) GO TO 60
         IV(1) = 65
         GO TO 999
 60   IF (IV(HC) .LE. 0 .AND. IV(RMAT) .LE. 0) GO TO 610
C
C  ***  COMPUTE  D**-1 * GRADIENT  ***
C
      DIG1 = IV(DIG)
      K = DIG1
      DO 70 I = 1, P
         V(K) = G(I) / D(I)
         K = K + 1
 70      CONTINUE
      V(DGNORM) =  V2NRM(P, V(DIG1))
C
      IF (IV(CNVCOD) .NE. 0) GO TO 510
      IF (IV(MODE) .EQ. 0) GO TO 440
      IV(MODE) = 0
      V(F0) = V(F)
      IF (IV(INITS) .LE. 2) GO TO 100
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 530
C
C  ***  COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S  ***
C
 80   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 90
         CALL  V2AXY(PP1O2, V(S1), NEGONE, V(HC1), V(H1))
         GO TO 100
 90   RMAT1 = IV(RMAT)
      CALL  L7SQR(PS, V(S1), V(RMAT1))
      CALL  V2AXY(PP1O2, V(S1), NEGONE, V(S1), V(H1))
 100  IV(1) = 2
C
C
C-----------------------------  MAIN LOOP  -----------------------------
C
C
C  ***  PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT  ***
C
 110  CALL  ITSUM(D, G, IV, LIV, LV, P, V, X)
 120  K = IV(NITER)
      IF (K .LT. IV(MXITER)) GO TO 130
         IV(1) = 10
         GO TO 999
 130  IV(NITER) = K + 1
C
C  ***  UPDATE RADIUS  ***
C
      IF (K .EQ. 0) GO TO 150
      STEP1 = IV(STEP)
      DO 140 I = 1, P
         V(STEP1) = D(I) * V(STEP1)
         STEP1 = STEP1 + 1
 140     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
 150  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
 160  IF (.NOT. STOPX(DUMMY)) GO TO 180
         IV(1) = 11
         GO TO 190
C
C     ***  COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX.
C
 170  IF (V(F) .GE. V(F0)) GO TO 180
         V(RADFAC) = ONE
         K = IV(NITER)
         GO TO 130
C
 180  IF (IV(NFCALL) .LT. IV(MXFCAL) + IV(NFCOV)) GO TO 200
         IV(1) = 9
 190     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 430
C
C. . . . . . . . . . . . .  COMPUTE CANDIDATE STEP  . . . . . . . . . .
C
 200  STEP1 = IV(STEP)
      W1 = IV(W)
      H1 = IV(H)
      T1 = ONE
      IF (IV(MODEL) .EQ. 2) GO TO 210
         T1 = ZERO
C
C        ***  COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE...
C
         RMAT1 = IV(RMAT)
         IF (RMAT1 .LE. 0) GO TO 210
         QTR1 = IV(QTR)
         IF (QTR1 .LE. 0) GO TO 210
         IPIV1 = IV(IPIVOT)
         CALL  L7MST(D, G, IV(IERR), IV(IPIV1), IV(KALM), P, V(QTR1),
     1               V(RMAT1), V(STEP1), V, V(W1))
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 260
C
 210  IF (H1 .GT. 0) GO TO 250
C
C     ***  SET H TO  D**-1 * (HC + T1*S) * D**-1.  ***
C
         H1 = -H1
         IV(H) = H1
         IV(FDH) = 0
         J = IV(HC)
         IF (J .GT. 0) GO TO 220
            J = H1
            RMAT1 = IV(RMAT)
            CALL  L7SQR(P, V(H1), V(RMAT1))
 220     S1 = IV(S)
         DO 240 I = 1, P
              T = ONE / D(I)
              DO 230 K = 1, I
                   V(H1) = T * (V(J) + T1*V(S1)) / D(K)
                   J = J + 1
                   H1 = H1 + 1
                   S1 = S1 + 1
 230               CONTINUE
 240          CONTINUE
         H1 = IV(H)
         IV(KAGQT) = -1
C
C  ***  COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP  ***
C
 250  DIG1 = IV(DIG)
      LMAT1 = IV(LMAT)
      CALL  G7QTS(D, V(DIG1), V(H1), IV(KAGQT), V(LMAT1), P, V(STEP1),
     1            V, V(W1))
      IF (IV(KALM) .GT. 0) IV(KALM) = 0
C
 260  IF (IV(IRC) .NE. 6) GO TO 270
         IF (IV(RESTOR) .NE. 2) GO TO 290
         RSTRST = 2
         GO TO 300
C
C  ***  CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE  ***
C
 270  IV(TOOBIG) = 0
      IF (V(DSTNRM) .LE. ZERO) GO TO 290
      IF (IV(IRC) .NE. 5) GO TO 280
      IF (V(RADFAC) .LE. ONE) GO TO 280
      IF (V(PREDUC) .GT. ONEP2 * V(FDIF)) GO TO 280
         IF (IV(RESTOR) .NE. 2) GO TO 290
         RSTRST = 0
         GO TO 300
C
C  ***  COMPUTE F(X0 + STEP)  ***
C
 280  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 999
C
C. . . . . . . . . . . . .  ASSESS CANDIDATE STEP  . . . . . . . . . . .
C
 290  RSTRST = 3
 300  X01 = IV(X0)
      V(RELDX) =  RLDST(P, D, X, V(X01))
      CALL  A7SST(IV, LIV, LV, V)
      STEP1 = IV(STEP)
      LSTGST = IV(STLSTG)
      I = IV(RESTOR) + 1
      GO TO (340, 310, 320, 330), I
 310  CALL  V7CPY(P, X, V(X01))
      GO TO 340
 320   CALL  V7CPY(P, V(LSTGST), V(STEP1))
       GO TO 340
 330     CALL  V7CPY(P, V(STEP1), V(LSTGST))
         CALL  V2AXY(P, X, ONE, V(STEP1), V(X01))
         V(RELDX) =  RLDST(P, D, X, V(X01))
         IV(RESTOR) = RSTRST
C
C  ***  IF NECESSARY, SWITCH MODELS  ***
C
 340  IF (IV(SWITCH) .EQ. 0) GO TO 350
         IV(H) = -IABS(IV(H))
         IV(SUSED) = IV(SUSED) + 2
         L = IV(VSAVE)
         CALL  V7CPY(NVSAVE, V, V(L))
 350  L = IV(IRC) - 4
      STPMOD = IV(MODEL)
      IF (L .GT. 0) GO TO (370,380,390,390,390,390,390,390,500,440), 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 360
C
C     ***  SWITCH MODELS  ***
C
         IV(MODEL) = 3 - IV(MODEL)
         IF (-2 .LT. L) GO TO 400
              IV(H) = -IABS(IV(H))
              IV(SUSED) = IV(SUSED) + 2
              L = IV(VSAVE)
              CALL  V7CPY(NVSAVE, V(L), V)
              GO TO 160
C
 360  IF (-3 .LT. L) GO TO 400
C
C  ***  RECOMPUTE STEP WITH NEW RADIUS  ***
C
 370  V(RADIUS) = V(RADFAC) * V(DSTNRM)
      GO TO 160
C
C  ***  COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST
C
 380  V(RADIUS) = V(LMAXS)
      GO TO 200
C
C  ***  CONVERGENCE OR FALSE CONVERGENCE  ***
C
 390  IV(CNVCOD) = L
      IF (V(F) .GE. V(F0)) GO TO 510
         IF (IV(XIRC) .EQ. 14) GO TO 510
              IV(XIRC) = 14
C
C. . . . . . . . . . . .  PROCESS ACCEPTABLE STEP  . . . . . . . . . . .
C
 400  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 430
         STEP1 = IV(STEP)
         TEMP1 = IV(STLSTG)
         TEMP2 = IV(W)
C
C     ***  SET  TEMP1 = HESSIAN * STEP  FOR USE IN GRADIENT TESTS  ***
C
         HC1 = IV(HC)
         IF (HC1 .LE. 0) GO TO 410
              CALL  S7LVM(P, V(TEMP1), V(HC1), V(STEP1))
              GO TO 420
 410     RMAT1 = IV(RMAT)
         CALL  L7TVM(P, V(TEMP1), V(RMAT1), V(STEP1))
         CALL  L7VML(P, V(TEMP1), V(RMAT1), V(TEMP1))
C
 420     IF (STPMOD .EQ. 1) GO TO 430
              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
 430  IV(NGCALL) = IV(NGCALL) + 1
      G01 = IV(W)
      CALL  V7CPY(P, V(G01), G)
      IV(1) = 2
      IV(TOOBIG) = 0
      GO TO 999
C
C  ***  INITIALIZATIONS -- G0 = G - G0, ETC.  ***
C
 440  G01 = IV(W)
      CALL  V2AXY(P, V(G01), NEGONE, V(G01), G)
      STEP1 = IV(STEP)
      TEMP1 = IV(STLSTG)
      TEMP2 = IV(W)
      IF (IV(IRC) .NE. 3) GO TO 470
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 450 I = 1, P
              V(K) = (V(K) - V(L)) / D(I)
              K = K + 1
              L = L + 1
 450          CONTINUE
C
C        ***  DO GRADIENT TESTS  ***
C
         IF ( V2NRM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4))  GO TO 460
              IF ( D7TPR(P, G, V(STEP1))
     1                  .GE. V(GTSTEP) * V(TUNER5))  GO TO 470
 460               V(RADFAC) = V(INCFAC)
C
C  ***  COMPUTE Y VECTOR NEEDED FOR UPDATING S  ***
C
 470  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 480
         CALL  S7LVM(PS, V(G01), V(HC1), V(STEP1))
         GO TO 490
C
 480  RMAT1 = IV(RMAT)
      CALL  L7TVM(PS, V(G01), V(RMAT1), V(STEP1))
      CALL  L7VML(PS, V(G01), V(RMAT1), V(G01))
C
 490  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 110
C
C. . . . . . . . . . . . . .  MISC. DETAILS  . . . . . . . . . . . . . .
C
C  ***  BAD PARAMETERS TO ASSESS  ***
C
 500  IV(1) = 64
      GO TO 999
C
C
C  ***  CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE  ***
C
 510  IF (IV(RDREQ) .EQ. 0) GO TO 600
      IF (IV(FDH) .NE. 0) GO TO 600
      IF (IV(CNVCOD) .GE. 7) GO TO 600
      IF (IV(REGD) .GT. 0) GO TO 600
      IF (IV(COVMAT) .GT. 0) GO TO 600
      IF (IABS(IV(COVREQ)) .GE. 3) GO TO 560
      IF (IV(RESTOR) .EQ. 0) IV(RESTOR) = 2
      GO TO 530
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE  ***
C
 520  IV(RESTOR) = 0
 530  CALL  F7HES(D, G, I, IV, LIV, LV, P, V, X)
      GO TO (540, 550, 580), I
 540  IV(NFCOV) = IV(NFCOV) + 1
      IV(NFCALL) = IV(NFCALL) + 1
      IV(1) = 1
      GO TO 999
C
 550  IV(NGCOV) = IV(NGCOV) + 1
      IV(NGCALL) = IV(NGCALL) + 1
      IV(NFGCAL) = IV(NFCALL) + IV(NGCOV)
      IV(1) = 2
      GO TO 999
C
 560  H1 = IABS(IV(H))
      IV(H) = -H1
      PP1O2 = P * (P + 1) / 2
      RMAT1 = IV(RMAT)
      IF (RMAT1 .LE. 0) GO TO 570
           LMAT1 = IV(LMAT)
           CALL  V7CPY(PP1O2, V(LMAT1), V(RMAT1))
           V(RCOND) = ZERO
           GO TO 590
 570  HC1 = IV(HC)
      IV(FDH) = H1
      CALL  V7CPY(P*(P+1)/2, V(H1), V(HC1))
C
C  ***  COMPUTE CHOLESKY FACTOR OF FINITE-DIFFERENCE HESSIAN
C  ***  FOR USE IN CALLER*S COVARIANCE CALCULATION...
C
 580  LMAT1 = IV(LMAT)
      H1 = IV(FDH)
      IF (H1 .LE. 0) GO TO 600
      IF (IV(CNVCOD) .EQ. 70) GO TO 80
      CALL  L7SRT(1, P, V(LMAT1), V(H1), I)
      IV(FDH) = -1
      V(RCOND) = ZERO
      IF (I .NE. 0) GO TO 600
C
 590  IV(FDH) = -1
      STEP1 = IV(STEP)
      T =  L7SVN(P, V(LMAT1), V(STEP1), V(STEP1))
      IF (T .LE. ZERO) GO TO 600
      T = T /  L7SVX(P, V(LMAT1), V(STEP1), V(STEP1))
      IF (T .GT.  R7MDC(4)) IV(FDH) = H1
      V(RCOND) = T
C
 600  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
 610  IV(1) = 1400
C
 999  RETURN
C
C  ***  LAST LINE OF  G7LIT FOLLOWS  ***
      END
      SUBROUTINE  L7NVR(N, LIN, L)
C
C  ***  COMPUTE  LIN = L**-1,  BOTH  N X N  LOWER TRIANG. STORED   ***
C  ***  COMPACTLY BY ROWS.  LIN AND L MAY SHARE THE SAME STORAGE.  ***
C
C  ***  PARAMETERS  ***
C
      INTEGER N
      REAL L(1), LIN(1)
C     DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IM1, JJ, J0, J1, K, K0, NP1
      REAL ONE, T, ZERO
      PARAMETER (ONE=1.E+0, ZERO=0.E+0)
C
C  ***  BODY  ***
C
      NP1 = N + 1
      J0 = N*(NP1)/2
      DO 30 II = 1, N
         I = NP1 - II
         LIN(J0) = ONE/L(J0)
         IF (I .LE. 1) GO TO 999
         J1 = J0
         IM1 = I - 1
         DO 20 JJ = 1, IM1
              T = ZERO
              J0 = J1
              K0 = J1 - JJ
              DO 10 K = 1, JJ
                   T = T - L(K0)*LIN(J0)
                   J0 = J0 - 1
                   K0 = K0 + K - I
 10                CONTINUE
              LIN(J0) = T/L(K0)
 20           CONTINUE
         J0 = J0 - 1
 30      CONTINUE
 999  RETURN
C  ***  LAST LINE OF  L7NVR FOLLOWS  ***
      END
      SUBROUTINE  L7TSQ(N, A, L)
C
C  ***  SET A TO LOWER TRIANGLE OF (L**T) * L  ***
C
C  ***  L = N X N LOWER TRIANG. MATRIX STORED ROWWISE.  ***
C  ***  A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L.  ***
C
      INTEGER N
      REAL A(1), L(1)
C     DIMENSION A(N*(N+1)/2), L(N*(N+1)/2)
C
      INTEGER I, II, IIM1, I1, J, K, M
      REAL LII, LJ
C
      II = 0
      DO 50 I = 1, N
         I1 = II + 1
         II = II + I
         M = 1
         IF (I .EQ. 1) GO TO 30
         IIM1 = II - 1
         DO 20 J = I1, IIM1
              LJ = L(J)
              DO 10 K = I1, J
                   A(M) = A(M) + LJ*L(K)
                   M = M + 1
 10                CONTINUE
 20           CONTINUE
 30      LII = L(II)
         DO 40 J = I1, II
 40           A(J) = LII * L(J)
 50      CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF  L7TSQ FOLLOWS  ***
      END
      SUBROUTINE  N3RDP(IV, LIV, LV, N, P, RD, RHOI, RHOR, V)
C
C  ***  PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 ***
C
      INTEGER LIV, LV, N, P
      INTEGER IV(LIV), RHOI(*)
      REAL RD(N), RHOR(*), V(LV)
C
C     ***  NOTE -- V IS PASSED FOR POSSIBLE USE BY REVISED VERSIONS OF
C     ***  THIS ROUTINE.
C
      INTEGER COV1, I, I1, I2, IEND, II, J, K, K1, NI, PU, PX, PX1, XNI
C
C  ***  IV AND V SUBSCRIPTS  ***
C
      INTEGER BS, BSSTR, COVMAT, COVPRT, COVREQ, LOO, NB, NEEDHD, NFCOV,
     1        NFIX, NGCOV, PRUNIT, RDREQ, REGD, RCOND, STATPR, XNOTI
C
      PARAMETER (BS=85, BSSTR=86, COVMAT=26, COVPRT=14, COVREQ=15,
     1           LOO=84, NB=87, NEEDHD=36, NFCOV=52, NFIX=83, NGCOV=53,
     2           PRUNIT=21, RDREQ=57, REGD=67, RCOND=53, STATPR=23,
     3           XNOTI=90)
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      PU = IV(PRUNIT)
      IF (PU .EQ. 0) GO TO 999
      IF (IV(STATPR) .EQ. 0) GO TO 30
         IF (IV(NFCOV) .GT. 0) WRITE(PU,10) IV(NFCOV)
 10      FORMAT(/1X,I4,50H EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOST
     1ICS.)
         IF (IV(NGCOV) .GT. 0) WRITE(PU,20) IV(NGCOV)
 20      FORMAT(1X,I4,50H EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTI
     1CS.)
         IF (IV(NFCOV) .GT. 0 .OR. IV(NGCOV) .GT. 0) IV(NEEDHD) = 1
C
 30   IF (IV(COVPRT) .LE. 0) GO TO 999
      COV1 = IV(COVMAT)
      IF (IV(REGD) .LE. 0 .AND. COV1 .LE. 0) GO TO 70
      IV(NEEDHD) = 1
      IF (IABS(IV(COVREQ)) .GT. 2) GO TO 50
C
      WRITE(PU,40) V(RCOND)
 40   FORMAT(/53H SQRT(RECIPROCAL CONDITION OF F.D. HESSIAN) = AT MOST,
     1       G10.2)
      GO TO 70
C
 50   WRITE(PU,60) V(RCOND)
 60   FORMAT(/54H SQRT(RECIPROCAL CONDITION OF (J**T)*RHO"*J) = AT MOST,
     1       G10.2)
C
 70   IF (MOD(IV(COVPRT),2) .EQ. 0) GO TO 210
      IV(NEEDHD) = 1
      IF (COV1) 80,110,130
 80   IF (-1 .EQ. COV1) WRITE(PU,90)
 90   FORMAT(/43H ++++++ INDEFINITE COVARIANCE MATRIX ++++++)
      IF (-2 .EQ. COV1) WRITE(PU,100)
 100  FORMAT(/52H ++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++)
      GO TO 999
C
 110  WRITE(PU,120 )
 120  FORMAT(/45H ++++++ COVARIANCE MATRIX NOT COMPUTED ++++++)
      GO TO 210
C
 130  IF (IABS(IV(COVREQ)) .LT. 3) GO TO 150
         WRITE(PU,140)
 140     FORMAT(/35H COVARIANCE = (J**T * RHO" * J)**-1/)
         GO TO 170
 150  WRITE(PU,160)
 160  FORMAT(/56H COVARIANCE = H**-1, WHERE H = FINITE-DIFFERENCE HESSIA
     1N/)
 170  II = COV1 - 1
      DO 180  I = 1, P
         I1 = II + 1
         I2 = II +  MIN(I, 5)
         II = II + I
         WRITE(PU,190) I, (V(J), J = I1, I2)
         IF (I .LE. 5) GO TO 180
         I2 = I2 + 1
         WRITE(PU,200) (V(J), J = I2, II)
 180     CONTINUE
 190  FORMAT(4H ROW,I3,2X,5G12.3)
 200  FORMAT(9X,5G12.3)
 210  IF (IV(COVPRT) .LT. 2) GO TO 999
      I = IV(REGD) + 4
      GO TO (230, 250, 270, 290, 310), I
      WRITE(PU,220) IV(REGD)
 220  FORMAT(/18H BUG... IV(REGD) =,I10)
      GO TO 999
 230  WRITE(PU,240) NB, IV(NB)
 240  FORMAT(/17H BAD IV(NB) = IV(,I2,3H) =,I10)
      GO TO 999
 250  WRITE(PU,260) NFIX, IV(NFIX)
 260  FORMAT(/19H BAD IV(NFIX) = IV(,I2,3H) =,I10)
      GO TO 999
 270  WRITE(PU,280) LOO, IV(LOO)
 280  FORMAT(/18H BAD IV(LOO) = IV(,I2,3H) =,I10)
      GO TO 999
 290  WRITE(PU,300)
 300  FORMAT(/42H REGRESSION DIAGNOSTIC VECTOR NOT COMPUTED)
      GO TO 999
 310  IV(NEEDHD) = 1
      XNI = 0
      I = MOD(IV(RDREQ)/2, 3) + 1
      GO TO (999, 330, 320), I
 320  XNI = IV(XNOTI)
      PX = P - IV(NFIX)
      PX1 = PX - 1
      IF (IV(LOO) .GT. 1) GO TO 400
 330  WRITE(PU,340)
 340  FORMAT (74H REGRESSION DIAGNOSTIC = 0.5 * G(I)**T * H(I)**-1 * H *
     1 H(I)**-1 * G(I)...)
      IF (XNI .LE. 0) GO TO 380
      WRITE(PU, 350)
 350  FORMAT(29H     I     RD(I)         X(I))
      DO 360 I = 1, N
         WRITE(PU, 370) I, RD(I), (RHOR(J), J = XNI, XNI+PX1)
         XNI = XNI + PX
 360     CONTINUE
 370  FORMAT(1X,I5,G13.3,4G15.6/(19X,4G15.6))
      GO TO 999
C
 380  WRITE(PU,390) RD
 390  FORMAT(6G12.3)
      GO TO 999
C
 400  WRITE(PU,410)
 410  FORMAT(/77H BLOCK REGRESSION DIAGNOSTIC = 0.5 * G(I)**T * H(I)**-1
     1 * H * H(I)**-1 * G(I))
      NI = IV(NB)
      K = IV(BS)
      K1 = IV(BSSTR)
      IEND = 0
      IF (XNI .GT. 0) GO TO 450
      WRITE(PU,420)
 420  FORMAT(28H BLOCK  FIRST  LAST    RD(I))
      DO 440 I = 1, NI
         I1 = IEND + 1
         IF (I1 .GT. N) GO TO 999
         IEND = IEND + RHOI(K)
         K = K + K1
         IF (IEND .GT. N) IEND = N
         WRITE(PU,430) I, I1, IEND, RD(I)
 430     FORMAT(I6,I7,I6,G12.3)
 440     CONTINUE
      GO TO 999
C
 450  WRITE(PU, 460)
 460  FORMAT(41H BLOCK  FIRST  LAST    RD(I)         X(I))
      DO 480 I = 1, NI
         I1 = IEND + 1
         IF (I1 .GT. N) GO TO 999
         IEND = IEND + RHOI(K)
         K = K + K1
         IF (IEND .GT. N) IEND = N
         WRITE(PU,470) I, I1, IEND, RD(I), (RHOR(J), J = XNI, XNI+PX1)
 470     FORMAT(I6,I7,I6,G12.3,3G15.6/(31X,3G15.6))
         XNI = XNI + PX
 480     CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF  N3RDP FOLLOWS  ***
      END

.