SUBROUTINE D7DOG(DIG, LV, N, NWTSTP, STEP, V)
C
C  ***  COMPUTE DOUBLE DOGLEG STEP  ***
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER LV, N
      REAL DIG(N), NWTSTP(N), STEP(N), V(LV)
C
C  ***  PURPOSE  ***
C
C        THIS SUBROUTINE COMPUTES A CANDIDATE STEP (FOR USE IN AN UNCON-
C     STRAINED MINIMIZATION CODE) BY THE DOUBLE DOGLEG ALGORITHM OF
C     DENNIS AND MEI (REF. 1), WHICH IS A VARIATION ON POWELL*S DOGLEG
C     SCHEME (REF. 2, P. 95).
C
C--------------------------  PARAMETER USAGE  --------------------------
C
C    DIG (INPUT) DIAG(D)**-2 * G -- SEE ALGORITHM NOTES.
C      G (INPUT) THE CURRENT GRADIENT VECTOR.
C     LV (INPUT) LENGTH OF V.
C      N (INPUT) NUMBER OF COMPONENTS IN  DIG, G, NWTSTP,  AND  STEP.
C NWTSTP (INPUT) NEGATIVE NEWTON STEP -- SEE ALGORITHM NOTES.
C   STEP (OUTPUT) THE COMPUTED STEP.
C      V (I/O) VALUES ARRAY, THE FOLLOWING COMPONENTS OF WHICH ARE
C             USED HERE...
C V(BIAS)   (INPUT) BIAS FOR RELAXED NEWTON STEP, WHICH IS V(BIAS) OF
C             THE WAY FROM THE FULL NEWTON TO THE FULLY RELAXED NEWTON
C             STEP.  RECOMMENDED VALUE = 0.8 .
C V(DGNORM) (INPUT) 2-NORM OF DIAG(D)**-1 * G -- SEE ALGORITHM NOTES.
C V(DSTNRM) (OUTPUT) 2-NORM OF DIAG(D) * STEP, WHICH IS V(RADIUS)
C             UNLESS V(STPPAR) = 0 -- SEE ALGORITHM NOTES.
C V(DST0) (INPUT) 2-NORM OF DIAG(D) * NWTSTP -- SEE ALGORITHM NOTES.
C V(GRDFAC) (OUTPUT) THE COEFFICIENT OF  DIG  IN THE STEP RETURNED --
C             STEP(I) = V(GRDFAC)*DIG(I) + V(NWTFAC)*NWTSTP(I).
C V(GTHG)   (INPUT) SQUARE-ROOT OF (DIG**T) * (HESSIAN) * DIG -- SEE
C             ALGORITHM NOTES.
C V(GTSTEP) (OUTPUT) INNER PRODUCT BETWEEN G AND STEP.
C V(NREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE FULL NEWTON
C             STEP.
C V(NWTFAC) (OUTPUT) THE COEFFICIENT OF  NWTSTP  IN THE STEP RETURNED --
C             SEE V(GRDFAC) ABOVE.
C V(PREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE STEP RETURNED.
C V(RADIUS) (INPUT) THE TRUST REGION RADIUS.  D TIMES THE STEP RETURNED
C             HAS 2-NORM V(RADIUS) UNLESS V(STPPAR) = 0.
C V(STPPAR) (OUTPUT) CODE TELLING HOW STEP WAS COMPUTED... 0 MEANS A
C             FULL NEWTON STEP.  BETWEEN 0 AND 1 MEANS V(STPPAR) OF THE
C             WAY FROM THE NEWTON TO THE RELAXED NEWTON STEP.  BETWEEN
C             1 AND 2 MEANS A TRUE DOUBLE DOGLEG STEP, V(STPPAR) - 1 OF
C             THE WAY FROM THE RELAXED NEWTON TO THE CAUCHY STEP.
C             GREATER THAN 2 MEANS 1 / (V(STPPAR) - 1) TIMES THE CAUCHY
C             STEP.
C
C-------------------------------  NOTES  -------------------------------
C
C  ***  ALGORITHM NOTES  ***
C
C        LET  G  AND  H  BE THE CURRENT GRADIENT AND HESSIAN APPROXIMA-
C     TION RESPECTIVELY AND LET D BE THE CURRENT SCALE VECTOR.  THIS
C     ROUTINE ASSUMES DIG = DIAG(D)**-2 * G  AND  NWTSTP = H**-1 * G.
C     THE STEP COMPUTED IS THE SAME ONE WOULD GET BY REPLACING G AND H
C     BY  DIAG(D)**-1 * G  AND  DIAG(D)**-1 * H * DIAG(D)**-1,
C     COMPUTING STEP, AND TRANSLATING STEP BACK TO THE ORIGINAL
C     VARIABLES, I.E., PREMULTIPLYING IT BY DIAG(D)**-1.
C
C  ***  REFERENCES  ***
C
C 1.  DENNIS, J.E., AND MEI, H.H.W. (1979), TWO NEW UNCONSTRAINED OPTI-
C             MIZATION ALGORITHMS WHICH USE FUNCTION AND GRADIENT
C             VALUES, J. OPTIM. THEORY APPLIC. 28, PP. 453-482.
C 2. POWELL, M.J.D. (1970), A HYBRID METHOD FOR NON-LINEAR EQUATIONS,
C             IN NUMERICAL METHODS FOR NON-LINEAR EQUATIONS, EDITED BY
C             P. RABINOWITZ, GORDON AND BREACH, LONDON.
C
C  ***  GENERAL  ***
C
C     CODED BY DAVID M. GAY.
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED
C     BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324 AND
C     MCS-7906671.
C
C------------------------  EXTERNAL QUANTITIES  ------------------------
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      REAL  SQRT
C/
C--------------------------  LOCAL VARIABLES  --------------------------
C
      INTEGER I
      REAL CFACT, CNORM, CTRNWT, GHINVG, FEMNSQ, GNORM,
     1                 NWTNRM, RELAX, RLAMBD, T, T1, T2
      REAL HALF, ONE, TWO, ZERO
C
C  ***  V SUBSCRIPTS  ***
C
      INTEGER BIAS, DGNORM, DSTNRM, DST0, GRDFAC, GTHG, GTSTEP,
     1        NREDUC, NWTFAC, PREDUC, RADIUS, STPPAR
C
C  ***  DATA INITIALIZATIONS  ***
C
C/6
C     DATA HALF/0.5E+0/, ONE/1.E+0/, TWO/2.E+0/, ZERO/0.E+0/
C/7
      PARAMETER (HALF=0.5E+0, ONE=1.E+0, TWO=2.E+0, ZERO=0.E+0)
C/
C
C/6
C     DATA BIAS/43/, DGNORM/1/, DSTNRM/2/, DST0/3/, GRDFAC/45/,
C    1     GTHG/44/, GTSTEP/4/, NREDUC/6/, NWTFAC/46/, PREDUC/7/,
C    2     RADIUS/8/, STPPAR/5/
C/7
      PARAMETER (BIAS=43, DGNORM=1, DSTNRM=2, DST0=3, GRDFAC=45,
     1           GTHG=44, GTSTEP=4, NREDUC=6, NWTFAC=46, PREDUC=7,
     2           RADIUS=8, STPPAR=5)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      NWTNRM = V(DST0)
      RLAMBD = ONE
      IF (NWTNRM .GT. ZERO) RLAMBD = V(RADIUS) / NWTNRM
      GNORM = V(DGNORM)
      GHINVG = TWO * V(NREDUC)
      V(GRDFAC) = ZERO
      V(NWTFAC) = ZERO
      IF (RLAMBD .LT. ONE) GO TO 30
C
C        ***  THE NEWTON STEP IS INSIDE THE TRUST REGION  ***
C
         V(STPPAR) = ZERO
         V(DSTNRM) = NWTNRM
         V(GTSTEP) = -GHINVG
         V(PREDUC) = V(NREDUC)
         V(NWTFAC) = -ONE
         DO 20 I = 1, N
 20           STEP(I) = -NWTSTP(I)
         GO TO 999
C
 30   V(DSTNRM) = V(RADIUS)
      CFACT = (GNORM / V(GTHG))**2
C     ***  CAUCHY STEP = -CFACT * G.
      CNORM = GNORM * CFACT
      RELAX = ONE - V(BIAS) * (ONE - GNORM*CNORM/GHINVG)
      IF (RLAMBD .LT. RELAX) GO TO 50
C
C        ***  STEP IS BETWEEN RELAXED NEWTON AND FULL NEWTON STEPS  ***
C
         V(STPPAR)  =  ONE  -  (RLAMBD - RELAX) / (ONE - RELAX)
         T = -RLAMBD
         V(GTSTEP) = T * GHINVG
         V(PREDUC) = RLAMBD * (ONE - HALF*RLAMBD) * GHINVG
         V(NWTFAC) = T
         DO 40 I = 1, N
 40           STEP(I) = T * NWTSTP(I)
         GO TO 999
C
 50   IF (CNORM .LT. V(RADIUS)) GO TO 70
C
C        ***  THE CAUCHY STEP LIES OUTSIDE THE TRUST REGION --
C        ***  STEP = SCALED CAUCHY STEP  ***
C
         T = -V(RADIUS) / GNORM
         V(GRDFAC) = T
         V(STPPAR) = ONE  +  CNORM / V(RADIUS)
         V(GTSTEP) = -V(RADIUS) * GNORM
      V(PREDUC) = V(RADIUS)*(GNORM - HALF*V(RADIUS)*(V(GTHG)/GNORM)**2)
         DO 60 I = 1, N
 60           STEP(I) = T * DIG(I)
         GO TO 999
C
C     ***  COMPUTE DOGLEG STEP BETWEEN CAUCHY AND RELAXED NEWTON  ***
C     ***  FEMUR = RELAXED NEWTON STEP MINUS CAUCHY STEP  ***
C
 70   CTRNWT = CFACT * RELAX * GHINVG / GNORM
C     *** CTRNWT = INNER PROD. OF CAUCHY AND RELAXED NEWTON STEPS,
C     *** SCALED BY GNORM**-1.
      T1 = CTRNWT - GNORM*CFACT**2
C     ***  T1 = INNER PROD. OF FEMUR AND CAUCHY STEP, SCALED BY
C     ***  GNORM**-1.
      T2 = V(RADIUS)*(V(RADIUS)/GNORM) - GNORM*CFACT**2
      T = RELAX * NWTNRM
      FEMNSQ = (T/GNORM)*T - CTRNWT - T1
C     ***  FEMNSQ = SQUARE OF 2-NORM OF FEMUR, SCALED BY GNORM**-1.
      T = T2 / (T1 +  SQRT(T1**2 + FEMNSQ*T2))
C     ***  DOGLEG STEP  =  CAUCHY STEP  +  T * FEMUR.
      T1 = (T - ONE) * CFACT
      V(GRDFAC) = T1
      T2 = -T * RELAX
      V(NWTFAC) = T2
      V(STPPAR) = TWO - T
      V(GTSTEP) = T1*GNORM**2 + T2*GHINVG
      V(PREDUC) = -T1*GNORM * ((T2 + ONE)*GNORM)
     1                 - T2 * (ONE + HALF*T2)*GHINVG
     2                  - HALF * (V(GTHG)*T1)**2
      DO 80 I = 1, N
 80      STEP(I) = T1*DIG(I) + T2*NWTSTP(I)
C
 999  RETURN
C  ***  LAST LINE OF D7DOG FOLLOWS  ***
      END

.