SUBROUTINE G7QTS(D, DIG, DIHDI, KA, L, P, STEP, V, W)
C
C  *** COMPUTE GOLDFELD-QUANDT-TROTTER STEP BY MORE-HEBDEN TECHNIQUE ***
C  ***  (NL2SOL VERSION 2.2), MODIFIED A LA MORE AND SORENSEN  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER KA, P
      REAL D(P), DIG(P), DIHDI(1), L(1), V(21), STEP(P),
     1                 W(1)
C     DIMENSION DIHDI(P*(P+1)/2), L(P*(P+1)/2), W(4*P+7)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  PURPOSE  ***
C
C        GIVEN THE (COMPACTLY STORED) LOWER TRIANGLE OF A SCALED
C     HESSIAN (APPROXIMATION) AND A NONZERO SCALED GRADIENT VECTOR,
C     THIS SUBROUTINE COMPUTES A GOLDFELD-QUANDT-TROTTER STEP OF
C     APPROXIMATE LENGTH V(RADIUS) BY THE MORE-HEBDEN TECHNIQUE.  IN
C     OTHER WORDS, STEP IS COMPUTED TO (APPROXIMATELY) MINIMIZE
C     PSI(STEP) = (G**T)*STEP + 0.5*(STEP**T)*H*STEP  SUCH THAT THE
C     2-NORM OF D*STEP IS AT MOST (APPROXIMATELY) V(RADIUS), WHERE
C     G  IS THE GRADIENT,  H  IS THE HESSIAN, AND  D  IS A DIAGONAL
C     SCALE MATRIX WHOSE DIAGONAL IS STORED IN THE PARAMETER D.
C     (G7QTS ASSUMES  DIG = D**-1 * G  AND  DIHDI = D**-1 * H * D**-1.)
C
C  ***  PARAMETER DESCRIPTION  ***
C
C     D (IN)  = THE SCALE VECTOR, I.E. THE DIAGONAL OF THE SCALE
C              MATRIX  D  MENTIONED ABOVE UNDER PURPOSE.
C   DIG (IN)  = THE SCALED GRADIENT VECTOR, D**-1 * G.  IF G = 0, THEN
C              STEP = 0  AND  V(STPPAR) = 0  ARE RETURNED.
C DIHDI (IN)  = LOWER TRIANGLE OF THE SCALED HESSIAN (APPROXIMATION),
C              I.E., D**-1 * H * D**-1, STORED COMPACTLY BY ROWS., I.E.,
C              IN THE ORDER (1,1), (2,1), (2,2), (3,1), (3,2), ETC.
C    KA (I/O) = THE NUMBER OF HEBDEN ITERATIONS (SO FAR) TAKEN TO DETER-
C              MINE STEP.  KA .LT. 0 ON INPUT MEANS THIS IS THE FIRST
C              ATTEMPT TO DETERMINE STEP (FOR THE PRESENT DIG AND DIHDI)
C              -- KA IS INITIALIZED TO 0 IN THIS CASE.  OUTPUT WITH
C              KA = 0  (OR V(STPPAR) = 0)  MEANS  STEP = -(H**-1)*G.
C     L (I/O) = WORKSPACE OF LENGTH P*(P+1)/2 FOR CHOLESKY FACTORS.
C     P (IN)  = NUMBER OF PARAMETERS -- THE HESSIAN IS A  P X P  MATRIX.
C  STEP (I/O) = THE STEP COMPUTED.
C     V (I/O) CONTAINS VARIOUS CONSTANTS AND VARIABLES DESCRIBED BELOW.
C     W (I/O) = WORKSPACE OF LENGTH 4*P + 6.
C
C  ***  ENTRIES IN V  ***
C
C V(DGNORM) (I/O) = 2-NORM OF (D**-1)*G.
C V(DSTNRM) (OUTPUT) = 2-NORM OF D*STEP.
C V(DST0)   (I/O) = 2-NORM OF D*(H**-1)*G (FOR POS. DEF. H ONLY), OR
C             OVERESTIMATE OF SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1).
C V(EPSLON) (IN)  = MAX. REL. ERROR ALLOWED FOR PSI(STEP).  FOR THE
C             STEP RETURNED, PSI(STEP) WILL EXCEED ITS OPTIMAL VALUE
C             BY LESS THAN -V(EPSLON)*PSI(STEP).  SUGGESTED VALUE = 0.1.
C V(GTSTEP) (OUT) = INNER PRODUCT BETWEEN G AND STEP.
C V(NREDUC) (OUT) = PSI(-(H**-1)*G) = PSI(NEWTON STEP)  (FOR POS. DEF.
C             H ONLY -- V(NREDUC) IS SET TO ZERO OTHERWISE).
C V(PHMNFC) (IN)  = TOL. (TOGETHER WITH V(PHMXFC)) FOR ACCEPTING STEP
C             (MORE*S SIGMA).  THE ERROR V(DSTNRM) - V(RADIUS) MUST LIE
C             BETWEEN V(PHMNFC)*V(RADIUS) AND V(PHMXFC)*V(RADIUS).
C V(PHMXFC) (IN)  (SEE V(PHMNFC).)
C             SUGGESTED VALUES -- V(PHMNFC) = -0.25, V(PHMXFC) = 0.5.
C V(PREDUC) (OUT) = PSI(STEP) = PREDICTED OBJ. FUNC. REDUCTION FOR STEP.
C V(RADIUS) (IN)  = RADIUS OF CURRENT (SCALED) TRUST REGION.
C V(RAD0)   (I/O) = VALUE OF V(RADIUS) FROM PREVIOUS CALL.
C V(STPPAR) (I/O) IS NORMALLY THE MARQUARDT PARAMETER, I.E. THE ALPHA
C             DESCRIBED BELOW UNDER ALGORITHM NOTES.  IF H + ALPHA*D**2
C             (SEE ALGORITHM NOTES) IS (NEARLY) SINGULAR, HOWEVER,
C             THEN V(STPPAR) = -ALPHA.
C
C  ***  USAGE NOTES  ***
C
C     IF IT IS DESIRED TO RECOMPUTE STEP USING A DIFFERENT VALUE OF
C     V(RADIUS), THEN THIS ROUTINE MAY BE RESTARTED BY CALLING IT
C     WITH ALL PARAMETERS UNCHANGED EXCEPT V(RADIUS).  (THIS EXPLAINS
C     WHY STEP AND W ARE LISTED AS I/O).  ON AN INITIAL CALL (ONE WITH
C     KA .LT. 0), STEP AND W NEED NOT BE INITIALIZED AND ONLY COMPO-
C     NENTS V(EPSLON), V(STPPAR), V(PHMNFC), V(PHMXFC), V(RADIUS), AND
C     V(RAD0) OF V MUST BE INITIALIZED.
C
C  ***  ALGORITHM NOTES  ***
C
C        THE DESIRED G-Q-T STEP (REF. 2, 3, 4, 6) SATISFIES
C     (H + ALPHA*D**2)*STEP = -G  FOR SOME NONNEGATIVE ALPHA SUCH THAT
C     H + ALPHA*D**2 IS POSITIVE SEMIDEFINITE.  ALPHA AND STEP ARE
C     COMPUTED BY A SCHEME ANALOGOUS TO THE ONE DESCRIBED IN REF. 5.
C     ESTIMATES OF THE SMALLEST AND LARGEST EIGENVALUES OF THE HESSIAN
C     ARE OBTAINED FROM THE GERSCHGORIN CIRCLE THEOREM ENHANCED BY A
C     SIMPLE FORM OF THE SCALING DESCRIBED IN REF. 7.  CASES IN WHICH
C     H + ALPHA*D**2 IS NEARLY (OR EXACTLY) SINGULAR ARE HANDLED BY
C     THE TECHNIQUE DISCUSSED IN REF. 2.  IN THESE CASES, A STEP OF
C     (EXACT) LENGTH V(RADIUS) IS RETURNED FOR WHICH PSI(STEP) EXCEEDS
C     ITS OPTIMAL VALUE BY LESS THAN -V(EPSLON)*PSI(STEP).  THE TEST
C     SUGGESTED IN REF. 6 FOR DETECTING THE SPECIAL CASE IS PERFORMED
C     ONCE TWO MATRIX FACTORIZATIONS HAVE BEEN DONE -- DOING SO SOONER
C     SEEMS TO DEGRADE THE PERFORMANCE OF OPTIMIZATION ROUTINES THAT
C     CALL THIS ROUTINE.
C
C  ***  FUNCTIONS AND SUBROUTINES CALLED  ***
C
C  D7TPR - RETURNS INNER PRODUCT OF TWO VECTORS.
C  L7ITV - APPLIES INVERSE-TRANSPOSE OF COMPACT LOWER TRIANG. MATRIX.
C  L7IVM - APPLIES INVERSE OF COMPACT LOWER TRIANG. MATRIX.
C L7SRT  - FINDS CHOLESKY FACTOR (OF COMPACTLY STORED LOWER TRIANG.).
C  L7SVN - RETURNS APPROX. TO MIN. SING. VALUE OF LOWER TRIANG. MATRIX.
C  R7MDC - RETURNS MACHINE-DEPENDENT CONSTANTS.
C  V2NRM - RETURNS 2-NORM OF A VECTOR.
C
C  ***  REFERENCES  ***
C
C 1.  DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE
C             NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH.
C             SOFTWARE, VOL. 7, NO. 3.
C 2.  GAY, D.M. (1981), COMPUTING OPTIMAL LOCALLY CONSTRAINED STEPS,
C             SIAM J. SCI. STATIST. COMPUTING, VOL. 2, NO. 2, PP.
C             186-197.
C 3.  GOLDFELD, S.M., QUANDT, R.E., AND TROTTER, H.F. (1966),
C             MAXIMIZATION BY QUADRATIC HILL-CLIMBING, ECONOMETRICA 34,
C             PP. 541-551.
C 4.  HEBDEN, M.D. (1973), AN ALGORITHM FOR MINIMIZATION USING EXACT
C             SECOND DERIVATIVES, REPORT T.P. 515, THEORETICAL PHYSICS
C             DIV., A.E.R.E. HARWELL, OXON., ENGLAND.
C 5.  MORE, J.J. (1978), THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMEN-
C             TATION AND THEORY, PP.105-116 OF SPRINGER LECTURE NOTES
C             IN MATHEMATICS NO. 630, EDITED BY G.A. WATSON, SPRINGER-
C             VERLAG, BERLIN AND NEW YORK.
C 6.  MORE, J.J., AND SORENSEN, D.C. (1981), COMPUTING A TRUST REGION
C             STEP, TECHNICAL REPORT ANL-81-83, ARGONNE NATIONAL LAB.
C 7.  VARGA, R.S. (1965), MINIMAL GERSCHGORIN SETS, PACIFIC J. MATH. 15,
C             PP. 719-729.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C     MCS-7906671.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL RESTRT
      INTEGER DGGDMX, DIAG, DIAG0, DSTSAV, EMAX, EMIN, I, IM1, INC, IRC,
     1        J, K, KALIM, KAMIN, K1, LK0, PHIPIN, Q, Q0, UK0, X
      REAL ALPHAK, AKI, AKK, DELTA, DST, EPS, GTSTA, LK,
     1                 OLDPHI, PHI, PHIMAX, PHIMIN, PSIFAC, RAD, RADSQ,
     2                 ROOT, SI, SK, SW, T, TWOPSI, T1, T2, UK, WI
C
C     ***  CONSTANTS  ***
      REAL BIG, DGXFAC, EPSFAC, FOUR, HALF, KAPPA, NEGONE,
     1                 ONE, P001, SIX, THREE, TWO, ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      REAL  SQRT
C/
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      REAL  D7TPR,  L7SVN,  R7MDC,  V2NRM
      EXTERNAL  D7TPR,  L7ITV,  L7IVM, L7SRT,  L7SVN,  R7MDC,  V2NRM
C
C  ***  SUBSCRIPTS FOR V  ***
C
      INTEGER DGNORM, DSTNRM, DST0, EPSLON, GTSTEP, STPPAR, NREDUC,
     1        PHMNFC, PHMXFC, PREDUC, RADIUS, RAD0
C/6
C     DATA DGNORM/1/, DSTNRM/2/, DST0/3/, EPSLON/19/, GTSTEP/4/,
C    1     NREDUC/6/, PHMNFC/20/, PHMXFC/21/, PREDUC/7/, RADIUS/8/,
C    2     RAD0/9/, STPPAR/5/
C/7
      PARAMETER (DGNORM=1, DSTNRM=2, DST0=3, EPSLON=19, GTSTEP=4,
     1           NREDUC=6, PHMNFC=20, PHMXFC=21, PREDUC=7, RADIUS=8,
     2           RAD0=9, STPPAR=5)
C/
C
C/6
C     DATA EPSFAC/50.0E+0/, FOUR/4.0E+0/, HALF/0.5E+0/,
C    1     KAPPA/2.0E+0/, NEGONE/-1.0E+0/, ONE/1.0E+0/, P001/1.0E-3/,
C    2     SIX/6.0E+0/, THREE/3.0E+0/, TWO/2.0E+0/, ZERO/0.0E+0/
C/7
      PARAMETER (EPSFAC=50.0E+0, FOUR=4.0E+0, HALF=0.5E+0,
     1     KAPPA=2.0E+0, NEGONE=-1.0E+0, ONE=1.0E+0, P001=1.0E-3,
     2     SIX=6.0E+0, THREE=3.0E+0, TWO=2.0E+0, ZERO=0.0E+0)
      SAVE DGXFAC
C/
      DATA BIG/0.E+0/, DGXFAC/0.E+0/
C
C  ***  BODY  ***
C
      IF (BIG .LE. ZERO) BIG =  R7MDC(6)
C
C     ***  STORE LARGEST ABS. ENTRY IN (D**-1)*H*(D**-1) AT W(DGGDMX).
      DGGDMX = P + 1
C     ***  STORE GERSCHGORIN OVER- AND UNDERESTIMATES OF THE LARGEST
C     ***  AND SMALLEST EIGENVALUES OF (D**-1)*H*(D**-1) AT W(EMAX)
C     ***  AND W(EMIN) RESPECTIVELY.
      EMAX = DGGDMX + 1
      EMIN = EMAX + 1
C     ***  FOR USE IN RECOMPUTING STEP, THE FINAL VALUES OF LK, UK, DST,
C     ***  AND THE INVERSE DERIVATIVE OF MORE*S PHI AT 0 (FOR POS. DEF.
C     ***  H) ARE STORED IN W(LK0), W(UK0), W(DSTSAV), AND W(PHIPIN)
C     ***  RESPECTIVELY.
      LK0 = EMIN + 1
      PHIPIN = LK0 + 1
      UK0 = PHIPIN + 1
      DSTSAV = UK0 + 1
C     ***  STORE DIAG OF (D**-1)*H*(D**-1) IN W(DIAG),...,W(DIAG0+P).
      DIAG0 = DSTSAV
      DIAG = DIAG0 + 1
C     ***  STORE -D*STEP IN W(Q),...,W(Q0+P).
      Q0 = DIAG0 + P
      Q = Q0 + 1
C     ***  ALLOCATE STORAGE FOR SCRATCH VECTOR X  ***
      X = Q + P
      RAD = V(RADIUS)
      RADSQ = RAD**2
C     ***  PHITOL = MAX. ERROR ALLOWED IN DST = V(DSTNRM) = 2-NORM OF
C     ***  D*STEP.
      PHIMAX = V(PHMXFC) * RAD
      PHIMIN = V(PHMNFC) * RAD
      PSIFAC = BIG
      T1 = TWO * V(EPSLON) / (THREE * (FOUR * (V(PHMNFC) + ONE) *
     1                       (KAPPA + ONE)  +  KAPPA  +  TWO) * RAD)
      IF (T1 .LT. BIG*AMIN1(RAD,ONE)) PSIFAC = T1 / RAD
C     ***  OLDPHI IS USED TO DETECT LIMITS OF NUMERICAL ACCURACY.  IF
C     ***  WE RECOMPUTE STEP AND IT DOES NOT CHANGE, THEN WE ACCEPT IT.
      OLDPHI = ZERO
      EPS = V(EPSLON)
      IRC = 0
      RESTRT = .FALSE.
      KALIM = KA + 50
C
C  ***  START OR RESTART, DEPENDING ON KA  ***
C
      IF (KA .GE. 0) GO TO 290
C
C  ***  FRESH START  ***
C
      K = 0
      UK = NEGONE
      KA = 0
      KALIM = 50
      V(DGNORM) =  V2NRM(P, DIG)
      V(NREDUC) = ZERO
      V(DST0) = ZERO
      KAMIN = 3
      IF (V(DGNORM) .EQ. ZERO) KAMIN = 0
C
C     ***  STORE DIAG(DIHDI) IN W(DIAG0+1),...,W(DIAG0+P)  ***
C
      J = 0
      DO 10 I = 1, P
         J = J + I
         K1 = DIAG0 + I
         W(K1) = DIHDI(J)
 10      CONTINUE
C
C     ***  DETERMINE W(DGGDMX), THE LARGEST ELEMENT OF DIHDI  ***
C
      T1 = ZERO
      J = P * (P + 1) / 2
      DO 20 I = 1, J
         T =  ABS(DIHDI(I))
         IF (T1 .LT. T) T1 = T
 20      CONTINUE
      W(DGGDMX) = T1
C
C  ***  TRY ALPHA = 0  ***
C
 30   CALL L7SRT(1, P, L, DIHDI, IRC)
      IF (IRC .EQ. 0) GO TO 50
C        ***  INDEF. H -- UNDERESTIMATE SMALLEST EIGENVALUE, USE THIS
C        ***  ESTIMATE TO INITIALIZE LOWER BOUND LK ON ALPHA.
         J = IRC*(IRC+1)/2
         T = L(J)
         L(J) = ONE
         DO 40 I = 1, IRC
 40           W(I) = ZERO
         W(IRC) = ONE
         CALL  L7ITV(IRC, W, L, W)
         T1 =  V2NRM(IRC, W)
         LK = -T / T1 / T1
         V(DST0) = -LK
         IF (RESTRT) GO TO 210
         GO TO 70
C
C     ***  POSITIVE DEFINITE H -- COMPUTE UNMODIFIED NEWTON STEP.  ***
 50   LK = ZERO
      T =  L7SVN(P, L, W(Q), W(Q))
      IF (T .GE. ONE) GO TO 60
         IF (V(DGNORM) .GE. T*T*BIG) GO TO 70
 60   CALL  L7IVM(P, W(Q), L, DIG)
      GTSTA =  D7TPR(P, W(Q), W(Q))
      V(NREDUC) = HALF * GTSTA
      CALL  L7ITV(P, W(Q), L, W(Q))
      DST =  V2NRM(P, W(Q))
      V(DST0) = DST
      PHI = DST - RAD
      IF (PHI .LE. PHIMAX) GO TO 260
      IF (RESTRT) GO TO 210
C
C  ***  PREPARE TO COMPUTE GERSCHGORIN ESTIMATES OF LARGEST (AND
C  ***  SMALLEST) EIGENVALUES.  ***
C
 70   K = 0
      DO 100 I = 1, P
         WI = ZERO
         IF (I .EQ. 1) GO TO 90
         IM1 = I - 1
         DO 80 J = 1, IM1
              K = K + 1
              T =  ABS(DIHDI(K))
              WI = WI + T
              W(J) = W(J) + T
 80           CONTINUE
 90      W(I) = WI
         K = K + 1
 100     CONTINUE
C
C  ***  (UNDER-)ESTIMATE SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1)  ***
C
      K = 1
      T1 = W(DIAG) - W(1)
      IF (P .LE. 1) GO TO 120
      DO 110 I = 2, P
         J = DIAG0 + I
         T = W(J) - W(I)
         IF (T .GE. T1) GO TO 110
              T1 = T
              K = I
 110     CONTINUE
C
 120  SK = W(K)
      J = DIAG0 + K
      AKK = W(J)
      K1 = K*(K-1)/2 + 1
      INC = 1
      T = ZERO
      DO 150 I = 1, P
         IF (I .EQ. K) GO TO 130
         AKI =  ABS(DIHDI(K1))
         SI = W(I)
         J = DIAG0 + I
         T1 = HALF * (AKK - W(J) + SI - AKI)
         T1 = T1 +  SQRT(T1*T1 + SK*AKI)
         IF (T .LT. T1) T = T1
         IF (I .LT. K) GO TO 140
 130     INC = I
 140     K1 = K1 + INC
 150     CONTINUE
C
      W(EMIN) = AKK - T
      UK = V(DGNORM)/RAD - W(EMIN)
      IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK
      IF (UK .LE. ZERO) UK = P001
C
C  ***  COMPUTE GERSCHGORIN (OVER-)ESTIMATE OF LARGEST EIGENVALUE  ***
C
      K = 1
      T1 = W(DIAG) + W(1)
      IF (P .LE. 1) GO TO 170
      DO 160 I = 2, P
         J = DIAG0 + I
         T = W(J) + W(I)
         IF (T .LE. T1) GO TO 160
              T1 = T
              K = I
 160     CONTINUE
C
 170  SK = W(K)
      J = DIAG0 + K
      AKK = W(J)
      K1 = K*(K-1)/2 + 1
      INC = 1
      T = ZERO
      DO 200 I = 1, P
         IF (I .EQ. K) GO TO 180
         AKI =  ABS(DIHDI(K1))
         SI = W(I)
         J = DIAG0 + I
         T1 = HALF * (W(J) + SI - AKI - AKK)
         T1 = T1 +  SQRT(T1*T1 + SK*AKI)
         IF (T .LT. T1) T = T1
         IF (I .LT. K) GO TO 190
 180     INC = I
 190     K1 = K1 + INC
 200     CONTINUE
C
      W(EMAX) = AKK + T
      LK = AMAX1(LK, V(DGNORM)/RAD - W(EMAX))
C
C     ***  ALPHAK = CURRENT VALUE OF ALPHA (SEE ALG. NOTES ABOVE).  WE
C     ***  USE MORE*S SCHEME FOR INITIALIZING IT.
      ALPHAK =  ABS(V(STPPAR)) * V(RAD0)/RAD
      ALPHAK = AMIN1(UK, AMAX1(ALPHAK, LK))
C
      IF (IRC .NE. 0) GO TO 210
C
C  ***  COMPUTE L0 FOR POSITIVE DEFINITE H  ***
C
      CALL  L7IVM(P, W, L, W(Q))
      T =  V2NRM(P, W)
      W(PHIPIN) = RAD / T / T
      LK = AMAX1(LK, PHI*W(PHIPIN))
C
C  ***  SAFEGUARD ALPHAK AND ADD ALPHAK*I TO (D**-1)*H*(D**-1)  ***
C
 210  KA = KA + 1
      IF (-V(DST0) .GE. ALPHAK .OR. ALPHAK .LT. LK .OR. ALPHAK .GE. UK)
     1                      ALPHAK = UK * AMAX1(P001,  SQRT(LK/UK))
      IF (ALPHAK .LE. ZERO) ALPHAK = HALF * UK
      IF (ALPHAK .LE. ZERO) ALPHAK = UK
      K = 0
      DO 220 I = 1, P
         K = K + I
         J = DIAG0 + I
         DIHDI(K) = W(J) + ALPHAK
 220     CONTINUE
C
C  ***  TRY COMPUTING CHOLESKY DECOMPOSITION  ***
C
      CALL L7SRT(1, P, L, DIHDI, IRC)
      IF (IRC .EQ. 0) GO TO 240
C
C  ***  (D**-1)*H*(D**-1) + ALPHAK*I  IS INDEFINITE -- OVERESTIMATE
C  ***  SMALLEST EIGENVALUE FOR USE IN UPDATING LK  ***
C
      J = (IRC*(IRC+1))/2
      T = L(J)
      L(J) = ONE
      DO 230 I = 1, IRC
 230     W(I) = ZERO
      W(IRC) = ONE
      CALL  L7ITV(IRC, W, L, W)
      T1 =  V2NRM(IRC, W)
      LK = ALPHAK - T/T1/T1
      V(DST0) = -LK
      IF (UK .LT. LK) UK = LK
      IF (ALPHAK .LT. LK) GO TO 210
C
C  ***  NASTY CASE -- EXACT GERSCHGORIN BOUNDS.  FUDGE LK, UK...
C
      T = P001 * ALPHAK
      IF (T .LE. ZERO) T = P001
      LK = ALPHAK + T
      IF (UK .LE. LK) UK = LK + T
      GO TO 210
C
C  ***  ALPHAK MAKES (D**-1)*H*(D**-1) POSITIVE DEFINITE.
C  ***  COMPUTE Q = -D*STEP, CHECK FOR CONVERGENCE.  ***
C
 240  CALL  L7IVM(P, W(Q), L, DIG)
      GTSTA =  D7TPR(P, W(Q), W(Q))
      CALL  L7ITV(P, W(Q), L, W(Q))
      DST =  V2NRM(P, W(Q))
      PHI = DST - RAD
      IF (PHI .LE. PHIMAX .AND. PHI .GE. PHIMIN) GO TO 270
      IF (PHI .EQ. OLDPHI) GO TO 270
      OLDPHI = PHI
      IF (PHI .LT. ZERO) GO TO 330
C
C  ***  UNACCEPTABLE ALPHAK -- UPDATE LK, UK, ALPHAK  ***
C
 250  IF (KA .GE. KALIM) GO TO 270
C     ***  THE FOLLOWING AMIN1 IS NECESSARY BECAUSE OF RESTARTS  ***
      IF (PHI .LT. ZERO) UK = AMIN1(UK, ALPHAK)
C     *** KAMIN = 0 ONLY IFF THE GRADIENT VANISHES  ***
      IF (KAMIN .EQ. 0) GO TO 210
      CALL  L7IVM(P, W, L, W(Q))
C     *** THE FOLLOWING, COMMENTED CALCULATION OF ALPHAK IS SOMETIMES
C     *** SAFER BUT WORSE IN PERFORMANCE...
C     T1 = DST /  V2NRM(P, W)
C     ALPHAK = ALPHAK  +  T1 * (PHI/RAD) * T1
      T1 =  V2NRM(P, W)
      ALPHAK = ALPHAK  +  (PHI/T1) * (DST/T1) * (DST/RAD)
      LK = AMAX1(LK, ALPHAK)
      ALPHAK = LK
      GO TO 210
C
C  ***  ACCEPTABLE STEP ON FIRST TRY  ***
C
 260  ALPHAK = ZERO
C
C  ***  SUCCESSFUL STEP IN GENERAL.  COMPUTE STEP = -(D**-1)*Q  ***
C
 270  DO 280 I = 1, P
         J = Q0 + I
         STEP(I) = -W(J)/D(I)
 280     CONTINUE
      V(GTSTEP) = -GTSTA
      V(PREDUC) = HALF * ( ABS(ALPHAK)*DST*DST + GTSTA)
      GO TO 410
C
C
C  ***  RESTART WITH NEW RADIUS  ***
C
 290  IF (V(DST0) .LE. ZERO .OR. V(DST0) - RAD .GT. PHIMAX) GO TO 310
C
C     ***  PREPARE TO RETURN NEWTON STEP  ***
C
         RESTRT = .TRUE.
         KA = KA + 1
         K = 0
         DO 300 I = 1, P
              K = K + I
              J = DIAG0 + I
              DIHDI(K) = W(J)
 300          CONTINUE
         UK = NEGONE
         GO TO 30
C
 310  KAMIN = KA + 3
      IF (V(DGNORM) .EQ. ZERO) KAMIN = 0
      IF (KA .EQ. 0) GO TO 50
C
      DST = W(DSTSAV)
      ALPHAK =  ABS(V(STPPAR))
      PHI = DST - RAD
      T = V(DGNORM)/RAD
      UK = T - W(EMIN)
      IF (V(DGNORM) .EQ. ZERO) UK = UK + P001 + P001*UK
      IF (UK .LE. ZERO) UK = P001
      IF (RAD .GT. V(RAD0)) GO TO 320
C
C        ***  SMALLER RADIUS  ***
         LK = ZERO
         IF (ALPHAK .GT. ZERO) LK = W(LK0)
         LK = AMAX1(LK, T - W(EMAX))
         IF (V(DST0) .GT. ZERO) LK = AMAX1(LK, (V(DST0)-RAD)*W(PHIPIN))
         GO TO 250
C
C     ***  BIGGER RADIUS  ***
 320  IF (ALPHAK .GT. ZERO) UK = AMIN1(UK, W(UK0))
      LK = AMAX1(ZERO, -V(DST0), T - W(EMAX))
      IF (V(DST0) .GT. ZERO) LK = AMAX1(LK, (V(DST0)-RAD)*W(PHIPIN))
      GO TO 250
C
C  ***  DECIDE WHETHER TO CHECK FOR SPECIAL CASE... IN PRACTICE (FROM
C  ***  THE STANDPOINT OF THE CALLING OPTIMIZATION CODE) IT SEEMS BEST
C  ***  NOT TO CHECK UNTIL A FEW ITERATIONS HAVE FAILED -- HENCE THE
C  ***  TEST ON KAMIN BELOW.
C
 330  DELTA = ALPHAK + AMIN1(ZERO, V(DST0))
      TWOPSI = ALPHAK*DST*DST + GTSTA
      IF (KA .GE. KAMIN) GO TO 340
C     *** IF THE TEST IN REF. 2 IS SATISFIED, FALL THROUGH TO HANDLE
C     *** THE SPECIAL CASE (AS SOON AS THE MORE-SORENSEN TEST DETECTS
C     *** IT).
      IF (PSIFAC .GE. BIG) GO TO 340
      IF (DELTA .GE. PSIFAC*TWOPSI) GO TO 370
C
C  ***  CHECK FOR THE SPECIAL CASE OF  H + ALPHA*D**2  (NEARLY)
C  ***  SINGULAR.  USE ONE STEP OF INVERSE POWER METHOD WITH START
C  ***  FROM  L7SVN TO OBTAIN APPROXIMATE EIGENVECTOR CORRESPONDING
C  ***  TO SMALLEST EIGENVALUE OF (D**-1)*H*(D**-1).   L7SVN RETURNS
C  ***  X AND W WITH  L*W = X.
C
 340  T =  L7SVN(P, L, W(X), W)
C
C     ***  NORMALIZE W  ***
      DO 350 I = 1, P
 350     W(I) = T*W(I)
C     ***  COMPLETE CURRENT INV. POWER ITER. -- REPLACE W BY (L**-T)*W.
      CALL  L7ITV(P, W, L, W)
      T2 = ONE/ V2NRM(P, W)
      DO 360 I = 1, P
 360     W(I) = T2*W(I)
      T = T2 * T
C
C  ***  NOW W IS THE DESIRED APPROXIMATE (UNIT) EIGENVECTOR AND
C  ***  T*X = ((D**-1)*H*(D**-1) + ALPHAK*I)*W.
C
      SW =  D7TPR(P, W(Q), W)
      T1 = (RAD + DST) * (RAD - DST)
      ROOT =  SQRT(SW*SW + T1)
      IF (SW .LT. ZERO) ROOT = -ROOT
      SI = T1 / (SW + ROOT)
C
C  ***  THE ACTUAL TEST FOR THE SPECIAL CASE...
C
      IF ((T2*SI)**2 .LE. EPS*(DST**2 + ALPHAK*RADSQ)) GO TO 380
C
C  ***  UPDATE UPPER BOUND ON SMALLEST EIGENVALUE (WHEN NOT POSITIVE)
C  ***  (AS RECOMMENDED BY MORE AND SORENSEN) AND CONTINUE...
C
      IF (V(DST0) .LE. ZERO) V(DST0) = AMIN1(V(DST0), T2**2 - ALPHAK)
      LK = AMAX1(LK, -V(DST0))
C
C  ***  CHECK WHETHER WE CAN HOPE TO DETECT THE SPECIAL CASE IN
C  ***  THE AVAILABLE ARITHMETIC.  ACCEPT STEP AS IT IS IF NOT.
C
C     ***  IF NOT YET AVAILABLE, OBTAIN MACHINE DEPENDENT VALUE DGXFAC.
 370  IF (DGXFAC .EQ. ZERO) DGXFAC = EPSFAC *  R7MDC(3)
C
      IF (DELTA .GT. DGXFAC*W(DGGDMX)) GO TO 250
         GO TO 270
C
C  ***  SPECIAL CASE DETECTED... NEGATE ALPHAK TO INDICATE SPECIAL CASE
C
 380  ALPHAK = -ALPHAK
      V(PREDUC) = HALF * TWOPSI
C
C  ***  ACCEPT CURRENT STEP IF ADDING SI*W WOULD LEAD TO A
C  ***  FURTHER RELATIVE REDUCTION IN PSI OF LESS THAN V(EPSLON)/3.
C
      T1 = ZERO
      T = SI*(ALPHAK*SW - HALF*SI*(ALPHAK + T* D7TPR(P,W(X),W)))
      IF (T .LT. EPS*TWOPSI/SIX) GO TO 390
         V(PREDUC) = V(PREDUC) + T
         DST = RAD
         T1 = -SI
 390  DO 400 I = 1, P
         J = Q0 + I
         W(J) = T1*W(I) - W(J)
         STEP(I) = W(J) / D(I)
 400     CONTINUE
      V(GTSTEP) =  D7TPR(P, DIG, W(Q))
C
C  ***  SAVE VALUES FOR USE IN A POSSIBLE RESTART  ***
C
 410  V(DSTNRM) = DST
      V(STPPAR) = ALPHAK
      W(LK0) = LK
      W(UK0) = UK
      V(RAD0) = RAD
      W(DSTSAV) = DST
C
C     ***  RESTORE DIAGONAL OF DIHDI  ***
C
      J = 0
      DO 420 I = 1, P
         J = J + I
         K = DIAG0 + I
         DIHDI(J) = W(K)
 420     CONTINUE
C
 999  RETURN
C
C  ***  LAST CARD OF G7QTS FOLLOWS  ***
      END

.